Collective dynamics effect transient subdiffusion of inert tracers in gel networks

Based on extensive Brownian dynamics simulations we study the thermally driven motion of a tracer bead in a cross-linked, dynamic gel network in the limit when the tracer bead's size is of the same size or even larger than the equilibrium mesh size of the gel. The analysis of long individual trajectories of the tracer bead demonstrates the existence of pronounced transient anomalous diffusion, accompanied by a drastic slow-down of the gel-bead relaxation dynamics. From the time averaged mean squared displacement and the van Hove cross-correlation function we elucidate the many-body origin of the non-Brownian tracer bead dynamics. Our results shed new light on the ongoing debate over the physical origin of sterical tracer interactions with structured environments.

Based on extensive Brownian dynamics simulations we study the thermally driven motion of a tracer bead in a cross-linked, dynamic gel network in the limit when the tracer bead's size is of the same size or even larger than the equilibrium mesh size of the gel. The analysis of long individual trajectories of the tracer bead demonstrates the existence of pronounced transient anomalous diffusion, accompanied by a drastic slow-down of the gel-bead relaxation dynamics. From the time averaged mean squared displacement and the van Hove cross-correlation function we elucidate the many-body origin of the non-Brownian tracer bead dynamics. Our results shed new light on the ongoing debate over the physical origin of sterical tracer interactions with structured environments. Modern single particle tracking technology [1] unveils the non-Brownian stochastic motion of submicron tracers or fluorescently labeled macromolecules in a variety of complex liquids. Anomalous diffusion of the form r 2 (t) ≃ t α with 0 < α < 1 [2] was observed in the cytoplasm and membrane [3][4][5] of living biological cells, as well as in dense polymer or protein solutions [6]. Causes for the observed subdiffusion may be the crowded state of the environment [7] or sticking effects between the tracer and the environment [8]. Another important origin for anomalous diffusion is the subtle interaction between the tracer particle and the dynamic confines of a structured, gel-like environment with a well defined mesh size [9].
The existence of a mesh-like, structured environment is a defining property for the transport in several systems. Thus, biological cells are internally equipped with a characteristic mechanical network consisting of actin and other biofilaments through which submicron particles diffuse or are actively transported [10]. Inside the eukaryotic nucleus, biomolecules diffuse through the complex chromatin network [11]. In cellular tissues the space between cells is filled with the mesh-like extracellular matrix [12]. Novel mobile clinical diagnosis tools are based on hydrogel films, into which pathogens such as viral particles need to diffuse [13]. Finally, particles in a biofilm move through a flexible, porous bacterial matrix [14].
How do particles diffuse through a flexible, thermally agitated mesh such as a hydrogel? The answer is trivial as long as the diffusing particles are either much smaller or much larger than the typical mesh size. In these cases they diffuse normally or are completely immobilized in the mesh [9]. The physically intriguing case is met when the size of the diffusing particles is comparable to the typical mesh size. This is the scenario we consider here (Fig. 1). From extensive Brownian dynamics simulations we elucidate the fundamental microand mesoscopic physical principles behind the (transient) anomalous tracer dynamics on the collective many-body level. The tracer's mean squared displacement (MSD) and the van Hove cross-correlation function demonstrate the occurrence of significant tracer subdiffusion and massive cooperative breathing effects of the mesh, allowing relatively large particles to move in the gel, albeit under massively reduced diffusive progress.
Tracer motion in gels. The motion of a tracer in a flexible mesh can be viewed as a convolution of the Brownian self-dynamics of the tracer and the thermal agitation of the interacting gel particles, important ingredients characteristic for the tracer-gel system being the sterical obstruction of the tracer by the mesh [15], tracer-gel interactions [16], and the elastic response of the network itself [17]. Remarkably, distinct anomalous tracer diffusion was observed in numerous experiments [9,16,[18][19][20] and simulations [21]. To explain this non-Brownian dynamics in the gel two scenarios are typically invoked: (i) In the first scenario trapping of the tracer occurs due to attractive tracer-gel interactions, that effect transient binding to the gel network [16,18,22,23], such that the tracer intermittently follows the gel chains' Rouse dynamics [18] or, alternatively, remains transiently immobilized for a given period [22]. (ii) In the second scenario sterical hindrance within a fractal organization of the accessible volume is suggested to cause the non-Brownian tracer dynamics [9,16,[19][20][21]24]. In these approaches the gel is represented by impenetrable beads (either static [21] or dynamic [24]), positioned on a lattice. For static networks simulations revealed that isolated, randomly positioned beads give rise to non-transient subdiffusion, whereas tracer diffusion in presence of locally clustered beads in the form of condensed cubes randomly distributed in space is only transiently anomalous [21]. The debate on the dominant mechanism in real gel systems remains open [16]. Notably, both mechanisms are only explained on a heuristic level and a better physical understanding is needed.
While lattice representations of gel structures offer a computationally convenient tool to study tracer diffusion therein [21,24], they suffer from the severe unphysical localization of the gel despite its inherent flexibility and exposure to a thermostat. Moreover, a static gel introduces spurious artifacts in the tracer dynamics, as quenched configurational disorder of obstacles itself introduces a strongly persistent memory in the motion of a Brownian particle at all obstacle densities [25]. As we will show here, the tracer motion is inherently related to the breathing of the gel, i.e., collective soft modes of local expansion and contraction of the thermalized mesh. Such breathing modes cannot be appropriately captured in the lattice models. We thus face a pressing need for more realistic continuum models with fully coupled, many-body tracer-gel dynamics to gain physical insight into the effective tracer motion in such systems.
Model. We consider a gel of 7 × 7 × 7 connected beads and a freely diffusing tracer bead ( Fig. 1) interacting with the repulsive part of the shifted Lennard-Jones potential where i and j for a specific pairwise interaction take on the values t for the tracer particle and g for gel beads. r ij ≡ |r i − r j | is the physical separation between i and j. σ ij are the separations at zero unshifted potential, i.e., σ ii is the diameter of species i. ǫ sets the respective energy scales, and Θ(x) is the step function. Nearest neighbor gel beads are connected by Morse springs with the potential with ζ = k e /(2ϕ d ). ϕ d is the potential depth, k e the force constant at the potential minimum, and d e ≡ r ij − r eq , where r eq denotes the equilibrium separation. r c is a cutoff distance and V s = V ij M (r c ). The particle positions are governed by the overdamped Langevin equations where D i is the diffusion coefficient of the tracer or gel bead, β −1 = k B T , and F ji = −∇ ri U ij (r ij ) is the force acting on particle i due to particle j, and equals the sum of the Lennard-Jones potential (1) and the Morse potential (2). ξ(t) is a delta-correlated Gaussian noise with zero mean and component-wise variance Since hydrodynamic interactions are screened in dense systems [26] we neglect them in our model. We use dimensionless units and express distances in units of the gel bead diameter σ gg , energies in units of β −1 and diffusivities in units of σ gg / √ ǫ gg , thus fixing the time unit t 1 = 1. We take ǫ gg = 1.2 corresponding to good solvent conditions for the corresponding linear chains of beads [27], and we set ǫ tg = ǫ gg . We choose r eq = 4, which also fixes the gel unit cell size L 0 = r eq . Finally, we consider various sizes σ tt of the tracer particle and spring stiffness k e . The free tracer diffusivity is given by the Stokes-Einstein relation, D ≃ 1/σ tt . We use combining rules such that σ tg = 1 2 (σ gg + σ tt ). The Langevin equations are solved with the Euler method with integration step 10 −4 under periodic boundary conditions. The initial configuration is chosen randomly to be slightly displaced from the energetically minimal configuration, not allowing significant overlap nor significant bond stretching. After extensive equilibration time of 10 2 , trajectories of length 10 4 are simulated and analyzed. In the following we drop the indices and use the remaining symbols for the tracer particle σ tt ≡ σ etc. The overall time to simulate and analyze a single trajectory of length T = 10 4 was 5-6 weeks of single-CPU time on a computer cluster, and our study thus reaches the current computational limit.
Analysis of tracer dynamics. We use the time averaged MSD at lag time ∆ over a trajectory of length T , to characterize the tracer dynamics. For a wide range of parameters we observe that δ 2 exhibits a distinct, transient subdiffusion between the initial free diffusion regime and the asymptotically normal diffusion regime at long lag times, see below. This behavior agrees with similar simulations [27].  The typical behavior of the MSD δ 2 is shown in Fig. 2. Since our model gel structure is comparatively open (r eq = 4), in agreement with intuition tracer particles of size comparable to the gel bead diameter almost do not feel the presence of the gel, and the transition from short to long time free diffusion is barely visible on the double-logarithmic scale. In contrast, as the size of the tracer grows, a pronounced subdiffusive power-law scaling δ 2 ≃ ∆ µ emerges for ∆ > 1, and the anomalous diffusion exponent µ delicately depends on the tracer size. We note that the displacement amplitude χ = |r(t+∆)−r(t)| on this time scale is several times larger than the unit cell of the gel, and we hence observe a significant tracer motion and not simply fluctuations around a localization. The latter will be shown explicitly below. After the transiently anomalous regime there follows the normally diffusive long time regime (not shown) [28]. These results demonstrate that the existence of spatial obstructions alone is indeed sufficient to effect transient anomalous diffusion, but they do not provide any details about the underlying dynamic mechanisms.
To gain insight beyond the ∆-scaling of the MSD δ 2 we compute the distribution of relative displacements of magnitude χ for given lag time ∆ along the trajectory, The results for this quantity for different tracer sizes are shown in Fig. 3. As expected from the behavior of δ 2 small tracer particles do not feel appreciable obstruction by the gel mesh, their relative displacement χ evolves smoothly as a single peak, which moves to progressively higher values while simultaneously broadening (Fig. 3a)). With growing tracer size, however, the major fraction of relative displacements represented by the dominant peak of P T spreads much more slowly and remains almost localized even for longer ∆, in particular, for σ = 3, as seen in Fig. 3b) and c). Moreover, P T slowly develops a second, larger displacement fraction at χ > 3.5 for σ = 2.5 and χ > 3 for σ = 3, respectively. We stress that the tracer is not localized in the same mesh cell on its trajectory, not even for the case σ = 3, but ventures into vicinal cells. We also point out that there exist an apparent time scale when most escape attempts across the mesh boundary to the next mesh cell become successful, corresponding to the spreading of the initial peak beyond the confines of a mesh cell. This time scale is of the order of ∆ ≈ 50 for tracers of size σ = 2.5 and significantly longer than ∆ ≈ 100 for σ = 3, suggesting that it is the very cell escape dynamics that bestows the subdiffusive nature on the tracer dynamics. A finite characteristic escape time restores the terminal normal diffusion.
Collective gel-tracer dynamics. Are the tracer passages from one mesh cell to the next simply governed by an averaged obstruction effect and thus decoupled from the breathing dynamics of the gel, or do they correspond to a much richer scenario involving many-body effects? We study this question via the time averaged van Hove crosscorrelation function (vHCF) is the position of the gel bead i at time t and N gel the overall number of gel beads. C −1 was chosen such that the probability density of observing a gel bead infinitely far away from the initial position of the tracer is unity at any lag time and corresponds to dividing the uncorrected G T (r, ∆) by the number density of beads, ρ = N gel /V . G T (r, 0) corresponds to the static pair cor-relation function. The vHCF effectively measures the memory loss of the environment about the instantaneous tracer location and thus reflects the correlations of the dynamics of the gel beads and the tracer particle.
The results for G T (r, ∆) for different tracer sizes are shown in Fig. 4. In part a) for the smallest tracer, gel beads quickly penetrate the tracer's initial position on the scale of ∆ ∼ 1, during which the tracer typically moved by χ ∼ 1, just about a bead diameter (compare Figs. 2 and 4a)). The sharper peak corresponding to nearest gel beads quickly decorrelates. On time-scales ∆ > 10 the small tracer typically moves a distance of (several) unit cells, and G T (0, ∆ 10) shows oscillations corresponding to vivid fluctuations of gel bead positions. On these and longer time scales the correlations are completely lost at distances beyond a very short cutoff, which is due to the fact that the gel motion is translationally confined by the elastic interactions. Hence, beyond a typical time scale of ∆ 10 the beads completely forget that the tracer particle was in their proximity.
The situation changes dramatically for larger beads. For σ = 2.5 we observe in Fig. 4b) that the position correlations are quite persistent, and even the position of second nearest neighbor beads is still appreciably correlated. We note that the tracer diameter easily fits a unit mesh cell. Even on a time scale of ∆ ∼ 10, on which the tracer has likely already escaped a unit mesh cell at least once (see Fig. 3b)), long range correlations persist and the gel bead positions have not decorrelated much. Looking at even longer time scales (∆ > 40), on which the tracer typically traverses several unit cells, the long range position correlations persist while the G T (0, t) region already shows oscillations. A similar but even more drastic picture is found in the case σ = 3, where, in contrast, the tracer typically remains inside a single unit cell while escaping it only occasionally.
From these results we are able to draw a full physical picture for the emergence of transient anomalous diffusion in flexible gels. For tracer particles of size comparable to the gel beads and for gel structures with mesh size of the order of several times the tracer size the observed diffusion is normal, with a a renormalized diffusion coefficient reflecting the average effect of tracer-gel collisions. While the tracer size grows the probability to collide with the gel increases strongly and the free diffusion coefficient is reduced, suggesting that on the time scale of the gel breathing modes the correlations in the tracer dynamics are not yet relaxed. The resulting coupling creates persistent long-range correlations on time scales on which the tracer already crosses several unit cells. Since the memory of the location of the tracer before an escape to an adjacent cell is not destroyed, the same holds true for the entrance to a new cell. That is, the beads of the cell into which the tracer is about to cross already experience effects of the vicinity of the tracer before the tracer actually enters. In turn, this demonstrates that the intermediate- Conclusion. We used extensive Brownian dynamics simulations to study the time dependent correlations for the motion of tracer particles in a flexible, thermally agitated gel. Going beyond the single trajectory analysis solely of the tracer dynamics we demonstrated that the distinct transient subdiffusion, that is frequently observed in experiments and simulations, is effected by the collective fluctuations of both the tracer and the vicinal gel beads. On the scale of the typical escape time of the tracer from a unit cell the dynamics for larger beads turn out to still be appreciably correlated. While memory effects are a common argument invoked in explaining anomalous diffusion phenomena, we here obtained concrete evidence identifying its physical origin. Namely, this memory arises from persistent many-body correlations in the entangled tracer-gel dynamics. These correlations for larger tracer particles clearly reach beyond the nearest neighbor gel beads. This new physical insight obtained from the van Hove cross correlation function and the displacement distribution sheds new light on the elusive problem of tracer-gel interaction in the experimentally important limit when the tracer size is comparable to the typical mesh size. We are confident that these results will inspire new experimental and theoretical in-vestigations of tracer motion in structured matrices.
Acknowledgment. AG acknowledges funding through an Alexander von Humboldt Fellowship. RM acknowledges funding from the Academy of Finland (FiDiPro scheme). The authors gratefully thank the National Institute of Chemistry, Slovenia for providing extensive access to the institutes computational facilities.