Self-subdiffusion in solutions of star-shaped crowders: non-monotonic effects of inter-particle interactions

We examine by extensive computer simulations the self-diffusion of anisotropic star like particles in crowded two-dimensional solutions. We investigate the implications of the area coverage fraction $\phi$ of the crowders and the crowder-crowder adhesion properties on the regime of transient anomalous diffusion. We systematically compute the mean squared displacement (MSD) of the particles, their time averaged MSD, as well as the effective diffusion coefficient. The diffusion appears ergodic in the limit of long traces, such that the time averaged MSD converges towards the ensemble averaged MSD and features a small residual amplitude spread of the time averaged MSD from individual trajectories. At intermediate time scales we quantify the anomalous diffusion in the system. Also, we show that the translational---but not rotational---diffusivity of the particles $D$ is a non-monotonic function of the attraction strength between them. Both diffusion coefficients decrease as $D(\phi)\sim (1-\phi/\phi^*)^2$ with the area fraction $\phi$ occupied by the crowders. Our results might be applicable to rationalising the experimental observations of non-Brownian diffusion for a number of standard macromolecular crowders used in vitro to mimic the cytoplasmic conditions of living cells.


Introduction
Over the recent years, deviations from the standard Brownian diffusion law [1] have been observed in a broad range of systems [2,3,4,5,6,7,8]. Depending on the physics of the system under consideration, various theoretical models are used to describe these deviations [2,3,4,5,6,7,8]. Such anomalous diffusion is typically characterised by the power-law growth of the mean squared displacement (MSD) of particles with time We distinguish subdiffusion for 0 < β < 1 and superdiffusion for 1 < β. Subdiffusion is an abundant phenomenon for passive motion in the world of live biological cells [4,5,6,7,8]. In the biological context subdiffusion was observed for particles ranging from small proteins [9] via messenger RNA molecules [10] in the cell cytoplasm to large chromosomal loci and telomeres in the nucleus [11] to sub-micron virus particles [12] as well as lipid granules [13]. The features of anomalous diffusion depend on the energy landscape and the physico-chemical interactions in the system of particles [14,15]. The advances of modern single particle tracking experiments [10,16,17] provide a wealth of high resolution experimental data to quantitatively compare the microscopic mechanisms of non-Brownian diffusion with known theoretical models. The latter include, inter alia, the continuous time random walk [18,19,20,21,22] or the equivalent formulation in terms of fractional diffusion equations [3,23], fractional Brownian motion [24], heterogeneous diffusion processes [25], scaled Brownian motion [26,27,28], as well as the fractional Langevin equation related to the viscoelasticity of the environment [29,30]. The cytoplasm of biological cells is a superdense [10] fluid consisting of proteins, nucleic acids, membranous structures, cellular machinery components, semiflexible filaments, etc. [31,32,33,34]. This macromolecular crowding (MMC) reaches volume occupancies of φ 30% [35]. In addition, the cytoskeletal meshwork [36] of eukaryotic cells impedes the diffusion of larger entities in cells, in particular, near the cell's plasma membrane. The cytoplasm in addition is highly heterogeneous both in prokaryotic and eukaryotic cells [37,38,39]. The anomalous diffusion of cell-related phenomena may represent a blend of more than one theoretical model representing the quality of the diffusion on different length and timescales [4,5,7,40,41,42,43,44,45].
A number of experimental [37,46], theoretical [47], and simulation [35,48,49,50,51,52,53,54,55,56,57] studies in recent years were devoted to tackling various aspects of particle diffusion in crowded environments. From the simulation perspective, for instance, the studies of tracer diffusion in non-inert [58], heterogeneously distributed and poly-disperse [59], restrictively mobile [48], squishy [47] and anisotropic [60,61] obstacles were performed. Despite the progress of analytical theories of crowded solutions some important diffusive characteristics can only be studied quantitatively by computer simulations. This is particularly true for crowders of the non-trivial of the Mercedes-Benz star like particles considered in the current paper (Fig. 1).
We here use computer simulations to unravel the implications of the particle shape and squishiness as well as the crowding fraction on the translational (D) and rotational (D r ) particle diffusivities in highly crowded solutions. Our main target is to gain insight into the physical behaviour of non-spherical crowders relevant for the situation in vitro where soft non-spherical and often non-inert crowders such as globular PEG and branched dextran polymers are routinely used to mimic the effects of MMC in living cells. Another important experimental example is the diffusivity of anisotropic lysozyme-like proteins studied by Brownian Dynamics simulations in crowded media [62]. It was demonstrated that-particularly in heavily crowded solutions-not only a transient subdiffusion of the protein centre of mass exists, but diffusion becomes also progressively anisotropic. This anisotropy of the translational diffusion pronounced on short-to-intermediate times disappears in the long time limit. The long time diffusivity values were shown to drop drastically with the protein concentration [62]. Moreover, the reduction of D r for Y-shaped proteins such as IgG γ-Globulin (molecular weight MW≈ 155 kDa) was shown to be stronger than for more spherical proteins such as Bovine serum albumin (MW≈66 kDa). These experimental observations based on fluorescence correlation spectroscopy measurements are supported by all atom Brownian Dynamics simulations [62]. The inclusion of hydrodynamic interactions revealed an additional reduction of D r of proteins [63]. The reader is also referred to the simulation study of Ref. [64] in which the self diffusion of star like polymers in the presence of hydrodynamic interactions [65] was examined in detail. The paper is organised as follows. In section 2 we introduce our simulation model, the physical observable we are interested in, and some details on the data analysis algorithms. We present the main findings of simulations in section 3. In section 4 the implications of our results for some cellular systems are discussed.

Simulation model and observables
We implement our computer code developed to simulate the particle diffusion of crowded solutions in which all particles are explicitly treated [52,53,54]. Here, we consider a two-dimensional system of Mercedes-Benz star shaped crowders, each consisting of four discs of diameter σ connected by elastic springs, see Fig. 1A. The elastic potential between the midpoint of the molecule and the centres of the outer monomers is where r c is the equilibrium distance and k s the spring constant. We also connect the outer monomers with springs of the force constant k s , namely, to mimic the softness of our triangular star like crowders. The equilibrium distances and constants are set to r c = 1.5σ, r o = 1.5 √ 3σ, and k s = 100k B T /σ 2 . The interaction between all beads is described by the 6-12 Lennard-Jones potential Here Θ(x) is Heaviside step function and C(r cut ) is a constant setting U LJ (r > r cut ) = 0. For a purely repulsive potential the standard cutoff distance r cut = 2 1/6 σ is used with the potential strength ǫ = k B T . For attractive interactions we set r cut = 2σ with varying adhesion strength ǫ = ǫ A between the monomers. This attraction acts between all the monomers of the stars. We use periodic boundary conditions within a square box of area L 2 . The packing fraction of N crowders in the system is defined as φ = NA/L 2 , where A = 4π(σ/2) 2 is the total area of the four monomers and N ∼ 10 2 is a typical number of stars used in our simulations. In most scenarios below the system size is L = 40σ and the total simulated trace length is ∼ 4 × 10 8 of elementary time steps. The dynamics of the two-dimensional position r i (t) of the ith monomer disc interacting with the other monomer discs is described by the Langevin equation Here ξ(t) represents Gaussian white noise with zero mean ξ(t) = 0 and correlator where k B is the Boltzmann constant, γ is the friction coefficient, and T the absolute temperature. In the following, we use σ and k B T as the basic units of length and energy, respectively. We simulate the system with the Verlet velocity algorithm with elementary time step ∆t = 0.005 for the total time T . The physical time scale in these simulations is the standard combination δτ = σ m/(k B T ) ≈ 0.36 ns [66], if we set the monomer diameter to σ = 6 nm and its mass to the average mass of cytoplasmatic crowders, MW≈ 68 kDa [51,67].
We track the positions of the centre monomers of all the crowder stars and their orientation with respect to the x-axis, denoted as θ i . From the trajectory of the ith crowder we calculate the time averaged translational (δ 2 i ) and rotational (δ 2 r,i ) MSDs as [7] and Here ∆ is the lag time along the trace. In addition to the individual MSDs δ 2 (∆) we compute the corresponding averages over the set of N individual trajectories, as well as their amplitude spread around this mean value.
The diffusion is called ergodic if the ensemble and time averaged MSDs coincide in the limit ∆/T → 0 and if the spread of δ 2 around the mean approaches the delta function in this limit [7]. A more accurate description of ergodicity can be achieved based on the so-called ergodicity breaking parameter EB. The latter is defined as the variance of the distribution of the dimensionless variable ξ = δ 2 / δ 2 , whose precise behaviour as a function of the lag time and the various model parameters, however, is beyond the computational scope of the current investigation.

Results
In Fig. 2 we present the behaviour of the translational and rotational MSDs for varying inter-particle attraction strength ǫ A at crowder packing fraction φ = 0.15. The initial crowder diffusion is ballistic, stemming from the simulation of inertial particles, see also Ref. [58]. At intermediate time scales of ∆ ∼ 0.1 . . . 10 we observe a non-monotonic behaviour of the time averaged MSD that we ascribe to the events of the first collision of a given crowder molecule with another crowder. We quantify the variation of the local scaling exponent in Fig. 6B below. In the long lag time limit the translational and rotational MSDs grow linearly with ∆ reflecting the Brownian behaviour of the crowder particles. In this limit the diffusion is ergodic, as we demonstrate in Fig. 3. This statement is not necessarily trivial: in many weakly non-ergodic systems the time averaged MSD turns out to be a linear function of the lag time ∆ while the ensemble averaged MSD scales as a power-law or logarithmically in time t. This phenomenon was observed in various experiments [13,42,43,44] and explained in terms of various stochastic processes [7,25,26,68,69,70,71,72,73]. Fig. 3 demonstrates both that to very good approximation ergodicity in the sense of the equality x 2 (∆) = δ 2 (∆) and that a very small amplitude scatter around the mean δ 2 (∆) exists and thus the time averages are reproducible quantities. We furthermore detail the dependence of the particle diffusivity in this Brownian limit versus the attraction strength and the filling fraction in Figs. 4 and 7, respectively. Let us be more specific. Fig. 3 illustrates the time averaged MSD for different lengths of the time series of the star diffusion as well as the superimposed ensemble averaged MSD shown as the bold black line. As can be seen from the figure, the amplitude scatter of single traces δ 2 around their mean remains small along the entire trajectory except when ∆ ∼ T , as expected. This growing spread as ∆ ∼ T is a standard feature of even canonical Brownian motion appearing due to progressively poorer statistics when taking the time average [7]. More importantly we observe that the amplitude spread of time averaged MSD at a fixed lag time ∆ decreases as the length T of the time traces increases. This property is ubiquitous for ergodic diffusion processes [7]. We note that the magnitude of the amplitude scatter that we observe for δ 2 for moderately adhering Mercedes-Benz stars are similar to that of a tracer in a network of sticky spherical obstacles, compare Fig. 3 above and Fig. 7 as well as the explanations in the text of Ref. [58]. Computing the magnitude of the mean time averaged MSD for varying trace length T we observe that its magnitude stays nearly unchanged with T (Fig. 3). In the short lag time regime the ballistic scaling is visible.   Given these observations our process is ergodic and thus fundamentally different from other anomalously diffusive systems such as those described by continuous time random walks or heterogeneous diffusion processes [25]. For the latter a pronounced scatter of the time averaged MSD trajectories around their mean and a clear dependence of the amplitude of δ 2 (∆) on T at fixed ∆ exist, that is, the system ages [7,74,75].
The particle diffusivities are defined as and obtained in the long time limit ∆ ≫ 1. Fig. 4 shows the values of D and D r obtained from a linear fit of the translational and rotational time averaged MSDs in the range ∆ = 10 3 . . . 10 4 . We find that while the rotational diffusivity D r decreases monotonically, the translational diffusivity remarkably exhibits a shallow yet significant maximum at ǫ * A ≈ 1k B T . This systematic trend persists for the variation of the crowder fraction in a quite broad range (Fig. 4). This implies that the self-diffusion of our star like crowders can be facilitated by a weak inter-particle attraction. This is one of the main conclusions of this study. One can rationalise this trend in the self-diffusion in terms of the concept of the effective crowder size that decreases for moderate attraction strengths ǫ A ≈ 1k B T . Fig. 4 illustrates that for progressively stronger star-star attraction their mutual diffusivity decreases eventually to zero due to aggregate formation, see also Fig. 5 and its discussion below.
We also detect a progressive aggregation of crowders at relatively large values of the crowder-crowder attraction strength ǫ A , as demonstrated in Fig. 1C and in the video files in the Supplementary Material. This is a well known phenomenon, for instance, in the glass transitions of dense suspensions of sticky hard spheres [76]. Accordingly, the decrease of the average diffusivity D of crowders as plotted in Fig. 4 emerges due to the averaging over an ensemble of particles that perform individual random motions. This average takes into account both particles forming transient aggregates as well as free particles. Roughly speaking the average diffusivity drops inversely proportionally to the number of particles in the cluster. The fraction of particles clustering in these aggregates increases with the mutual attraction strength. The average diffusivity therefore progressively decreases with ǫ A due to a larger fraction of particles in the transient aggregates. At large attraction strength the majority of particles belong to big clusters.
At a fixed cohesiveness ǫ A of our Mercedes-Benz stars vicinal crowders create a rough energy landscape for the self-diffusion and the hopping of a given crowder particle. As the MMC fraction φ increases, the binding events give rise to more a prolonged particle aggregation and reduced self-diffusivity. Above a critical MMC fraction the barrier height exceeds the thermal energy unit k B T thus increasing the lifetime of crowder aggregates significantly. For stronger star-star attraction the formation of essentially permanent aggregates sets in for less crowded systems leading to an inhomogeneous, phase-separated spatial distribution, see Fig. 1C.
A similar non-monotonicity of the translational diffusivity D at similar strengths of the particle-crowder attraction was found in Ref. [51] for the tracer diffusion in dense suspensions of spherical Brownian particles. While we here detect that the attraction strength yielding the highest value of D is a function of the crowding fraction φ of the stars for the spherical particles, the stickiness facilitating the particle diffusivity was almost φ-independent in Ref. [51]. The non-monotonic D(ǫ A ) dependence was interpreted in Ref. [51] in terms of the roughness of the free energy landscape for the tracer diffusion using the concept of the chemical potential. Interestingly, the tracer diffusivity was also non-monotonic in φ in a static regular array of sticky obstacles, as quantified in Ref. [58].
We checked the universality of the observed dependencies for D(ǫ A ) and D(φ) also for spherical particles. Namely, we simulated just a single monomer of our Mercedes-Benz stars with the given adhesive properties. The diffusivity was indeed found to reveal a maximum at ǫ opt A ∼ 0.5 . . . 1k B T (not shown), indicating some universality of this a priori counter-intuitive faster diffusion for a weak inter-particle attraction [51]. Note also that for a polymer chain diffusing in an array of sticky obstacles a weak chain-obstacle attraction can also substantially enhance the polymer diffusivity [51,77].
To rationalise the observed behaviour of D(ǫ A ) we calculate in Fig. 5 the potential of the mean force between two crowders as In this reconstructed free energy ρ(r) is the average radial distribution function of the centre monomer of the crowders in the steady state long time limit. As the mutual attraction strength ǫ A increases we observe that the potential well at the separation r ≈ σ becomes deeper, see the first well in Fig. 5. Concurrently, the distance at which F (r) sharply increases becomes shorter for larger ǫ A . For a stronger star-star attraction the crowders feature a more organised appearance, resulting in measurable oscillations of ρ(r) and F (r), as evidenced in Fig. 5. These trends indicate that the effective crowder radius gets smaller with increasing ǫ A , and at an optimal value ǫ opt A the crowders approach one another more closely yet without sticking. This in turn might result in a faster average diffusivity D at ǫ A ≈ ǫ opt A , as we indeed observe. An effective reduction of the crowder size at optimal attraction strength is one important causealbeit possibly not the only one-for this facilitated diffusion. In the current system, the equilibrium distance of the outer monomers from the cenral monomer is reduced by about 2% for the inter-monomer attraction strength of 2k B T . Even higher values of ǫ A give rise to the formation of large clusters of crowders, see the Supplementary Material. As shown in Fig. 4 the corresponding diffusivity of an average particle drops dramatically.
In Fig. 6A we show the translational and rotational MSD for varying packing fraction of crowders φ. As expected-from a linear fit to the long time time averaged MSD-the diffusivity is a monotonically decreasing function of φ, as evidenced by Fig. 7. For more crowded systems the tracer diffusion gets more obstructed and the magnitude of the corresponding mean time averaged MSD δ 2 decreases. To elucidate these effects further we evaluate from the time averaged MSD traces of Fig. 6A for the translational motion the local diffusion exponent [7,78] For the rotational motion the exponent β r (t) is defined analogously. We observe a ballistic regime with β = 2 in the particle diffusion at short times, (Fig. 6B). This ballistic regime is followed by a decrease and further increase of the scaling exponent at ∆ ∼ 1. These non-monotonic trends are also clearly visible from the behaviour of the time averaged MSD traces themselves as a function of the lag time ∆, see Fig. 6A. We find that the variations of the scaling exponent for translational and rotational motions of the star like crowders appear correlated, indicating a coupling of these diffusion modes [16]. In the plots for the scaling exponent β(t) in Fig. 6B the significant spike-like signal at ∆ ∼ 1 is interpreted as an effect of the first collision of particles and the resulting onset of an effective confinement. We note that even in effective one-particle theories pronounced oscillations occur at the crossover point between initial ballistic and the overdamped regime [79,80].
With increasing crowder fraction φ we also observe a more pronounced range of anomalous diffusion for lag times of the order of ∆ ∼ 1 . . . 100. This range appears strongly correlated between rotational and translational particle motion, as shown in   6B. For rotational diffusion the scaling exponent drops practically to zero for times longer than those of the initial ballistic growth, and the corresponding mean time averaged MSD trace δ 2 r exhibits a short plateau (Fig. 6B). In the long lag time limit the exponent becomes Brownian β ≈ 1. Such a transient subdiffusion was observed for a number of systems [5,6,7], see also the Introduction. Especially in dense colloidal systems close to the glass transition φ = φ * this subdiffusion is accompanied by an exponential growth of the solution viscosity η = η(φ) which is divergent at φ → φ * [81]. The colloidal glasses also exhibit progressive particle localisation effects as discussed in Refs. [81,82].
Remarkably, the relative variation of the translational and rotational diffusivities with the crowding fraction of stars is quite similar. For comparison, we plot in Fig. 7 the theoretical prediction for dense suspensions of hard spheres [83,84] with the critical packing fraction for our system of φ * (ǫ A = 1k B T ) ≈ 0.52. Above this value φ * both translational and rotational diffusivities of the crowders essentially vanish. At this critical crowding fraction the inter-particle attraction becomes so strong that the self-diffusion is almost completely localised and the motion of particles corresponds more to a very restricted wiggling and jiggling. As expected, when the star-star interactions become stronger aggregate formation sets in for less crowded systems, and thus the critical value φ * is diminished (not shown). Albeit this theory in Refs. [83,84] is developed for three dimensional suspensions in the presence of hydrodynamic interactions, it agrees remarkably well with our results, as shown in Fig. 7. The reader is also referred to Ref. [85] for experimental data of the crowding dependent diffusivity of colloidal particles and alternative theoretical predictions for the diffusivity D(φ). We note that Ref. [86] suggest exponential rather than power-law forms for the particle diffusivity in crowded solutions.

Discussion and conclusions
We performed extensive computer simulations and theoretical data analysis of the diffusion of crowders with a branched structure. A simple example of such spiky but responsive crowders in two dimensions are deformable Mercedes-Benz like stars employed here. Their outer monomers are inter-connected by an elastic potential bestowing upon it a certain degree of responsiveness-an important characteristics for many polymeric crowders [47]. We also incorporated in the simulations an inter-particle attraction strength which represents another realistic feature of solutions of non-ideal crowders in vitro.
We found that the diffusion of our Mercedes-Benz star like crowders is ergodic and, within accuracy, Brownian in the long time limit. We examined the behaviour of the ensemble averaged MSD and the time averaged MSDs of the crowders in a wide range of the crowder fraction φ and the inter-crowder attraction strength ǫ A . As a function of the crowding fraction we demonstrated that both translational D and rotational D r diffusivities follow the analytical decrease (13) of D(φ) predicted for suspensions of hard spheres. The dependence of the star-star attraction strength is more remarkable. Namely, the translational diffusivity shows a weak yet systematic nonmonotonic dependence on ǫ A for the solutions at all crowding fractions studied herein. The rotational diffusivity, in contrast, is a monotonically decreasing function of the inter-particle attraction strength ǫ A . Thus, a relatively weak inter-monomer attraction can facilitate the lateral diffusion and also induce a certain degree of clustering and spatial heterogeneities in crowded solutions of non-inert particles. These effects will impact the diffusion of a tracer particle in crowded solutions-such as those of PEG, dextran, or Ficoll-used in vitro to mimic the crowded conditions in living cells [47,49]. In addition to the proof of the ergodic long lag time diffusion shown in Fig. 3  Of course, our planar triangle-like stars still represent a quite primitive system to mimic the non-ideal shape of real crowders in experimentally relevant setups. Future investigations including a three dimensional pyramid-like shape of crowders with longer polymeric arms will further elucidate the physical consequences of non-spherical and squishy crowders, and potentially exhibit additional unexpected behaviour. Moreover, not only the self-diffusion is to be studied but also the diffusion of tracer particles of various sizes and shapes in such crowded suspensions [56] as well as poly-disperse mixtures of crowders should be investigated. Some asymmetry may also be incorporated in the crowder shape. Recently, single particle tracking measurements allowed one to rationalise the translational and rotational diffusivities of micron-size symmetric and asymmetric boomerang-shaped particles in two dimensions [87]. It was observed that the regimes of Brownian diffusion exist at short and long times while a coupling of D and D r gave rise to subdiffusion at intermediate times.