## Abstract

We present a model that assesses the different elastic responses of a semiflexible network, which either (i) is constrained to deform in an affine way or (ii) is permitted to thermally fluctuate and deviate from affine response. The thermal, non-affine response of the network is achieved using a Metropolis Monte Carlo algorithm with dynamic step size. We find that non-affine deformations soften the network dramatically at low strains and make the eventual nonlinear strain stiffening far more pronounced. We show that the effect of these non-affine deformations are very sensitive to the degree variation in the lengths of filaments connecting cross-links. Where there is high variation, non-affine deformations allow internal stresses to relax, giving rise to a smaller range of tensile forces in filaments and a dramatic reduction of network stiffness. This highlights that non-affine deformations are crucial in small strain response of stiff polymer networks.

## 1. Introduction

Semiflexible polymer networks are collections of stiff filaments that are physically connected to one another via chemical cross-links or entanglements. The elasticity of these networks has wide-ranging implications, particularly in biology—the cytoskeleton of eukaryotic cells (Alberts *et al.* 1994; Fernandez & Ott 2008), the extracellular matrix (Rosenbloom *et al.* 1993; Wen *et al.* 2007) and the fibrin network in blood clots (Weisel 2004; Guthold *et al.* 2007; Wen *et al.* 2007)—are all semiflexible networks. The elasticity of these systems is complex, the response of a single filament depends on the response of all filaments connected to it, hence there are many interactions. These interactions most often mean that one must invoke assumptions about how the network responds to an imposed strain to make analytical progress. Models of conventional elastomer networks (e.g. rubbers) assume that the strain field is affine, which is to say that it is everywhere homogeneous and the same as that applied to the bulk. This assumption is remarkably successful for rubbers, correctly predicting the elasticity for a large range of strains, see Flory (1969), Deam & Edwards (1976) and Treloar (2005).

The affine strain assumption may not however be a physically sensible assumption for semiflexible networks as highlighted by Head *et al.* (2003), Onck *et al.* (2005) and Heussinger *et al.* (2007). Semiflexible networks are composed of much stiffer polymers, which typically bend only weakly owing to thermal fluctuations, hence such networks resemble scaffolds rather than the melts typical of rubbers (see, for example, Gardel *et al.* 2004*a* for microscopy images). This bending rigidity can be cast in terms of a length scale (the persistence length *l*_{p}), which measures the correlation length of tangent vectors along the chain backbone. For semiflexible polymers, this length *l*_{p}≈*L*, where *L* is the typical length of filaments. The interplay of these two length scales produces many of the non-trivial properties of semiflexible chains, in particular a highly nonlinear free energy as a function of extension, as in Fixman & Kovac (1973), Wilhelm & Frey (1996) and Blundell & Terentjev (2009). This means that given a distribution of initial filament lengths, the shorter ones will cost much more in free energy to stretch by the same factor compared with longer ones. Under these circumstances, one might wonder whether the affine strain will not vastly overestimate the stiffness of the network.

There are two further physical concerns associated with the affine strain assumption. Firstly, it only allows one to consider the stretching of a filament between cross-links and does not allow the possibility that filaments can bend at a cross-link. Such a concern has indeed been the subject of some recent research (Head *et al.* 2003; Onck *et al.* 2005; Heussinger *et al.* 2007; Huisman *et al.* 2008), where the effect of a bending deformation over stretching is shown to be crucial to understanding small strain elasticity. Secondly, under the affine strain assumption, the cross-links connecting filaments are not allowed to fluctuate thermally, an effect that recently has shown to be important in understanding the large strain properties of conventional elastomers (Xing *et al.* 2007). One might legitimately ask whether these effects will also be important in semiflexible networks. The importance of non-affine deformations is not limited to cross-linked semiflexible filaments either. It has been shown, for example, in Lieleg *et al.* (2007) that non-affine deformations dominate the mechanical response of bundled actin filament networks too, and so an investigation of the effect of non-affinity is timely.

In this paper, we develop a Metropolis Monte Carlo algorithm that calculates the elasticity of a model network allowed to find its true thermal equilibrium in three dimensions, given an imposed strain. That is, the nodes in the network (representing the cross-links) are permitted to move in a non-affine way, and allowed to thermally fluctuate. Our model network is of course a simplification, but it explicitly allows us to independently vary the energy cost of stretching filaments and the energy cost to bend around cross-links. It therefore enables us to examine how these factors affect the influence of non-affine deformations on network elasticity. We point out that this is still very much a theoretical approach and that true experimental systems have important complications, such as the kinetic unbinding of cross-linkers (Lieleg & Bausch 2007) and finite cross-linker stretch modulus (Broedersz *et al.* 2008), which we do not treat here. Analytical work on non-affine deformations in semiflexible networks has been attempted successfully, but is often confined to two dimensions (Heussinger *et al.* 2007; Mizuno *et al.* 2007) and, while simulations of networks in three dimensions have also been studied, none of the works examine the effect of thermally fluctuating cross-links. We use this model to investigate the effect both non-affine strains and thermally fluctuating cross-links have on the elasticity of the network by comparing with the elasticity of the same network that has been assumed to deform homogeneously. We find that both effects soften the network at small strains, but make the eventual nonlinear response more dramatic at larger strains.

## 2. Model network

The network is modelled as a collection of *N*^{3} nodes. The nodes are initially placed on a regular *N*×*N*×*N* cubic lattice (figure 1*a*). Each node is uniquely identified by a triad of integers (*i*,*j*,*k*), which is the position vector of the lattice point at which it is initially placed. In general, the position of the (*i*,*j*,*k*) node will *not*, however, be (*i*,*j*,*k*) since the node is allowed to move. The (*i*,*j*,*k*) acts simply as a convenient label.

Connections in the network (representing the physical filament connections) are formed between the (*i*,*j*,*k*) node and all of its six nearest neighbours (*i*+1,*j*,*k*), (*i*−1,*j*,*k*), etc., with the exception that if the node initially lies on one, two or three of the exposed faces of the lattice then the number of nearest neighbours will be fewer—namely five, four or three, respectively. This means that an *N*×*N*×*N* network will have 3(*N*^{3}−*N*^{2}) connections (each of the *N*^{3} nodes has three connections; however, since the three faces at *i*=*N*, *j*=*N* and *k*=*N* will have counted *N*^{2} connections, which do not belong to the network, we subtract these). Each one of the connections in the network—which we label with the index *q*—will have a separation *R*(*q*). While on the lattice points, all *R*(*q*)=1, so that we have taken our unit of length to be the lattice spacing. This completely defines the topology of the network.

The chosen case of *z*=6 whereby internal nodes connect to six others represents the marginally rigid case in that the Maxwell counting condition (Maxwell 1864*a*,*b*) for rigidity yields zero free degrees of freedom. Although few physical networks possess such coordination, this is an interesting test case since it lies on the line between floppy and rigid and under shear is predicted to have a zero shear modulus at zero strain (Thorpe & Duxbury 1999).

The positions of all the nodes, excluding those in the top and bottom layers (*k*=0 and *k*=*N*) are then shifted from the lattice sites by adding a Gaussian smear about their original lattice site positions (*δi*,*δj*,*δk*), which are random variables drawn from a Gaussian of zero mean and variance *α*^{2}:
2.1These *r*_{ijk} define the initial, zero energy state of the network. The separations between nodes in the network after this process will also be Gaussian distributed according to
2.2where the variance in the spans is twice that in the positions because spans are formed from the difference in two positions. The mean separation between nodes is always unity.

In this initialized state, we label the spans *R*_{0}(*q*) to emphasize that these separations are the ones that give a minimum free energy (taken to be zero). The top and bottom layers of the network are left unchanged because throughout the entire algorithm these surfaces are where the strain is applied. Strictly speaking, therefore, the connections in the layers *k*=0 and *k*=*N* are distributed with a variance of *α*^{2} as opposed to 2*α*^{2} as in the bulk of the network. Physically one can think of the nodes in the top and bottom layers as being rigidly grafted to two external plates which impose the deformation. The parameter *α* is a measure of the quenched disorder in the initial separations between nodes in the network (figure 1*a*–*c*).

The square lattice is anisotropic, and this anisotropy will remain even in the networks with non-zero *α* since there is still long range correlation in positions. However, we are considering only one form of deformation in this work (shear deformations), so the reader should bear in mind that what applies for shear in this paper could potentially be different if we consider other modes of deformation (e.g. extension) since the network is inherently anisotropic.

### (a) Network energy

With the initial zero energy state of the system defined, we need a way of associating energies with any other configuration of the network. This is done by associating a free energy with every connection in the network and adding up over all connections. We use the free energy expression from Blundell & Terentjev (2009), which is valid for stiff filaments and which has the benefit of being simple and hence computationally fast to evaluate. As an example, the separation between nodes (*i*,*j*,*k*) and (*i*+1,*j*,*k*) is *R*=|*r*_{i+1jk}−*r*_{ijk}| and the free energy associated with that connection would be:
2.3where the subscript on the free energy is to remind that it comes purely from considering entropic stretching, and where, as usual, *l*_{p}=*A*/*k*_{B}*T* and *L* is the contour length of the filament connecting nodes (*i*,*j*,*k*) to (*i*+1,*j*,*k*).

The contour length can be calculated from the initial separation between the nodes *R*_{0} by demanding that the tension in the connection (*i*,*j*,*k*) to (*i*+1,*j*,*k*) was initially zero. This is equivalent to demanding that at upon formation of the network, the free energy *in each individual connection* is minimum. This does not allow for initial internal stresses as would be allowed under a so-called tensegrity structure (see, for example, Brangwynne *et al.* 2006). Although it has been shown experimentally that cross-linked networks often possess pre-stressed regions (Schmoller *et al.* 2008; Lieleg *et al.* 2009), considering such effects would inevitably lead to having to assume an initial pre-stress distribution, which would depend strongly on material preparation. To our minds, it is therefore more physically sensible and simpler to assume a completely relaxed network. Furthermore, requiring zero tension in each connection yields an expression for *L* in terms of *R*_{0}
2.4This method of calculating *L* given a separation *R*_{0} breaks down when *R*_{0}>*π*^{3/2}*l*_{p}/4≈1.39*l*_{p}. The reason being that beyond this point the minimum free energy occurs for *R*_{0}=0, which makes *L* undefined. We therefore restrict ourselves to stiff filaments whereby *l*_{p}≫*R*_{0}, which is self-consistent because the formula (2.3) is only valid in this limit, as discussed in Blundell & Terentjev (2009). Given a set of initial node positions, is it therefore possible to calculate the contour lengths of all connections in the network this way?

One issue with equation (2.3) is that it diverges as . This difficulty has been discussed in detail in Kierfeld *et al.* (2004) and Blundell & Terentjev (2009), and the resolution adopted here is based on a similar argument. In short, the divergence is unphysical for filaments with because at some point there is a cross-over to a region whereby it becomes energetically favourable to mechanically stretch the filament and the assumption of rigid in-extensibility must break down. To account for this effect, we define a sharp cutoff extension *x*_{c}=*R*_{c}/*L*, which is the extension at which the entropic stretching stiffness *K*_{e}(*x*)=2*k*_{B}*T*/(*l*_{p}*L*(1−*x*)^{3}) becomes higher than the mechanical stiffness of the filament (here, the integral is over the cross-sectional area of the filament and *E* is Young’s modulus, see, for example, Landau & Lifshitz (1986)). For *R*/*L*>*x*_{c}, the energy is simply that of pure mechanical stretch:
2.5and
2.6where *R*_{h} is defined by matching the force across the boundary at *R*=*R*_{c}. That is, at the point of cross-over *R*_{c}, both the force and the gradient of the force (stiffness) are continuous. Given the contour lengths of all connections in the network, the free energy of the entire network with conformation {*r*_{ijk}} is calculated by working out the set of all separations between nodes *R*(*q*), calculating the energy associated with each connection using equation (2.6) and summing over all *q* connections.
2.7where *q* labels all connections in the network and we draw the reader’s attention to the notational difference between *F* the free energy of a single filament, and , the free energy of a network. By subtracting the energy of the initial state, we ensure that the energy in the initial, undeformed state of the network is zero. It is important to point out here that equation (2.7) assumes no interaction between connections—for example, does not (yet) account for a bending energy across a node. We examine the effect of including bending energies in §7.

Calculating the free energy of the network as a function of imposed strain, under the assumption of affine strain, is now straightforward. We deform all nodes in the network by the same deformation gradient ** Λ** so that
2.8for all connections

*q*. The resulting free energy of the network will be given by 2.9which is easily calculated.

The true free energy of the network that is allowed to thermally fluctuate and deform non-affinely is obtained using the following Metropolis algorithm, as in Metropolis *et al.* (1953), combined with dynamic updating of the trial step size to achieve faster convergence as discussed in Frenkel & Smit (1996). We choose to calculate the true free energy including thermal fluctuations for two reasons: firstly, the effect of thermally fluctuating cross-links has been shown to be important, and, secondly, because the energy landscape of the network is rugged, hence a simple energy minimization algorithm may not in fact find an equilibrium area of phase space.

Deform all nodes under an affine strain so that

*r*_{ijk}()=*Λ**Λ**r*_{ijk}for all nodes (*i*,*j*,*k*).For all nodes, except those in the top and bottom layers (

*k*=0 and*k*=*N*, which are held fixed), we form a trial state by making a small Gaussian smear about the current positions′*r*_{ijk}(*Λ*)=*r*_{ijk}(*Λ*)+*δ**r*_{ijk}, where*δ**r*_{ijk}is a vector where each component is drawn from a Gaussian distribution with zero mean and a variance*σ*^{2}(which we dynamically update, see, for example, Frenkel & Smit (1996)).Calculate the energy of both the original state and the energy of the trial state

Accept the trial configuration with probability given by the ratio of the statistical weights of each state , where if then

*p*=1.Record the statistical weight of the current state and add it to the total statistical weight of all trials so far . This forms an estimate for the partition function.

Repeat steps 2–5.

The free energy for a given imposed strain is then calculated via 2.10

As with all Metropolis Monte Carlo algorithms, we must be careful to discard trials before the algorithm converges on the desired region of phase space. The minimum number of trials one must discard, , before equilibrium statistics recorded is estimated by looking at a physical quantity, in this case, the average energy (figure 2) and looking for the number of trials after which there is no change in the physical quantity with increasing trial number. In this case, . Beyond this point, we are in equilibrium and so can begin recording statistics. The variance of the trial steps is updated every 30 trials by 5 per cent so as to drive the algorithm towards an acceptance probability of *p*≈0.5 as discussed in, for example, Frenkel & Smit (1996).

## 3. Affine deformation

We now examine the effect quenched disorder in initial filament separations *α* has on the elasticity of our semiflexible filament network. We restrict ourselves initially to considering shear deformations, although the algorithm can handle any form of deformation in principle. In this case, the deformation gradient tensor is
3.1We initially consider the free energy and modulus of the network as a function of applied strain assuming that the strain field is affine, with the view to comparing it to non-affine deformations, which we do in §4.

In figure 3, we plot the free energy density and in figure 4, the modulus as a function of imposed affine strain. The units used for the free energy per unit volume and hence the modulus are *k*_{B}*T*/*L*^{3}. For a network composed of rod of length 1 μm at room temperature 1(*k*_{B}*T*/*L*^{3})≈0.1 Pa. The modulus and free energy are plotted for the three networks shown in figure 1 with *N*=5, composed of filaments with *l*_{p}/*L*=5 (here *L* is the mean length of connections) with *α*=0.0 (filled squares), *α*=0.5 (open circles) and *α*=1.0 (open squares).

For *α*=0 at small strains, the free energy scales like . This agrees with what we would expect. Under the affine shear strain applied, for zero disorder only vertical connections are stretched (those connecting adjacent *k* points like (*i*,*j*,*k*) to (*i*,*j*,*k*+1)). For small shear strains, a vertically oriented filament will be stretched by *R*−*R*_{0}≈1/2*ϵ*^{2}. Since the free energy is a minimum at *R*=*R*_{0}, Taylor expanding the energy about *R*_{0} gives a lowest order term, which scales like . This means that the modulus (figure 4) goes to zero at zero strain. This is what we expect since for the chosen cubic geometry with free hinges, the shear deformation is a floppy mode (completely soft up to first order in the deformation; Thorpe & Duxbury (1999)).

As one increases the strain, one observes at *ϵ*≈0.2 a region of dramatic nonlinear behaviour whereby modulus increases by two orders of magnitude over strains of 10 per cent. This nonlinearity is a consequence of the individual nonlinearity in connections and makes quantitative sense. The average filament length is unity, so if *l*_{p}/*L*=5, using equation (2.4) the average filament will have an initial extension of *x*_{0}≈0.965 and so the divergence of the free energy would occur at strains where *R*−*R*_{0}≈1/2*ϵ*^{2}≈0.035, which corresponds to strains in the region of 30 per cent, consistent with figure 3. Beyond these strains, we see the cross-over to mechanical stretching of the filament whereby the modulus converges to the same value for all networks.

Increasing the quenched disorder to *α*=0.5 alters the physical behaviour of the network in three important ways. Firstly, the network becomes dramatically stiffer. Figure 4 shows that the modulus at zero strain is finite, and consistently approximately two orders of magnitude larger at small strains. This stiffening of the network is a feature of the affine strain assumption. Increasing the quenched disorder means that instead of filaments being the same length, some are longer and some are shorter. Since the free energy function (2.6) is a nonlinear (convex) function of extension, the shorter filaments are in fact much stiffer than the longer ones. Assuming all filaments stretch by the same factor owing to affine deformation, therefore, means that the shorter filaments contribute much more to the free energy than the longer ones. This results in a stiffening of the network.

The second important feature is that the free energy has a weaker functional dependence on strain in the small strain region as illustrated by the small-dashed line in figure 3. This quadratic dependence is linear elasticity—the shear modulus as shown in figure 4 being constant at small strains. This dependence is understandable because as one increases the disorder in node positions, some connections will be oriented close to the principal stretch direction of the applied strain—for the shear deformation applied this is the (1,0,1) direction. These connections will be stretched by a factor Δ*R*∼*ϵ* and since the free energy of each filament remains a quadratic function about *R*_{0}, the free energy scales like . We also note that the nature of the stiffening is not as dramatic as in the case of zero disorder, which we can attribute to the fact that if we have distributed lengths but at the same time assume all these filaments deform by the same factor, they will stiffen at different values of the applied strain *ϵ* hence the nonlinear elasticity is spread over a wider value of strains and so stiffening is smoothed out.

These results extend for a further increase in quenched disorder to *α*=1.0. The network stiffens by an order of magnitude and the nonlinear stiffening becomes less dramatic. It is comparatively straightforward to see why increasing the disorder in the network results in stiffer networks under the affine strain assumption. What is less clear is whether this is a physically sensible assumption when one has some disorder in the network.

## 4. Non-affine deformation

For the same three values of the disorder parameter *α*, we examine the shear modulus as a function of strain for the same network as in §3 that is now permitted to deform non-affinely and thermally fluctuate. The important point here is that all nodes except those in the top and bottom layers can fluctuate away from the affine positions and are accepted or rejected based on the Metropolis *et al.* (1953) algorithm. This means that, given enough trials (), the positions of the nodes that the algorithm samples will be those of a truly thermally equilibrated network. Figure 5 plots the shear modulus for the three thermally equilibrated networks of increasing disorder and, for comparison, the shear modulus predicted by affine response from §3. For ease of explanation, we refer to the moduli of networks allowed to deform non-affinely by *G*_{NA}, and the moduli of affinely deformed networks by *G*_{A}. The modulus (also in units of *k*_{B}*T*/*L*^{3}) is defined as the discrete second derivative of the free energy
4.1for some finite difference *δϵ*, which in the cases here were taken to be 0.01. The effect of quenched disorder in filament positions has directly the opposite effect on the network that is permitted to deform non-affinely compared with ones constrained to deform affinely. For the case of the perfectly ordered network *α*=0.0, *G*_{NA} is indistinguishable from *G*_{A}. This result is reassuring since in such an ordered system where all connections have the same length, the minimum in free energy for the network must be a state where all vertical filaments have stretched by equal amounts. If this were not the case, one connection stretching more than the average must imply another connection, which stretches less—since the vertical filaments are fixed at the boundaries. Since all filaments are initiated in a state of minimum free energy and possess a free energy function which is convex, in the absence of any quenched variation in filament lengths, a deformation in which one filament stretches more than the average and another less will always be higher in energy than the affine strain deformation. This however is not true for a network with quenched disorder.

As we introduce a variation in filament lengths (beginning with *α*=0.5), we see two significant results. Firstly, contrasting affine and non-affine responses we see that, for strains up to ≈ 30 per cent, *G*_{NA} is two or three orders of magnitude softer than *G*_{A}. For strains below about 20 per cent, *G*_{NA} shows only a weak dependence on strain—changing by at most an order of magnitude—though this is hard to see from the data because errors in *G*_{NA} are large when compared with true variations. Beyond this soft region, up to ≈ 20 per cent strains, is a region of dramatic stiffening, with *G*_{NA} increasing by three orders of magnitude between strains of 20 and 40 per cent. The affine response in contrast stiffens more weakly over a larger range of strains. Reassuringly, both moduli *G*_{NA} and *G*_{A} tend to the same value at large strains as one would expect since in this region the dominant contribution to the modulus comes from mechanical stretching of filaments. The other significant observation is that *G*_{NA} for *α*=0.5 is consistently smaller than *G*_{NA} for *α*=0.0. That is, even when we compare two non-affinely deforming networks, the one with the bigger quenched disorder in positions is softer—the reverse of the trend for the affine case.

The same trend is seen on a further increase in disorder to *α*=1.0. *G*_{NA} and *G*_{A} show a wider disparity and, again, the network is softer than the equivalent networks with lower values of *α*.

The influence of non-affine deformations on the modulus–strain relationship *G*(*ϵ*) striking and is a large departure from what one would expect from affine response, particularly at lower strains. Many experiments on semiflexible networks report modulus–stress relationships *G*(*σ*) since it enables one to disentangle the true elastic response from the viscous response (see Fernandez & Ott 2008 for a discussion). We plot modulus–stress relationships for the three networks of increasing disorder in figure 6.

We observe that although the quenched disorder in the network *α* has a dramatic effect on *G*(*ϵ*), it has almost no effect on *G*(*σ*). All three datasets follow an approximate scaling relationship at low stresses of *G*∼*σ*^{ν}, with the value of *ν* lying between 1 and 1.5. This is a significant result since these scaling relationships have been observed in single fibroblasts (Fernandez *et al.* 2006) and *in vitro* networks of actin filaments (Gardel *et al.* 2004*a*,*b*), where in the latter the results we interpreted as an indication that an affine theory captures the correct elasticity of semiflexible networks. Here, we observe that in fact non-affine responses are also consistent with these experimental findings and that to distinguish between affine and non-affine, one should focus on *G*(*ϵ*) relations.

## 5. A scalar measure of non-affinity

We would like to try and understand why non-affine deformations lead to dramatically softer networks, when there is a quenched disorder in node positions. Clearly, the idea that strain is distributed homogeneously though the network does not work well for these stiff filament networks. We hypothesize that perhaps a more physically sensible picture is one in which the *stress* is homogeneous. The picture however is complicated by being in three dimensions. For equilibrium, the only condition on the stress field is that it should be divergenceless (∇_{i}*σ*_{ij}=0)—not necessarily constant. Nevertheless, appealing again to a simple one-dimensional model, we can understand why constant stress could be more physical.

If one models the network as a series of springs connected in series then the force along this line must be constant. For *α*=0.0, all filaments have the same stiffness, hence constant force also means constant extension and the connections will deform affinely. An increase in the variation of filament lengths (say to *α*=0.5) also increases the variation in filament stiffnesses in a dramatic way since the stiffness is a highly nonlinear function of length. That is, increasing *α* means the network is composed of some much longer (softer) filaments and an equal number of shorter (stiffer) filaments. In a simple picture where forces must be constant this would result in the longer, softer filaments stretching more than the shorter, stiffer ones, and the overall stiffness would be dominated by the softer filaments as they add in reciprocals. This would explain the observed softening of the network with increasing *α*, but it would also predict two important trends. It would firstly predict that increasing the quenched disorder in the network should result in larger non-affine deformations. Secondly, even if stress is not perfectly constant in the non-affine network, it would predict that the variation in forces observed in the non-affine case should be far less than for the affine case.

We test both of these predictions. Firstly, we plot a scalar measure of non-affinity Δ both as a function of increasing length variation *α* (figure 7), and as a function of strain (figure 8). The non-affinity measure Δ is defined as
5.1where the angled braces subscripted by *T* refer to ensemble average and subscripted by *q* means averaged over all network connections. This measure is therefore the modulus of the ensemble average deviation of a connection span from affine response which is then averaged over all connections. It is zero if the deformation is affine or if the deformation is uniformly distributed around an affine deformation. Non-zero values of Δ therefore imply truly non-affine strains.

We do indeed see that as one increases the disorder in the network, the measure of non-affinity increases. We also see that if there is some quenched disorder in the positions, then measure of non-affinity increases as one applies larger strains. This strongly suggests that for networks where there is a quenched disorder in filament lengths and positions, non-affine deformations are a crucial factor in the elasticity and enable the network to deform with a far smaller cost in free energy than would be predicted by the affine theory. Notice that the non-affinity measure Δ is non-zero even for the zero-disorder case *α*=0. This is seemingly in conflict with the previous section where zero disorder was shown to follow affine response closely. This effect is in fact owing to the non-zero temperature at which the simulation is run as we discuss in the following section.

Figure 9 plots histograms showing the distribution of forces throughout the network. We plot the frequency of occurrence *n* of the logarithm of modulus of the forces in the 300 network connections for a network with *α*=1.0, at increasing values of strain. Note that the force is measured on a log scale, i.e. if all 300 connections were under equal tension *f* then there would be a single peak of height 300 at . On each histogram, we plot the distribution of forces for affine response (grey) and non-affine/thermal response (white).

The most important observation by contrasting non-affine to affine is that the variance in the forces is consistently lower for the non-affine case at all strains, except zero strain. The forces in connections that undergo affine response are, as we would predict, spread over a wider range of values, and in particular have a longer tail in the high-force limit (a few connections that are subject to very high forces). Forces in connections are allowed to deform non-affinely are more strongly peaked at intermediate forces and do not have a long tail with a small proportion of the network connections carrying very large forces. This reinforces the point that affine response will overestimate the stiffness of the network because it assumes all connections deform the same way leading to some being under large tensions. The non-affine response allows at least some of the filaments to relax by moving away from affine positions and this shows up in the narrower distribution in forces. The point we want to make is that although the stress is far from constant—the peak is not a delta function—it is far narrower than the distribution of the forces seen for the affine case. Since these distributions are in log-space, the typical reductions in forces observed are on the order of 10^{2}−10^{3}(*A*/*L*^{2}) for some filaments (see histogram corresponding to *ϵ*=0.2, for example), which is a very significant reduction. It does not seem unreasonable therefore that semiflexible networks might be better modelled by a constant stress rather than a constant strain approach.

## 6. The effect of temperature

One of the most striking observations from figures 7 and 8 is that the non-affinity measure Δ is not zero for the case *α*=0, for small strains. In this perfectly cubic lattice where each of the connections is the same length, one would expect there to be no non-affine deformations, since they all have a higher free energy. How then is it possible that we observe non-zero values of Δ? The resolution of this comes from the definition of Δ
Although the thermal average (that is, averaging over samples drawn from the Metropolis algorithm) should yield zero for samples distributed symmetrically about the affine position, for a finite number of trials, this is never going to be the case and one is left with a small non-zero values of Δ. We verify that this non-zero value of Δ is indeed owing to the effect of a finite temperature by considering the temperature dependence of Δ^{2} for the cases *α*=0 and *ϵ*=0.

Moving a node away from its zero energy position weakly stretches some filaments and compresses others. Provided that the deviation in position is small, the penalty associated with this will be quadratic in the position because to first order the potential of semiflexible polymers is quadratic. Hence, for small ** R**(

*q*)−

*Λ*

*R*_{0}(

*q*), the energetic penalty associated with

**(**

*R**q*)−

*Λ*

*R*_{0}(

*q*) is

*E*∼(

**(**

*R**q*)−

*Λ*

*R*_{0}(

*q*))

^{2}. By equipartition, we should therefore observe that 〈(

**(**

*R**q*)−

*Λ*

*R*_{0}(

*q*))

^{2}〉=Δ

^{2}∼

*k*

_{B}

*T*. Figure 10 plots Δ

^{2}against

*T*on logarithmic scales and we do indeed see a scaling Δ

^{2}∼

*T*as we would expect. This verifies that the non-zero values of Δ observed for small strains in figures 7 and 8 are due to thermal fluctuations.

## 7. Bending across nodes

One of the obvious criticisms one could level against our model is that the nodes are essentially free hinges. In reality, the cross-links in a semiflexible network simply link two filaments together, and while this cross-link itself may allow the filaments to rotate with respect to one another, there will still be an energy penalty for bending a single filament through the cross-link (figure 11). The effect of bending across a cross-link has been considered in some two-dimensional models of semiflexible networks before by Head *et al.* (2003) and by Heussinger *et al.* (2007) and one might wonder therefore if it has an influence in our model too. It is clear that whether or not bending will be important depends on the imposed deformation. If one imposes compression, for example, then a freely hinged regular structure will be floppy, while in contrast the same structure would not be floppy to shear deformations.

To check the influence of bending across a cross-link, we examined the shear modulus and measure of non-affinity Δ by including an energy that penalizes bending across a link. By this we mean the following: if we have a node labelled by (*i*,*j*,*k*) then the connections and are assumed to be the same filament as are connections that link between *j*’s and *k*’s (figure 11 shows the connections belonging to the same filament in the same colour).

Referring to figure 11, we mean that there should be an energy penalty associated with the angle *θ*, but not *ϕ*. Therefore, we modify the energy function (2.9) to include an additional term 1/2*k*_{b}*θ*^{2}, where *k*_{b} is a measure of bending rigidity. A rod with bending modulus *A* separated into links of length *b*, bent through an angle *θ*, the associated energy is 1/2(*A*/*b*)*θ*^{2}, so *k*_{b}=*A*/*b*. Since we have taken the lattice spacing *b*=1 for the network, this means we use the same value for *k*_{b} as for *A* in the free energy of the individual connections. The shear modulus (figure 12) and non-affinity measure (inset) are not affected by including this bending energy term. This is of relevance to biological networks of polymers that are cross-linked by molecules that allow free rotation.

On the other hand, if cross-links are much stiffer (for example, the ARP2/3 complex that forms the branching in dendritic actin Blanchoin *et al.* (2000)), the picture changes. In this case, it is no longer true that the effect of an additional rod reduces the number of degrees of freedom by 1. In fact, for a rigid linker, each rod reduces the degrees of freedom in the network by 3. In such a case, the network is always rigid independent of the network coordination and shows a much weaker stiffening on deformation as shown in figure 13. This highlights that the softening of networks owing to non-affine deformations will only be observed if the underlying network has low-energy deformations that can be explored. Put another way, affine theory places additional constraints on the way in which vertices in the network deform. If these constraints prohibit the network from exploring low-energy deformations (as in the case of freely hinged networks) then the affine theory will be a poor approximation and non-affine deformations will be crucial to network elasticity. If, on the other hand, other effects (such as rigid cross-links) are placing additional constraints on the network, low-energy deformations may no longer exist and the effect of non-affine deformations will be less dramatic.

We stress that this discussion of the importance of bending and non-affine deformations is very much dependent on the type of deformation applied. Although a weak bending energy at vertices is not relevant in shear geometry, it will certainly be important under simple compression. Shear necessarily stretches some filaments and because stretching energies are highly nonlinear functions of extension these dominate the network response. Under compression, it is not necessary to stretch connections and one would therefore expect bending energy to play a far more important role.

It is also crucial to point out that the influence of bending stiffness on mechanics and the extent to which non-affine deformations soften a network relative to affine deformations will depend strongly on the coordination of the network. In the case of the marginally rigid networks considered here, the effects are dramatic. For over-constrained networks composed of many cross-linked filaments this may no longer hold (Head *et al.* 2003) and for under-constrained networks one would expect the influence of bending stiffness to be crucial (Heussinger *et al.* 2007). For networks that lie in the middle between these extreme cases, however, our work shows that the effect of disorder and non-affine deformations are of great importance.

## 8. Conclusions

Semiflexible networks are complicated interacting systems and to construct models of elasticity one usually needs to make assumptions about how the imposed deformation is distributed in the network. The usual approach is often to assume that the strain field is affine—locally being the same everywhere—which works well for conventional elastomer networks. For semiflexible networks, however, an affine theory vastly overestimates the modulus and the approximation gets worse at non-zero strains. That is, if one wants to model the modulus–strain relationship for a semiflexible network, non-affine deformations must be accounted for. These non-affine deformations are able to relax some filaments such that the distribution of forces in the network is far narrower for a network allowed to deform non-affinely compared with one constrained to deform affinely. We propose therefore that modelling semiflexible networks in constant stress regimes might be a more physical approach than modelling them in the constant strain regime. For example, one line of potential enquiry might be to investigate whether some properties of semiflexible networks can be explained using simple one-dimensional models of filaments in a constant force ensemble whereby one obtains the average strain of the network as a function of the applied stress by inverting the expression in equation (2.3).

The influence of quenched disorder in a semiflexible network is also important. Since semiflexible filaments have such nonlinear force–extension relationships, filaments of different lengths have dramatically different stiffnesses. This means that, in general, the network would like to find deformation fields that stretch the softer, longer filaments at the expense of the shorter, stiffer ones since such deformations must be lower in energy.

We find that, for the shear geometry considered, including the effect of bending energies across cross-links for physical values of the bending stiffness, does not appreciably alter the elasticity of the networks.

## Acknowledgements

We would like to thank John Biggins, Robert Lee, Panayotis Benetatos and Alessio Zaccone for helpful discussions. We also acknowledge that the comments of the reviewers have helped clarify certain points. This work was supported by the TCM EPSRC partnership.

- Received November 16, 2010.
- Accepted January 25, 2011.

- This journal is © 2011 The Royal Society