Inelastic mechanics of sticky biopolymer networks

We propose a physical model for the nonlinear inelastic mechanics of sticky biopolymer networks with potential applications to inelastic cell mechanics. It consists in a minimal extension of the glassy wormlike chain (GWLC) model, which has recently been highly successful as a quantitative mathematical description of the viscoelastic properties of biopolymer networks and cells. To extend its scope to nonequilibrium situations, where the thermodynamic state variables may evolve dynamically, the GWLC is furnished with an explicit representation of the kinetics of breaking and reforming sticky bonds. In spite of its simplicity the model exhibits many experimentally established non-trivial features such as power-law rheology, stress stiffening, fluidization, and cyclic softening effects.


Introduction
In many studies of cell mechanics and dynamics, the cell is characterized as a viscoelastic body [1][2][3]. It is an intriguing question to what extent such mechanical behaviour can be rationalized in terms of comparatively simple polymer physics models. In this respect, the comparison of cell rheological data and minimalistic in vitro reconstituted constructs of the cytoskeleton, such as pure actin solutions [4][5][6][7][8] or crosslinked actin networks [9][10][11][12][13][14], has recently provided many new insights. Important progress has also been achieved in the development of phenomenological mathematical descriptions. This includes approaches related to the tube model [15][16][17][18], tensegrity-based approaches [19][20][21], effective-medium models [22][23][24], and some others [25,26]. In particular, the glassy wormlike chain (Gwlc) model [27], a phenomenological extension of the standard Wlc model of semiflexible polymers [28] has been successful in describing a plethora of rheological data for polymer solutions [7,29] and living cells [30] over many decades in time with a minimum of parameters. However, all these studies were primarily concerned with viscoelastic behaviour, while the latest investigations have underscored the glassy [31][32][33] fragile [34,35], and inelastic [34][35][36] character of the mechanical response of living cells. Even for biopolymer networks in vitro, experiments operating in the nonlinear regime had so far to resort to special protocols that minimize plastic flow [8,37] in order to make contact with dedicated theoretical models.
The aim of the present contribution is to overcome this restriction by extending the Gwlc to situations involving inelastic deformations. As a first step, we concentrate onto reversible inelastic behaviour, where the deformation does not alter the microscopic ground state. The protocol applied by Trepat et al. [34] provides a paradigmatic example. Cells are subjected to a transient stretch such that, after some additional waiting time in the unstretched state, the (linear) material properties of the initial state are recovered. The simplification for the theoretical modelling results from the assumption that not only the macro-state but also the micro-state of the system may to a good approximation be treated as reversible under such conditions; i.e., we assume that the complete conformation of the polymer network, including the transiently broken bonds between adjacent polymers, is constrained to eventually return to its original equilibrium state. For the time-delayed hysteretic response of the network to such protocols one could thus still speak of a viscoelastic ("anelastic") response in an operational sense, but we refrain from doing so in view of the fundamentally inelastic nature of the underlying stochastic process -in contrast to the reversible softening effects observed in [38], for example. Indeed, by simply allowing bonds to reform in new conformational states, the model developed below can readily be extended to arbitrary irreversible plastic deformations, as will be demonstrated elsewhere [45]. Before entering the discussion of our model, we would also like to point out that the proposed (inelastic) extension of the Gwlc is strongly constrained by identifying the newly introduced parameters with those of the original (viscoelastic) Gwlc model, where possible. Despite its increased complexity, the extended model will therefore enable us to subject the underlying physical picture to a more stringent test than hitherto possible by comparing its predictions to dedicated experiments. Moreover, unlike current state-of-the-art simulation studies [25] it is not limited to rod networks but is firmly routed in a faithful mathematical description of the underlying Brownian polymer dynamics. This paper is organized as follows. First, we review some basic facts about the Gwlc in section 2.1. Next, in section 2.2, we introduce our extended reversible inelastic version, which we formulate using the notion of an effective interaction potential as in the original construction of the Gwlc in [27]. (A preliminary account of the basic procedure and some of its cell-biological motivation including reversible bond-breaking kinetics has recently been given in a conference proceedings [39].) Sections 3.1 and 3.2 explain the physical mechanism underlying the mechanical response under pulsed and periodically pulsed loading, while section 3.3 illustrates its phenomenology. We demonstrate that the model exhibits the hallmarks of nonlinear cell mechanics: strain/stress stiffening, fluidization, and cyclic softening [34][35][36]. Section 3.4 investigates the relevance of the lately quantified structural heterogeneities in networks of semiflexible polymers [40] for the mechanical properties, before we conclude and close with a brief outlook.

The glassy wormlike chain
The glassy wormlike chain (Gwlc) is a phenomenological extension of the wormlike chain (Wlc) model, the well-established standard model of semiflexible polymers. A broad overview over Wlc and Gwlc dynamics can be found elsewhere [28]. The Wlc describes the mechanics of an isolated semiflexible polymer in an isothermal viscous solvent. In the weakly bending rod approximation, a solution of the stochastic differential equations of motion for the Wlc is possible via a mode decomposition ansatz for the transverse displacement of the polymer contour from the straight ground state. The individual modes labelled by an index n are independent of each other and decay exponentially with rates τ −1 n . For convenience, we set the thermal energy k B T = 1, so that the bending rigidity can be identified with the persistence length p , in the following. Using this convention, the Wlc expression for the transverse susceptibility of a polymer of contour length L (with respect to a point force) reads [27] α Wlc (ω) = L 3 Here, f is an optional backbone tension, and f L ≡ π 2 p /L 2 the Euler buckling force for a hinged rod of length L. The different powers of n in the denominator give notice of the competition of bending and stretching forces. In the Gwlc, interactions of the polymer with the surrounding network (e.g. excluded volume interactions or sticky contacts) are reflected in an altered Wlc mode relaxation spectrum [27,29]. The intuitive -albeit not fully microscopic -picture underlying the formulation of the Gwlc is illustrated in figure 1, which depicts a test polymer reversibly bound to the background network via potential wells at all topological contacts, the so-called entanglement points. Consider now a generic point somewhere along the polymer backbone (not an entanglement point). It can relax freely until the constrictions are being felt, which slow down the contributions from long wavelength bending modes. The Gwlc translates this intuition into a simple prescription for the mode spectrum: the short-wavelength modes are directly taken over from the wormlike chain model, while modes of wavelength λ n longer than the typical contour length Λ between adjacent bonds are modified to account for the slowdown. Motivated by the physical picture illustrated in figure 1, the slowdown of the relaxation of a wavelength λ n is expressed by an Arrhenius factor exp [E (λ n /Λ − 1)] for breaking (λ n /Λ − 1) potential energy barriers of height E simultaneously. Accordingly, the phenomenological recipe to turn a Wlc into a Gwlc reads: Upon inserting this into (1) one obtains the Gwlc susceptibility α Gwlc (ω). The microscopic "modulus" for transverse point excitations of a generic backbone element on a test polymer is then defined as the inverse susceptibility, g Gwlc (ω) = 1/α Gwlc (ω). An approximate expression for the macroscopic shear modulus is obtained along similar lines [22,27].
In the original equilibrium Gwlc theory, Λ was assumed to be a constant on the order of the entanglement length of the network, Λ L e . Note, however, that L e is a geometric quantity (which is determined by the polymer concentration and the persistence length) while the contour length Λ between closed bonds clearly depends on the state of the bond network. One therefore has to allow for an increase of Λ with the number of broken bonds in non-equilibrium applications. This issue is explored in the following.

Effective interaction potential and bond kinetics
All mechanical quantities calculated within the Gwlc model crucially depend on the interaction length Λ. In previous applications of the model [7,27,30] it was assumed that Λ remains constant -equal to its equilibrium value and unaffected by the deformation of the sample. In other words, the equilibrium theory allowed for statistical bond fluctuations but not for a dynamical evolution of the parameters characterizing the thermodynamic state of the bond network. An obvious starting point for generalizations of the model to non-equilibrium situations is therefore to consider the number of closed bonds, and therefore also Λ, as dynamic variables, dependent on the strain-and stresshistory of the network. We now provide a mean-field description to account for such a dynamically evolving bond network. For clarity, we return to the intuitive picture underlying the Gwlc where the (possibly crosslinker-or molecular-motor-mediated) complex interactions between the polymers are summarized into an effective interaction potential for a test segment against the background, as sketched in figure 2. The same idea has also been used before in many related situations (e.g. [26,41,42]). The generic potential exhibits three essential features: a well-defined bound state, a Sketch of the schematic effective interaction potential for a test segment against the background. A polymer segment can either be "bound" by a transient crosslinker or a sticky patch on the backbone of another polymer (left potential well) or "unbound" and merely confined by the surrounding network to a tube like cage (right well). When considering an ensemble of contacts, the fraction of closed bonds depends on the free energy difference U between both states. The transition rates k − and k + between the bound and unbound state depend on the barrier heights E and (E − U ), respectively. An externally applied force can tilt the potential and favour binding or unbinding events.
well-defined unbound state, and a barrier in between (figure 2). The confining well corresponding to the unbound state represents the tube-like caging of the test polymer within the surrounding network [15,17,43]. For ease of notation, we further introduce the dimensionless fraction of closed bonds or, bond fraction which is simply the ensemble-averaged fraction of sticky contacts that are actually in the bound state. The minimum bond distance Λ 0 is typically on the order of L e , but may be somewhat larger in situations where the bonds are mediated by crosslinker molecules or by partially sterically inaccessible sticky patches (as e.g. for helical molecules [12,44]). A quantitatively fully consistent way of calculating the dynamics of ν would involve solving the full Fokker-Planck equation for a Wlc-Wlc contact -a formidable program to be pursued elsewhere [45]. Here, for the sake of our qualitative purposes, we chose to concentrate onto the physical essence of the discussion, and keep the mathematical structure as simple and transparent as possible. We therefore approximate the dynamics by a simple exponential relaxation as familiar from the standard example of reacting Brownian particles, conventionally described by Kramers theory with a Belllike force dependence [46,47]. Using this simplification and assuming a schematic interaction potential as depicted in figure 2, the value of ν is determined by a competition of bond breaking and bond formation with reaction rates k − and k + , respectively. Both rates are represented in the usual adiabatic approximation (meaning that the equilibration inside the wells is much faster than the barrier crossing and external perturbations) by [46] where τ 0 is an intrinsic characteristic Brownian time scale for bond breaking and formation ‡, and f is the force pulling on the bond. Noting that the fraction of open bonds is 1 − ν, we can then write down the following rate equation for the fraction of closed bondsν The time dependence of ν(t) leads via (3) to an implicit time dependence of the Gwlc-parameter Λ(t) and thereby of all observables derived from the Gwlc. The time-dependent force f (t) in (5) may be thought to result from an externally imposed stress protocol or from internal dynamical elements such as molecular motors setting the network under dynamic stress. Hence, via an appropriately chosen set of slowly changing state parameters f (t), U (t), E(t), . . . the model can in principle accommodate for the active biological remodelling in living cells and tissues [39]. (For a discussion of the relation between the microscopic f and the macroscopic stress, see [29].) Note that for constant force, the stationary value of ν, obtained by settingν ≡ 0, does not depend on E, as it should be (the steady state is independent of the transition state).

Viscoelastic properties
Whenever the bond kinetics can be disregarded (ν = 0), the viscoelastic properties are simply those of the bare Gwlc [7,27,29,30], which can basically be characterized as Wlc short time dynamics followed by a slow, highly stretched logarithmic relaxation resembling power-law rheology with a non-universal exponent within the typical experimental time windows. Strong perturbations (stress or strain) then give rise to a pronounced stiffening due to the steeply nonlinear response of semiflexible polymers under elongation [2]. This well established behaviour can be understood as the canvas ‡ Strictly speaking the fluctuations in the potential wells at x b and x u are characterized by different Brownian time scales depending on the width of the potential wells. At the present stage, we do not bother to distinguish these time scales nor to fine-tune their numerical values, and identify them, for the sake of simplicity, with the entanglement time scale τ Le in our numerical calculations. against which we aim to discern the characteristic mechanical signatures of the bond kinetics.
An observable that will be of particular interest in the following is the complex microscopic "modulus" g ≡ g Gwlc (ω 0 , Λ 0 /ν, f ), introduced in the previous section. It is used to determine the linear as well as (via a nonlinear update scheme, see Appendix A) the nonlinear force response of the system. To understand how the time-dependence (5) of ν affects this important quantity, it is instructive to first examine the dependence of g on ν. To this end, we formally consider ν temporarily as an independent variable instead of determining it from (5). Note that this approach nevertheless makes sense as for a fixed set of parameters, ν can still take any value between zero and one. This is due to the freedom of choosing an initial state, which can be imagined to have evolved from the prior deformation history.
For fixed other parameters, both the real part g (figure 3a) and the imaginary part g of g increase monotonically with ν. We emphasize that g is not simply proportional to g and therefore the loss angle δ = arctan(g /g ) also depends on ν (figure 3b). For small ν, the loss angle is large, corresponding to fluid-like behaviour. With increasing ν, the system becomes more solid-like. Increasing the barrier height E makes the dependence of g on ν more pronounced ( figure 3). Conversely, as can be expected, the dependence on the barrier height E vanishes with decreasing bond fraction ν. Note that due to the boundary conditions the limit of a Newtonian fluid (δ = π/2) is not recovered when formally taking the limit of vanishing bond fraction (ν → 0, Λ → ∞). Finally, the potential difference U solely determines the dynamics of ν via (5), and therefore influences the modulus g solely through the equilibrium value for ν.
After these general considerations, we now concentrate on the nonlinear and nonequilibrium dynamic response of the extended inelastic Gwlc model, which results from the coupled relaxation of the viscoelastic polymer network and the transient bond network.

Stress-stiffening versus rate-dependent yielding
A convenient way to characterize the mechanical properties of an inelastic material is a force-displacement curve. For a perfectly linear elastic (Hookean) body, it would simply consist of a straight line, whereas a perfectly plastic body would feature as a rectangle delineated by a yield force f y and an arbitrary plastic strain value. For our qualitative purpose, we identify the average transverse displacement of the test polymer segment (which is used to determine the force response, see Appendix A) with the reaction coordinate x of the schematic potential sketched in figure 2; and we identify the force f entering the reaction rates with the backbone tension of the test polymer. As a characteristic length scale for the transverse displacements we use the width (x u − x t ) of the effective confinement potential, which is a measure of a typical mean-square displacement of the polymer in the unbound state. Using this convention, we now consider the effect of a time-symmetric displacement pulse on the evolution of the force f (t) and bond fraction ν(t). For technical convenience, we use a Gaussian shape for the displacement pulses, but the qualitative conclusions to be drawn are largely independent of the precise protocol. The duration of the displacement pulse, which sets the time scale for the dynamic response, is used as the unit of time in the following. Here, we are not interested in short-time tension-propagation and single-polymer dynamics [48], hence a lower bound for pertinent pulse durations is provided by the interaction time scale τ Λ 0 τ 0 . On the other hand, for pulse durations longer than τ 0 e E the system deformation would be quasistatic so that no genuinely dynamic effects of the bond kinetics could be observed.
For a small Gaussian displacement pulse of relative amplitude 0.73 (c.f. figure 4a and Appendix A), the shape of the curve predicted by our inelastic Gwlc model shows all features of a viscoelastic medium (figure 4b, blue dashed curve). It starts with a nearly linear regime for relative displacements 0.4, where a weak stiffening sets in. Due to dissipation in the medium, the path back to the initial state takes its course at lower force, which causes a weak hysteresis. This is essentially the viscoelastic response that already the bare Gwlc model would have predicted. For a larger relative deformation amplitude of 2.2, however, the predictions of the bare Gwlc and the inelastic Gwlc diverge. The bare Gwlc predicts a strong stiffening (figure 4b, red dotted curve), whereas the full model exhibits an initial stiffening regime followed by a pronounced softening (figure 4b, solid red curve). A very interesting observation is that "softening" in this case does not only mean a decrease of the slope of the force-extension curve, but that the slope actually becomes negative over an extended parameter region, once an operationally defined threshold force f y is exceeded. This "flow state" bears more resemblance to plastic flow than to viscoelastic relaxation. The reason for this qualitatively new phenomenon is the force-induced bond breaking. Theoretically this is best exemplified by the time-dependence of the (experimentally not easily accessible) x u x t f Ν min d Figure 4. Transition from a stiffening to a softening response at a rate-dependent, operationally defined yield force f y : Gaussian deformation pulses (a) with relative displacement amplitudes of 0.73 (dashed) and 2.2 (solid); corresponding forcedisplacement curves and bond fraction evolutions (b, c; dotted curve in b: response of a Gwlc without bond breaking); order-of-magnitude estimate for the rate-dependence of the yield force, (7) with const.= 2.2, versus the numerically obtained maximum forces of the force-displacement curves for various mean force rates (d); the minimum bond fraction (d, inset) depends non-monotonically on the force rate; E = 10.5, U = 2, fraction ν(t) of closed bonds (figure 4c). A glance at the bond fraction during the weak deformation scenario (figure 4c, blue dashed curve) reveals how the limit of a bare Gwlc is recovered, namely whenever the deformation is not sufficiently violent and persistent to significantly decrease the fraction of closed bonds. For the large deformation scenario, in contrast, the bond fraction stays only initially constant (figure 4c, solid red curve). As a consequence of the large strain the force rises steeply, as can be seen from the strong stiffening in figure 4b (solid red curve). At the yield force, the bond fraction suddenly decreases to nearly half of its initial value during a very short time. The decrease in the bond fraction is accompanied by a somewhat slower drop in the force, which is reflected in the sudden softening of the force-displacement curve. The bond fraction eventually recovers on a much slower time scale, which is roughly given by τ 0 exp(E −U ). Note that even though the effects of the inelastic response dominate the stress-strain curve, the viscoelastic relaxations from the underlying Gwlc model are still present. They could hardly be disentangled from the inelastic contributions, though, without an underlying faithful model of the viscoelastic response at hand.
In summary, we observe a competition between force-amplitude dependent stressstiffening and force-rate dependent yielding events. If the backbone force stays much smaller than the yield force, the adaptation of the bond network requires an adiabatically long time, on the order of the bond lifetime τ 0 e E , so that plain viscoelastic Gwlc behaviour is observed on the pulse time scale. In contrast, if the backbone tension f reaches the yield force f y , the bond fraction declines sharply, thereby switching the response from the Wlc/Gwlc-typical stiffening to a pronounced softening.
The rate dependence of the yield force f y can be estimated by setting the time it takes to reach the yield force, f y /ḟ , equal to the force-dependent time scale of bond opening, k −1 − (f ), from (4), (7) Here, LW(x) is the positive real branch of the Lambert W function. In figure 4d, this estimate is compared with results from numerical evaluations for Gaussian displacement pulses at different average rates. Apart from the numerical errors, the slight deviations from the estimate may be attributed to the fact that the force rate is not constant for the Gaussian protocol. They can be mostly eliminated by using force ramp protocols of various slopes instead of the Gaussian displacement pulses. For not too low rates, the rate dependence can be approximated by the even simpler relation where the force rate has to be normalized by a suitable force scale. The minimum bond fraction reached during the application of the pulse, which we interpret as a measure of the degree of fluidization, depends non-monotonically on the rate (figure 4d, inset). Qualitatively, this can be understood by noticing that two different factors influence the fluidization, namely the maximum force attained during deformation and the time over which the force is applied. The maximum force is simply the yield force f y , with the rate dependence in (7) and figure 4d, while the reciprocal rate itself sets the time scale. For high rates (x u − x t )ḟ 1 (in our example) the rate-dependence of the maximum force wins and the minimum ν reached decreases with increasing rate. For low rates, where the rate-dependence of the maximum force is much weaker (figure 4d) the minimum ν decreases with increasing pulse duration, viz. decreasing rate. Note that for slow pulses (low rates), the bond fraction after the pulse may be quite different from the minimum ν, as significant recovery may already take place during the pulse.

Transient remodelling and cyclic softening
The substantial changes in the material properties accompanying bond breaking can be exemplified by monitoring the linear elastic modulus g (ω 0 ) (measured at a fixed frequency ω 0 ) in response to a strain pulse (figure 5a). Apart from the usual Wlc/Gwlc-typical stress-stiffening below the threshold force f y , the modulus is seen to drop below the value it had before the deformation pulse, where it apparently saturates (remember that the deformation vanishes for t > 1 and that the recovery takes roughly a time τ 0 e E−U 1). We thus observe what we call a "passive", "physical" remodelling of the bond network as opposed to the "active", "biological" remodelling of the cytoskeleton of living cells in response to external stimuli. These passive remodelling effects have recently been observed for human airway smooth muscle cells [34]. (A more quantitative discussion of the relation to the experimental observations will be given elsewhere.) The fact that the deformation pulse leads to a decrease in g and an increase in the loss angle δ (figure 5c) suggests the notion of fluidization [34].  Figure 5. The internal force f and the bond fraction ν largely determine the storage modulus g and the loss angle δ: normalized storage modulus g n (a) and loss angle δ (c) during displacement pulses of large and small amplitudes, 2.2 (solid), 0.73 (dashed); the stronger deformation pulse causes a strong and persistent overshoot of the loss angle, indicating fluidization; slow recovery after the pulse, (b, d); parameters as in figure 4, ω 0 = 3.
The passive remodelling described here is a transient phenomenon, because, after cessation of the external perturbations, the bond fraction will eventually recover its equilibrium value. This indicates that also the change in the system properties is a transient effect, which is demonstrated by the recovery from fluidization in figure 5 b & d. Thus, while the fluidization resembles a plastic process on short time scales, on long times scales, the phenomenology is more similar to pseudo-and superelasticity as observed in shape-memory alloys [49].
To demonstrate that the transient passive remodelling also affects the nonlinear material properties, we consider a series of three pulses, next (figure 6). The force response to such a protocol is depicted in figure 6b. The force-displacement curve for the second stretching (second left branch of the solid curve) is less steep than for the first stretching (very left branch of the solid curve) and the yield force is lower. This indicates a cyclic softening, or viscoelastic shake-down effect, closely related to the fluidization of the network by strain. The strength of the shake-down depends on the fraction of transiently broken bonds, and hence on the rate as well as on the amplitude of the applied deformation (c.f. figure 6c, solid red curve). Upon repeated application of deformation pulses the force-displacement curves settle on a limit-curve that is essentially preconditioned by the initial deformation pulse and the inelastic work it performed on the sample. For more gentle protocols that only break a smaller fraction of bonds, the initial fluidization would be not as pronounced and one would obtain a gradual shake-down, which converges to a limit curve only after many cycles.  Figure 6. Protocol (a) and response (b, c) for three subsequent Gaussian strain pulses; force-displacement curves (b) for a single minimum interaction distance (Λ 0 = 22.3, solid) and for a distribution of Λ 0 according to [40] with mean Λ 0 = 22.3 (dashed); the bond fraction ν (c) for Λ 0 = 0.5 Λ 0 (dash-dotted), Λ 0 (solid), and 2Λ 0 (dashed); other parameters as in figure 4.

Introducing multiple length scales
So far, we have assumed that there is one well-defined characteristic length scale Λ 0 for the polymer interactions, which is on the order of the entanglement length L e and interpreted as the minimum average contour distance between adjacent bonds of the test polymer with the background network. This is clearly a mean-field assumption. Recent combined experimental and theoretical studies have established that the local tube diameter and entanglement length in pure semiflexible polymer solutions actually exhibit a skewed leptocurtic distribution with broad tails [40]. To include this information in the above analysis is not straightforward, since the elastic interactions between regions of different entanglement length are not known, a priori. For the sake of a first qualitative estimate, the simplest procedure seems to be to average the above results over a distribution of Λ 0 , corresponding to a parallel connection of independent entanglement elements. Qualitatively, this renders the abrupt transition from stiffening to softening somewhat smoother (c.f. figure 6b, dashed curve). Moreover, the initial stiffening is slightly more pronounced. This is due to the contributions of Λ 0 < Λ 0 , which exhibit both stronger stiffening and fluidization. Nevertheless, all qualitative features -like stiffening and remodelling effects and a "flow state" indicating fluidization or inelastic shakedown -can still be well discerned. This is consistent with the assumption that pure semiflexible polymer solutions are at least qualitatively well described by a single entanglement length [17,18,22,27,50].

Conclusions and outlook
We have presented a theoretical framework for a polymer-based description of the inelastic nonlinear mechanical behaviour of sticky biopolymer networks. We represented the polymer network on a mean-field level by a test polymer described by the phenomenologically highly successful Gwlc model. Beyond the equilibrium statistical fluctuations captured by the original model, we additionally allowed for a dynamically evolving thermodynamic state variable Λ(t) characterizing the network of thermoreversible sticky bonds, which is updated dynamically upon bond breaking according to an appropriate rate equation. Within our approximate treatment we could obtain a number of robust and qualitatively interesting results, which are not sensitive to the precise parameter values chosen, nor to the choice of the transverse rather than the longitudinal susceptibility in (1). For sufficiently strong deformations, changes in the mean fraction of closed bonds, reflected in Λ(t), can influence material properties to a degree at which the material behaves qualitatively different from the usual viscoelastic paradigm. In particular, our theory predicts a pronounced fluidization response, which develops upon strong deformations on top of the intrinsic viscoelastic stiffening response provided by the individual polymers. After cessation of loading, the system slowly recovers its initial state. These observations are in qualitative agreement with recently published experimental results obtained for living cells [34], and a more quantitative comparison with dedicated measurement results for cells and in vitro biopolymer networks will be the subject of future work. We found the yield force for the onset of the fluidization to be sensitive to the deformation rate. Moreover, the fluidization response was shown to be accompanied by a cyclic softening or shake-down effect. Taking into account the spatial heterogeneity of biopolymer solutions by a distribution of entanglement lengths leads to a smoothing of the force-displacement curves without affecting their qualitative characteristics. One may still expect qualitatively new effects in situations with unusually broad distributions of entanglement lengths, such as for strongly heterogeneous (e.g. phase separating) networks.
A parameter that was found to be very important for cells [51] but has not been discussed much in this contribution is the prestressing force f 0 (see appendix Appendix A), which is present in adhering cells even in the absence of external driving. Experimentally, it was established that increasing prestress is correlated with a higher stiffness of adhering cells [51] which in turn is correlated with a higher stiffness of the substrate [52]. When naively treating the prestressing force just as an additive contribution to the overall force, our model predicts the opposite: the additional force breaks more bonds and the network becomes softer and more fluid-like. To reconcile these apparently contradicting trends, one can appeal to the notion of homeostasis, which essentially amounts to postulating that the cell actively adapts its structure such that a certain set of thermodynamic state variables remain in a physiologically sensible range [39]. This basically implies that the cell will respond to stiff substrates or persistent external stresses by a biological remodelling that corresponds to an increase of E and U , and possibly to a decrease of Λ 0 . By adapting also the internal prestress f 0 , the cell may then avoid, apart from the structural collapse, an equally undesirable loss of flexibility. We note, however, that internal stresses are indeed observed to disrupt the cytoskeleton, as implied by the simple physical theory presented here, if they are not permanently balanced by a substrate [53].
Even though the nonlinear effects presented in this contribution depend crucially on the bond kinetics, it should not be overlooked that the Gwlc is an equally indispensable ingredient of the complete theory. First, it provides the constitutive relation that gives the force in response to an infinitesimal displacement, which in turn governs the bond kinetics. Secondly, the time-dependent linear susceptibility strongly filters the dynamic remodelling of the bond network in practical applications, in which one rarely has access to the microstate of the underlying bond network. In summary, the extended model thus integrates experimentally confirmed features of the Gwlc response such as slow relaxation, power-law rheology, viscoelastic hysteresis, and shear stiffening with a simple bond kinetics scheme, which allows less intuitive complex nonlinear physical remodelling effects like fluidization and inelastic shake-down to be addressed.
Finally, we want to point out that it is straightforward to extend the above analysis in a natural way to account for irreversible plastic deformations [45]. To this end, one has to allow for the possibility that the transiently broken bonds reform somewhere else than at their original sites (as always assumed in the foregoing discussion) after strong deformations with finite residual displacement. This can be included by accommodating an additional term acting as a "slip rate" in (5) [45].
We now face the problem that (A.2) is an implicit equation for f Gwlc (t), which depends on ν(t), which in turn depends via f (t) on f Gwlc (t). We therefore use a two-step Euler scheme: As initial values for f and ν, we choose the prestressing force, f (t = 0) = f 0 , and the steady-state value under the prestressing force, ν(t = 0) = ν(t → ∞|f = f 0 ), respectively. We then choose a sufficiently small time step ∆t and apply the following iterative rule: f (k∆t) ≈ f 0 + f GW LC (k∆t|f ((k − 1)∆t, ν((k − 1)∆t)) (A.5) ν(k∆t) ≈ν(k∆t|f (k∆t)), (A. 6) where k ∈ N + . In the limit ∆t → 0, this procedure converges to the exact solution.