Elasticity of 3D networks with rigid filaments and compliant crosslinks

Disordered filamentous networks with compliant crosslinks exhibit a low linear elastic shear modulus at small strains, but stiffen dramatically at high strains. Experiments have shown that the elastic modulus can increase by up to three orders of magnitude while the networks withstand relatively large stresses without rupturing. Here, we perform an analytical and numerical study on model networks in three dimensions. Our model consists of a collection of randomly oriented rigid filaments connected by flexible crosslinks that are modeled as wormlike chains. Due to zero probability of filament intersection in three dimensions, our model networks are by construction prestressed in terms of initial tension in the crosslinks. We demonstrate how the linear elastic modulus can be related to the prestress in these network. Under the assumption of affine deformations in the limit of infinite crosslink density, we show analytically that the nonlinear elastic regime in 1- and 2-dimensional networks is characterized by power-law scaling of the elastic modulus with the stress. In contrast, 3-dimensional networks show an exponential dependence of the modulus on stress. Independent of dimensionality, if the crosslink density is finite, we show that the only persistent scaling exponent is that of the single wormlike chain. We further show that there is no qualitative change in the stiffening behavior of filamentous networks even if the filaments are bending-compliant. Consequently, unlike suggested in prior work, the model system studied here cannot provide an explanation for the experimentally observed linear scaling of the modulus with the stress in filamentous networks.


Introduction
The mechanical properties of biological cells are governed by the cytoskeleton, a viscoelastic composite consisting of three main types of linear protein polymers: actin, microtubules, and intermediate laments. These lamentous polymers are crosslinked by various binding proteins and constitute a dynamic complex network that maintains the structural integrity of the cell with the capacity for dynamic reorganization needed for active processes. Many in vitro studies have focused on reconstituted networks with rigid crosslinks. [1][2][3][4][5][6][7][8][9][10][11][12] In the cytoskeleton, however, many of the crosslinks are themselves extended and highly compliant proteins. Such exible crosslinks can strongly affect the macroscopic network elasticity. [13][14][15][16][17][18][19][20][21] Indeed, experimental studies show that composite networks can have a linear modulus as low as $1 Pa, while being able to stiffen by up to a factor of 1000. 11,14 Here we analyze 3-dimensional (3D) composite networks theoretically, and we offer physical simulations thereof. Our networks are composed of randomly oriented rigid laments that are connected by highly exible crosslinks, each of which is modeled as a wormlike chain (WLC), 22,23 which has been shown to accurately describe exible crosslinkers, such as lamin. 24,25 In our approach we assume that the laments are much more rigid than the crosslinks, meaning that the network elasticity is dominated by the entropic stretching resistance of the crosslinks.
In our theoretical analysis we adopt the widely employed assumption of affine deformations. 16,19,26 Under this premise, the network is assumed to deform affinely on the length scale of the laments, which in turn is assumed to be much longer than the contour length of the crosslinks. Using a single lament description in the limit of a continuous distribution of crosslinks along the lament, we obtain the asymptotic scaling behavior of the elastic modulus with the stress in the nonlinear regime. We show that in 1-dimensional (1D) networks, the elastic modulus scales with the second power of the stress, whereas it scales with the third power in 2-dimensional (2D) networks. Remarkably, there is no power law scaling in 3D-in fact, the elastic modulus of a 3D composite network increases exponentially with the stress. Numerical evaluation of the affine theory at nite crosslink densities-as opposed to a continuous distribution of crosslinks-shows that (i) the only asymptotic scaling is the modulus scaling with an exponent 3/2 with the stress and that (ii) the dependence on dimensionality of the system is limited to an intermediate-stress regime. These ndings are in agreement with our extensive physical simulations of 3D composite networks. For all cases, the elastic modulus diverges at a nite strain.
Our theoretical analysis is inspired by the mean-eld model proposed by Broedersz et al. 16,26 In sharp contrast to our theoretical analysis and to the results of our physical simulations, however, these authors predict linear scaling of the elastic modulus with applied stress. In particular, for any nite strain, the elastic modulus remains nite in their model. While this linear scaling of the elastic modulus is in accordance with what has been observed experimentally, 13,20,21 we here argue that this model does not adequately capture the elastic response of networks with rigid laments and permanent (i.e., non rupturing or rebinding) crosslinks of nite length.
In ref. 19, the authors ruled out that the experimentally observed approximately linear scaling of the modulus with the stress might be due to enthalpic (linear) stretching compliance of the crosslinks or laments. Here, we complement their analysis by physical simulations that take into account bending of laments. Our results empirically show that the inclusion of bending rigidity does not impact the nonlinear stiffening behavior of composite networks either. We therefore conclude that the theoretical explanation for the linear scaling of the modulus with stress in experiments remains a challenging open problem.
By physical simulations, we also study the role of prestress. We show that in contrast to 1D and 2D networks, 3D networks experience an initial tension due to non-intersecting laments resulting in initially stretched crosslinks, and are therefore prestressed. The modulus in the linear deformation regime is then governed by this prestress; indeed, it is higher than the modulus expected from the affine theory. Our simulations additionally indicate that if the network is allowed to relax initial tension by unbinding and rebinding of crosslinks, the impact of prestress on the elastic modulus in the linear regime becomes insignicant, although the prestress does not relax all the way to zero.
The remainder of the article is organized as follows. First, we present the affine theory of composite networks in Section 2. Under the assumption that deformations of the network are affine on the length scale of the laments, we derive expressions for the stress and modulus in 1D, 2D, and 3D. We then present our physical simulation model and describe our network generation procedure in Section 3. We expand on the implications of our 3D simulation procedure in Section 4; in particular, we explain the emergence of prestress. We then discuss the results of our simulations in the linear deformation regime in Section 5 and indicate which results from the affine theory are still valid. Finally, we analyze the simulation results in the nonlinear regime in Section 6.

Theory
In this section we analytically derive the stress and modulus of a composite network under the assumption of affine deformations on the length scale of the laments. We consider a collection of N rigid laments of length L that are permanently connected by nN/2 exible crosslinks of contour length l 0 , where n is referred to as the crosslink density, i.e., the number of crosslinks per lament. The laments are assumed to be perfectly rigid, i.e., they neither bend nor stretch, and the crosslinks are modeled via the WLC interpolation formula 23 where k B T is the thermal energy, l p is the persistence length and u $ 0 is the end-to-end distance of the crosslink. Assuming l 0 [ l p this force-extension relation implements a crosslink restlength of zero and shows a characteristic stiffening with divergence of force as u / l 0 . Eqn (1) can be integrated to yield the energy † (up to a constant) Imposing affine deformations on the lament level fully determines the deformation eld u on the sublament level. In the following analysis, we consider a single representative lament subject to an extensional strain of the surrounding medium that it is embedded in and crosslinked to.

1D network calculation
We start with a one dimensional system, i.e., 1D extensional strain 3, and advance in dimensionality by converting an applied shear strain g into the orientation dependent extensional strain 3(g) felt by the lament.
In the rest frame of the lament, the end-to-end distance of a crosslink at distance x from the center of the lament is given by |u(x,3)| ¼ |3x| (see Fig. 1(a)). For notational convenience, we consider positive 3 only. Under the assumption that the crosslink density is high enough that one can consider the associated distribution as uniformly continuous, the total energy of a lament in 1D is given by Substituting eqn (2) into eqn (3), this expression can be integrated analytically (see Appendix A.1). Following the described approach for the linear regime of the WLC force-extension relation, i.e., for u ( l 0 , the linear modulus may be extracted as where E/V is the energy per unit volume V in the network and 3 is a small strain. 27 For a 1D system this yields G aff 0 ¼ k B T l p l 0 being the † More precisely, it is a free energy, which includes both energetic (bending) and entropic terms for the crosslinks (not for the laments).
linear spring constant of a crosslink and r: ¼ NL/V the total length of laments per unit volume. The same holds for the modulus in 2D and 3D, but with different numerical prefactors: 1/96 and 1/192, respectively. 16,19,26 Next we show that one can extract a functional relation between the nonlinear modulus and stress in the nonlinear regime, based on simple asymptotic scaling analysis. It follows from above that there is a strain 3 d : ¼ l 0 /(L/2) at which the outer most crosslink (at x ¼ L/2) reaches maximum extension. For 3 / 3 d the energy diverges as follows with '$' dened via E $ f 5 E/f / const. The upper index 'div' always indicates that we are only taking into account the diverging part of the 1D lament energy. We express stress and differential elastic modulus via s ¼ 2 , respectively, in order to obtain s 1D $ 1/(1 À 3/3 d ), and K 1D $ 1/ (1 À 3/3 d ) 2 . We arrive at the asymptotic scaling relation This scaling relation between the modulus and stress in 1D has also been derived in previous work. 19 Next we consider scaling relations in 2D and 3D.

2D network calculation
We perform similar calculations as in 1D, while taking into account that the extensional strain 3, which results from a shear strain g on a 2D system, depends on the orientation of the lament under consideration. In the small-strain limit one thus obtains where q˛[0, p] is the angle between the lament and the shear direction (see Fig. 1(b)). Substituting this expression into eqn (4) and averaging over all orientations leads to ðg=2Þ sin 2q dq; where we assume g $ 0 for notational convenience; the upper integration limit is reduced to p/2 because |sin 2q| is p/2-periodic. Note that we do not take into account a redistribution of lament orientations under the shear transformation. This approach, as well as the small-strain approximation for 3(g, q), are justied if L [ l 0 , since then the strain g d : ¼ 4l 0 /L at which the integrand diverges is small. Differentiating eqn (7) with respect to g and neglecting the weaker (logarithmically) diverging part of the integrand lead to an expression for the stress, as g / g d : ; The divergence of the stress is of the form s 2D $ 1/(1 À (g/g d )) 1/2 and hence K 2D $ 1/(1 À g/g d ) 3/2 . Therefore, the asymptotic scaling behavior for the differential modulus in two dimensions is given by Note the difference of the scaling relations to the ones in the 1D scenario. Stress shows a weaker divergence with strain than in 1D but a stronger dependence on the differential modulus.
Integration of the diverging part of the stress further shows that the energy at maximum strain is nite-in contrast to the 1D setting, where the energy diverges at the critical strain. This is an effect introduced by orientational averaging only.

3D network calculation
For a 3D network, the extensional strain on a lament in the small-strain limit is given by Fig. 1 Sketch of the assumptions of the affine theory: (a) 1D: a filament (green) of length L is connected to its surrounding through n crosslinks (blue) that have zero extension at zero strain. The surrounding of the filament is subject to a uniform extensional strain 3 ¼ DL/L. Since the filament itself is assumed to be perfectly rigid, all deformation goes into the crosslinks (drawn in the y-direction for better visualization). The deformation of a crosslink at distance x from the center of the filament is given by u ¼ 3x (deformation field depicted by the horizontal gray arrows). (b) For a 2D system, the extensional strain on a filament at angle q with the axis in the shear direction is given by 3 z (g/2) sin 2q, for a small shear strain g ¼ Dx=h ¼ tan w.
in the usual spherical coordinates. In direct analogy to the 2D case (see eqn (9)), the averaged stress close to g d ¼ 4l 0 /L can be written as with g $ 0; the upper integration limit for the f integration is reduced to p/2 because |cos f| is p-periodic and symmetric about p/2 on [0,p]. If we carry out the f integral and expand the integrand around q ¼ p/4, in order to integrate over q (see Appendix A.2 for details), we obtain s 3D $ Àln(1 À g/g d ) and hence K $ 1/(1 À g/g d ).
Consequently, K does not scale with s as a power law; instead, one obtains with a real constant c. The absence of asymptotic power law scaling sets 3D networks apart from 1D and 2D networks. In 3D, we observe the weakest (logarithmic) divergence of stress with strain. Integrating the diverging part of the stress shows that the energy again remains nite for g / g d . Finite crosslink density. By considering the limit of innite crosslink density, we have derived theoretical scaling relations for strain stiffening by integrating along a lament's backbone (see eqn (3)). For any real system, however, the crosslink density is nite and eqn (3) turns into a sum where {x i } are the crosslink binding sites along the lament. Fig. 2 shows numerical results for the behavior of the corresponding differential modulus K for nite n, obtained by numerical evaluation of eqn (14) and appropriate orientational averaging. Note that the asymptotic scaling behavior of K in the limit of innite crosslink density inuences a nite network's behavior in the intermediate-stress regime (see the inset of Fig. 2); however, near the critical strain, the differential modulus scales as K $ s 3/2 , i.e., like the response of a single WLC. Furthermore, for 1D and 2D systems the theoretical scaling exponents in the limit of innite crosslink densities can (in the intermediate regime) indeed be approached by increasing n. In contrast, as shown above, in 3D the theoretically derived scaling of K is exponential in s. Such an exponential increase is quantied by an (in principle) indenitely increasing maximal slope with increasing n in the ln K versus ln s plots; e.g., for n ¼ 60 the maximal slope is 3.49, for n ¼ 3000 it is 5.82. However, for any nite n, eventually there is always universal scaling of K $ s 3/2 , resulting from the single WLC force-extension relation, independent of the dimensionality of the system. Indeed, for any given n, the integral representation eqn (3) becomes invalid close to g ¼ g d due to the divergence of the WLC energy. The numerical results in Fig. 2 have been obtained without the small-strain approximation for the extension of the laments. However, the redistribution of the lament orientations under shear has not been taken into account in Fig. 2. Calculations including this effect show that it may both decrease and increase the maximum intermediate slope in the ln K versus ln s plot and shi the peak to larger stress values depending on the maximum strain g d . In any case, the asymptotic scaling regime remains unchanged.
In the next section we introduce the simulation framework that we use to study 3D networks consisting of many laments and crosslinks, relaxing the assumption of affine deformations.

Simulation model
We perform quasistatic simulations of 3D networks that consist of N rigid laments of length L, permanently crosslinked by a collection of nN/2 crosslinks of length l 0 . All lengths are measured in units of the side length of the cubic periodic simulation box. A typical set of parameters is N ¼ 3000, Each lament is modeled as perfectly rigid, implying that its conguration can be described by its two endpoints only, which are constraint to stay at distance L. The exible crosslinks are modeled as a central force acting between the two binding sites. In particular, we use the WLC interpolation formula (eqn (1)) and the corresponding energy (eqn (2)). In all data that are presented, forces are measured in units of (k B T)/l p . There are no additional degrees of freedom introduced through the crosslinks, since their conguration is represented via the endpoints of the laments, in terms of barycentric coordinates.
In order to generate an initial network conguration we proceed as follows. We generate N randomly distributed laments by rst randomly choosing their centers of mass in our simulation box and by then picking a random orientation for each lament. For crosslinking we apply the following iterative procedure. We randomly choose two points on distinct laments and insert a crosslink if the corresponding point-to-point distance is shorter than a certain threshold al 0 . Here a˛[0,1) serves as a parameter to vary the initially allowed crosslink lengths in the system. This procedure is repeated until the desired number of crosslinks is reached; see Fig. 3 for an illustration of the nal conguration. Since we perform quasistatic simulations, the system must be at static equilibrium at all times. As practically all crosslinks will be stretched beyond their rest-length aer the initial network generation, we minimize the energy (of the crosslinks) before subjecting the simulation box to any deformation. ‡ For energy minimization we use the freely available external library IPOPT, 28 which requires the gradient and the Hessian of the system's energy function. It might happen during the optimization process that individual crosslinks reach extensions u larger than their contour length l 0 . Acceptance of these solutions is prohibited by setting the energy to innity (10 19 ) for u $ l 0 in eqn (2); without this modication it would become negative in that regime. The length constraints for the laments are realized via Lagrange multipliers.
In order to extract elastic properties of the network we perform quasistatic shearing by applying an affine incremental shear strain dg to the network, with subsequent rescaling of laments to length L (see Fig. 1). We apply Lees-Edwards shearing periodic boundary conditions. 29 The magnitude of dg is determined by calculating the maximum affine shear that leaves all crosslinks below their contour length. Due to the rescaling of lament lengths, a nonaffine deformation component is introduced. This nonaffinity may lead to crosslinks being overstretched aer all. In this case, we iteratively halve the shear strain until the length of all crosslinks remains below their contour length. Aer each shear increment, the energy is minimized. We apply a xed upper bound of 1% strain on dg in order to stay reasonably close to the previous solution. This increases numerical efficiency and accelerates convergence because it allows us to use a warm-start procedure that reuses Lagrange multipliers from one minimization as initial guesses for the next one. Moreover, the application of small shear steps reduces the likelihood of discontinuously jumping between local energy minima.
We stop shearing when the achievable increment in shear strain becomes smaller than a chosen threshold due to crosslinks that are very close to their maximum extension. During the entire simulation process, we record network parameters in the equilibrated states-in particular, the energy E as a function of shear strain g. This allows us to extract the shear stress as well as the differential shear elastic modulus Derivatives are taken by rst interpolating E(g) with a cubic spline. We dene the linear shear elastic modulus as In the following section we discuss the implications of our specic simulation model, in particular with respect to the network structure, and contrast it with previous studies that have been carried out mostly in 2D.

Initial tension and prestress
As mentioned in Section 3, our network generation results in a non-zero initial energy E 0 at zero strain. Indeed, by randomly placing (zero-diameter) laments in a 3D container, laments have zero probability to intersect; thus, crosslinks have nite initial extension with probability one. This is different from 2D, where randomly placed laments mutually intersect with a probability approaching one as their number increases. Indeed, so-called Mikado models, 19,[30][31][32] where laments are crosslinked at their intersection sites only, exhibit no forces at zero strain.
In contrast, the initial stretching of crosslinks in our networks results in an initial tension before any deformation. For a quantitative analysis we measure a global variant of this effect by what we call total prestress s 0 , which measures the Fig. 3 Example of an initially generated network that has not been relaxed into static equilibrium yet. Rigid filaments are shown in green, and flexible crosslinks in blue. Short crosslink or filament fragments correspond to filaments/crosslinks that cross the periodic boundaries of the simulation box. For the sake of visual appearance, the network is much sparser than the systems that are studied in the remainder of this article, and the ratio of filament to crosslink length is much smaller, N ¼ 300, n ¼ 10, normal stress § component orthogonal to the shear planes.{ More precisely, we measure the single sided (e.g., upward) normal component of the force that is acting on a given shear plane, by summing up the normal components of the forces exerted by each crosslink and lament passing through the given shear plane, see Fig. 4(a). The normal stress is then given by dividing by the surface area of the shear plane. Note that s 0 does not depend on the choice of a particular shear plane; indeed, if the total stress was changing during vertical movement of a shear plane, then this would immediately contradict force balance in the system. Intuitively, one might expect negative normal stresses (pulling down on the upper face of the simulation box), since crosslinks are contractile. However, since laments withstand compression, it is possible to construct systems that exhibit positive normal stress. This suggests the existence of congurations with zero normal stress.k Indeed, so-called tensegrity structures, 34 which are in static equilibrium in the absence of boundary conditions satisfy this criterion-while still being able to store arbitrary amounts of energy (see Fig. 4(b)). Empirically, our simulations show that the random networks generated by the procedure described in Section 3 exhibit negative initial normal stresses throughout. Their integrity is provided through the application of periodic boundary conditions. Note, in particular, that our setup enforces conservation of the volume of the simulation box. In general, it would be possible to relax the prestress by letting the volume of the simulation box change. However, we did not follow this approach in the study presented here, in order to ensure that the lament length remains signicantly smaller than the size of the simulation box.
In the following, we relate total prestress to the linear elastic response of our networks.

Linear regime
In previous work, 16,19,26 an expression for the linear modulus in 3D was derived under the assumption of affine deformations and in the absence of any initial tension in the network. Our simulations show that the linear elastic modulus depends on the initial tension in the network.
One scenario that clearly demonstrates the dependence of the linear modulus G 0 (dened in eqn (15)) on the initial tension is shown in Fig. 5 where the admissible maximum initial crosslink length was varied.
For a more quantitative analysis we have designed a method that allows us to change initial tension for a network with a xed set of simulation parameters. We rst randomly generate a network as described above and let it relax into static equilibrium. We then remove a given amount (5%) of the most-stretched crosslinks in the system. Then we reconnect those crosslinks randomly again, and let the network relax. This procedure is repeated N rel times. Thereby, we successively decrease the system's initial tension, and therefore also its total energy, see the inset of Fig. 6. Not only does the total energy decrease, we also observe a change in the distribution of forces (see Fig. 6). As long as one performs the crosslink binding-unbinding procedure over a small enough fraction of crosslinks, the network remains nearly isotropic. Fig. 4 (a) Measuring the total prestress s 0 by extracting the normal component of the total force acting on a shear plane. We sum up all the forces acting on one side of the plane exerted by (i) the crosslinks passing through (here f 2 and f 4 ) and (ii) the filaments passing through (here f 1 + f 3 )-then we project onto the normal vector n. (b) A tensegrity structure (here: Snelson's X 33 ) remains in static equilibrium without the application of boundary conditions. The forces acting on any plane add up to zero, i.e. no plane carries any total prestress although it is under tension locally. Fig. 5 Differential elastic modulus K as a function of strain g for different levels of initial tension. The initial tension in the network is varied by changing the initially admissible maximal crosslink length al 0 . The linear modulus G 0 ¼ K| g¼0 increases with the initial tension in the network (initial tension increases with a). It is also evident that the divergence of K occurs at a strain g d that decreases with increasing a. Here: N ¼ 3000, n ¼ 60, L ¼ 0.3, l 0 ¼ 0.03. § Note that our notion of prestress is not to be confused with the constant prestress externally applied in bulk rheology experiments, which is a shear stress in general. { Although we could in principle dene total prestress as the normal component of the stress acting on any plane in our system we prefer to use the shear plane as this simplies the forthcoming analysis. k Note that individual crosslinks are still under tension; however, the total normal force acting on the shear plane vanishes.
It is apparent from the inset of Fig. 7 that the linear elastic modulus is reduced by increasing the number of relaxation steps, as expected. Fig. 7 also shows the dependence of linear modulus G 0 on the total prestress s 0 , which has been introduced in Section 4. We varied s 0 via the above described procedure, and measured G 0 with the shearing protocol described in Section 3. Aer a certain number of relaxation steps the empirical value for G 0 equals the value G aff 0 expected from affine theory (see Section 2.1). Relaxing initial tension further, we reach moduli even below G aff 0 . This is possible because the network can rearrange nonaffinely, thereby soening its response. Over a certain range of total prestresses, we observe linear scaling of G 0 with s 0 , a phenomenon, which has been discussed in other contexts before (see for example ref. 35). We explain the linear regime as follows. For small strains the normal component s N of the stress acting on shear planes is close in magnitude to the total prestress s 0 , i.e., s N z s 0 . For small strains given by shear angles w z 0, total forces acting on the shear planes make an angle 4 with the direction normal to the shear planes (see Fig. 8). Our simulations show that tan 4 f tan w and that the constant of proportionality remains unchanged in the linear scaling regime. Therefore, shear satises where s S is the component of the stress acting on shear planes in the shear direction, see Fig. 8. Hence, the linear elastic shear modulus G 0 dened via s S ¼ G 0 g is proportional to the total prestress s 0 via eqn (16). However, for very small total prestresses, i.e., aer many relaxation steps, the modulus shows a steeper than linear dependence on s 0 . Indeed, in this regime the aforementioned constant of proportionality becomes larger.
This effect might be attributed to the fact that for small s 0 , tensegrity type elements (see Fig. 4(b)), which do not contribute to the total prestress but carry energy, contribute signicantly to the measured shear stress, thereby increasing 4 (see Fig. 8). Furthermore, affine theory predicts linear scaling of the modulus G 0 with crosslink density n. Fig. 9 shows that this linear scaling is indeed reproduced in our simulations, independent of the prestress. Moreover, by changing the prestress via our relaxation procedure it is possible to reach comparable slopes to what is predicted by the affine theory.
The next section deals with the nonlinear elastic response of the simulated networks, and relates it to the theoretical results that were derived in Section 2. Fig. 6 Distribution of forces in crosslinks for a system without or with N rel ¼ 100 relaxation steps. The relaxation procedure cuts the large force tail of the initial distribution and establishes a sharper peak at small forces. The inset shows the total energy E in the system, normalized by the initial energy E 0 , as a function of the number of relaxation steps N rel . Fig. 7 Linear elastic modulus G 0 normalized by the affine prediction G aff 0 as a function of total prestress s 0 normalized by the total prestress s * 0 immediately after initial network generation. The total prestress is reduced via the procedure described in Section 5. For small total prestress, G 0 exhibits superlinear dependence on s 0 . Up to s 0 ¼ s * 0 , we observe linear scaling G 0 f s 0 , as predicted by the model. The straight line is drawn as a guide to the eye, representing linear scaling. Parameters: N ¼ 3000, n ¼ 60, L ¼ 0.3, l 0 ¼ 0.06, and a ¼ 0.5. The inset shows differential elastic modulus K versus shear strain g for systems with the varying number of relaxation steps N rel˛{ 0, 50, 100, 150}. G 0 goes down with increasing N rel . Parameters: N ¼ 3000, n ¼ 60, L ¼ 0.3, l 0 ¼ 0.03, and a ¼ 0.5. Fig. 8 The initial network carries a total prestress s 0 . After a small shear g ¼ tan w has been applied it exhibits shear stress s S and normal stress s N , with tan 4 ¼ s S /s N .

Critical strain
The networks that we study are inherently nonlinear because crosslinks are WLCs with nite length l 0 (see eqn (1)), resulting in pronounced strain stiffening at a critical strain g c . Stress diverges at a higher strain g d . In our simulations, we dene the critical strain g c to be the strain where K/G 0 z 3. In the affine theory, g d and g c scale linearly with the ratio of crosslink to lament length l 0 /L. In our simulations, we cannot conclusively report on this dependence because the accessible ranges for l 0 and L are quite limited. On the one hand, there exists an upper limit for L (therefore also for l 0 , since l 0 /L ( 1 should hold) to be signicantly smaller than the simulation box. On the other hand, L and l 0 are bounded from below due to computational limitations-this is because we need to increase the number of laments in order to keep networks homogenous.
For ranges that are accessible to our simulations, we obtain the following results. If we x l 0 , then we observe linear scaling g c f 1/L for systems where no relaxation procedure has been applied (see Fig. 10(a)). Relaxed systems, however, sometimes show a less than linear dependence. This effect might be due to anisotropies induced by the relaxation procedure. If we x L, then the dependence of g c on l 0 is slightly less than linear (see Fig. 10(b)).

Differential modulus
It remains to discuss the dependence of the differential modulus on stress, the affine theory of which has been derived in Section 2. For nite crosslink densities, the only persistent scaling behavior is K $ s 3/2 , as g approaches g d -due to the fact that eventually the single WLC response dominates. In an intermediate regime, above the critical stress s c ¼ s(g c ), we observe slopes (d ln K/d ln s) > 3/2. The majority of the simulations shows intermediate slopes around 2 or slightly above, mostly independent of simulation parameters, but there are also realizations that show maximum slopes up to 3.5 (see Fig. 11). These higher slopes and the nal scaling K $ s 3/2 are in accordance with the predictions of affine theory. Indeed, a slope of 3.5 is the maximum slope predicted by the affine theory when using the same crosslink density as in the simulation (Fig. 2). There are, however, differences between theory and simulation in terms of slope proles since various assumptions are made by the theory that do not hold in the simulations: a randomly generated network does not have a uniform crosslink density along the laments, these systems are prestressed, and there is   no perfect isotropy. Moreover, the networks do not deform perfectly affinely.

Nonaffinity
In order to study to what extent simulation results deviate from affine theory, apart from prestress, nonuniform crosslink density, and anisotropy, we analyze the nonaffinity of the network deformation under shear. For a single lament, we dene its differential nonaffinity with respect to the center of mass by where dr aff and dr are the 3D coordinates of lament's center of mass aer applying an incremental shear strain dg without and with relaxation, respectively. We let dG denote the average of the differential nonaffinities over all laments. Affine approximations imply dG ¼ 0. Fig. 12 shows that center of mass deformations are mostly affine for small strains. However, the differential nonaffinity increases starting at a strain around g c and eventually diverges as g / g d . This can be understood, since the networks are strain stiffening, such that small incremental strain can induce a large increase in the forces of individual crosslinks, thereby inducing large local rearrangements during energy minimization.
While increasing shear strain, there are force chains 36-38 developing in the network, which carry most of the tension, and which cannot reduce their strain due to the fact that they span the entire system (see the inset of Fig. 13). We quantify this effect by considering tension proles along laments. The tension s at position x along a lament is given via sðxÞ¼ X |xi| . |x| f cl ðu i Þ, where {x i } are the crosslink binding sites and {u i } are their extensions (u i ¼ 3x i in affine theory). Fig. 13 shows tension proles averaged over all laments for both, theoretical and simulated systems at various strains. In the simulations there is non-zero tension at zero strain due to prestress. With increasing g, the simulations resemble the proles expected from affine theory. However, when approaching the maximum strain g d , the emergence of selective paths (force chains) that carry most of the tension becomes evident. The highly stretched crosslinks dominate the averaged tension proles and therefore lead to jumps in the tension curves at the respective binding sites along the lament (green solid curve in Fig. 13).

Bending
Thus far we have restricted our theory and simulations to rigid laments that can neither bend nor stretch. In ref. 19, the authors considered nite stretching compliance of laments, while bending compliance was assumed to be zero. They report that nite stretching stiffness does not impact the nonlinear stiffening regime of a composite network apart from the expected convergence (to some constant value) of the modulus at high strains. Here we complement this analysis by considering laments that have nite bending but no stretching compliance. We performed simulations on a 2D network because of the relative computational ease compared to the 3D case. In addition to the energy stored in the crosslinks, we consider bending energy of the form E b ¼ kq 2 /(2l av ), where k is the bending rigidity, q is the angle through which the laments bend locally, and l av ¼ (l 1 + l 2 )/2 is the average distance between two adjacent pairs of crosslinks (Fig. 14). We show the results in Fig. 12 Differential nonaffinity dG as a function of scaled shear strain g/ g c for a system with N ¼ 3000, n ¼ 60, L ¼ 0.3, l 0 ¼ 0.06, and a ¼ 0.05. Fig. 13 Average tension s as a function of position x along the filament for various strain values. Tension s is normalized by its maximum absolute value s 0 . Dashed curves correspond to theoretical results for n ¼ 60 at g ¼ g c (blue), g x g d (green). Solid curves show simulation data, with N ¼ 3000, n ¼ 60, L ¼ 0.3, l 0 ¼ 0.06, and a ¼ 0.5. The inset shows a snapshot of the same system at maximum strain g d x 0.6 where only the 15 most stretched crosslinks and the corresponding filaments are shown. They form singular paths that span the whole system, thereby preventing further stress reduction via nonaffine rearrangements in these finite systems. Fig. 15. The range of bending rigidity was chosen such that the linear modulus was still determined by the so stretching modes of the crosslinks, so that bending did not impact the linear regime. As can be seen from these plots, bending compliance does not impact the nonlinear stiffening regime either-since bending modes are geometrically prohibited for large strains.
Thus, in isolation, neither bending nor stretching compliance of laments impacts the nonlinear stiffening regime of composite networks. These ndings suggest that the theoretical models at present cannot explain the K $ s scaling observed in experiments.

Conclusions
We have studied the elastic properties of composite crosslinked lamentous networks in 3D analytically and numerically. We modeled such networks as a collection of rigid laments connected by WLC crosslinks.
Based on the affine theory introduced in ref. 19 we derived asymptotic power law scaling exponents for the differential elastic modulus with stress, in the limit of innite crosslink density. In this case, the scaling exponents depend on the dimensionality of the system. In particular, we showed that 3D systems no longer exhibit a power law. Furthermore, we showed that for nite crosslink densities, the only persistent regime (over several orders of magnitude of stress) is the s 3/2 scaling, as it is derived from the single WLC force-extension relation, eqn (1). This is in sharp contrast with the model proposed in ref. 16 and 26, where linear scaling was suggested, independent of the dimensionality of the system. This model implies nite stress at any strain and therefore does not apply to composite networks of rigid laments with exible crosslinks of nite length.
We further developed a simulation framework that allows us to measure the elastic response of random lamentous networks with WLC crosslinks. One important property of these 3D networks is that, by construction, they are prestressed due to initial extensions of the crosslinks. In addition to geometrical constraints, active elements such as motors can induce prestress as well. 39 We showed that the prestress in a network can dominate the linear response and might therefore be a feature that is worthwhile analyzing in experimental systems.
Regarding the nonlinear response, we observed divergence of stress (and differential modulus) at nite strain. Close to this strain we measured power law scaling of the differential modulus with stress, with an exponent 3/2, just as expected for a single WLC. In an intermediate-stress regime we observed local exponents that span the entire range of theoretically derived values for systems of differing dimensionality. The fact that our simulation results do not always resemble the predictions of a 3D affine theory, in this intermediate regime, may be attributed to nonaffine deformations. Extracting the exact set of assumptions-such as uniform crosslink density, isotropy, or zero prestress-that are responsible for these discrepancies is le for future investigation.
Experiments (see, e.g., ref. 13, 20, and 21) have shown that in the nonlinear regime the differential modulus scales approximately linearly with the shear stress. We did not nd such a regime in our simulations-neither when working with rigid laments nor when incorporating nite bending stiffness (or enthalpic stretching as done in ref. 19). Therefore, we argue that none of the currently available theories can adequately explain the linear scaling of the differential modulus observed experimentally. It could possibly be that the WLC model does not accurately describe the elastic response of a single crosslink throughout the whole experimentally accessible regime. We speculate, however, that the linear scaling might be due to thermal uctuations of the laments, which have not been considered so far.
Fig. 14 Sketch of the local bending geometry of a filament (green) with crosslinks attached (blue). The local bending energy is given by E b ¼ kq 2 /(2l av ), with k being the bending rigidity and l av ¼ (l 1 + l 2 )/2. Fig. 15 (a) Differential modulus K as a function of shear stress s, rescaled by linear modulus G 0 and critical stress s c ¼ s(g c ), respectively, for various bending rigidities k. The solid straight line indicates power law scaling K $ s 3/2 . (b) Differential modulus K as a function of shear strain g, rescaled by linear modulus G 0 and critical strain g c , respectively. Parameters: N ¼ 800, L ¼ 1, l 0 ¼ 0.1, and system-size L x ¼ L y ¼ 6.