The diffusion of doxorubicin drug molecules in silica nanoslits is non-Gaussian, intermittent and anticorrelated

This paper reports the results of all-atom Molecular Dynamics simulations of a doxorubicin molecule diffusing in water and confined to a nanoscale silica channel. Our statistical analysis reveals anomalous diffusion, due to strong interactions with the surfaces of the channel, and exhibiting a mixture of characteristic features. These include non-Gaussian and anticorrelated displacements, non-ergodic diffusion, and pronounced ageing. The results would be relevant to modelling a system of drug-delivery via nanoporous silica particles. See Ralf Metzler, Michael W. Finnis et al ., Phys . Chem . Chem . Phys ., 2020, 22 , 27955. PAPER Yasushi Ohga, Manabu Abe et al . Dynamic solvent effects in radical–radical coupling reactions: an almost bottleable localised singlet diradical ISSN 1463-9076 rsc.li/pccp PCCP


Introduction
In recent years considerable progress has been made in discovering and characterising physical and biological systems that feature non-Brownian diffusion. [1][2][3][4][5][6][7][8][9] Both single-molecule tracking experiments and computer simulations provide powerful means to study the diffusion dynamics in a variety of environments. In 2017, 3D single-molecule tracking experiments of large organic molecules diffusing on a functionalised silica surface clearly showed the presence of a hopping mechanism whereby the molecules undergo large displacements when momentarily detached and otherwise remain in states of low mobility. 10 Previous 2D tracking experiments with smaller organic molecules indicated a similar behaviour. 11 A recent study by Sarfati and Schwartz, 12 which tracked rhodamine molecules in a thin film of water condensed on a silica substrate, is probably the most directly related to our present theoretical work. They observe anticorrelated diffusion near the surface and clearly identify non-Brownian, non-Gaussian motion. The molecule's dynamics was described by fractional Brownian motion (FBM), to account for the predominance of anticorrelated displacements, additionally interrupted by intermittent surface-adsorption events described by a continuous-time random walk (CTRW). In mathematical language, this corresponds to a subordination of the antipersistent motion to the adsorption times. Previous MD studies of diffusion of organic molecules in a solution confined between silica surfaces 13,14 produced trajectories that were too short (o50 ns) to enable a detailed analysis within the framework of anomalous diffusion.
Instead of a linear growth rate for the ensemble-averaged mean-squared displacement (EAMSD) of the particles, non-Brownian systems exhibit a nonlinear power-law scaling of the EAMSD, namely [15][16][17][18][19][20][21][22][23] Here, P(r,t) denotes the probability density function of particle displacements, a is the anomalous scaling exponent, and K a is the generalised diffusion coefficient. A given diffusion process is called sub-and superdiffusive when 0 o a o 1 and a 4 1, respectively, while motion is Brownian for a = 1 and ballistic for a = 2. Subdiffusivity can have different physical origins: [15][16][17][18][19][20][21][22][23] when it arises due to a predominance of anticorrelated displacements (as occurs, for example, when a particle is diffusing in a viscoelastic environment) the dynamics is described by FBM, while when the subdiffusivity is caused by the intermittent interruption of the particle's movement such that the theoretical sojourn time is infinite, the model used to described the motion is a CTRW. Brownian motion and FBM produce Gaussian displacement distributions whose standard deviation spreads as t a/2 , while CTRWs exhibit non-Gaussian distributions. 22 Non-Gaussian Brownian motion of single-particle trajectories was first reported in 2009 24 and attributed to a diffusing object that intermittently undergoes different Brownian processes. 25 In 2017 a system was reported where, instead of a superposition of Brownian motion processes, one might be observing a superposition of FBM processes. 4 Non-Gaussianity is further discussed in ref. 26.
In this work we study the diffusion of a doxorubicin molecule (DOX) in a nanofilm of water confined between two silica surfaces. DOX is considered one of the most effective anticancer drugs ever developed, 27 being routinely used in treatment. 28 However, its toxicity can lead to cardiac complications and failure, 28 limiting the benefits of the therapy. One of the primary prevention strategies is to use liposomal DOX, whereby DOX is encapsulated within nanospheres. 28 Commercialised in 1996, one such formulation was the first successful application of nanotechnology for drug delivery. 27 Since mesoporous silica was first suggested as a suitable material to carry drug molecules in drug-delivery applications, 29 a number of research groups have explored this option, as reviewed in ref. 30 and 31, partly because of the ease with which the silica surface can be functionalised and the porosity controlled. For example, the release rate of DOX from hollow mesoporous silica spheres can be manipulated by tuning the pore sizes, while a lower pH, such as that found in cancer cells, enhances the rate of release. 32 Controlling the release rate is important in order to keep drug concentrations in the serum at an optimal level, thus improving the efficiency of the therapy by avoiding high-and low-concentration peaks. 30 It is also important to achieve low release rates in the bloodstream and enhanced release at the target (exploiting, for example, the lower pH in cancer cells), in the absence of a more sophisticated design enabling other mechanisms to trigger release, such as, for example, that in ref. 33.
A central question in such applications is how efficiently DOX exits the silica carrier by diffusion. Though single-molecule tracking experiments of DOX diffusing inside mesoporous silica do exist, 34 the data resolution is insufficient to provide a decisive insight regarding the underlying physical diffusion mechanism. Such information is important to understand better the effects of the molecular interactions and also to scale the observed dynamics to longer times and greater lengthscales. The motivation of our present work is therefore to better understand the diffusion of DOX in nanoconfinement in order to address this question. In our simulations, confinement is only between opposite planar surfaces, however this enables us to discover how nanoconfinement affects the diffusion process. The separation between the two surfaces was chosen to match the pore diameters of the nanoparticles reported in ref. 35, based on the synthesis method developed in ref. 36. The insights gained are interesting and unexpected: they could inform anomalous-diffuion studies of other systems as well as studies of DOX diffusing in cylindrical pores.
The paper is organised as follows. In Section 2 we present the details of the atomistic MD simulations. In Section 3 we present the main results and statistical data analysis. We consider the EAMSD and time-averaged MSD (TAMSD), the ageing properties of the underlying random walk, the distribution of particle displacements, the distribution of waiting times, and the position autocorrelation function. We also discuss possible mathematical and physical models of diffusion compatible with the data. In Section 4 we summarise our findings from the perspective of physical transport mechanisms and discuss some implications.

Details of computer simulations
Standard Gromacs 2018.2 was used to perform the MD simulations. The force field parameters describing the DOX molecule were defined using the Antechamber package 37 and the general Amber force field (GAFF) 38 for organic molecules. Discrete partial atomic charges were assigned using the constrained electrostatic potential energy methodology, 39 consistent with the GAFF parameterisation strategy, and the quantum-mechanical electrostatic potential calculated with the Gaussian09 software. 40 To achieve this, the structure of DOX was extracted from Protein Data Bank ID: 4DX7 (which is correctly protonated at physiological pH, see Fig. 1a) and the geometry optimised with the hybrid B3LYP functional and the 6-31G* basis set. 41,42 Subsequently, the quantum-mechanical electrostatic potential of DOX was calculated using the Hartree-Fock formulation and the 6-31G* basis set, which is known to reproduce the relative free energies of solvation of organic molecules. 43 Finally, a least squares procedure is used to fit a charge to each atomic centre in the molecule. 39 To simulate the silica surface we used the INTERFACE force field, 44,45 known to be suitable for simulations of biomolecules in solution in contact with, for instance, metals and ceramics. The force field is compatible with TIP3P 46 and flexible SPC water models; the possibility to combine the silica surface with a flexible SPC water model is an advantage compared to the other commonly used silica force field for simulating organic molecules near silica, 47 which is only compatible with TIP3P water.
One of these water models SPC/Fw 48 yields the self-diffusion coefficient of water similar to that measured experimentally. 48,49 The TIP3P water model, in contrast, predicts a diffusion coefficient more than twice the experimental value. 48,50 This enhanced water diffusivity is also known to affect the mobility of organic molecules: for example, the diffusivity of glucose in TIP3P water confined in a silica channel was overestimated by 30%. 13 Another study of protein diffusivity in TIP3P and SPC/Fw water models further demonstrated the advantage of SPC/Fw 50 and the SPC/Fw model was also used recently to study the interactions of amino acids with a calcium-carbonate surface. 51 The structure of the amorphous silica surface was taken from the Surface Model Database created by the developers of the INTERFACE force field. 45 The surface has 4.7 silanol groups per nm 2 of which 0.45 silanols per nm 2 are unprotonated; instead, the charge is compensated by a sodium ion. This type of surface is found in most porous silica immersed in a slightly acidic solution (above pH 5). 45 The silica slab in our simulations has periodic boundary conditions, and its parallel surfaces are separated by a slab of water in the {x,y}-plane, see The velocity Verlet scheme was used to integrate the equations of motion. The stretching of the covalent bonds involving hydrogen atoms in the system was constrained to be zero using the LINCS method. 52 The correct treatment of this motion would require a quantum-mechanical treatment of the protons, while within the present classical treatment a fixed hydrogen-oxygen bond length provides a better representation than an oscillating bond. 53 Without the need to resolve ultrafast oscillations, the time step of the simulation was set to 2 fs.
We first performed a short equilibration of the system using a Parinello-Rahman barostat in order to adjust the water density and remove a thin layer of vacuum introduced between the water and the silica. The next phase of equilibration was done separately for each simulation. Namely, the COM of the DOX molecule was tethered to the centre of the channel using a harmonic potential, in order to avoid surface adsorption before the start of the data collection. Velocities were randomly assigned to each atom of the system corresponding to a Maxwell distribution at 298 K. The temperature was controlled using a Nose-Hoover thermostat. 54 After 5 ns of simulation time, the force restraining the COM was removed and the data recording started. This point marks the beginning of each simulated trajectory.
Velocity-scaling thermostats, such as the Nose-Hoover thermostat, are most suitable for approximating the real diffusive dynamics. 55 However, when a rescaling thermostat is applied, the kinetic energy from the high-frequency modes is gradually transferred to low-or zero-frequency modes (such as the COM velocity). As a result, the system can ''freeze'' while it acquires a large translational velocity, a phenomenon called the flying ice cube effect. 56 To avoid this, the COM velocities of the liquid, including DOX, Na + , Cl À and H 2 O species, and the atoms of the silica surface, were removed at every time step in the simulations. The fifty simulations were performed in single nodes of 32 CPUs, with a daily yield of 75 ns per trajectory.

Surface adsorption leads to intermittent diffusion
The molecule was seen to intermittently adsorb to the surface and adopt a large number of configurations with respect to the surface, remaining in these positions of low mobility for long periods of time, often lasting hundreds of nanoseconds and in some cases for the entire one microsecond of simulation time. A sample trajectory is shown in Fig. 3. In Fig. 4 we depict four sample configurations in which the DOX molecule is close to the silica surfaces. In some cases the molecule's backbone is perpendicular to the surface (with some atoms locked in a  pocket of the silica surface), see Fig. 4d, whereas in other situations it can be positioned flat with respect to the surface ( Fig. 4a and b). When adsorbed the amplitude of the molecule's COM ''oscillations'' is often less than 1 Å along the {x,y}-plane, with sporadic larger oscillations, see inset in Fig. 5c, corresponding to the adsorbed molecule in Fig. 4a. The exact pattern of the oscillations depends on the adsorption configuration.
Surface-detachment events are rare and often accompanied by quick lateral diffusion of the molecule in the channel, see Fig. 5 and 6. This two-state process with alternating periods of adsorbed motion and much less restricted bulk-like diffusion in the channel is physically similar to the picture of CTRWs: 22,57 these feature distinct sojourn periods (here the bound states) and jump events. A similar phenomenon of intermittent, almost fully immobile states was observed, for instance, for the motion of Ka channels in the membranes of living human kidney cells. [58][59][60] The density of water perpendicular to the channel axis varies non-monotonically as the silica surface is approached, first increasing slightly before decreasing steadily into the surface, which indicates a degree of water structuring typical for liquidsolid interfaces, see Fig. 2. When adsorbed, DOX displaces the water layer from the silica surface (notice the decrease in the water density in Fig. 2).

Mean-squared displacements show anomalous diffusion
To quantify the above observations we plot the ensemble-and time-averaged MSD (EAMSD and TAMSD respectively) of DOX, see Fig. 7. The EAMSD is computed from    TAMSD is defined as a sliding-window average along the trajectory where T is the total length of the trajectory over which the measurement is made and Dt the time step by which the window slides. The mean TAMSD is the average over all trajectories, The EAMSD starts linearly and after E1 ns turns into a subdiffusive behaviour with exponent a E 0.4. In stark contrast, the mean TAMSD starts with exponent a E 0.5 at short lag times and has an extended region of nearly normal diffusion at intermediate lag times, with a E 0.9, see Fig. 7a. This disparity between EAMSD and TAMSD is often called weak ergodicity breaking 22,61,62 and is not equivalent to non-ergodicty in the stronger mixing sense. 63 Even a well defined process such as Brownian motion will display stochastic variations of the TAMSD from one trajectory to the next when measured over a finite time period. 22,62 Similar fluctuations can also be seen in the single-trajectory power spectrum. 64,65 For non-Brownian motion the distribution of such amplitude fluctuations are often characteristic to identify the specific underlying motion, e.g., to distinguish FBM from a CTRW process. 22,66 The spread of the individual TAMSD is here quantified in Fig. 7b based on We observe that at very short lag times (D = 10 ps) the distribution is approximately bell-shaped but the maximum occurs just above the average. At lag time D = 1 ns the distribution has a clear dip in the probability at x = 1. At later lag times, at D = 100 ns, the difference between TAMSDs becomes more pronounced, with the emergence of traces with x 4 3 and a large fraction of realisations with x close to zero. This broad distribution suggests subdiffusive CTRWs as a possible model for the observed diffusion. For subdiffusive CTRWs, however, theoretically the TAMSD should grow strictly linearly with lag time and the MSD should have a sublinear scaling with time, [17][18][19]22 in contrast to our observations presented in Fig. 7. In fact the combination of CTRW with other processes, such as FBM and diffusion on a fractal, was observed in the single-particle tracking data from ref. 58 and 66.
Ageing is a another characteristic of subdiffusive CTRWs: the longer the time, T, allowed for the particle to move, the greater the probability it gets trapped for a long period of time (of the order of the recorded trajectory). In Fig. 5c, for example, the molecule displaces by more than 6 nm along the x-axis during the first 300 ns (intermittently undergoing short adsorption events) and then adsorbs in a very stable configuration (see Fig. 4a) for the remaining 700 ns. As these long periods of immobility are included in the sliding-window average, the TAMSD decreases. 62 This is observed in Fig. 8a. A good fit by equation which describes the ageing behaviour of subdiffusive CTRWs, 62 is achieved for a E 0.66, see Fig. 8b. While only over a decade in time, the agreement with the scaling law DT aÀ1 is quite remarkable. Another manifestation of ageing is observed when an initial segment of the trajectory (termed the ageing time t a ) is discarded and as a result the TAMSD decreases. In Fig. 8c, we take a segment of T = 400 ns and calculate the mean TAMSD when this segment starts at t = 0 (as in Fig. 8a) and at t = t a . As t a increases we observe a decrease in the mean TAMSD. The results after normalising the magnitude of the aged TAMSD to that of the non-aged for D = 10 ns and different values of T were plotted in Fig. 8d as a function of the universal rescaled variable t a /T and fitted to a well defined relation characteristic of subdiffusive CTRWs, 67-69 with a = 0.65, If we look at the percentage of time the molecule spends at different heights across the channel, Fig. 9, we see that during the first 100 ns the molecule is more likely (by an order of magnitude) to be in the central region of the channel compared to the last 100 ns of the trajectory. This is another illustration of the phenomenon discussed above: as the system evolves, the likelihood of finding the DOX molecule adsorbed in a very stable configuration increases. This is also clear from the sharp increase in the peak closest to the left surface and the overall increase in the probability of finding DOX close to the right surface.

Distribution of DOX displacements
The distribution of DOX displacements turns from an almost Gaussian shape at short times (1 ps) to a clearly non-Gaussian distribution at longer times in Fig. 10. A stretched Gaussian, which has been previously used to characterise non-Gaussian distributions of displacements 70,71 Pðx; tÞ $ exp À 1 d xðtÞ sðtÞ d " # was fitted to our data. The results of the fit depend on the exact binning used for the respective histograms, see ref. 6, as well as on whether the regions of small or large displacements are in the focus of the analysis. Fitting the region of small displacements, which means 1D displacements of this system, which is isotropic in the {x,y}-plane, we find that after time intervals of 10, 50, 100, 200, and 500 ps the best-fit exponents decrease as d E 1.69, 1.47, 1.38, 1.26, and 1.14, respectively. Fitting the tails of the displacements distributions after the same time intervals we find d E 1.2, 0.85, 0.75, 0.6, and 0.5, respectively. We note a rather dramatic decrease of d in both situations. Concurrently, the standard deviation s(t), defined by eqn (6), grows with diffusion time for these two fitting scenarios in a way inconsistent with the MSD growth. We thus do not use this fitting approach.
As an alternative, we use the sum of four Gaussians with different amplitudes and widths. In order to obtain a fit independent of the bin width, the Gaussians were integrated over consecutive bin widths, producing discrete probability points, which were then fitted to the data. The probability of measuring a displacement of size, x, in the interval (a,b) is therefore given by where w i is the weight of the i Gaussian curve and s i (t) its standard deviation at time t.
The four Gaussian curves produce an excellent fit, as illustrated in Fig. 11. While the argument of Occam's razor would favour fewer adjustable parameters, the following discussion develops our reasoning to apply a description based on four Gaussians. Two Gaussian curves were insufficient to describe the tails of the distribution, three curves produced a good fit at   short lag times but could not describe the tails correctly at larger times, while four produced an excellent fit for all computed distributions. Importantly, this choice yields a scaling of the standard deviation, encoded in this displacement distribution, consistent with the TAMSD growth, see Fig. 12a. Moreover, the fractional contribution of each Gaussian to the distribution remains approximately constant, see inset of Fig. 12a (compare  ref. 3 and 6). Namely, the average weights of the fitted curves (in the order of increasing mobility) are % w i = 44.8, 42.2, 11.8, and 1.2%. The narrowest Gaussian stops widening after E 100 ps. The time evolution of the remaining three Gaussian curves could be fitted with the relation s i (t) B t gi where g i E 0.18, 0.32, and 0.43 respectively, see Fig. 12a. The spreading dynamics is therefore slower than Brownian motion (g = 0.5) for all cases. The different spreading rates would explain the change in the distribution's shape observed in Fig. 10.
From the macroscopic symmetry of the system, in the limit of infinitely long time and large simulation cell, the molecule would spend equal times in regions equally distant from the midplane of the channel, and we see in Fig. 12b that this symmetry of the distribution versus height in the channel is not exactly satisfied in the simulations. In Fig. 12b the channel is divided into four sections: a central region (green) and three regions (yellow, red and brown), each one closer to the surface than the previous one. The boundaries of these sections are chosen to match features in the distribution of time versus distance: for example, the regions closest to the surface end at the maxima of the distribution. This subdivision of the channel in four sections serves to illustrate our physical interpretation of the four fitted Gaussians by comparing their weights in the fit with the probability of finding the COM of the molecule at certain heights in the channel: % w 1 = 1.2% (DOX spends 1.2% of the time in green channel region), % w 2 = 12.9% (11.8% in the yellow region), % w 3 = 42.2% (40.3% in the red region) and % w 4 = 44.8% (45.6% in the brown region).
In the region highlighted in green, the molecule's COM is at least 7.8 Å from the surface; since DOX is E1.5 nm long, in this region the molecule is not in direct contact with the surface. The two Gaussians with the highest weights and slowest spreading dynamics correspond to the distribution peaks close to each surface.
Taken literally, this model would correspond to a division of the channel into layers 1, 2, 3 and 4, exhibiting four different mobilities of the diffusing molecule, between which the molecule hops at random intervals. Layer 1 would be a sub-channel centred on the mid-plane of the channel, where the mobility is highest, and layer 4 would be adjacent to the walls of the  channel, where the mobility is most retarded by attachments. Layers 2 and 3 would fill the space between layers 1 and 4. Although there is no such clear discrete layering of the channel in our simulation, the physical picture we have is very similar if we interpret the four Gaussians to be a discrete representation of the continuously increasing mobility along the channel as a function of the distance of the molecule's COM from either surface.

Distribution of attachment times
A molecule is considered immobile if it stays for a given time t esc within a region of the {x,y}-plane with the escape radius r esc . For these two parameters the values of t esc = 0.5, 1 ns and r esc = 1, 2, 3, 5, and 10 Å were used in the analysis.
We obtain long tails Bt À(1+g) (leading to infinite theoretical sojourn times for g o 1), see Fig. 13. Such scale-free distributions of immobilisation times are observed, e.g., in ref. 58 and 66. However the exponent g is shown to decrease from E0.65 for r esc = 3 Å to E0.1 for r esc = 5 Å (both computed for t esc = 0.5 ns). Although functionally similar to the CTRW-process, the scaling exponent g varies strongly, complicating a CTRW-based interpretation of this data, see below for more details.

Displacement-autocorrelation function
A sensitive quantifier of the underlying diffusion process is the displacement autocorrelation function (ACF) of the drug molecules defined as where e is the time interval over which the displacement is measured, t represents the timelapse over which the correlation of the displacements is measured, t 0 represents the beginning of each trace fragment used to calculate the average and N is the number of these fragments. In Fig. 14 we demonstrate that the displacement-autocorrelation function reveals antipersistent features, corresponding to the negative part of the correlation function. Physically this means that a displacement in one direction is likely followed by one in the opposite direction. This is consistent with the slower spreading dynamics of the Gaussian curves fitted to the displacement distributions when compared with Brownian motion. The data for e = 10 ps is compared in Fig. 14 to a subdiffusive FBM process with a scaling exponent of a E 0.47, producing a good fit. Note that a subdiffusive unconfined CTRW does not exhibit negative values in its displacement-autocorrelation function. 17,22 The depth of the minimum decreases for larger e values approaching a constant negative value asymptotically, see the inset of Fig. 14. In contrast, the FBM model predicts a minimum value independent of e. Similar behaviour was recently reported for the diffusion of rhodamine molecules confined in a water layer between a silica surface and a vapour phase. 12 The mininum value, in the FBM description, also approaches asymptotically a constant negative value, which is interpreted as an indication of true anticorrelation. 12 Previously a dependence of the minimum value on e had also been reported for the diffusion of nanoparticles in crowded dextran solutions, see ref. 72. Here, however, the minimum value approached zero asymptotically, which the authors attributed to a transient FBM process evolving towards Brownian motion. 72,73 We measure the frequency with which n consecutive correlated displacements (or steps) of the drug molecule occur. For this calculation we consider the x and y directions separately. The results should be independent of the direction along which the steps were measured, so for optimal statistics an average was taken over both axes. The occurrence of each sequence of correlated steps is normalised by the number of times the steps change from being positively correlated to negatively correlated or vice versa. For example if we have the following six displacements: (+ + + À À À) we have two positively correlated steps (+ + +), followed by one negatively correlated step (+ À), followed by two Fig. 13 Normalised frequency of escape events from the circle of radius r esc in the {x,y}-plane for r esc = 3 Å (black, left axis) and r esc = 5 Å (blue, right axis), with t esc = 0.5 ns for both panels. Both distributions are shown together in the inset (sharing the same y-axis). positively correlated steps (À À À). The correlation therefore changed twice. In Fig. 15a we see that the decay in the correlation is independent of the time interval used to define the displacement, as is the case of the FBM model.
The asymmetry of the curve, which exhibits a clear predominance of negatively correlated steps, highlights the already observed deviation from Brownian motion. The data coincides well with an FBM diffusion process with a = 0.30, see Fig. 15b, although there are some discrepancies at the tails. This could be due to the fact that we have a superposition of FBM processes or due to additional effects from the complex DOX-wall interactions.

Discussion
We studied in detail the motion of a single DOX molecule in a nanochannel with silica surfaces. Our simulations shed light on previously inaccessible time-and length-scales, providing new insight into diffusion in nanoconfinement.
We observed that the DOX molecule adsorbs for prolonged periods on the silica surfaces and diffuses rapidly when detached. Once it adsorbs in a stable configuration it often remains in that position for the rest of the trajectory. This leads to so-called ageing, where on average the TAMSD is larger if measured at the beginning of the trajectory, (the molecule is less likely to have found a stable adsorption configuration), than when it is calculated from a segment later in the trajectory (when it is more likely to have adsorbed for a prolonged period of time). The motion of DOX, with its prolonged interruptions, therefore exhibits subdiffusive CTRW traits. The specific characteristics of the oscillatory motion induced by the presence of the surface could be causing the FBM-like character (with pronounced negative displacement correlations when adsorbed and almost vanishing correlations in the central region of the channel). The simulated data also exhibit properties that are not consistent with either of these two stochastic models. For example, while the overall shape of the autocorrelation function matches the FBM prediction, its minimum value depends on e. A summary of the properties observed in our data compared to those of the FBM and CTRW models is shown in Table 1. Features of diffusion arising, as ours do, from different physical mechanisms within a system, have been observed previously. 12,58,66,74 We fitted the diffusion of DOX by a time-independent superposition of four Gaussian processes, evolving at a rate slower than that associated with Brownian motion. The fit by four Gaussian curves is suggestive of the model discussed by Wang et al., 25 a particle intermittently undergoing different Gaussian processes leading to an overall non-Gaussian displacement distribution. In reality the molecule is moving through a continuum of Gaussian diffusion regimes parallel to the surfaces as it explores the width of the channel.
It has been demonstrated theoretically that ageing (specifically: a decrease in the TAMSD with increasing T as d 2 ðDÞ can arise from the superposition of two Gaussian Brownian processes when a significant difference between the diffusion coefficients exists and the theoretical average time the diffusing object spends in the slow mode is infinite. 75 Our system suggests a variation of this scenario, where we have a superposition of subdiffusive-FBM-like processes. At short lag times (in the order of several picoseconds), the difference in the displacements associated with these different processes is not pronounced, resulting in a FBM-like nonlinear TAMSD. As the lag times become longer the differences become larger and the a exponent approaches the CTRW value of 1. In fact, the trajectories (see Fig. 6), might appear similar to those of the noisy CTRW, 76 where the movements during periods of attachment are described by Ornstein-Uhlenbeck noise.
In this system we might be observing a more complex motion due to the size of the molecule, with its 69 atoms sampling different parts of the channel simultaneously. Eliazar et al., 77 describe a scenario in which a superposition of Ornstein-Uhlenbeck processes can lead to FBM. A variation of this description could explain the origin of the observed FBM-type dynamics. Alternatively, as suggested in ref. 12, this FBM motion could be an indication that the water near the silica surface is viscoelastic.  Experimental studies measuring the percentage of DOX released from hollow mesoporous silica nanoparticles over a period of 30 hours showed that release was faster as the pore size was increased from 3.6 nm to 12.6 nm. 32 As the pore diameter increases, the central region associated with fast diffusion increases leading to an overall increased diffusion rate. Furthermore, in a sufficiently large channel, Brownian motion would be recovered in the central region.

Outlook
In order to facilitate the design of new drug-carriers with optimal drug-release rates, it is useful to be able to extract the diffusion coefficient of DOX (or other drug molecules) through a variety of silica membranes (or other porous materials) either from singlemolecule tracking experiments or computer simulations. The ultimate goal would be to produce mathematical expressions that predict the drug release rate at timescales that are relevant for medical applications and can be compared with experimental data.
This study provides insight into the complexity of the diffusion process even in the simplest scenario, and shows that an analysis in which the diffusion coefficient is extracted simply from the classical Brownian motion equation, as in the single-molecule tracking experiments of DOX in silica nanopores, 34 is not valid. Furthermore, this anomalous diffusion approach might prove advantageous compared to models that try to approach the problem by considering classical diffusion combined with an adsorption/desorption process (such as for example in ref. 78 and 79). Within the anomalous diffusion approach, all adsorption/desorption information is contained within the distribution of displacements and attachment times, the MSD and correlation analysis, thus avoiding the challenge of modelling the adsorption process explicitly.
We finally remark here that diffusion was modelled with a finite set of different Brownian motions with different diffusivities in the context of NMR signals, 80 and an analogous model has recently been proposed for single-particle dynamics. 81,82

Conflicts of interest
There are no conflicts to declare.