Non-Markovian modeling of non-equilibrium fluctuations and dissipation in active viscoelastic biomatter

Based on a Hamiltonian that incorporates the elastic coupling between a tracer and active particles, we derive a generalized Langevin model for the non-equilibrium mechanical response of active viscoelastic biomatter. Our model accounts for the power-law viscoelastic response of the embedding polymeric network as well as for the non-equilibrium energy transfer between active and tracer particles. Our analytical expressions for the frequency-dependent response function and the positional autocorrelation function agree nicely with experimental data for red blood cells and actomyosin networks with and without ATP. The fitted effective active-particle temperature, elastic constants and effective friction coefficients of our model allow straightforward physical interpretation.

Based on a Hamiltonian that incorporates the elastic coupling between a tracer and active particles, we derive a generalized Langevin model for the non-equilibrium mechanical response of active viscoelastic biomatter.Our model accounts for the power-law viscoelastic response of the embedding polymeric network as well as for the non-equilibrium energy transfer between active and tracer particles.Our analytical expressions for the frequency-dependent response function and the positional autocorrelation function agree nicely with experimental data for red blood cells and actomyosin networks with and without ATP.The fitted effective active-particle temperature, elastic constants and effective friction coefficients of our model allow straightforward physical interpretation.
Viscoelastic gels such as permanently or transiently cross-linked networks of semiflexible polymers are important soft biological materials [1][2][3].The polymeric nature of such gels is responsible for their salient rheological properties, including their frequency-dependent response to external forces.The cell cytoskeleton is a prime example of a viscoelastic gel [4,5] and determines cell shape, motility and division [6,7].Even though the cytoskeleton has a complex structure, its mechanical properties at small deformations are determined mainly by the underlying F-actin network [8,9].Even in equilibrium, the dynamic properties of F-actin are complex and interesting [10].In vivo, constant F-actin (de)polymerization and interactions with myosin motors consume naturally occurring adenosine-tri-phosphate (ATP) and drive the viscoelastic network into a non-equilibrium (NEQ) state [11][12][13], where the motors act as active cross-links and produce local tension [2].
Another illustrious example of active viscoelastic biomatter are red blood cells (RBCs), whose structural components include a filamentous protein (spectrin) scaffold underneath the enclosing cell membrane [14][15][16][17][18].The mechanical response and spectra of membrane fluctuations (RBC flickering) display active NEQ signatures that are attributed to the membrane ion-pump activity and the ATP-induced reorganization of the spectrin network [19][20][21].
Supported by high-resolution experimental data [13,18,20,[22][23][24][25][26], a key line of inquiry has been to examine violations of the fluctuation-dissipation theorem (FDT) in these NEQ systems.The FDT holds for nearequilibrium (near-EQ) dynamical processes, but not for active NEQ processes, as indeed corroborated by microrheology experiments [13,20].Typical experiments involve a micron-sized tracer bead immersed in the active medium whose motion is tracked and which is either allowed to move spontaneously or is driven by a laser beam (passive versus active microrheology).Active microrheology allows to simultaneously measure the mechanicalresponse and the positional autocorrelation function and thereby to check the validity of the FDT.
Theoretical studies of FDT violation focused on externally forced systems [27][28][29], glasses [30] and active systems [13,20,31].In a generalized approach based on the Langevin equation, NEQ fluctuations have been modeled by an additive athermal noise acting on particles in a viscous fluid [21,32].A similar strategy has been used to model active dynamical undulations of elastic biological membranes [21,[33][34][35] and several computational models have explicitly incorporated active pumps [20,36].A non-equilibrium multicomponent elastic model for the interaction of motors and a tracer particle was recently shown to describe the frequency-dependence of the experimentally observed FDT violation [13] rather well [37].
In this paper, we model the non-equilibrium mechanical response in active biomatter by taking into account the medium viscoelasticity and the elastic coupling between tracer and the active particles.We derive analytical expressions for mechanical-response and spatial correlation functions, by simultaneously fitting these functions to experimental microrheological data for RBCs [20] and reconstituted actomyosin networks [13] we extract all parameters describing the viscoelastic medium and the elastic coupling between tracer and active particles.We find that power-law viscoelasticity describes the experimental data better than the standard Maxwell and Kelvin-Voigt viscoelastic models.The extracted viscoelastic exponent is close to κ = 3/4 for actomyosin networks and close to unity for RBC, in agreement with previous studies.Interestingly, in the presence of ATP the effective temperature ratio between active and tracer particles is above 20 for both systems, demonstrating that these systems are very far from equilibrium.Our analysis and the quantitative agreement established with experiments demonstrates the generic capability of our multicomponent elastic model and our hybrid approach that combines relevant particle-based and hydrodynamic aspects of the problem into a single framework and enables closed-form predictions.
Inspired by typical microrheology experiments [13, 20, arXiv:2207.11307v1[cond-mat.soft]22 Jul 2022 24], we assume a single tracer bead of mass M , embedded in an isotropic viscoelastic network and coupled elastically via harmonic springs of strength k with n active particles of mass m.The tracer is trapped in an external harmonic potential of elastic constant K.The Hamiltonian reads where X is the position of the tracer and x i are the activeparticle positions (i = 1, • • • , n), Ẋ and ẋi are the tracer and active-particle velocities.We assume that the system is isotropic and therefore use one-dimensional displacement variables, the extension to multidimensional variables is straightforward [37,38].The system dynamics is assumed to follow a generalized Langevin equation (GLE), which in addition to the deterministic (Hamiltonian) forces, contains thermal and active random forces as well as time-dependent friction forces that describe the polymeric network viscoelasticity.The latter is modeled using a general memory kernel K(t) [39,40], yielding the GLEs Here, γ and γ a are drag coefficients and F (t) and f i (t) are random forces acting on the tracer and active particles, respectively, with zero mean and second moments given by In NEQ the temperatures characterizing the random forces acting on the active and tracer particles, T a and T , are unequal and the NEQ parameter α = T a /T − 1 becomes non-zero [37].Guided by previous results for biological membranes and polymeric networks [13,[41][42][43] and theoretical works [44,45], we use a normalized power-law memory kernel, where Γ(x) is the Gamma function.Note that by construction, the tracer particle and the active particles themselves obey the FDT since the random force correlation equals the friction kernel K + (t), NEQ is therefore introduced via the coupling of tracer and active particles.Also, the friction kernel of active and tracer particles is assumed to be identical, which reflects that both particles are assumed to be coupled to the same polymeric network.
By summing over the active particles and introducing collective variables x(t) = i x i (t)/n and f (t) = i f i (t)/n, we obtain two coupled GLEs for the tracer and the collective active-particle coordinate.The collective random noise has zero mean and strength By Fourier transforming, Eq.( 3) can be solved for X(ω) = dt X(t) e ιωt to yield with the effective friction kernel and the effective noise given by From Eq. ( 4) we can directly calculate the positional autocorrelation function, C(ω) = dt X(0)X(t) e ιωt .From the linear-response relation , the imaginary part of the response function ˜χ (ω) follows, see Supplemental Material [46] for explicit derivations.The FDT violation can be quantified by the spectral function [37] where e ιωt is the spectral noise autocorrelation and K tot (ω) the real part of the Fourier-transformed memory kernel [46].For an equilibrium system characterized by α = 0, i.e. when the active and tracer particles have the same temperature, the FDT is recovered [47] and Ξ(ω) = 0. Assuming overdamped motion with M → 0 and m → 0, we obtain explicitly  [20] and our model predictions for the positional autocorrelation function ω C(ω)/(2kBT ) in Eq. ( 9) (solid lines) and the imaginary part of the response function χ (ω) Eq. ( 8) (broken lines) for (a) the EQ scenario, i.e.ATP-depleted RBCs, and (b) the NEQ scenario, i.e.RBCs in the presence of ATP.In the insets, the spectral function Ξ(ω) from the experiments (symbols) is compared with our analytical prediction Eq. ( 10).(c) Comparison of the real and imaginary parts of the response function χ (ω) and χ (ω) for the NEQ and EQ scenarios, indicating a larger response (i.e. a softening and less dissipation) in the NEQ case in particular at low frequencies.The dotted lines in all subfigures indicate the model prediction for the real part of the response function χ (ω) (see SI).
where the active-particle number n only appears implic-itly via γa = nγ a and k = nk, so we are left with six parameters, the non-equilibrium parameter α, the powerlaw exponent κ, the tracer and active-particle drag coefficients γ and γa , the tracer confinement strength K and the tracer-active-particle elastic coupling strength k.We consider two recent microrheology experiments on RBC flickering [20] and cross-linked F-actin networks with myosin II molecular motors [13], both reporting FDT violation by independently measuring χ (ω) and C(ω).For both experimental systems, it was shown that the FDT remains valid at short times or high frequencies, while strong FDT deviations are observed at long times or low frequencies.We digitize the experimental data reported for χ (ω) and ω C(ω)/(2k B T ) and first fit their high-frequency limits, which from Eq. ( 8) and Eq. ( 9) are predicted as χ (ω) sin( κπ 2 )/(γ ω κ ) and ω C(ω)/(2k B T ) sin( κπ 2 )/(γ ω κ ) and obtain κ and γ for each data set.Then, we simultaneous fit Eqs. ( 8)-( 9) in the entire frequency range to find the remaining four parameters, the fitting parameters are shown in Table I, see Supplemental Material [46] for details on the fitting procedure and a discussion of the uniqueness of the resulting fit parameters.
We start with the discussion of the RBC data, where the flickering dynamics of the RBC membrane was investigated by tracking the motion of a tracer bead attached to the RBC membrane in an equatorial position.
In the experiments both the in vivo (NEQ) and the ATPdepleted scenarios were considered [20], we denote the latter scenario as EQ although a finite residual ATP concentration is present, the consequences of which will be discussed further below.In Fig. 1a,b we find excellent agreement between the experimental data and the fits for χ (ω) and ω C(ω)/(2k B T ) for both EQ and NEQ scenarios in the entire frequency range.We also show the real (storage) part of the response function χ (ω) from our model as dotted lines in Fig. 1a,b,c, see Supplemental Material [46] for the explicit expression.In the EQ case χ (ω) and ω C(ω)/(2k B T ) agree rather nicely and display a weak minimum around ω/(2π) = 2 s −1 , that is caused by the competing elastic and friction effects, see Supplemental Material [46] for an asymptotic analysis of the response for low and high frequencies.In contrast, in the NEQ case the scaled autocorrelation ω C(ω)/(2k B T ) is much larger than the response function χ (ω), indicating strong FDT violation.The spectral function Ξ(ω) in the insets illustrates that FDT violation occurs predominantly at low frequencies.The comparison in Fig. 1c demonstrates that both real and imaginary parts of the response function increase for NEQ at low frequencies, i.e. the effective stiffness of the network decreases at NEQ (see Supplemental Material [46] for a discussion of the low-/high-frequency asymptotic response).The viscoelastic power-law exponents for NEQ and EQ RBC experiments turn out to be κ = 0.936 and κ = 0.926 (see Table I), respectively, which fall in the previously reported range κ = 0.92−1.07for fresh and aged RBCs [43].
The fitted trap stiffness in the NEQ case, K = 1.66±0.03pN/µm, agrees with the experimentally reported value 1.56±0.3pN/µm [20], for the EQ scenario our fitted stiffness is substantially higher, which presumably reflects added stiffness due to the viscoelastic EQ environment of the tracer particle.Our fit results for the coupling strength between tracer and active particles k is of the same order of K. Also, the fitted values of the tracer bead friction coefficient γ in Table I are close to the experimentally reported value for the bare tracer bead γ = 0.085 pN s/µm [20] for both NEQ and EQ scenarios, the fit results for the active-particle friction coefficients γa are of the same order.The good agreement between the fit parameters and the known experimental parameters shows that the model and the fitting procedure are robust.The non-equilibrium parameter in the NEQ case turns out to be α = 36 and thus indicates strong departures from EQ behavior.In the EQ case we obtain a non-zero value α 0.985, which indicates weak NEQ due to incomplete ATP depletion in the experiment, as reflected by weak deviation between χ (ω) and ω C(ω)/(2k B T ) in Fig. 1a.
Mizuno et al. [13] have studied the in vitro mechanics of active cross-linked actomyosin networks in the presence and absence of ATP, again referred to as NEQ and EQ cases.
The myosin-II motor proteins perform stepwise motion along the actin filaments in the presence of ATP, which causes sliding motion between two filaments connected by motors and thereby produces contractile active forces in the network.As shown in Fig. 2a, b, we again find excellent agreement between the experimental data and our fits for χ (ω) and ω C(ω)/(2k B T ) for both the EQ and the NEQ scenarios in the entire frequency range.In the EQ case χ (ω) and ω C(ω)/(2k B T ) are practically identical and show (similar to RBC) two weak maxima that are well produced by the model, see Supplemental Material [46] for an asymptotic analysis.In the NEQ case the scaled autocorrelation ω C(ω)/(2k B T ) is much larger than the response function χ (ω), indicating strong FDT violation with a magnitude rather similar to the RBC case.For the EQ case the value α = −0.035indicates equilibrium in the absence of ATP.In both EQ and NEQ cases the fitted power-law exponent in Table I is in excellent agreement with κ = 3/4 predicted by theory [1,10,48] and experiments [13].The comparison in Fig. 1c demonstrates that both real and imaginary parts of the response function are at low frequencies for the NEQ case smaller than for the EQ case, i.e. the effective stiffness of the network increases in the NEQ case, in contrast to the RBC case.This behavior is primarily caused by the elastic model parameters: we see that K, which describes the stiffness of the tracer bead confinement, and k, which describes the coupling strength between the tracer and the active particles, increase for actomyosin networks when going from the EQ to the NEQ scenario, whereas the opposite trend is observed for the RBC flickering data, see Table I.
In summary, we introduce a simple model for viscoelastic active biomatter that incorporates the elastic coupling between a tracer bead and active particles and the medium viscoelasticity through a memory kernel and yields closed-form analytical predictions for the autocorrelation and response function of the tracer bead.These predictions are in excellent agreement with experimental data for RBCs and actomyosin networks and in particular allow to quantify the departure from EQ by the NEQ parameter α, which corresponds to the effective temperature ratio of the reservoirs coupled to the active particles and the tracer bead in our model.We find α ≈ 23 and 36 in the presence of ATP and α ≈ 0 and 1 in the absence of ATP in the two experimental systems we considered, clearly signaling the strong NEQ character induced by ATP.Note that the effective temperature of the active particles does not correspond to the actual physical temperature, but rather characterizes the strength of the stochastic force the motor proteins produce and the power they dissipate.Our model has a wide range of applicability for active biomaterials and experimental setups and with suitable modifications can for example predict cross-correlations in two-particle microrheology experiments [49].
We acknowledge Hélène Berthoumieux for useful comments on manuscript.We can split the response function given in Eq. (S.8) into its real and imaginary parts as where A is the real part of the denominator given by Now, for the dissipative part of the response function, we have The Fourier transform of response function reads (S.17) where χ(t) = 0 for t < 0 due to causality.Since χ(t) is a real function one has χ * (ω) = χ(−ω).Thus, Eq. (S.16) yields and from Eq. (S.13) we have Combining Eqs.S.18 and S.19 we arrive at Using the fluctuation-dissipation theorem in the time domain, where θ(t) denotes the Heavyside function, in the Fourier domain we have Comparing Eqs.S.20 and S.21 we realize that in equilibrium we have Ξ(ω) = 0.In the next step, we show that departure from equilibrium happens when α deviates from zero.Using Eqs.S.4, S.11 and S.20, after some intermediate steps we find 22) where we have used the definition T a ≡ T (α + 1).Thus to recover the equilibrium state, i.e.Ξ(ω) = 0, the NEQ parameter α should be equal to zero.

II. RESPONSE FUNCTIONS, POSITIONAL AUTOCORRELATIONS AND SPEC-TRAL FUNCTIONS FOR DIFFERENT MODELS
Using memory kernels for different viscoelastic models (see Table S.1), we can derive different expressions for the response and positional autocorrelation functions.

A. Derivation of memory kernels in Fourier space
Here we briefly review the derivation of the Fourier transform of the memory kernels used in this paper.Power-law model:

Newtonian fluid:
, (S.27) , (S.28) , (S.29) (S.34) , (S.37) One can see that in the limiting case τ M → 0, the Maxwell model becomes equivalent to the Newtonian fluid model.Likewise, the same limit is obtained for the Kelvin-Voigt model in the limits k KV → 0 and k KV,a → 0 and for the power-law model in the limit κ → 1.

III. ASYMPTOTIC ANALYSIS
Defining = ω κ , we expand all functions in the low-frequency limit in powers of according to  Results for the NEQ case.Here, there is a clear deviation between the asymptotic expressions for χ (ω) (dotted black lines) and ω C(ω)/(2k B T ) (dashed magenta lines) in the low-frequency regime.
(c) The comparison of χ (ω) in the EQ and NEW cases shows the NEQ stiffening effects.In all cases, the agreement between the asymptotic expressions and the full model predictions is good in the asymptotic frequency regimes.
In the high-frequency limit, we expand in inverse powers of according to From these asymptotic expansions, we obtain that in the low-frequency regime, χ (ω) displays a maximum located at The physical origin of these maxima is easily obtained from the analytical expressions.In fact, the maximum in the low-frequency regime involves γa , while the maximum in the high-frequency regime does not.This shows that the low-frequency maximum is determined by the active-particle friction, while the high-frequency maximum is dominated by the probe particle friction.This is also obvious from the frequency ω * 1 , at which the positional correlation function displays a maximum in the low-frequency regime, . (S.56) Setting α = 0 gives the same value as the low-frequency maximum of the response function, which implies that the difference in the position and value of the maximum of the positional autocorrelation function and the low-frequency maximum of the response function is imposed by the non-equilibrium parameter only.Note that the high-frequency expression for the maximum of χ (ω) coincides with that of ω C(ω) 2k B T .But, clearly, for NEQ there is no maximum in the experimental data for ω C(ω) 2k B T in the high-frequency regime, which implies that the first two terms of the asymptotic high-frequency expansion are not sufficient for describing the positional autocorrelation function at high frequencies in the NEQ case.

NELS
In this section, we briefly discuss fitting results to the NEQ experimental data (both actomyosin network and RBC flickering cases) for different viscoelastic memory kernels,    S.2.

V. FITTING PROCEDURE
A. Fitting models to experimental data For both experiments discussed in this paper, we have digitized the EQ and NEQ data from Refs.[1,2] using WebPlotDigitizer [3].All fits have been performed in the log-log plane.
To determine κ and γ in the power-law model and γ in the Newtonian fluid and Kelvin-Voigt models, we use the high-frequency limit expression of the dissipative part of the A.A. and R.R.N. acknowledge funding from the European Research Council (ERC) under the European Union's Horizon 2020 Research and Innovation Program under Grant Agreement No. 835117.A.A. acknowledges the hospitality by the Institute for Research in Fundamental Sciences (IPM), Tehran, during the initial stages of this work.
t)e ιωt = γ (S.25)We treat the elastic term of the Kelvin-Voigt model as a potential contribution, which stays outside the memory kernel.Therefore, the memory contribution of the Kelvin-Voigt model is the same as that of the Newtonian fluid.
changed variables as ζ = −ιωt.We note that the integral in Eq. (S.26) converges only for κ < 1.Now, exploiting the memory kernels introduced above, we obtain expressions for the real part of the response function χ (ω), the imaginary part of the response function χ (ω), the positional correlation function C(ω) and the spectral function Ξ(ω).

FIG. S. 1 .
FIG. S.1.Actomyosin network: Asymptotic analysis.The same data is shown as in Fig. 2 in the main text, except that here we include the results from the asymptotic analysis.(a) For the EQ case the asymptotic expressions for χ (ω) (dotted black lines) are indistinguishable from the asymptotic expressions for ω C(ω)/(2k B T ) (dashed magenta lines) in both low-and high-frequency regimes.(b) 1 ) = k (γ + γa ) 2 tan πκ 2 8K K γ2 a + k (γ + γa ) 2 .(S.52)In the high-frequency regime, χ (ω) displays a second maximum located at these maxima are in very good agreement with the experimental data, as shown in Fig. S.1.

FIG. S. 2 .
FIG. S.2.Fitting different viscoelastic models to in vitro actomyosin experiment data: (a-d) the data for the dissipative part of the response function (dark blue triangles), χ (ω), and the normalized positional autocorrelation function (red circles), ω C(ω)/(2k B T ), extracted directly from experimental data for F-actin network reported in Ref [1], are shown.In the inset the spectral function is shown (black circles), Ξ(ω).Continuous curves for ω C(ω)/(2k B T ), χ (ω) and χ (ω) denote the model predictions, given in Sec.II.Fitting parameters are given in Table S.2.

3 .
FIG. S.3.Fitting different viscoelastic models to RBC flickering experiment data: (a-d) the data for the dissipative part of the response function (dark blue triangles), χ (ω), and the normalized positional autocorrelation function (red circles), ω C(ω)/(2k B T ), extracted directly from experimental data for F-actin network reported in Ref [1], are shown.In the inset, the spectral function is shown (black circles), Ξ(ω).Continuous curves of ω C(ω)/(2k B T ), χ (ω) and χ (ω) denote the model predictions, given in Sec.II.Fitting parameters are given in Table S.2.

FIG. S. 4 .
FIG. S.4.Fitting the high-frequency limit of the imaginary part of the response function and rescaled correlation function in the log-log plane.Filled circles and triangles are digitized experimental data for χ (ω) and ω C(ω)/(2k B T ), respectively, broken lines are the asymptotic fits.Fitting parameters for each case are given in the main text.
FIG. S.5.Actomyosin network in EQ (a-f): contour plots of the cost function Eq.S.57 versus different fitting parameters.Tilde symbols denote rescaled quantities with respect to the optimal values extracted from the fitting procedure.White filled circles show the fitting parameter values obtained by our fitting procedure.

TABLE I .
Model parameters obtained from fits to the experimental data FIG. 1. RBC flickering.Comparison of experimental data (symbols)

TABLE S .
1. Memory kernels for different viscoelastic models and their Fourier transform.To obtain the memory kernel for active particles, one has to replace γ → γ a everywhere in this table.Note that for the Kelvin-Voigt model one also has to replace k KV → k KV,a , see Sec.II A for more

TABLE S .
2. Actomyosin network experiment (NEQ case) and RBC flickering experiment (NEQ case): Fitting parameters extracted and