Diffusion is a basic transport mechanism that exists widely in nature, and plays an important role in maintaining life activities . The mean squared displacement 〈r2(t)〉 (MSD) of a tracer is a central measured quantity in both experiments and simulations. How 〈r2(t)〉 varies with the time t determines three different situations . When the exponent α of the power-law relationship feature 〈r2(t)〉 ∝ tα is larger than 1, the motion of a tracer is called superdiffusion. The diffusion is normal as α = 1, while it is anomalous when α < 1.
Compared to the diffusion of solutes in porous medium, biomacromolecules diffuse in complex and crowded intracellular environments. For example, the cytoplasm is a soft core colloid made up of biomacromolecules occupying a volume fraction up to 30%-35%, such as ribosomes, proteins, and RNA [3-5]. It has shown that biomembranes are exceptionally crowded with proteins [6-8]. Therefore, the diffusive properties of a tracer in crowded environments have been a hot research topic for decades [6-42]. Benefitting from the significant advances in experimental techniques, the diffusion of molecules in living cells could be measured directly [11-24]. Computer simulations have been proven to be a powerful tool in studying the diffusion of molecules in complex and crowded environments [6, 25-42].
The diffusion rate of proteins was found to be reduced in crowded environments since the molecular crowding resulted in an increasing effective viscosity of the medium . Furthermore, it has been shown that the molecular crowding also leads to the anomalous diffusion of tracers [7-10]. A typical example is the work of Banks et al. , in which the streptavidin diffuses anomalously in highly concentrated solutions containing dextran acting as crowders. However, Dix et al.  reported that both the diffusion of the proteins in in vitro solutions and in living cells are normally Brownian with a reduced diffusion rate. The complexity of intracellular environments and the fact that the subdiffusive behaviors observed in in vivo experiments can not always be reproduced by in vitro ones motivates us to clarify whether the molecular crowding itself could give rise to the anomalous diffusion of tracers.
Besides this basic question, how the size of tracers and crowders and the concentration of crowders affect the diffusive properties of a tracer is another issue. Banks et al.  showed that the subdiffusion exponent of streptavidin decreased as the size of dextran (crowders) got larger. However, a recent theoretical work by Sharp  suggested that large crowders were less effective on causing crowding than smaller ones. Then, there comes a question. That is, if it is the molecular crowding that causes the anomalous diffusion of tracers, how the subdiffusion exponent changes with the increasing size of crowders with other conditions being the same. In fact, we will show later that it is the relative size of tracers to crowders rather than their absolute sizes that determines the magnitude of the subdiffusion exponent. By using 2D Langevin dynamics simulations in this work, we attempt to get answers of these doubts proposed above.Ⅱ. MODEL AND METHODS
We consider a two-dimensional geometry as a square box with a side length L=100σ0. Here σ0 is the length scale of the system. A certain area fraction of crowders move inside the box, and are treated as mobile and impenetrable disks. The size of the tracer and crowders in unit of σ0 is denoted as σt and σc, respectively.
The purely repulsive interactions between the tracer-crowder and the crowder-crowder are described by a truncated Lennard-Jones (LJ) potential
Here ε is the depth of the potential well. The LJ radius used to calculate the tracer-crowder and the crowder-crowder interactions is σLJ=(σt+σc)/2 and σLJ=σc, respectively.
In the Langevin dynamics simulations, the motion of each particle is governed by the Newton's laws of motion 
where −∇ULJ, ξivi and FiR is the conservative, frictional, and random forces, respectively. mi, ξi, and vi is the mass, friction coefficient and velocity of the ith particle, respectively. Note that ξi could be time-dependent in a underdamped system . However, we considered an overdamped system in the present work so that ξi keeps constant for a specified particle and does not vary with the time. The mass of a tracer is mt=m0(σt/σ0)2 with m0 being the mass of a particle with a size of σ0. Certainly, mc=m0(σc/σ0)2 for crowders. For the tracer, ξi=ξ0(σt/σ0) with ξ0 being the friction coefficient of a particle of size σ0. For crowders, ξc=ξ0(σc/σ0) as it should be. The Langevin equation is then integrated in time by the method described by Ermak and Buckholz . The periodic boundary conditions are applied to both directions.
The parameters ε and m0 fix the energy and mass scales of the system, respectively. ε and m0 together with σ0 determine the corresponding time scale tLJ=(m0σ02/ε)1/2 and force scale f=ε/σ0. In our simulations, the thermal energy is set as kBT =1.0ε with kB being the Boltzmann constant and T being the absolute temperature.
The tracer is initially placed at the center of the box. Then, crowders are randomly added into the box. The tracer and crowders go through thermal collisions described by the Langevin thermostat to obtain an equilibrium configuration, which is taken as the initial configuration of our simulations. Then, 1000 independent trajectories of the tracer in 1×108 time steps are obtained with the time step △t=0.005tLJ.
The subdiffusive behavior of the tracer can be distinguished from a normal one by its nonlinear relationship between the mean squared displacement (MSD) and the diffusion time t. The MSD is defined as
where ri(0) and ri(t) is the initial position and the position at time t of the tracer in the ith trajectory, respectively. In our simulations, 〈r2(t)〉 is obtained by averaging over n=1000 independent trajectories. Namely, the MSDs of the tracer in the simulations are ensemble averaged.Ⅲ. RESULTS AND DISCUSSION
For the case that no crowder is in the simulation box, the solution of the Ornstein-Uhlenbeck process described in Eq.(2) for the tracer without LJ potential is well known as 
where d is the spatial dimensionality. mt/ξt is the characteristic time. At short time scales (t ≪ mt/ξt), Eq.(4) can be rewritten in an initial ballistic form,
and when t ≫ mt/ξt, Eq.(4) can be rewritten in the form of overdamped Brownian motion.
In our 2D system, the area fraction of crowders
Here Nc is the number of crowders. Considered that there is only a tracer in the simulation box, the average distance between the center of two neighboring crowders is estimated as L/Nc1/2. Thus, the average spacing between neighboring crowders with respect to their surfaces R is
By substituting Eq.(7) into Eq.(8), we should have
It is clearly shown by Eq.(9) that R increases as the size of crowders σc gets larger at a fixed
In the following, we investigate the anomalous diffusion of a tracer at intermediate time scales and its normal diffusion at long time scales in crowded environments. Specifically, the friction coefficient of the medium ξ0, the area fraction of crowders
As plotted in Fig. 1(a), the MSDs of a tracer scaled by the time 〈r2(t)/t display three time-dependent regimes. At short time scales, ballistic diffusions of the tracer with 〈r2(t)〉 ∝ t2 are observed. At long time scales, normal diffusions of the tracer with 〈r2(t)〉 ∝ t are recovered. These two regimes replicate the crossover behavior of 〈r2(t)〉 predicted by Eq.(5) and Eq.(6). Note that this time-dependent feature in a 〈r2(t)/t vs. t curve has also been observed in previous studies [12, 24, 32-34]. Interestingly, at intermediate time scales, the transient anomalous diffusions of the tracer appear clearly, characterized by the negative slopes β < 0 of 〈r2(t/t vs. t curves. According to the values of β, the subdiffusion exponent α could be obtained since α=1+β.
It has been reported very recently by Sentjabrskaja et al.  that the slow mobility of the matrix is crucial for the observation of the anomalous diffusion of a tracer, where the system approaches the glass transition with the volume fraction of large particles
It is indicated by our simulation results here that the anomalous diffusion of a tracer can be induced by the crowding effect alone. But what is responsible for this phenomenon in physics? To clarify this question, we have calculated the probability to find a tracer of a diameter σt=4.0 having traveled a distance r in a sea of crowders of a diameter σc=3.0 at different time steps, ρ(r, t'). As shown in Fig. 2(a), peaks in ρ(r, t') appear with the evolution of time. Notably, the peak at t'=1000 time steps is broadened and decays slower compared with that at t'=100 time steps. It indicates that the tracer spends more time in a local trap. Furthermore, ρ(r, t') curves show a non-Gaussian feature, which is consistent with a very recent work of Ghosh et al. . To have a direct impression of this behavior, a typical trajectory of the tracer in 1000 time steps is presented in Fig. 2(b). It can be seen clearly that the motion of the tracer is highly localized at most of its time accompanied with a jump event between the two local traps. In the following, this behavior is referred to the "cage effect". With this, we can come to a conclusion that the trapping of the tracer in one location for varying and widely distributed periods combined with the infrequent "jumps" between locations gives rise to the anomalous diffusion of the tracer. As the matrix mobility gets slower caused by increasing ξ0 or
Next, we turn to investigate the effect of the sizes of the tracer and crowders on the anomalous diffusion behavior. Here, the size ratio of the tracer to crowders defined as δ=σt/σc is adopted as the independent variable. Obviously, δ characterizes the size disparity between the tracer and crowders. Figure 3 shows clearly that the subdiffusion exponent α has a minimum αmin=0.75 as a function of δ. The emergence of αmin could be understood by considering the following two limiting cases. For δ ≪ 1, the tracer should not be trapped by crowders since it can diffuse normally in the void spaces unoccupied by them. While for δ ≫ 1, the crowders in the system act as a role similar to solvent molecules, and thus the tracer should display a normal diffusive behavior.
Furthermore, α=0.91 observed at the smallest δ=1/3 we have studied indicates only slightly anomalous diffusion. Generally, we can certainly make the claim that the diffusion of a tracer is anomalous as α ≤ 0.90. Figure 3 shows that there exist two critical size ratios: one is the lower critical size ratio δc, 1 ≈ 0.37, and the other is the upper critical size ratio δc, 2 ≈ 4.0. If δc; 1 ≤ δ ≤ δc; 2, the tracer should exhibit a subdiffusive behavior in crowded environments. Our observation of δc, 1 here is in good agreement with the work of Sentjabrskaja et al. , in which they reported the existence of a lower critical size ratio δc ≈ 0.35 obtained from the simulations. They also got a slightly smaller experimental value of δc ≈ 0.28, which is attributed to the polydispersity of tracers used in experiments . The value of δc, 1 could be estimated in theory as follows. For the trapping of a tracer by surrounding crowders, the size of tracer σt should be larger than the average spacing between neighboring crowders with respect to their surfaces R,
Thus, we get an area fraction of crowders-dependent
The normal diffusion of a tracer is expected if δ gets large further. We noted that in an experimental study of Banks et al. , when the streptavidin with a hydrodynamic radius of 4.9 nm diffuses in a solution with different concentrations of dextrans with a radius of gyration Rg=0.82 nm, its subdiffusion exponent α are always larger than 0.95. At this point, the size ratio of the tracer (streptavidin) to crowders (dextran) is δ ≈ 5.98. Obviously, the diffusion of streptavidin could be deemed as almost normal in this case, which is a good proof of our simulation results.
Interestingly, the minimum αmin=0.75 is achieved at δ=1. Namely, the subdiffusion of a tracer is the most pronounced as its size equals to that of crowders. Through collisions with surrounding crowders, a tracer breaks the trapping cage formed by them to diffuse out. Assuming that the momentum is conserved at the moment of a collision between the tracer and a crowder, the energy loss caused by the collision would be the largest at δ=1. In the worst case, the tracer is instantaneously stationary after the collision. Clearly, this is very unfavorable for the breaking of the trapping cage for the tracer, and it is the reason why αmin appears at this time. Moreover, the minimum αmin=0.75 is in accord with the experimental observation by Banks et al. , where the subdiffusion exponent has a limit value α=0.74 ± 0.02.B. Normal diffusion at long time scale
As stated above, a tracer undergoes the transition from a subdiffusive behavior at intermediate time scales to a diffusive one at long time scales in crowed environments. The recovery of the normal diffusion of a tracer is signalled by the platform in the 〈r2(t)/t vs. t curves, as plotted in Fig. 1. The diffusion coefficient of a tracer in crowed environments D can be calculated from the plateau value.
There are two length scales that are relevant for the diffusion of a tracer in crowed environments, i.e., the size of tracer σt and the average spacing between neighboring crowders with respect to their surfaces R. The dimensionless quantity R/σt characterizing the confinement condition of a tracer in crowded environments is referred to the confinement parameter. A smaller R/σt suggests a more pronounced confining effect. The reduced diffusion coefficient D/D0 as a function of R/σt is plotted in Fig. 4(a) with σc=4.0 and σt varying from 1.0 to 4.0. Here, D0 is the diffusion coefficient of the tracer in the case that no crowders are in the system. D/D0 almost keeps at a constant near 1 for R/σt ≥ 3, indicating the diffusion of the tracer is almost not affected. However, for R/σt < 3, a significant decrease in D/D0 is observed. The decrease is more rapid as R/σt gets smaller due to the more and more pronounced crowding and confining effect of the systems. Similar observations have been reported in polymer nanocomposite systems [51-54].
By plotting D/D0 as a function of ln (R/σt), we find that D/D0 increases linearly with ln (R/σt) for ln (R/σt) < 1, see Fig. 4(b). In addition, the data points fall on a master curve as δ < 1, which is also indicated by Fig. 4(a). Considered that the area fraction of crowders ϕ decreases with δ as ϕ ≈ (δ+1) -2, a larger δ means a lighter degree of crowding at the same value of ln (R/σt). The collapse of data points at δ ≤ 1 indicates that the relative mobility of tracers with σt < σc does not rely on the degree of crowding. In contrast, the relative mobility of tracers with σt > σc is sensitive to the degree of crowding.Ⅳ. CONCLUSION
To summarize, we have investigated the anomalous diffusion of a tracer at intermediate time scales and its normal diffusion at long time scales in crowded environments by using 2D Langevin dynamics simulations. The emergence of subdiffusion of a tracer depends on the size ratio of the tracer to crowders δ. If δ falls between a lower critical size ratio δc, 1 and a upper one δc, 2, the anomalous diffusion appears purely due to the molecular crowding. Further analysis indicates that the physical origin of subdiffusion is the "cage effect". Moreover, the subdiffusion exponent α decreases with the increasing medium viscosity and the degree of crowding, and gets a minimum αmin=0.75 at δ=1. As to the normal diffusion of a tracer at long time scales, we find that the diffusion coefficient D keeps almost the same as D0 for R/σt ≥ 3 followed by a significant decrease for R/σt < 3. More importantly, the relative mobility of a tracer in crowded environments relies on the magnitude of δ. For δ ≤ 1, the relative mobility is independent of the degree of crowding. In contrast, it is sensitive to the degree of crowding for δ > 1. These results are helpful in deepening the understanding of the diffusive properties of biomacromolecules that lie within crowded intracellular environments, such as proteins, DNA and ribosomes.Ⅴ. ACKNOWLEDGMENTS
This work is supported by the National Natural Science Foundation of China (No.21225421 and No.21474099), the National Basic Research Program of China (No.2014CB845605).
|||R. Phillips, J. Kondev, and J. Theriot, Physical Biology of the Cell, London: Taylor and Francis Group (2008).|
|||B. D. Hughes, Random Walks and Random Environ-ments. Oxford, UK: Clarendon Press (1995).|
|||S.B. Zimmerman, and S.O. Trach, J. Mol. Biol. 222 , 599 (1991). DOI:10.1016/0022-2836(91)90499-V|
|||R.J. Ellis, and A.P. Minton, Nature 425 , 27 (2003). DOI:10.1038/425027a|
|||A. P. Minton, J. Cell Sci. 119 , 2863 (2006). DOI:10.1242/jcs.03063|
|||J. H. Jeon, M. Javanainen, H. Martinez-Seara, R. Metzler, and I. Vattulainen, Phys. Rev X6 , 021006 (2016).|
|||R. Metzler, J. H. Jeon, and A. G. Cherstvy, Biochem. Biophys. Acta BBA-Biomembr. 1858 , 2451 (2016). DOI:10.1016/j.bbamem.2016.01.022|
|||A. Godec, and R. Metzler, Phys. Rev E92 , 010701(R) (2015).|
|||R. Metzler, and J. Klafter, Phys. Rep. 339 , 1 (2000). DOI:10.1016/S0370-1573(00)00070-3|
|||F. Höfling, and T. Franosch, Rep. Prog. Phys. 76 , 046602 (2013). DOI:10.1088/0034-4885/76/4/046602|
|||T. J. Feder, I. Brust-Mascher, J. P. Slattery, B. Baird, and W. W. Webb, Biophys. J. 70 , 2767 (1996). DOI:10.1016/S0006-3495(96)79846-6|
|||D.S. Banks, and C. Fradin, Biophys. J. 89 , 2960 (2005). DOI:10.1529/biophysj.104.051078|
|||J. Dix, and A.S. Verkman, Annu. Rev. Biophys. Biomol. Struct. 37 , 247 (2008). DOI:10.1146/annurev.biophys.37.032807.125824|
|||O. Seksek, J. Biwersi, and A. S. Verkman, Biophys. J. 75 , 557 (1998). DOI:10.1016/S0006-3495(98)77545-9|
|||N. Periasamy, and A.S. Verkman, J. Cell Biol. 37 , 6316 (2009).|
|||M. Ario-Dupont, G. Foucault, M. Vacher, F. Devaux, and S. Cribier, Biophys. J. 78 , 901 (2000). DOI:10.1016/S0006-3495(00)76647-1|
|||M. Platani, I. Goldberg, J. R. Swedlow, and A. I. Lamond, J. Cell Biol. 151 , 1561 (2000). DOI:10.1083/jcb.151.7.1561|
|||A. S. Verkman, Science 27 , 27 (2002).|
|||E. O. Potma, W.P. de Boeij, L. Bosgraaf, J. Roelofs, P.J.M. Van Haastert, and D. A. Wiersma, Biophys. J. 81 , 2010 (2001). DOI:10.1016/S0006-3495(01)75851-1|
|||Y. Cheng, R.K. Prud'homme, and T. L. Thomas, Macromolecules 35 , 8111 (2002). DOI:10.1021/ma0107758|
|||M. Platani, I. Goldberg, J. R. Swedlow, and A. I. Lamond, Nat. Cell Biol. 4 , 502 (2002). DOI:10.1038/ncb809|
|||M. Wachsmuth, T. Weidemann, G. MÄuller, M.W. Hoffman-Rohrer, T. A. Knoch, W. Waldeck, and J. Langowski, Biophys. J. 84 , 3353 (2003). DOI:10.1016/S0006-3495(03)70059-9|
|||N. Fatin-Rouge, K. Starchev, and J. Buffle, Biophys. J. 86 , 2710 (2004). DOI:10.1016/S0006-3495(04)74325-8|
|||E. Vilaseca, I. Pastor, A. Isvoran, S. Madurga, J. L. Garces, and F. Mas, Theor. Chem. Acc. 128 , 795 (2011). DOI:10.1007/s00214-010-0840-5|
|||M. Weiss, M. Elsner, F. Kartberg, and T. Nilsson, Biophys. J. 87 , 3518 (2004). DOI:10.1529/biophysj.104.044263|
|||J.P. Bouchaud, and A. Georges, Phys. Rep. 195 , 127 (1990). DOI:10.1016/0370-1573(90)90099-N|
|||A. Y. Grosberg, Phys. Rev. Lett. 85 , 3858 (2000). DOI:10.1103/PhysRevLett.85.3858|
|||J. H. Jeon, V. Tejedor, S. Burov, E. Barkai, C. Selhuber Unkel, K. Berg-Sørensen, L. Oddershede, and R. Metzler, Phys. Rev. Lett. 106 , 048103 (2011). DOI:10.1103/PhysRevLett.106.048103|
|||K. Nakazato, and K. Kitahara, Prog. Theor. Phys. 64 , 2261 (1980). DOI:10.1143/PTP.64.2261|
|||Y. He, S. Burov, R. Metzler, and E. Barkai, Phys. Rev. Lett. 101 , 058101 (2008). DOI:10.1103/PhysRevLett.101.058101|
|||E. Barkai, Y. Garini, and R. Metzler, Phys. Today 65 , 29 (2012).|
|||J. H. Jeon, H. Martinez-Seara Monne, M. Javanainen, and R. Metzler, Phys. Rev. Lett. 109 , 188103 (2012). DOI:10.1103/PhysRevLett.109.188103|
|||S. K. Ghosh, A. G. Cherstvy, and R. Metzler, Phys. Chem. Chem. Phys. 17 , 1847 (2015). DOI:10.1039/C4CP03599B|
|||(a) M. J. Saxton, Biophys. J. 52, 989 (1987); (b) M. J. Saxton, Biophys. J. 58, 1303 (1990); (c) M. J. Saxton, Biophys. J. 64, 1053 (1993); (d) M. J. Saxton, Biophys. J. 66, 394 (1994); (d) M. J. Saxton, Biophys. J. 81, 2226 (2001).|
|||E. Vilaseca, A. Isvoran, S. Madurga, I. Pastor, J. L. Garces, and F. Mas, Phys. Chem. Chem. Phys. 13 , 7396 (2011). DOI:10.1039/c0cp01218a|
|||S. Burov, J. H. Jeon, R. Metzler, and E. Barkai, Phys. Chem. Chem. Phys. 13 , 1800 (2011). DOI:10.1039/c0cp01879a|
|||R. Metzler, J. H. Jeon, A. G. Cherstvy, and E. Barkai, Phys. Chem. Chem. Phys. 16 , 24128 (2014). DOI:10.1039/C4CP03465A|
|||T. Miyaguchi, and T. Akimoto, Phys. Rev E91 , 010102(R) (2015).|
|||H. Berry, and H. Chaté, Phys. Rev E89 , 022708 (2014).|
|||H. J. Kim, Phys. Rev E90 , 012103 (2014).|
|||T. Neusius, I. M. Sokolov, and J. C. Smith, Phys. Rev E80 , 011109 (2009).|
|||F.D. A. Aarão Reis, and D. di Caprio, Phys. Rev E89 , 062126 (2014).|
|||D. S. Banks, C. Tressler, R. D. Peters, F. Höfling, and C. Fradin, Soft Matter 12 , 4190 (2016). DOI:10.1039/C5SM01213A|
|||K. A. Sharp, Proc. Natl. Acad. Sci. USA 112 , 7990 (2015). DOI:10.1073/pnas.1505396112|
|||M.P. Allen Allen, and D.J. Tildesley, Computer Simulation of Liquids. New York: Oxford University (1987).|
|||A. S. Bodrova, A. V. Chechkin, A. G. Cherstvy, H. Safdari, I. M. Sokolov, and R. Metzler, Sci. Rep. 6 , 30520 (2016). DOI:10.1038/srep30520|
|||D.L. Ermak, and H. Buckholz, J. Comput. Phys. 35 , 169 (1980). DOI:10.1016/0021-9991(80)90084-4|
|||W. T. Coffey and Y. P. Kalmykov, The Langevin Equa-tion, Singapore: World Scientific, (2012).|
|||S. K. Ghosh, A. G. Cherstvy, D. S. Grebenkov, and R. Metzler, New J. Phys. 18 , 013027 (2016). DOI:10.1088/1367-2630/18/1/013027|
|||T. Sentjabrskaja, E. Zaccarelli, C. De Michele, F. Sciortino, P. Tartaglia, T. Voigtmann, S. U. Egelhaaf, and M. Laurati, Nat. Commun. 7 , 11133 (2016). DOI:10.1038/ncomms11133|
|||S. Gam, J. S. Meth, S. G. Zane, C. Chi, B.A. Wood, M. E. Seitz, K. I. Winey, N. Clarke, and R. J. Composto, Macromolecules 44 , 3494 (2011). DOI:10.1021/ma102463q|
|||S. Gam, J. S. Meth, S. G. Zane, C. Chi, B. A. Wood, K. I. Winey, N. Clarke, and R. J. Composto, Soft Matter 8 , 6512 (2012). DOI:10.1039/c2sm25269d|
|||C. Lin, S. Gam, J. S. Meth, N. Clarke, K. I. Winey, and R. J. Composto, Macromolecules 46 , 4502 (2013). DOI:10.1021/ma4000557|
|||J. Choi, M.J. A. Hore, N. Clarke, K. I. Winey, and R. J. Composto, Macromolecules 47 , 2404 (2014). DOI:10.1021/ma500235v|