Nonlinear elasticity of semiflexible filament networks

We develop a continuum theory for equilibrium elasticity of a network of crosslinked semiflexible filaments, spanning the full range between flexible entropy-driven chains to stiff athermal rods. We choose the 3-chain constitutive model over several plausible candidates, and derive analytical expressions for the elastic energy at arbitrary strain, with the corresponding stress-strain relationship. The theory fits well to a wide range of experimental data on simple shear in different filament networks, quantitatively matching the differential shear modulus variation with stress, with only two adjustable parameters (which represent the filament stiffness and pre-tension in the network, respectively). The general theory accurately describes the crossover between the positive and negative Poynting effect (normal stress on imposed shear) on increasing the stiffness of filaments forming the network.


Introduction
Networks of semiflexible filaments and fibers are common in biological systems, where dynamic structures of tunable strength, elasticity, and adjustable response are required. [1][2][3] From the microscopic scale of cytoskeleton to the macroscopic scale of fibrous tissue, such networks utilize stiff or semiflexible filaments, allowing for rich mechanical response, 4 dynamic remodelling, 5,6 and controlled structural failure, 7 while remaining an open structure allowing easy through access for solvents and solutes. 8 Inspired by these observations, filament networks are produced in industry, and now stand as a promising area in manufacturing functional materials. 9,10 The physics of a single semiflexible chain or filament is well understood with the aid of worm-like chain model, in various implementations and approximations. [11][12][13][14][15] The full theory is capable of accurately describing the force-extension relationship, as well as the magnitude of transverse fluctuations, for a range of filament stiffness spanning from very high (a rigid elastic rod) to very flexible (a classical polymer coil) and a range of end-to-end extensions from zero up to the full stretch where the entropic force has a characteristic divergence. When such filaments form a macroscopic crosslinked network, the main scaling features of its nonlinear mechanical response are still determined by the single filament properties, as reviewed in ref. 16 and 17. There have been many key contributions to the elastic theory of stiff and semiflexible networks, including ones highlighting the issues of strain non-affinity. 18,19 In order to obtain the constitutive relationship of a semiflexible network, a widely used approach follows the following procedure (explained in detail in chapters 3 and 4 of ref. 20): the shear stress on ij plane is calculated by summing the contributions of the tension along the i axis of all the chains crossing the j plane. This approach has been very successful in deriving viscoelastic properties of polymer solutions and melts under shear flow; 20,21 Storm et al. 22 have applied the same stress-strain relation to a crosslinked network by assuming the affine elastic strain acting on each filament and the probability distribution of semiflexible filament length of Wilhelm and Frey. 23 This model qualitatively describes how the network behaves when being deformed; however, it can only be used numerically (since no closed expressions are possible) and it omits the pre-stress acting on the chains, as in this model the tension is assumed as zero when there is no deformation. Wilhelm and Frey also produced a numerical simulation of filament network elasticity on their own, 24 using the Mikado model of connectivity and athermal rigid rods as elastic elements. Although important issues of percolation rigidity threshold are exposed, this work cannot be used to describe most experimentally relevant (and most biological) filaments. More recently, Palmer and Boyce 25 proposed a closed analytical form of elastic free energy, with the corresponding constitutive relation, using the same force-extension filament relationship as Storm et al. They applied the '8-chain model' originally introduced in the context of ordinary rubber elasticity 25 with an approximate expression for individual filament elasticity applied for each strand. This might be the best attempt in formulating the nonlinear elasticity of semiflexible network to date. On the other hand, Unterberger et al. 26 developed a '1-chain' model for a network, by assuming filaments have a homogenous orientational distribution in the equilibrium system, and the average (imposed) stretch of the network l is a p-root average of deformations of individual filaments: where l*(y,f) is the deformation of an individual filament, p is the averaging parameter and A the area of the unit sphere. This procedure allows the local stretch to fluctuate around its average value, and may partly account for the non-affinity. This model could only be applied numerically to reproduce several key mechanical features (shear and normal stress). Most good models discussed here make a successful fitting of shear stress data of different filament networks, which captures the generic stress-stiffening effect originating from the characteristic divergence of the force-extension curve of an individual filament. Usually, the authors employ the versions of celebrated Marco-Siggia interpolation formula, 12 which gives the correct response for filaments near full-extension (but much less so for more flexible filaments). In fact, the well-known MacKintosh scaling of differential shear modulus with stress, K p s 3/2 , in the stress-stiffening regime is entirely based on the mentioned divergence of an individual filament near full extension 27 (it holds even for an athermal network of undulating filaments, as long as the tensile force is proportional to the inverse-square of the compression, as is the case near full extension 28,29 ). Therefore, a mere agreement (good fitting) of shear stressstrain curves is not a sufficient test of different theories. In particular, they must simultaneously describe the effect of negative normal stress in the network of stiffer filaments, with a positive normal stress (also known as the Weissenberg effect) for networks made of more flexible chains. 30,31 It turns out that none of the mentioned theories, including the 8-chain model of Palmer and Boyce, are able to produce the correct normal stress, or address the network stability increasing with filament pre-tension. 32 In this paper we develop and discuss a continuum theory of semiflexible network (in closed analytical form) that addresses these issues, while retaining accurate fitting of a wide range of shear stress data.

Network theory
Before introducing free energy of a semiflexible network, we need to discuss the properties of a single semiflexible filament. A filament connecting two neighboring crosslinks in a network is sketched in Fig. 1(a), with coordinates r(s), where 0 r s r L c is the arc-length coordinate along the chain. Different approaches to the chain (in)extensibility have been tested over the years, 33 from the strict constraint to the requirement that the length of filament remains constant on average (while small local fluctuations are allowed), 15,23 to the models that explicitly include filament stretching. 22,34 It turns out that in the regime of high extension, when there are no foldbacks (hairpins) on the chain, all length-constraining models give the same divergence of the force-extension, f p 1/(1 À x) 2 , with the end-to-end ratio x = x/L c , where x is the end-to-end length of the filament, Fig. 1(a). This was used by the famous interpolation formula of Marko and Siggia. 12 More recently, a complete theory of filament entropic elasticity has been developed, which spans the full range of extensions and the full range of bending modulus. 15 In the worm-like chain model, 11,35 the bending energy can be expressed by 1 2 k Ð Lc 0 ds @ 2 rðsÞ @s 2 2 , where k is the bending rigidity. The key physical quantity for describing the stiffness of a polymer chain is the persistence length l p = k/k B T. When the contour length of the filament, L c , is comparable with its persistence length l p , the chain is regarded as ''semiflexible''.
Combining the effects of enthalpy arising from bending and entropy of conformation fluctuations, the closed form of the single chain free energy 15 can be expressed as a function of its end-to-end factor, x = x/L c : where c = k/2k B TL c is a dimensionless stiffness parameter reflecting the competition between bending and thermal energy; in our notation c = l p /2L c . In ref. 15 one can find the comparison with several notable models of semiflexible filament (Marko-Siggia and Ha-Thirumalai) and where they deviate from the accurate expression (1). As shown in Fig. 1(b), if the value of c is smaller than a critical value c* = p À3/2 E 0.18, the minimum of the free energy (or the force-free natural length) will be at x = 0; such a chain can be regarded as flexible. When c { c*, one recovers the Gaussian entropic spring form: F chain E (2k B T/pl p L c )x 2 . The Marko-Siggia (in fact, Fixman-Kovac) limit commonly used for quick fitting of force-extension curves is reached for flexible chains at high extension c { c*, x -1: F chain E (k B T/pl p L c )/(1 À x). When c c c*, the chain can be regarded as a stiff rod with a natural , and its energy is dominated by elastic bending: F chain E (p 2 k/L c ) [1 À x]. Under compression x o x 0 such a filament undergoes Euler buckling instability. Note that on forming a crosslinked network, x does not need to be equal to x 0 for each strand in equilibrium, meaning there can be pre-tension in the network. Fig. 1 (a) A sketch of semiflexible network mesh unit, with a bent filament connecting two crosslinks with its contour length L c , and end-to-end length x being the mesh size. (b) The relationship 15 between the tensile force, f (x)/k B T, and end-to-end factor, x = x/L c , plotted for semiflexible chains with different stiffness c, and marking the position x 0 for a chain to stay at the force-free state.
We now construct the continuum elastic free energy of a network of such filaments using the methodology that was successfully developed in rubber elasticity. 36,37 If the sample shape is changed with a deformation tensor E, then the corresponding Cauchy-Green tensor 38 is C = EE T , with eigenvalues: l 1 2 , l 2 2 and l 3 2 . The values l 1 , l 2 , l 3 can be interpreted as stretching ratios along the principal directions of deformation. A class of simple, yet powerful theories of rubber elasticity is formulated in terms of invariants of the C-G tensor, . Famous examples of this class are neo-Hookean (Gaussian), Mooney-Rivlin, and Gent models (see ref. 39 for review). Note that polymer networks are frequently treated as incompressible, so I 3 = 1.
For a polymer network, the chains are distributed and crosslinked randomly throughout the material, which is the main difficulty in obtaining an analytical expression of the total free energy of the network (unless the individual chain is Gaussian). To incorporate a more complicated chain, such as eqn (1), into a constitutive framework, it is necessary to have a model that relates the chain deformation to the applied affine strain. This can be accomplished by representing an element of volume in the average network as a ''mesh cell'', which is then symmetrically multiplied to fill the volume. Here we need to distinguish ''affinity'' and ''non-affinity'' in a filament network constructed by the unit-cell method. A ''mesh cell'' is deformed differently from the material, as these cells are constructed in orientation aligned with the principal stretching direction of the material, which is referred to as ''non-affinity'' by Palmer and Boyce. 25 All mesh cells deform in the same way, with repeated deformed structure. However, the real non-affinity can arise from the different responses among different ''mesh cells'', due to local force relaxation. In this work we refer to the concepts of ''affinity'' or ''non-affinity'' in the local context of individual filament junctions and mesh cells, rather than the global one between the orientation of a mesh cell and the macroscopic deformation of the material. In assuming a symmetry of repeated mesh cells, we automatically discard the effects of the local non-affine deformations, 40,41 i.e. take all mesh cells deforming uniformly, which we know must be the case at least for stiff filaments. 18,19 Since there is no quantitative way of assessing the degree of the error introduced by the affine approximation (in our interpretation), we will have to look at the fits to the experimental data for validation. Within a mesh cell, chains or crosslinks have several possible arrangements, reflecting what one assumes about the topology of the network mesh, see In 1-chain model, 15,26,36 one end of a polymer chain is fixed at the center of a sphere, while the other end is on the sphere surface at an arbitrary orientation (y,j), distributed isotropically. When deformed by E, with stretching ratios along principal directions l 1 , l 2 , and l 3 , the lengths of the three semi-axes will change from x, to xl 1 , xl 2 and xl 3 , respectively. The elastic energy density of the network can be expressed as an orientational average of a deformed filament: where n is the density of crosslinked chains andlðy; jÞ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi sin 2 y cos 2 jl 1 2 þ sin 2 jl 2 2 À Á þ cos 2 yl 2 3 q . If one takes the 'entropic spring' limit of the Gaussian chain, the average F 1c reduces to the classical neo-Hookean rubber-elastic expression However, in the general case the free energy in 1-chain model cannot be expressed in analytical form as a function of strain invariants, which renders it less convenient. In the following, we compare the 3-and 8-chain models, and decide on the 3-chain model preference, in particular due to the failure of 8-chain model in reproducing the normal stress.
In the 3-chain model, a primitive cubic is constructed with lattice points representing the crosslinking sites, and the edges are aligned along the principle directions of deformation tensor E. Three chains are linked with their end-to-end vectors along the edges and the equilibrium mesh size x. On deformation, the lengths of three perpendicular edges at one lattice point become l 1 x, l 2 x and l 3 x, respectively. Then the free energy density of a semiflexible network can be expressed as Eqn (3) can be rearranged as a function of the strain invariants: where c = l p /2L c and x = x/L c as used in eqn (1). The body-centered cubic cell of 8-chain model is constructed with eight filaments connected from the center point to all eight lattice points. 25 The edges of the cell are x as we define the mesh size, while the chains are shorter by factor ffiffi ffi 3 p =2. The important feature of the high-symmetry 8-chain model is that on deformation all chains change their distance by exactly the same amount, from ffiffi ffi 3 p x=2 to simply x ffiffiffiffi I 1 p 2, since the mesh cell is aligned along the principal axes of deformation tensor E. The free energy density of the network in this case is given by the single-chain expression directly: When the elastic energy is expressed as a function of strain invariants, the stress tensor of an incompressible material (with I 3 = 1) can be obtained as: 38,39 where C ij is the Cauchy strain, and P the Lagrangian multiplier for incompressibility, the value of which determined by the boundary conditions.

Simple shear deformation
Let us consider the simple shear deformation such that in Cartesian coordinates a point (x, y, z) in an original material will change to (x + gz,y,z) after being sheared; g is the shear strain. The incompressibility is satisfied automatically, and the remaining strain invariants are: I 1 = I 2 = 3 + g 2 . The shear stress in 3-chain model can be obtained from the general constitutive relation (6) as: while the corresponding expression in 8-chain model is: One must distinguish the nominal shear modulus G(g) = s xz (g)/g, from the differential shear modulus K(g) = qs xz /qg, 17,22 while the linear shear modulus of the network G 0 is equal to the differential modulus K 0 at g -0.
Take actin and fibrin network under simple shear deformation as an example. Both 3-chain and 8-chain models fit the shear-experiment data from ref. 22 and 27 equally well, in factperfectly, as shown in Fig. 3(a). The stiffness, c, and the initial end-to-end factor, x = x/L c , obtained by fitting with 3-chain model are a bit smaller than those in 8-chain model, but both in the semiflexible regime. This is because the eight chains in a BCC are stretched equally (in principal axes) and share the deformation, while the three chains in a PC are stretched differently and the most stretched one contributes the most to the nonlinear elastic energy. We believe the intrinsic heterogeneity in the 3-chain model, and the fact that the chain lying along the maximum principle stretch direction dominates the response of the whole cubic in semiflexible networks, is closer to the realistic case of filament network. In addition, we shall see below that the 8-chain model cannot explain negative normal stress. Hence we apply the 3-chain model in the following. Fig. 3(b) shows fits to experimental data for a wide variety of semiflexible filaments. In this figure we plot the shear modulus G(g) instead of stress, because the fitting convergence is much better when the initial data section is close to 1. Fitted parameters (c,x) are listed in Table 1 below, along with other parameters for each material that we list from the literature.
First of all, one may be surprised that the persistence length of collagen fibers is quoted as 20 mm, when there is a large body of literature that would claim that collagen fibers are stiff athermal rods with persistence length of centimeters. This is all to do with the way a sample is prepared, and we use/quote the data 42 where the collagen was apparently less aggregated than in a typical extra-cellular matrix. The same ambiguity will apply to actin networks as well, below. Another point to note about the fitted values in Table 1 is about neurofilament, which has c o c*, that is, rather flexible chains. Taking the fitted value for x = 0.15 = x/L c , this means that L c in this network was B1.6 mm, i.e. about 3.6 times longer than l p . By calculating for different filaments listed in Table 1, we see that the fitted x is smaller than x 0 , indicating all of the filaments in Fig. 3(b) are pre-compressed, rather than pre-stretched in the equilibrium state of the network.
As Table 1 shows, no matter what the effective stiffness of the examined biofilaments, the ratio of fitted parameters c/x = l p /2x remains close to 1. Later, when we discuss the network stability (Fig. 6), it will become clear that all the networks we examined here lie very near the stability boundary. It is not clear to us whether the fact that the mesh size is close to the filament persistence length is an unintended result of different crosslinking density in experiments, 22 or is a relevant and universal biological feature.  22 for sheared actin, collagen, vimentin, fibrin and neurofilament networks are optimally obtained for the scaled ratio G(g)/G 0 , which has a limit of 1 at g -0; the fitting parameters for 3-chain model are given in Table 1. Table 1 Fitting parameters (c,x) for collagen, actin, vimentin, fibrin and neurofilament data obtained from in ref. 22. Also shown are the linear modulus G 0 extracted from the original data and used for scaling in Fig. 3, the literature values of l p , and the calculated mesh size x = l p (x/2c) It is clear that networks of biological filaments have a great variety even within the same substance. Fig. 4 shows the published data for in vitro crosslinked actin networks reported by different groups, all performing the simple-shear experiment (the data is digitized from references given in the plot). The stress-strain plotted in log-log format allows a clear identification of the linear regime s = G 0 g (the modulus varying between 95 Pa and 1 Pa for different sets), and the subsequent stiffening at higher shear. All curves in Fig. 4(a) are fitted by the same eqn (7) 46 As a result the network at 'cycle 7' has lower G 0 but higher stiffening.
To find the onset of non-linearity in our theory, the differential shear modulus K(g) can be expanded in powers of shear strain g: The crossover point when these two terms become comparable with each other, in a more stiff network with x -1, can be approximated as g t $ 1 À x 2 À Á ffiffi ffi 6 p . The stress at this point is i . This crossover stress s t increases with the temperature as (k B T) 2 /k, closely matching experimental observations. 47 At the end of range, when the shear strain approaches the point of divergence in stress s xz , (g -1/x À x in eqn (7)), the stress and the differential modulus can be approximated as: KðgÞ Both expression diverge due to the same vanishing denominator, while maintaining the obvious scaling relation K B s 3/2 , which has been reported by different theories and experiments. 27,29,[47][48][49][50] We see this limit exposed clearly in Fig. 4(b) for very different actin networks. In this plot, the data points are the same as in Fig. 4(a), and we also plot the predicted curves of K(s) obtained by differentiating eqn (7), using parameters G 0 , c and x from the fitting in Fig. 4(a), for each data set. One has to make a comment here, in the context of K B s 3/2 scaling and different experiments. In many cases, such as in Unterberger et al., 26 the actin network was crosslinked by HMM and apparently retains some transient activity, producing the stress-softening and plasticity at higher stress (this is also the case in ref. 46 before the actin filaments formed bundles, or with F-actin crosslinked by filamin). Of course, our theory is not intended to deal with network plasticity (we assumed all crosslinks permanent) and therefore we only retained the experimental data points in the early stress-stiffening regime to see the 3/2 scaling.
On the other hand, just by examining the actual data for the actin network of Storm et al., 22 one might conclude that it strongly and systematically deviates from the 3/2 scaling. In fact, the authors of ref. 22 develop their own theory invoking various additional factors (e.g. filament extensibility) to account for this deviation. However, our basic theory, assuming permanent crosslinks, bulk incompressibility and inextensible filaments, evidently fits both s(g) and K(s) data very well. The fact that the data (and the predicted curve) do not appear to follow the 3/2 scaling is due to the fact that, for these values of c and x, the final crossover to this characteristic scaling regime would occur at an even higher stress (at which point the network would probably not survive in practice).

Normal stress
Though simple shear deformation was helpful in this analysis, the important issue of normal stress remains controversial. This is mainly because of the uncertain boundary conditions, see ref. 51 for detail. Since the experiments reporting normal stress measurements are most commonly conducted in rotating cylindrical geometry of a standard rheometer, 30,31 we will consider this geometry and realistic boundary conditions to describe the response of the material in its normal direction when a shear is applied, as shown in Fig. 5(a). Suppose the height and the radius of the undeformed cylinder are h 0 and R 0 , respectively; on deformation they may become h = l h h 0 and R = l R R 0 . Incompressibility maintains l h l R 2 = 1. In the (r, y, z) coordinate system, after rotating the top plate by the angle Y, the coordinates change as: r ! r ffiffiffiffi ffi l h p ; y ! y þ Yz=h 0 ; z ! l h z: where g(r) = rY/h 0 denotes the shear strain, which is a function of the radial position in this parallel-plate geometry. The strain invariants become: I 1 = l h 2 + [2 + g 2 ]/l h and I 2 = 2l h + [1 + g 2 ]/l h 2 . Given the shear strain at the outermost surface, g 0 = g(R 0 ), the total free energy then becomes a function of the stretching ratio l h along the z axis, after integration over radius: The equilibrium l h can be obtained by minimizing this free energy. When the cylinder radius becomes smaller (or larger) upon shear, with its height becoming correspondingly larger (or smaller) -this phenomenon is called the positive (or negative) Poynting effect. 52,53 This geometric effect with stress-free top plate corresponds to the positive (negative) normal stress required to maintain the fixed plate separation in a more common rheometry experiment. 30,31 The 8-chain model fails in obtaining the correct normal stress in a sheared network: due to its core assumption that eight chains in a cell are identically deformed (stretched) upon shear, the elongation ratio of the network along the stretching direction is always larger than 1, making the normal stress always positive. Fig. 5(b) shows the results for equilibrium l h , presented as a response to an oscillating imposed shear g 0 , for two model materials with different filament stiffness and pre-tension. A flexible network (c = 0.1, x = 0.1) has a positive Poynting effect, or positive normal stress is required to counter the expansion along the height (in other setting, this is called the Weissenberg effect). In a flexible network (c o c*), when the mesh size is much smaller than the contour length of the subchains connecting the neighboring crosslinks, i.e., x 0 { 1, the entropic energy plays a main role. Though chains are stretched along the principal extension direction in the shear geometry illustrated in Fig. 5(a), the chains in other two principle directions are more likely being compressed, leading to the material contraction along the radial and circumferential directions. However, if the mesh size is comparable with the contour length of the subchains, 0 { x 0 o 1, the material can behave in a negative Poynting effect manner: the height of the cylindrical sample contracts, or negative normal stress is required to counter that and maintain the fixed height. This is because the force acting on the chains directed along the principal extension in Fig. 5(a) is close to a divergence if x 0 -1, and it causes less energy when the material is compressed in the longitudinal direction, rather than stretched. Similar reason works also for a network of more stiff filaments (see blue curve for c = 0.3 4 c*, x = 0.6).
We did not specifically calculate the normal stress here (this would be a cumbersome process involving the full tensor form of eqn (6), first fixing the pressure P from the condition of zero radial stress on the free outer surface). However, the magnitude of normal stress can be easily estimated from the linear relationship s h B 3G 0 (l h À 1), which uses the fact that the normal strain is quite small (i.e. the linear regime is justified) and the Young modulus is 3G 0 . Taking the fibrin values in Table 1 as an example, the normal stress is about 11.5 Pa under a shear strain of g 0 = 0.5. These are very close to the observations in experiments. 30,31 The magnitude of s h C 20 Pa for the same shear also accurately matches the actin results obtained in ref. 26. Fig. 6 gives the full 'phase diagram' of the stiffness-tension parameter space (c,x) with phase boundaries separating positive/negative normal stress regions. In order to generate this diagram, for each parameter set (c,x) we have calculated the value of l h by minimizing F(l h ;g 0 ), under the given shear. In this way the boundary separating the positive (l h 4 1, s h 4 0) and the negative (l h o 1, s h o 0) normal stress regions in Fig. 6 is calculated. There is also a weak dependence of this boundary on the magnitude of the applied shear strain, which is due to the inherent non-linearity of stress-strain response; however, we are not showing this in the figure to avoid clutter. Fig. 6 shows that a loose flexible network usually has a positive Poynting effect, while a network of more stiff filaments has a negative Poynting effect, especially when the filaments are crosslinked with increasing pre-tension.

Network stability
One can see from Fig. 3(b) and 4(b) that the linear regime when s = G 0 g persists for a different range of strain in different Fig. 5 (a) A cylinder-shaped sample under oscillating shear deformation; (b) the relationship between imposed oscillating shear strain g at the outer surface of the cylinder (black curve) and the normal strain l h À 1 for a flexible network (red curve) with c = 0.1, x = 0.1, and a more stiff network (blue curve) with c = 0.3, x = 0.6. Fig. 6 The phase diagram of the network response at different stiffness, c, and filament pre-tension, x, showing the boundaries of positive/negative Poynting effect and the boundary of network stability (the dashed line of the neutral filament, x 0 , was defined in Fig. 1). Three 'softer' filament networks from Table 1  filaments; the onset of hyperelastic regime is determined by how much pre-tension is in the filaments, i.e. how x compares with x 0 (c) from eqn (1). The linear shear modulus G 0 can be easily obtained from our theory: in 3-chain model it is given by the first term in eqn (9). The marginal rigidity condition G 0 Z 0 determines the strand pre-tension that is required for achieving a stable network. The stability criterion here is: Note that the expression in the right hand side is always greater than c* = p À3/2 as defined in Section 2. Flexible chains with c o c* are therefore always stable in the rubbery network, that is, G 0 4 0 always. On the other hand, stiff filaments have to be crosslinked with x/L c exceeding the pre-tension threshold given by eqn (14), which turns out to be slightly lower than x 0 (c) defined in Fig. 1(b). In other words, there have to be tensile forces acting on the crosslinked filaments in the network in order for it to be mechanically stable with a non-zero shear modulus. This notion is familiar from the ''tensegrity'' concept in biology and engineering. 54 For stiff athermal filament network, the window of pre-tension between the linear modulus G 0 = 0 at x E 1 À 1/p(2c) 2/3 , and G 0 -N at x -1 (we assume inextensible chains) is very narrow. The stability boundary of the network is plotted in Fig. 6. The condition for G 0 = 0 gives the equilibrium case labelled as g = 0 in the phase diagram. However, the full analysis shows that the magnitude of the shear strain modifies the stability condition, as represented by the coloured curves for g = 0.2 and 0.5. This shift means that the region of mechanical stability of filament network expands on increasing deformation, which matches exactly what the recent paper by Sharma et al. states. 32 To summarize, in this work we develop a continuum elastic theory of a network of semiflexible filaments, by implementing the general free energy of one semiflexible chain into that of a disordered network. On reflection, we choose the 3-chain model as most closely matching the realistic system -and achieve various quite stringent fits for a number of different experimental data sets over the full range of nonlinear stress-strain range. We demonstrate that the general theory produces the conditions for positive/negative Poynting effect (normal stress) in a network under imposed shear. The greatest weakness of this model is omitting the effects of local strain non-affinity, which might play an important role when crosslink density of the network is small and filaments stiff. However, the very broad agreement of our affine theory with so many different experiments is reassuring. The effect of tensegrity, or linking of filament pre-tension to the network stability and the magnitude of the shear modulus G is an unexpected result of this theory. We believe that the presented analytical continuum model can provide as an efficient and portable tool for studying the mechanical properties of semiflexible networks.