Non-universal tracer diffusion in crowded media of non-inert obstacles

We study the diffusion of a tracer particle, which moves in continuum space between a lattice of excluded volume, immobile non-inert obstacles. In particular, we analyse how the strength of the tracer-obstacle interactions and the volume occupancy of the crowders alter the diffusive motion of the tracer. From the details of the partitioning of the tracer diffusion modes between trapping states when bound to obstacles and bulk diffusion, we examine the degree of localisation of the tracer in the lattice of crowders. We study the properties of the tracer diffusion in terms of the ensemble and time averaged mean squared displacements, the trapping time distributions, the amplitude variation of the time averaged mean squared displacements, and the non-Gaussianity parameter of the diffusing tracer. We conclude that tracer-obstacle adsorption and binding triggers a transient anomalous diffusion. From a very narrow spread of recorded individual time averaged trajectories we exclude continuous type random walk processes as the underlying physical model of the tracer diffusion in our system. For moderate tracer-crowder attraction the motion is found to be fully ergodic, while at stronger attraction strength a transient disparity between ensemble and time averaged mean squared displacements occurs. We also put our results into perspective with findings from experimental single-particle tracking and simulations of the diffusion of tagged tracers in dense crowded suspensions. Our results have implications for the diffusion, transport, and spreading of chemical components in highly crowded environments inside living cells and other structured liquids.


I. INTRODUCTION
Macromolecular crowding (MMC) abounds in living biological cells, with up to φ ≈ 30 . . . 35% of the volume of the cytoplasmic liquid being occupied by large biopolymers such as proteins, nucleic acids, ribosomes, as well as membranous structures, and other complexes [1][2][3][4][5]. These volume-excluding and often non-inert obstacles alter the diffusion behaviour of cellular components and the rates of biochemical reactions taking place in this highly complex liquid [6][7][8][9][10][11]. These changes occur both due to an enhanced solution viscosity [12] and the sheer physical obstruction imposed on particle diffusion due to the presence of the obstacles.
The mean squared displacement (MSD) of a tracer particle in such crowded solutions often becomes anomalous [13][14][15][16][17][18] where D β is the anomalous diffusion coefficient of dimension cm 2 /sec β and β the anomalous diffusion exponent. Its typical range 0 < β < 1 indicates slower-than-Brownian, subdiffusive motion [13,14]. Often, the subdiffusion (1) turns out to be transient and at long times the MSD crosses over to Brownian Motion [18,19]. A wide range of scaling exponents β ∼ 0.4 . . . 0.9 has been reported for the obstructed diffusion of tracers of various sizes and surface properties in crowded solutions inside cells. Examples include the motion of small proteins [20], mRNA molecules [21], telomeric chromosomal loci [22], Cajal bodies [23], lipid and insulin granules [24], and natural virus particles [25,26] in the cytoplasm of living cells. Further examples are membrane lipids and membrane-bound proteins [27][28][29], water molecules associated to membranes [30], hair bundles in ears [31], as well as several examples for the motion of tracers such as proteins in complex liquids [32][33][34]. Subdiffusive regimes in crowded systems have also been observed and modelled for actively driven particles [35], and the dependence of the effective diffusivity D(φ) reveals a minimum as function of the MMC volume fraction φ [36].
Anomalous diffusion of the form (1) is modelled in terms of a wide range of stochastic processes [13][14][15][16][17][18][19]. These include continuous time random walks (CTRW) [13,14,37], fractional Brownian motion [16,17,38,39] and the closely related fractional Langevin equation motion [16,17,[39][40][41], as well as diffusion processes with space [42][43][44][45] and time [45][46][47] dependent diffusion coefficients. CTRW models are closely related to trap models, in which the tracer is successively immobilised [13,[48][49][50]. Despite advances in simulations [29,51] and theoretical approaches [16,17], there is no consensus on the physical understanding of the subdiffusion of passive tracer particles in crowded solutions, that would be directly applicable to the cytoplasm of living cells [15,18]. In particular, it is likely that different physical origins dominate for different tracer sizes and shapes, as well as length and time scales of the diffusion. The lack of a consensus picture for diffusion in MMC environments suggests that the observed anomalous diffusion is not universal but depends on specific parameters.
To address such specific origins for deviations from Brownian motion, extensive computer simulations of passive tracer diffusion were performed by several groups. The recent approaches of Refs. [52,53], for instance, considering the tracer motion in lattices of immobile, ran-domly positioned obstacles or regularly ordered obstacles jiggling in a confining potential demonstrated that the anomalous diffusion regime is governed by the obstacle volume occupancy, with significant subdiffusive motion observed at higher crowder volume fraction φ. This subdiffusion is transient and can be quantified by the dependence of the local anomalous diffusion exponent β(φ, t) = d log r 2 (φ, t) d log(t) (2) and the effective diffusion coefficient D(φ, t), where we included the explicit dependencies on the MMC volume fraction φ and time t. In denser obstacle lattices, diffusion is more localised and the value of β(t) smaller [53]. Remarkably, the transient subdiffusion regime was shown to disappear at higher obstacle mobility [52,53]. Moreover, the distribution of particle trapping times in dynamical cages formed between the crowders observed in simulations for mobile, confined, and static obstacles was shown to be inconsistent with CTRW models [52,54] Recent experimental advances in the field of obstructed obstacle-mediated tracer diffusion, including the regimes of transient anomalous diffusion, are presented in Refs. [55][56][57].
The study of obstructed diffusion by means of simulations was pioneered by Saxton in a series of papers [58][59][60]. In particular, the effects of the fraction of lattice sites occupied by crowders and of their diffusivity were examined. Specifically, effects of tracer-obstacle binding on the anomalous diffusion properties were studied and connected to a binding energy landscape for immobile point-like obstacles positioned at a fixed concentration on the lattice [59]. Obstructed diffusion of point-like tracers in a lattice of randomly positioned, static obstacles was investigated in Ref. [61] and shown to give rise to a reduction of the tracer diffusivity D with the obstacle concentration. Nonzero values of the diffusivity D even for very densely packed obstacles appear in such a model due to the existence of a percolation structure. True long-time subdiffusion can only be realised at the percolation threshold in such lattices of randomly distributed obstacles [18]. On a cubic lattice, the critical percolation thresholds corresponds to 31% [18]. For a random walker on the infinite incipient cluster, the scaling exponent of the MSD is β ≈ 0.697 [62]. The physical reason is the formation of a labyrinthine-like environment [19], in which the tracer needs to escape dead ends and cross narrow causeways present on all scales. For lattices below the percolation transition as well as for regularly positioned obstacles, as in the current study, the anomalous diffusion regime is transient.
An interesting alternative to the modelling of transient anomalous diffusion are Lorentz gas-based models which were developed to exploit the localisation transitions on a percolation network of overlapping spherical obstacles [63,64]. The scaling relations for the suppression of the tracer diffusivity as the system approaches the critical percolation densityφ was determined, namely D(φ) ∼ [(φ −φ)/φ] µ , where the percolation exponent is µ ≈ 2.88 [63]. At the percolation threshold, persistent anomalous diffusion with exponent β = 2/6.25 ≈ 0.32 was found, for even denser systems the particles are eventually localised [63].
Here, we extend the class of systems considered by Saxton [59] based on transient binding of tracer particles to physical obstacles. We perform extensive Monte-Carlo simulations of tracer diffusion on 3D lattices of sticky spherical obstacles of varying radius R. In this obstruction-binding diffusion model we examine the trapping time distributions of the tracers, the time averaged MSD -which is a more relevant observable when compared to experimental situations than the ensemble averaged MSD (1) -and the effective tracer diffusivity D. The model parameters are systematically varied, including the crowder radius R and thus the volume fraction of crowders φ and the tracer-crowder binding energy ǫ A . A schematic of the system is shown in Fig. 1 along with sample trajectories of a tracer particle. These novel features substantially extend the known simulations results for obstructed tracer diffusion on 2D lattices of reflecting spherical [52,53] and cylindrical [65] obstacles. Despite the difference of mobile polymer obstacles to our scenario of ordered reactive crowders, our results show interesting similarities with the diffusion of tracer particles in dense solutions of non-inert polymer chains recently reported by the Holm group [66]. We will discuss the consequences of this similarity below.
The paper is organised as follows. In Section II we introduce the basic notations and the quantities to be analysed. We outline the computational scheme and theoretical concepts. In Section III we report the main simulations results and support them by theoretical scaling arguments. We analyse the effects of the MMC volume fraction and the strength of the obstacle-tracer binding. Moreover, we compute the ensemble and time averaged particle displacements as well as the distributions of particle trapping times to the sticky obstacles. To rationalise the stochastic behaviour and determine the concrete underlying effective diffusion model, we also systematically compute the non-Gaussianity parameter G. In Section IV the conclusions are drawn and possible applications of our results to some experimental systems discussed.

II. SIMULATION MODEL AND APPROXIMATIONS
To mimic the conditions of a crowded environment, we consider a primitive cubic lattice every site of which is occupied by a spherical obstacle, as shown in Fig. 1. The maximal size of the obstacle R max for the conditions of close packing is R max = a/2 where a is the lattice constant. In the following we will use a = 2 in dimensionless units. The maximal volume occupancy by obstacles on such a static cubic lattice is φ max = π/6 ≈ 0.524, com- Figure 1: Schematic of the three-dimensional tracer-obstacle system used in our simulations, for the obstacle radius R = 0.6 and tracer-obstacle binding strength ǫA = 2, 6, and 10 (from left to right). The tracer trajectories, as obtained directly from the simulations, are of the same length in all the three panels. Note that the particle traces are rendered with a finite thickness although the tracer particle is point-like. The fraction of time that the particle spends in the surface-bound diffusion mode grows with ǫA. pared to φ max = π/ √ 18 ≈ 0.740 for the densest packing of spheres in 3D [67]. The obstacles are considered immobile in our simulations.
In the simulations presented below, the point-like tracer starts in the centre of a cage, at the maximal distance from the eight surrounding obstacles. At the very start the tracer particle thus performs free motion, until it encounters the surface of a crowding particle to which it can subsequently bind. The length scale of the spatial heterogeneity l ⋆ in the system is of the order of the free path of the tracer between neighbouring obstacles, As we will show such heterogeneities effect subdiffusion at intermediate time scales t ⋆ , while at much longer time-scales the diffusion becomes Brownian, as demonstrated in Fig. 3. As many binding-unbinding events take place during the length of the recorded traces in our simulations, the initial particle position does not affect the long time dynamics. The tracer particle becomes adsorbed onto the obstacle surface with the binding energy ǫ A and stays in the bound state for the average adsorption time t ads,i . While bound, the tracer diffuses along the spherical obstacle surface with the same diffusion coefficient as in the free unbound state, that is, it moves along the surface of a crowder sphere with D ads = D 0 . The tracer is considered unbound once it separates from the obstacle for more than the distance 0.1R max , see also the definition of the interaction potential below. [87] In our simulations we place a single tracer on the crowder lattice and then average over many individual traces. The repeated binding and unbinding events separating the particle motion between surface and bulk diffusion lead to an effective distribution between the modes of tracer motion, which is reflected in the tracer particle MSD and other diffusive characteristics such as average trapping times, see below.
The attractive interactions between the mobile tracer particle and immobile crowder spheres are modelled in terms of the Lennard-Jones 6-12 potential, whose attract-ive branch is cut off at the distance r cutoff , namely (3) The parameter σ is connected to the obstacle radius by R = 2 1/6 σ corresponding to the minimum of U LJ (r). The attractive potential U LJ (r) is truncated at the critical distance r cutoff = R + 0.1R max thus mimicking an attractive shell of a fixed width 0.1R max around the obstacles. Thus, the thickness of the attractive layer around obstacles of different sizes is the same. The vertical energy shift ǫ cutoff is a constant which sets U LJ (r cutoff ) = 0.
We simulate the motion of the point-like tracer of mass m with coordinate r(t) in the presence of friction based on the Langevin equation Here γ is the friction coefficient coupled to the strength of the Gaussian noise through where δ(·) denotes the Dirac delta-function and k B T represents the thermal energy. The noise has zero mean, and has vanishing correlations in the different Cartesian directions. The sum in Eq. (4) runs over the positions of all crowding particles R j . The fact that we consider a point-like particle is not a severe restriction, as a finite size of the tracer particle would correspond to a re-normalisation of the crowder radius, compare also Ref. [61]. Concurrently the surface diffusivity of the tracer along the crowders would need to be adjusted. In our simulations we neglect tracer-obstacle hydrodynamic interactions, which can, in principle, affect the long-time behaviour of the system. In particular, because of their long-range 1/r-nature [10], the diffusing particles can feel the obstacles at a finite distance without direct collisions, see e.g. Ref. [36]. Note however that hydrodynamic interactions have recently also been demonstrated to affect the short-time tracer diffusion dynamics in fluids [69,70].
In free space the solution of the Ornstein-Uhlenbeck process (4) without the last term is well known [71], describing the crossover from initial ballistic motion before the characteristic time m/γ to overdamped, Brownian motion This solution is reproduced by our simulations, see Fig. 2.
The local scaling exponent β(t) is computed as discretised logarithmic derivative from the MSD traces obtained from the simulations (i.e., we compute the local derivative of log(MSD) with respect to the logarithmically sampled time, see also Eq. (3) of Ref. [57]). The particle mass is m = 1 throughout the paper (we made sure that the code works fine for varying particle mass and solution friction, as shown in Fig. 2). In all figures below, time t is shown in units of the simulation step δt = 0.001 of the Verlet velocity integration scheme, the displacements appear in units of the lattice constant a. The obstacle size below is given in terms of the maximal geometrically allowed radius on the square lattice. Note that even zero-sized crowders R = 0 have an attractive shell of finite width around them, as determined by the specific nature of the attractive potential (3). In the text, however, when we talk about "free diffusion", no crowders are included in the simulations at all. In the presence of the sticky, excluded volume obstacles an exact solution is not known, and we thus analyse this case by simulations. We find that the fraction of time that the tracer particle spends in the surface-bound mode increases with the volume occupancy by obstacles and with the tracer-obstacle affinity ǫ A . The effect of MMC on the long-time particle diffusivity D(φ) reveals a non-trivial dependence at larger ǫ A , as shown below.

III. RESULTS
We discuss the simulations results with respect to three main quantities. In section III A we study the ensemble averaged MSD and the associated effective diffusivity. We study the MSD r 2 (t) of the tracer particles at varying volume occupancy φ of obstacles and the tracerobstacle affinity ǫ A . The main results are shown in Figs. 3 and 4. Generally, we observe that for small interaction strengths ǫ A the anomalous diffusion exponent β(t) varies along the time evolution of the MSD r 2 (t) . At short times it starts out with the underdamped ballistic motion (7) of the above Ornstein-Uhlenbeck process, as demonstrated in Fig. 2. Such a short time superdiffusive behaviour is due to inertial effects and was indeed observed experimentally, for instance, in single particle tracking studies of fluorescent beads in sucrose solutions [57]; see also the detailed studies of Refs. [69,70] of inertial effects for the particle diffusion in a fluid. Subsequently, Fig. 3 demonstrates that the MSD crosses over to a transient subdiffusive regime with 0 < β(t) < 1. In the long time limit the tracer particle performs Brownian motion with a linear scaling of the MSD and β(t) = 1. Concurrently this Brownian motion regime is characterised by a reduced diffusivity D(φ) as compared to unrestricted Brownian motion of the tracer. The dependence of the effective diffusivity D(φ) on the MMC volume fraction φ is shown for different interaction strengths in Fig. 4. For weakly adhesive obstacles, the diffusion becomes monotonically slower for larger crowders positioned on the lattice, see Fig. 4.
Let us look at these behaviours in some more detail. The variation of the scaling exponent β(t) with the crowding fraction and tracer-crowder attraction strength is shown in Fig. 3. We observe that as both ǫ A and φ increase, the transient subdiffusion regime progressively extends over a larger time window. Remarkably, this anomalous diffusion spans up to two decades of time for substantial tracer-obstacle attraction strengths, see the plot of β(t) for ǫ A = 6k B T in Fig. 3. For relatively weak tracer-obstacle attraction, at ǫ A = 2k B T , only marginally anomalous tracer diffusion was detected, with β ∼ 0.85 . . . 1. Overall, the subdiffusion regime extends over one to two decades of time, similar to the results for inert static, randomly-positioned obstacles [61] or for the motion of tracers in a random channel with sticky surfaces [72].
A physically similar renormalisation of the particle diffusivity was discovered in Ref. [73] for glassy states in sticky-particle systems at relatively large volume fractions φ. The phase diagram of the hard sphere mixture with a short range inter-particle attraction as well as a self-diffusive MSD dependence were examined, for instance, by simulations in Ref. [73]. The implications of the square-well sphere-sphere interaction potential ǫ A and volume fraction φ of crowders were rationalised in detail. A progressive slowing down of the particle selfdiffusion in the attractive hard-sphere mixtures as functions of φ and ǫ A observed in Figs. 4 and 5 of Ref. [73]  is similar to the properties of the tracer diffusion on the lattice of moderately-sticky crowders examined here. For non-attracting obstacles corresponding to ǫ A = 0 the ratio D(φ)/D 0 of the diffusing coefficient of the tracer particle versus its value in an un-obstructed environment as function of the obstacle volume occupancy φ in a simple mean field approach is predicted to scale as [49] based on a volume exclusion argument in section 8.3 of Ref. [49]. Here x = 2R/a is the relative obstacle size with respect to the lattice spacing as used in simulations. The reader is also referred to Ref. [74] for the next-order corrections in the dependence of D(φ) on the volume fraction of crowders. Indeed, in absence of an attraction between the tracer and the obstacle surfaces, the behaviour of D(φ)/D 0 as function of R shown in Fig. 4 is in quite good agreement with the prediction (10). For moderate attraction, ǫ A = 2k B T , the decrease of D(φ)/D 0 with R becomes less pronounced at R values larger than R ≈ 0.6. For even stronger tracer-obstacle interaction, ǫ A ≥ 4k B T , remarkably the long-time diffusivity D(φ) becomes non-monotonic with the crowder size R, as shown in Fig. 4. Physically, at higher obstacle volume fraction φ the available space for tracer diffusion becomes effectively reduced from the three dimensional volume to a lower dimensional space. This creates pathways or channels between the "cages" created by the obstacle and can effectively speed up the exploration of space by the tracer particle at higher φ fractions. Note that this effect would be modified when the surface diffusivity were considerably smaller than the volume diffusivity. However, as shown in the discussion of the trapping times below, another contribution to this speed up-effect could be that for larger crowder radius R the tracer is in a limbo between vicinal attractive surfaces and thus manoeuvres between obstacles without binding to them.
We note that in Ref. [65] the diffusion coefficient of a tracer on an array of cylindrical obstacles on a static, two dimensional square lattice was analysed in terms of the generalised Fick-Jacobs equation and by Brownian dynamics simulations. For the relative diffusivity as function of the crowding fraction φ the analogous behaviour D(φ)/D 0 = 1 − φ = 1 − π(R/R max ) 2 /4 was found, where π/4 is the maximal surface coverage in this situation [65].

B. Statistics of tracer trapping times
The non-monotonic dependence of the long-time tracer particle diffusivity on the crowding fraction is also manifested in the non-monotonic variation of the times that the tracer spends in the obstacle-adsorbed state. In Fig. 5 we present the statistics of individual adsorption times to the crowding particles along a very long trajectory containing many binding-unbinding events. In the main graph of Fig. 5 we observe the distribution of tracerobstacle adsorption times features a peak at short t ads , while the tails of the histograms indicates the expected exponential decay. The corresponding mean adsorption time t ads is the average over all the binding events encountered in our simulations for a particular set of the model parameters. As we can see, for larger crowders these exponential tails become progressively longer with increasing obstacle size up to some R < 0.6×a/2, that is, the duration of binding events can become significantly longer. In the inset of Fig. 5 we use a logarithmic abscissa to pronounce the initial peak of the histograms. From this plot it becomes clear that several extremely long binding events can shift the mean t ads of the binding time significantly with respect to the most likely value, compare the relative positioning of the maximum of the histograms and the mean values indicated by the coloured dots.
Thus, a small number of extremely long binding events govern the corresponding mean adsorption time t ads . Here, we mention the related study of Ref. [75] in which the modes of surface versus bulk diffusion of a tracer in spherical domains were investigated. Also note that a different tracer diffusivity in the obstacle-bound mode as compared to the bulk diffusion can give rise to new interesting effects. In particular, an optimisation of the overall passage times of a tracer in the target-search problems on a lattice of trapping sites (obstacles) with a likely slower diffusivity should be analysed in the future.
Let us now turn to the total time of adsorption t A experience by the tracer particle during a trajectory of duration T generated. We thus sum up all the adsorption times experienced by the tracer, Adding to this quantity the excursion times in the bulk between the obstacles, we have T = t A + t B (the subscripts A and B denote the quantities in the adsorbed and bulk state, respectively). As expected, the total adsorption time t A obtained from the simulations is an increasing function of the obstacle size R, compare Fig. 6. The initial increase of t A with the crowder size R can be understood in terms of the concept of a "phantom sphere'. Namely, the surface of a crowding sphere available for surface diffusion by the crowding particles is S(R) = 4πR 2 . If the remaining volume a 3 − 4πR 3 /3 in a cubic unit cell framed by eight crowding spheres were converted into a ("phantom") sphere, the surface of that sphere would amount to Following this crude argument to simply divide up the time spent on the surface of the crowders and on the phantom sphere surface corresponding to the free volume of a unit cell in the obstacle lattice, the ratio of the adsorption time to the total time T of a sufficiently long trajectory can then be written as where we again used the dimensionless variable x = R/R max = 2R/a. The maximal volume of the phantom sphere is given by the volume a 3 of the unit cell, corresponding to an effective radius of the phantom sphere of R max,ph = (a/2)(6/π) 1/3 . The minimal effective radius of the phantom sphere corresponds to obstacles, which just touch each other on the square lattice, namely R min,ph = (a/2)(6/π − 1) 1/3 . The ratio t A /T given by Eq. (14) is quadratic in the obstacle size R for small volume occupancy by obstacles. This result is in agreement with our simulations results for moderate tracer-obstacle binding, as shown in Fig. 6. In fact, for the attraction strength ǫ A = 2k B T the agreement with result (14) is quite good given the simplicity of our model. In the shown double-logarithmic scale of Fig. 6 over the range R/R max = 0.1 . . . 0.9 the law (14) in fact only weakly deviates from the quadratic scaling t A /(t A + t B ) ∼ (π/6) 2/3 x 2 . Naturally, the data for increasing binding affinity progressively deviate from the formula (14), for the strongest binding strength shown the particles move almost exclusively on the obstacle surface, independently of the obstacle size.
We observe that for small obstacle sizes and weak tracer-obstacle attraction strengths the total adsorption time t A at fixed R grows exponentially with ǫ A , t A ≃ exp (|ǫ A |/[k B T ]), corresponding to the Boltzmann activation in an equilibrium system. As demonstrated in Fig. 6 on the right, for strong tracer-crowder attraction a saturation in t A (ǫ A ) is reached, such that the activation curve for the ratio t A /(t A + t B ) is analogous to the expression for a simple two level system, For larger crowders, when the tracers are confined predominantly to the obstacle surface, the saturation effect at larger attraction strengths is more pronounced, see Fig. 6.
In agreement with the results of Fig. 3 for the MSD, at moderate tracer-obstacle binding the adsorption time progressively increases for more voluminous obstacles on the lattice. In contrast, at strong tracer-crowder adsorption the adsorption time initially increases with the obstacle size R = a[3φ/(4π)] 1/3 (15) but above a critical volume fraction φ the value of t A decreases, as seen in Fig. 6 for the values ǫ A = 8k B T and 10k B T , as well as for the mean values shown in Fig. 5. This decrease of the mean adsorption time and thus a stronger contribution of bulk excursions is consistent with the non-monotonic dependence of D(φ)/D 0 at large tracer-obstacle binding strengths ǫ A and with the enhanced diffusivity for larger strongly adhesive obstacles, see the curve for ǫ A = 8k B T in Fig. 4. We ascribe this small yet somewhat counterintuitive effect to a competition of binding to neighbouring surfaces due to which the tracer particle is in a limbo in the bulk, possibly in conjunction with the reduced effective dimensionality of the environment mentioned above.

C. Time averaged MSD and non-Gaussianity parameter
To obtain more insight into the characteristics of the diffusive motion of the tracer particle in the crowded environment, we compute the time averaged MSD [15][16][17][18][19] from the trajectory r(t) of the tracer of range t = 0, . . . , T . In Eq. (16) ∆ is the so-called lag time, which defines the size of a window slid along the trajectory r(t).
The length of the trajectory T is also referred to as measurement time. In addition to the individual time traces δ 2 (∆) we also consider the average over N individual trajectories. When measured time series are not long enough, this quantity provides smoother curves if sufficiently many trajectories N are available.
In Fig. 7 we show the time averaged MSD (16) along with the MSD r 2 (t) of N = 100 individual trajectories. Let us first focus on the case of a moderate tracerobstacle attraction strength, ǫ A = 2k B T . We observe that the individual time traces δ 2 (∆) initially grow ballistically and then cross over to normal diffusion, as indicated by the slopes. The results for individual time traces δ 2 (∆) show only minute amplitude variations at different lag times, quantitatively similar to the spread of time traces of regular Brownian particles. Of course, when ∆ approaches the measurement time T , the statistics of the time average defining δ 2 (∆) worsen and some amplitude scatter occurs. The average (17) almost perfectly coincides with the ensemble average MSD r 2 (t) , compare the blue and green curves in Fig. 7. The latter observation corroborates the ergodic nature of the tracer motion, that is, the equivalence of ensemble and long time average of physical observables, here r 2 (∆) = δ 2 (∆) [15][16][17]. The linear long time scaling of the MSD defines the effective diffusion coefficient D(φ) we shown in Fig. 4. Note that prolonged adsorption periods of the tracer on obstacle spheres correspond to effective trapping and delays the growth of either x 2 (∆) or δ 2 (∆). For more information on the violation of the equivalence between time and ensemble averaged physical observables in anomalous-diffusive stochastic processes we refer to the recent review in Ref. [17].
As shown for the binding time statistics above, the trapping times are exponentially distributed and thus the long time motion converges to regular Brownian motion with a reduced, effective diffusivity D(φ). This scenario is therefore fundamentally different from subdiffusive CTRWs [14,37], in which the characteristic trapping time diverges. Our retarded Brownian motion-based physical rationale is consistent with experimental observations of protein diffusion in dense dextran solutions and with Monte-Carlo simulations of tracer diffusion on lattices of immobile inert obstacles as reported in Ref. [56]. In addition, we checked that the average time averaged MSD features no dependence on the trace length T : the values of δ 2 (∆) almost perfectly overlap for different trace-lengths T , as demonstrated in Fig. 8.
A different situation is encountered when we consider strong tracer-obstacle attraction, ǫ A = 6k B T in Fig. 7. We immediately observe that up to t = ∆ ≈ 10 2 the time averaged MSD δ 2 (∆) significantly differs from the corresponding ensemble average r 2 (t) . That is, on these time scales the systems exhibits the disparity r 2 (∆) = δ 2 (∆). For times exceeding t = ∆ ≈ 10 2 the agreement between both quantities becomes excellent, the system is asymptotically ergodic. Individual curves δ 2 (∆) for single time traces show a somewhat increased spread around their mean δ 2 (∆) , however, this is still within the range expected for (asymptotically) ergodic processes [76], and is significantly different from weakly non-ergodic processes such as subdiffusive CTRWs [15-17, 19, 77] or heterogeneous diffusion processes [17,42]. For diffusion processes of the subdiffusive CTRW type the spread of individual time averaged MSD traces is expected to be finite even for vanishing lag times [16,17]. This kind of behaviour is definitely not observed in our simulations. Here a very narrow spread of δ 2 in the whole range of tracer-obstacle affinities and obstacle sizes is observed, see, e.g., Fig. 7. The transient non-ergodic features observed here imply that the relaxation time towards ergodic behaviour is increased for longer trapping times when the tracer-obstacle attraction is more pronounced. The fact that the disparity r 2 (∆) = δ 2 (∆) is most pronounced around the turnover from initial ballistic to terminal Brownian motion is consistent with observations of confined stochastic processes driven by correlated Gaussian noise [17,79]. We now address a quantity that is based on the fourth order moment of the time trace δ 2 (∆). This non-Gaussianity parameter G(∆) was shown to be a sensitive experimental indicator of the type of effective stochastic process driving the tracer particle in crowded complex fluids [57]. The non-Gaussianity parameter in three dimensions is defined by [18] G(∆) = 3 5 For diffusion processes with a stationary Gaussian distribution of increments such as Brownian motion and fractional Brownian motion we have G = 0, while the parameter G becomes non-zero for processes with nonstationary increments and/or non-Gaussian distributions such as subdiffusive CTRWs or heterogeneous diffusion processes [44,57].
We observe that for the simulated diffusion process in the presence of attractive obstacles the non-Gaussianity parameter shown in Fig. 9 is close to zero for almost the entire length of the traces, apart from the initial regime of the motion including inertial effects. These short time deviations from G ≈ 0 are particularly pronounced for relatively large obstacles with strong tracer-obstacle attraction, as shown for the different parameters analysed in Fig. 9. The fact that G ≈ 0 together with the equivalence of the ensemble and time averaged MSDs at sufficiently long times are, of course, a mirror for the Brownian nature of the observed motion. The smaller non-zero values of G detected in the long-time limit in Fig. 9 are due to small discrepancies between the ensemble and time averaged MSDs.
At short to moderate times the non-Gaussianity parameter substantially deviates from the zero value characteristic of Brownian motion for strongly adhesive and relatively large crowders, see, e.g., the blue curve in Fig. 9. The deviations of G(∆) occur at time scales of t ∼ 1 . . . 30 when the tracer diffusion is inherently non-ergodic, as we see from the right panel of Fig. 7 plotted for the same parameters (ǫ A = 6k B T and R/R max = 0.6). At very short times, t ≪ 1, on which the MSD and time averaged MSD coincide, the non-Gaussianity parameter respectively assumes very small values. Otherwise, the deviations from G = 0 we observe are within the range typically measured in experiments [57] for asymptotically ergodic stationary-increment processes.

IV. CONCLUSIONS
We studied the passive motion of a tracer particle in an ordered, stationary array of attractive spherical crowders. Based on a truncated attractive Lennard Jones interaction between tracer and crowding obstacles, we simulate long individual trajectories of the tracer based on the Langevin equation. The resulting motion features a transition from initial ballistic flights corresponding to the inertial particle motion to a long time Brownian diffusion behaviour. The magnitude of the effective diffusion coefficient of this terminal Brownian motion is reduced for increasing obstacle density and tracer-obstacle attraction. However, a distinct non-monotonicity in this behaviour for large obstacle radii and high attraction strength is observed, likely due to competing attraction from multiple obstacles for which the tracer particle is in a limbo in the bulk. The same effect also leads to a decrease of the time t A (φ) spent adsorbed to the obstacle surfaces during a long trajectory. The long time Brownian dynamics was consistently shown to be associated with approximately vanishing non-Gaussianity parameter. These dynamic features point out the crucial role of varying tracerobstacle binding strengths in the analysis of crowded systems, as performed here. To put these qualitative statements on a more physical foundation, additional simulations will be necessary. In particular, we will study the time averaged van Hove cross-correlation functions [78] for the tracer motion.
At intermediate times the tracer particle motion is anomalous, with a distinctly time-dependent scaling exponent β(t). In this regime the trapping to the obstacles becomes the dominant mechanism. According to our results the time and ensemble averaged MSDs are equivalent for moderate tracer-obstacle attraction strength, however, a transient non-ergodic disparity between the two observables is observed over a range of some 2.5 orders of magnitude for stronger tracer-crowder binding. This transient form of weakly non-ergodic behaviour was previously observed for correlated Gaussian processes under confinement [34,79] of the otherwise ergodic process [39]. This behaviour should be kept in mind when precise diffusive properties are to be analysed from measured or simulated time traces of our system.
Extensions of the current model should consider a range of surface diffusivities of tracer particles bound to obstacles. Moreover, the arrangement of crowders should release the static, ordered arrangement on a lattice. Thus, off-lattice simulations could include the motion of crowders around their equilibrium positions, similarly to the analysis in Ref. [53]. Differences in the sizes of individual crowders and a certain randomness in the tracer-obstacle affinity would be additional relevant generalisations of the current system. Mobile obstacles were shown to profoundly shrink the time range of the transient subdiffusive motion compared to immobile crowders [53]. However, the generality of these findings needs to be explored in a broader parameter range. Currently, it remains elusive to arrive at realistic models capturing the richness of real MMC in living biological cells with their wide variety of crowder shapes, surface properties, persistence length and degree of branching as well as a poly-disperse size distribution, in addition to cellular structural elements, as well as charge effects [5]. Finally, active processes such as energy-consuming transport in living cells [35,80,81] needs to be added to achieve a closed picture of all facets of cellular dynamics.
Our study complements several other recent analyses of tracer motion in crowded environments. Thus, tracer diffusion in a system of relaxed and stretched polymer chains in the presence of tracer-polymer attraction was studied by Langevin dynamics simulations with the Espresso package [82]. We observe that such obstructed diffusion with sticky obstacles resembles our current results. For instance, the evolution of the time-dependent MSD scaling exponent β (see Fig. 4 in Ref. [82]) shows a transition from the initial ballistic regime to a subdiffusive regime at intermediate times, and further to Brownian motion in the long-time limit. The anomalous diffusion regime was shown to span a larger time window for the relaxed as compared to the stretched self-avoiding chains [82]. Higher polymer densities and stronger tracerpolymer interactions yield wider regions of tracer subdiffusion [82]. A continuation of this study for tracer diffusion in a cylindrical pore with surface-grafted polymer brushes of varying density showed that the subdiffusive regime is more pronounced for weak-to-moderate tracerpolymer interaction [66].
We also mention an experimental study of impeded colloidal diffusion in transient polymer networks with varying colloid-polymer binding interactions [83]. For the diffusion of an inert tracer in a responsive elastic network system, when the tracer size is of the same order as the unit cell of this gel, transient subdiffusion was reported and shown to involve characteristic collective dynamics of tracer and gel [78]. The transient subdiffusion in this study is in contrast to the experimentally observed longtailed distribution of trapping times of sub-micron tracers in semi-flexible, inherently dynamic networks of crosslinked actin [84] and thus further underlines the nonuniversal character of the dynamics in crowded systems. We finally mention that diffusion in two-dimensional, oriented fibrous networks in the presence of repulsive and attractive particle-obstacle interactions was in fact studied experimentally in connection with hydrogel-like structures of the extracellular matrix [85].
Is crowding in cells merely an effect of cramming a rich multitude of different bio-molecules into a minimal volume, or does it have an evolutionary purpose in giving rise to dynamic phenomenon such as (transi-ent) subdiffusion [21,86]? The combination of advanced single particle technology and other experimental methods along with improved in silico studies will lead to significant advances in the understanding of these still elusive questions.