Diffusion in crowded colloids of particles cyclically changing their shapes

A simple model of an active colloid consisting of dumbbell-shaped particles that cyclically change their length without propelling themselves is proposed and analyzed. At nanoscales, it represents an idealization for bacterial cytoplasm or for a biomembrane with active protein inclusions. Our numerical simulations demonstrate that non-equilibrium conformational activity of particles can strongly affect diffusion and structural relaxation: while a passive colloid behaves as a glass, it gets progressively fluidized when the activity is turned on. Qualitatively, this agrees with experimental results on optical tracking of probe particles in bacterial and yeast cells where metabolism-induced fluidization of cytoplasm was observed.

turnover cycles are only weakly non-reciprocal and the resulting propulsion forces are too small [7,8]. Therefore, such systems constitute a special class of active colloids, where individual particles are repeatedly changing their shapes, but do not swim.
Even in absence of self-propulsion, persistent energydriven conformational changes in macromolecules create non-thermal fluctuating hydrodynamical flows around them. Because passive tracer particles are advected by such fluctuating flows, their diffusion can become enhanced [9][10][11]. However, the analysis has so far been limited to dilute systems and therefore its results are not directly applicable to crowded colloids.
At sufficiently large volume ratios, passive colloids are known to have glass-like properties, manifested by slow relaxation, subdiffusion, and non-ergodicity (see, e.g., review [12]). In vivo experiments [13] on optical tracking of particles inside bacterial or yeast cells have been performed -with a surprising result that, at least for the probe particles with relatively large sizes, such glass behavior is only characteristic under starvation conditions 40003-p1 or in absence of the chemical fuel, such as ATP. The metabolism fluidizes the cytoplasm, leading to the recovery of classical transport properties within it.
In this letter, we propose an idealized model of oscillatory active colloids where the activity level and the rate of energy supply can be gradually controlled. The individual particles forming the considered colloid are active dimers, or dumbbells [9,14]. It is known that such active dimers reproduce, in an approximate way, the mechanochemical activity in real enzymes and protein machines [15]. Moreover, they have been already employed [16] in large-scale hydrodynamic simulations of the colloids (but still under less crowded conditions where glass effects are not yet seen). In contrast to the previous publications [9,11,15,16], our dumbbells are however further simplified: we assume that the natural length of the elastic spring that connects two beads in a dumbbell is periodically varied with time. The dumbbells interact via a soft parabolic repulsive potential and stochastic Langevin dynamics is assumed. In the present version, hydrodynamic interactions between the particles are dropped.
By running numerical simulations for relatively small 2D systems, we demonstrate that diffusion of dumbbells gets strongly enhanced when shape oscillations of its constituting particles are turned on. Further statistical investigations reveal that, while the passive colloid behaves as a glass, classical diffusion properties become recovered as conformational activity is increased.
The model. -Enzymes are single-protein catalysts that convert substrate(s) into product(s) in each turnover cycle. In most enzymes, the cycles are accompanied by internal mechanochemical motions, i.e., by repeated changes in the conformation of a protein. Biological molecular motors and machines are also enzymes, with the only difference that mechanochemical motions are employed by them to manipulate other macromolecules or to produce mechanical work. To do this, chemical energy supplied with the substrate (often, ATP) is used. The bacterial cytoplasm or a cellular biomembrane are essentially densely packed colloids of active proteins that repeatedly change their shapes.
Typically, mechanochemically active enzymes have a domain structure and their ligand-induced internal motions consist of the changes in distances between the domains and in the mutual orientation of them. Principal properties can be already reproduced in a simple active dimer model of an enzyme (see review [15]).
In the model, an enzyme protein is viewed as consisting of two beads (domains) connected by an elastic spring. The natural length of this spring depends on the ligand state of the enzyme: it is longer in absence of a ligand, but gets shorter after the substrate-enzyme and productenzyme complexes are formed, returning to the original longer length after the product release. Thus, an enzyme behaves as a mechanochemical oscillator that undergoes an elongation and a contraction in each turnover cycle [15]. This model has been successfully used to investigate collective hydrodynamic effects in numerical simulations for large populations of catalytically active enzymes [16]. In our study, it will be however further simplified: we will not explicitly consider the ligand states and chemical transitions between them. Instead, it shall be just assumed that the natural length of the spring connecting two beads in a dimer is periodically (with some additional random drift) changing with time. Physically, the modulation period should be considered as corresponding to the turnover time of an enzyme.
Explicitly, we assume that the natural length i of dumbbell i oscillates with time as where the oscillation phase satisfies the equation Here, ω i is the mean oscillation frequency and ζ i (t) is the internal noise that takes into account stochastic variations in cycle times. This Gaussian noise is delta-correlated in time and independent for different dumbbells, The parameter η controls the characteristic coherence time for oscillations of the shape. In our numerical simulations, we assume that all dumbbells are identical and have the same oscillation frequency ω. The model can however be readily extended to allow for random variation of these parameters. Thus, the time-dependent elastic energy of dumbbell i is where k is the stiffness of the internal spring and r (i) 12 = |r (1) i −r (2) i | is the distance between the first and the second beads in the dumbbell i.
The considered colloid consists of N dumbbells located within a volume of linear length L (periodic boundary conditions are used). Beads belonging to different dumbbells interact via a soft repulsive potential where r is the half-distance between the beads. The parameter u 0 specifies the repulsion strength and 2R is the distance between the particles below which the repulsion starts. Note that there is no repulsion between the beads in the same dumbbell. In molecular dynamics (MD) simulations of colloids, either full Newton dynamics or reduced stochastic dynamics can be used. However, at long times, both descriptions become equivalent [17]. Below, stochastic Langevin dynamics will be employed. Furthermore, in the present version, possible hydrodynamic interactions between the particles 40003-p2 are dropped, so that the effects of direct collisions are more clearly seen. They can however be introduced in the future into the model either by explicitly including the water molecules or in the framework of the multiparticle collision dynamics approximation (see, e.g., [16]).
Under the stochastic Langevin dynamics, the velocity of the bead n = 1, 2 in dumbbell i is Here μ is the mobility of the beads and represents the potential experienced at position r where T is the temperature, k B is the Boltzmann constant, It is convenient to measure all lengths in units of the repulsion interaction radius R and the time in units of the relaxation time τ = (μk) −1 of the dumbbell. On the considered molecular scales, a convenient unit of energy is the thermal energy k B T .
After rescaling, the model is characterized by a set of dimensionless parameters: Through periodic modulation of the natural length of the spring connecting two beads, energy is persistently supplied to a dumbbell. To estimate its mean supply rate, we can notice that, in the steady state, it should be equal to the mean rate at which energy is dissipated by the dumbbell. A simple calculation yields that the energy ΔE supplied to an active dumbbell per one oscillation period is Additionally, an important parameter of the model is the fraction φ of the total volume occupied by the dumbbells. Figure 1 shows an example of the simulated active colloids. Three supplementary videos of active colloids, based on the simulations, are available in (see S1.mpg, S2.mpg and S3.mpg).
The choice of parameters. -Depending on the parameter values, our model of an oscillatory colloid can describe various systems, including, for example, populations of biological microorganisms that cyclically change their shapes. The focus in the present study is however on the nanoscale phenomena within single living cells. Therefore, the parameter values typical for active proteins in biological cells shall be used.
Protein domains, corresponding to dumbbell beads in our model, typically have the size of about 10-20 nm. They are so stiff that one domain cannot deform another and penetrate inside it. Therefore, a repulsive hard-core potential could have been a good candidate to describe repulsive interactions between them. There are however also soft electrostatic interactions between proteins that extend over a few nanometers. Additionally, a protein is surrounded by a layer of hydrated water that is soft [18]. Recently, non-contact effective interactions between proteins in water were analyzed by direct molecular dynamics simulations [19]. Therefore, the hard repulsion core is effectively surrounded by a soft interaction shell. In our simple model, a single interaction potential (5) is used. Nonetheless, one can approximately interpret this potential as having a hard core and a soft outside shell. The boundary between them can be set at a radius σ at which the repulsion potential becomes much larger than the thermal energy k B T . Below, we take u 0 = 100 k B T /R 2 in eq. (5) and define the bead "hard-core" radius by the condition that u(r = σ) = 40k B T , which yields approximately σ = 0.68R (see fig. 2). Note however that such radius σ is sensitive to the choice of the threshold potential value. It should be stressed that the model does not have a true hard-core repulsion potential and the beads can penetrate one into another to some extent.
If particles of linear size σ are randomly distributed at the mean distance of d between them, their volume fraction is about φ 3D = (σ/d) 3 or the area fraction about φ 2D = (σ/d) 2 in the 2D case, where the numerical prefactors that depend on particle shapes are dropped.
In bacteria, about 30 percent of cytosol is typically occupied by proteins [20]. Assuming that all proteins have the same size σ and using the above estimate, this yields the relative distance of about d/σ = 0.3 −1/3 1.5 between them. If, for example, σ = 0.68R, this corresponds to d R ≈ 1.47σ, so that the neighbours of a protein would 40003-p3 typically be located within a soft interaction shell from it. A similar situation is characteristic for lipid membranes in biological cells. Typically, protein inclusions make up about 40 percent of the membrane mass and most of such inclusions (like, e.g., ion pumps) are active and cyclically change their shapes. Note that, in contrast to cytosol, the biomembrane represents a 2D colloid. In this letter, simulations for a two-dimensional system are performed. We shall have 246 dumbbells within the area with the linear size of L = 40 R. Assuming that the radius of a bead is σ = 0.68 R, this yields the area fraction of φ 2D = 0.45.
The relaxation time τ of the dumbbell should be about the characteristic slow conformational relaxation time in a protein, which is of the order of milliseconds. On the other hand, the turnover cycle time of an enzyme, corresponding in our model to the modulation period 2π/ω, typically takes tens of milliseconds. Therefore, we can choose Ω = 0.1. The initial oscillation phases will be randomly chosen for different dumbbells and noise with Y = 10 −1 will be moreover included into the phase evolution equation (2).
When choosing the parameters a 0 and a 1 , it is important that, even at the maximal length of a dumbbell, an additional bead cannot enter into the space between the beads within it, i.e., that the condition (a 0 + a 1 ) < 4(σ/R) = 2.72 is satisfied. Below, we choose a 0 = 1.5 and vary a 1 between 0 and 1, so that this condition holds even at the largest oscillation amplitude.
The dimensionless stiffness of the spring connecting the beads is κ = 100. The energy supplied to our model enzyme per single cycle is ΔE = 7.8 k B T at a 1 = 0.5 and ΔE = 31.1 k B T at a 1 = 1. For comparison, the chemical energy supplied in the reaction of ATP hydrolysis, often powering protein machines, is about 20 k B T .
Numerical simulations. -In our two-dimensional simulations, periodic boundary conditions have been used. To prepare the initial configuration, the following procedure was employed: We started with a system where all dumbbells had a zero natural length, 0 = 0, and they were regularly distributed forming a two-dimensional grid. Then, a short numerical simulation of this system over the time interval of 100τ was performed. During this simulation, the system's temperature was raised 5-fold and, moreover, the natural length of dumbbell springs 0 , together with the modulation amplitude 1 , was gradually increased from zero to 0 = 1.5R and 1 = a 1 R. The equations were then integrated over further 100τ . Configurations of particles established at the end of such preparatory simulations were used as initial conditions for the actual simulations with the duration of 50000τ .
To explore the diffusion phenomena, we traced motion of the centers of mass, ρ i = (r (1) i + r (2) i )/2, for all dumbbells i. Thus, trajectories were obtained that could be further analyzed. Averaging was always performed over all 246 different dumbbells in one simulation.
Note that at very short times, corresponding to displacements of mass centers much less than the mean distance between the particles, interactions between the dumbbells do not play a role. In this short-time regime, free Brownian motion of an isolated dumbbell, described by the Langevin equation (6) without the interaction terms, should be observed. For the considered system, the diffusion coefficient for an isolated dumbbell is D 0 =0.005 R 2 /τ . At such very short times, the mean-square displacement (MSD) of a dumbbell from its initial position should be Δρ 2 (t) = 4D 0 t.
The mean-square displacements of dumbbells, determined in our simulations, are shown as functions of time for different activity levels in fig. 3.
It can be noticed that the dependences in this figure have the form characteristic for colloidal glasses [12,21]. As is usually done, they can be analyzed by fitting to power laws t β . As is seen in fig. 3, one gets different exponents β − and β + in the intermediate-and long-time regimes. The crossover between the regimes occurs at times t cross .
Both at short and at long times, classical diffusion is observed (exponents β + vary between 0.97 at a 1 = 0 and 0.99 at a 1 = 1). In the intermediate regime, subdiffusion with the exponents β − ranging between 0.55 and 0.75 takes place. The dependences of the subdiffusion exponent β − and the crossover time t cross on the activity level a 1 are shown in figs. S1 and S2 in the Supplementary Material Supplementarymaterial.pdf (SM).
While classical diffusion is again recovered at long times, the diffusion coefficient D at such times is however smaller than D 0 , indicating that diffusion in the colloid becomes suppressed.
Both the suppression of diffusion at long times and the observed subdiffusion at intermediate times are typical glass effects. They suggest that caging of particles takes place [12]. The free diffusive motion of a particle is blocked by other particles surrounding it and forming a cage. Displacements over large distances are only possible if, by a rare fluctuation, the particle was able to escape from its cage. The dependence of the long-time diffusion coefficient of particles on the activity level a 1 is shown in fig. 4. We see that diffusion becomes enhanced when non-equilibrium conformational activity of dumbbells is introduced and gradually increased. Alternatively, it can be said that suppression of diffusion with respect to that for free particles becomes then less strong. Furthermore, as is seen in figs. S1 and S2 in the SM, the subdiffusion exponent β − gets larger at higher activity levels and the classical diffusion regime sets on earlier (i.e., at the shorter cross-over times t cross ) with an increase in the parameter a 1 .

40003-p4
Similar behavior is observed in equilibrium colloids when the volume or area fractions of particles are decreased; it corresponds to a transition from the glass to the fluid phase [12,21]. The above results suggest that effective fluidization of a colloid can also take place because of the non-equilibrium conformational activity of the particles forming it.
To further explore this conjecture, structural relaxation in the model has been numerically analyzed. To do this, we have determined the scattering function Because of the isotropy of the system, this function depends only on k = |k|. The determined functions F 2 (k max , t) are shown for three different activity levels in fig. 5(a). Here, k max = 2π/b, where b = L/N 1/2 is the mean distance between the dumbbells. The structural relaxation time τ α is defined [12] by the equation The dependence of τ α on the activity level a 1 is displayed in fig. 5(b). As we can see, structural relaxation gets much faster when non-equilibrium conformational activity in the particles takes place. The structural relaxation time decreases by an order of magnitude in comparison to the passive colloid (a 1 = 0) when a 1 = 1; it becomes then close to the oscillation period T = 2π/Ω = 62.8 of active dumbbells.
To further explore statistical properties of trajectories, we determined statistical distributions of displacements |Δρ| within a given time Δt. For classical diffusion, such distributions should be 40003-p5  Therefore, the ratio p(|Δρ|)/|Δρ| has then the Gaussian form.
In fig. 6(a), we show the normalized histograms where the frequencies of displacements at three activity levels within time Δt = 1000τ are divided by |Δρ|. Additionally, this figure shows the fits of such distributions to the Gaussian form. It can be noticed that, for a passive colloid, the distribution deviates from the Gaussian dependence for large displacements Δρ. When non-equilibrium conformational activity is switched on, the deviations become smaller and they practically disappear at a 1 = 1.
To quantitatively characterize the deviations, the non-Gaussianity coefficient is introduced, where ζ d = 1 + 2/d and d is the dimensionality of the system. This coefficient vanishes in the Gaussian case.
The computed non-Gaussianity coefficients at different activity levels are shown in fig. 6(b). One can see that such coefficients decrease at higher activity levels, thus further supporting our conclusion that fluidization of the system takes place.
To check for possible finite-size effects, some simulations have been repeated for a system of the double linear size. Figure S3 in the SM demonstrates that the computed time dependences of MSD for the single-and double-size systems practically coincide.
In another test, we have taken as an initial condition a snapshot from the oscillatory system with a 1 = 1. Then, the simulation was continued, but with the natural lengths i of dumbbells frozen at their values in such initial snapshot. Thus, an equilibrium system was constructed with exactly the same, but static, size distribution of dumbbells as in the oscillatory case (see supplementary video S3.mpg).
The dependence of MSD on time in the frozen system is shown in fig. S4 in the SM. From this dependence, we could find that the long-time diffusion constant is D frozen /D 0 = 0.072. This value is closer to the diffusion constant D/D 0 = 0.045 in the passive colloid and much less than the diffusion constant D/D 0 = 0.45 in the respective oscillatory case. Hence, we have demonstrated that active oscillations, not the dispersion of sizes or possible overlaps, are responsible for the observed diffusion enhancement and for the fluidization of a colloidal glass.
Discussion. -In this letter, we have proposed a simple model of a non-equilibrium colloid of active particles reciprocally changing their shapes. The parameters of the model could be chosen in such a way that the conditions typical for metabolically active bacterial cytoplasm and for cellular membranes with active protein inclusions were roughly reproduced.
Colloids of circular particles, whose radii periodically changed, were previously considered [22] as a model for epithelial cellular tissue. This model is not applicable to protein colloids since, due to incompressibility of water and lipid bilayers, the volume of a particle should be conserved. Instead of the glass behavior, an analog of the yielding transition in amorphous solids was observed [22].
In contrast, our simulations have demonstrated that passive colloids of dumbbells behave as glasses, but the glass properties fade away when conformational oscillations are introduced. They were performed for 2D systems, but, as typical for crowded colloids [12,21], similar behavior can be expected in the 3D case too. Thus, we conclude that, under sufficiently strong conformational activity of the particles, effective fluidization of a glass-like colloid may take place. The fluidization becomes possible at the rate of energy supply of about 10 k B T per a particle and a turnover cycle.
Previously, a similar conclusion has been drawn in the study [13] where positions of fluorescent genetically engineered particles with the sizes varying from 50 nm to 150 nm were optically tracked in E. coli and in yeast cells. Classical diffusion was characteristic in presence of metabolism. However, the diffusion was much slowed down and pronounced glass properties were observed under starvation conditions or when ATP, i.e., the chemical fuel, was depleted in the cells.
The experimental effects depended however strongly on the size of the probe particles, getting less pronounced for the tracers of a smaller size [13]. In our study, additional probe particles were not introduced and only trajectories of small dumbbells themselves were tracked. Therefore, our results cannot be directly compared with them. Nonetheless, some comments can be made.
According to ref. [23], cages in colloidal glasses are characterized by a hierarchical onion-like structure, with the smaller ones enclosed within the larger ones. If a smaller cage is destroyed and a particle escapes from it, it may still stay confined within a larger cage. Moreover, characteristic times of diffusion processes get progressively increased with the cage size, following a self-similarity law [23]. It seems natural to assume that a large probe particle can be only confined within the cages of the respective large size. Moreover, one might also expect that diffusion for probe particles would be qualitatively the same as for the smaller particles, though scaled up in time. The crossover from subdiffusion to classical diffusion would occur at MCD about the probe particle size. In the experiments [13], only displacements at relatively long times could be resolved. For the smallest probe particles, the observation times could have been shorter than t cross , so that the glass-like effects indeed remained weak.
For glass-like colloids, it is known that their statistical single-particle properties can be already reproduced using small systems, whereas many-particle correlations should have strong size effects [21] (see also [12,23]). Thus, we could reliably investigate diffusion of individual particles in the present study with a small system, but further investigations aimed, e.g., at exploring dynamic heterogeneity effects, shall require working with the systems of a larger size. In the future, simulations accounting for hydrodynamic interactions can also be performed.
Based on our results, strong diffusion enhancement can be expected not only in bacterial cytoplasm, but also for biological membranes with protein inclusions, provided that chemical energy needed to maintain shape oscillations in membrane proteins is persistently supplied. It would be interesting to experimentally check this.
The observed fluidization of a colloidal glass can be interpreted as an effect of an additional non-thermal noise caused by active conformational changes in the particles constituting it. Hence, our study supports the conjecture [24] (see also [25]) that, in presence of metabolism, active intracellular noise might prevail over thermal noise and determine transport phenomena in the cells.
At the end, we want to stress that, although the discussion in this letter was centered on the nanoscale processes within biological cells, the proposed model is universal; it can be applied to other, natural or artificial, systems at different length and time scales as well. Indeed, artificial non-equilibrium colloids of small oscillating dumbbellshaped particles can be readily designed. * * * Stimulating discussions with T. Ooshida and S. Komura are gratefully acknowledged. This work was supported by JSPS KAKENHI Grants No. JP19K03765 and No. JP19J00365 in Japan.