Nonlinear and heterogeneous elasticity of multiply-crosslinked biopolymer networks

In the graph representation, the vertices correspond to the crosslinks, while the edges reflect the polymer contours (figure 1). To each vertex i, we assign a 3D vector ${{\bf{r}}}_{i}$ to denote the position of the crosslink; to each edge linking two vertices, say i and j, we assign a length ${{\ell }}_{{ij}}$ to denote the length of the polymer contour between these two vertices. For a given graph topology, the free-energy (to be specified later) is a function of the coordinates ${{\bf{r}}}_{i}$ and the lengths ${{\ell }}_{{ij}}$.

Zoom In Zoom Out Reset image size

The assignment of the crosslink coordinates ${{\bf{r}}}_{i}$ and the contour lengths ${{\ell }}_{{ij}}$ is performed as follows. We begin with a $L\times L\times L$ periodic simulation box. Next, ${N}_{\mathrm{poly}}$ polymer chains, each of length ${{\ell }}_{{\rm{c}}}$, are placed at random inside this box. These polymers are modeled as discretized WLCs, with persistence length ${{\ell }}_{{\rm{p}}}$; the discretization step-size (or bond length) is denoted by b, so each chain contains ${{\ell }}_{{\rm{c}}}/b+1$ vertices linked by the bonds. We ignore excluded-volume interactions between bonds, and so it is fairly straightforward to 'grow' the polymers. For any polymer, the first bond ${{\bf{b}}}_{12}$ (between vertices 1 and 2) is placed randomly in the box with random orientation. The orientation of the second bond (between vertices 2 and 3) is sampled from the probability distribution appropriate for the WLC model, $P({\hat{{\bf{b}}}}_{23})\propto \mathrm{exp}[({{\ell }}_{{\rm{p}}}/b)\;{\hat{{\bf{b}}}}_{12}\cdot {\hat{{\bf{b}}}}_{23}]$, $| {\hat{{\bf{b}}}}_{{ij}}| =1$, and so forth until the entire chain has been grown. Since for stiff chains (${{\ell }}_{{\rm{p}}}/b\gg 1$), the latter probability distribution can be very narrow, we use the integral-transform method, rather than rejection-sampling, for maximum efficiency. This process is repeated for all polymer chains, after which we end up with a 'gas' of non-interacting WLC polymers (but still without crosslinks).

To introduce the crosslinks, in principle, we prepare a list of all vertex pairs. The list is sorted according to the pair separation, so that pairs with the smallest distance between the two vertices appear first. Next, the first $i=1,\ldots ,{N}_{\times }$ pairs satisfying certain criteria (soon to be prescribed) are selected from the list, and the crosslink i is placed at position ${{\bf{r}}}_{i}$ defined to be the pair's geometric center. The set of coordinates ${{\bf{r}}}_{i}$ thus obtained define the graph vertices. The polymer chains between the vertices define the edges of the graph, which are stored in a connectivity table. Following this procedure, each crosslink is, by means of the edges, connected to 2–4 other crosslinks (we excise all dangling polymer ends). Note that when a crosslink has only two neighbors, this implies that its connecting edges do not belong to the same filament. Indeed, the order in which the edges originating from a crosslink are stored in the connectivity table is important, as we would like to distinguish between those edges that belong to one polymer, ${i}_{1}\leftrightarrow i\leftrightarrow {i}_{2}$, and the other edges that belong to the other polymer, ${i}_{3}\leftrightarrow i\leftrightarrow {i}_{4}$. This distinction is needed later on to define the three-body contributions to the free-energy. Lastly, for each edge, we store the length ${{\ell }}_{{ij}}$ of the connecting polymer chain (that is, by counting the number of vertices). Once these lengths have been stored, the graph representation is complete, after which the vertex coordinates are discarded.

The criteria, alluded to above, for the inclusion of any vertex pair as a crosslink, are as follows: (1) the vertices must belong to different polymer chains, that is, a single polymer chain is never crosslinked to itself. (2) The length of the polymer contour between two connected crosslinks may not become arbitrarily small: ${{\ell }}_{{ij}}\gt {{\ell }}_{\mathrm{cutoff}}$. (3) Each polymer must be crosslinked at least twice. (4) No disjoint sets of vertices in the network may remain. Note that the last two criteria imply that ${N}_{\times }$ cannot be chosen arbitrarily small.

The total free-energy of the network is given by

$$\begin{eqnarray}F & = & \displaystyle \sum _{[{ij}]}{f}_{2}({r}_{{ij}},{{\ell }}_{{ij}}) & & +\displaystyle \sum _{[{ijk}]}{f}_{3}({r}_{{ij}},{r}_{{jk}},{\theta }_{{ijk}},{{\ell }}_{{ij}},{{\ell }}_{{jk}}),\end{eqnarray} \tag{ 1 }$$

where ${r}_{{ij}}=| {{\bf{r}}}_{{ij}}|$, ${{\bf{r}}}_{{ij}}={{\bf{r}}}_{j}-{{\bf{r}}}_{i}$, ${\theta }_{{ijk}}={\mathrm{cos}}^{-1}[{{\bf{r}}}_{{ij}}\cdot {{\bf{r}}}_{{jk}}/({r}_{{ij}}{r}_{{jk}})]$.

The first term describes the two-body interactions, or pair-potentials, with the sum over all edges $i\leftrightarrow j$ in the network graph. The pair-potential is

$$\begin{eqnarray}&&{f}_{2}({r}_{{ij}},{{\ell }}_{{ij}})=-{k}_{{\rm{B}}}T\mathrm{log}P({{\bf{r}}}_{{ij}},{{\ell }}_{{ij}},{{\ell }}_{{\rm{p}}}),\end{eqnarray} \tag{ 2 }$$

where $P({{\bf{r}}}_{{ij}},{{\ell }}_{{ij}},{{\ell }}_{{\rm{p}}})$ is the radial end-to-end separation probability density of the WLC strand of length ${{\ell }}_{{ij}}$ and persistence length ${{\ell }}_{{\rm{p}}}$ connecting crosslinks i to j integrated over all orientations of both its end-tangents. Its functional form is taken from [30], which provides a good description for the tensile, compressional, buckling, and post-buckling strain regimes of a semiflexible WLC (see figure 2) within a biologically relevant range of stiffnesses.

Zoom In Zoom Out Reset image size

The second term in equation (1) describes the three-body interaction. In this case, the sum is over all triplets $i\leftrightarrow j\leftrightarrow k$, each representing a single polymer connecting crosslink i to j to k consecutively (with no other crosslinks in-between), and ${\theta }_{{ijk}}$ is the angle between ${{\bf{r}}}_{{ij}}$ and ${{\bf{r}}}_{{jk}}$. The idea is that the polymer backbone should persist through crosslink j and therefore penalize three-point bending between these three crosslinks. Its functional form is taken from [31]:

$$\begin{eqnarray}&&{f}_{3}({r}_{{ij}},{r}_{{jk}},{\theta }_{{ijk}},{{\ell }}_{{ij}},{{\ell }}_{{jk}})=\displaystyle \frac{3}{2}\displaystyle \frac{{k}_{{\rm{B}}}T{{\ell }}_{{\rm{p}}}{\theta }_{{ijk}}^{2}}{({{\ell }}_{{ij}}+{{\ell }}_{{jk}})}.\end{eqnarray} \tag{ 3 }$$

In a later publication, we will elaborate further on equations (1)–(3) and their interactions, particularly its underlying assumptions and the phenomenology that they neglect.

In what follows, we prepared our networks to mimick in vitro F-actin networks, for which the persistence length ${{\ell }}_{{\rm{p}}}\approx 17\;\mu {\rm{m}}$ at physiological temperatures, and the filament length is taken to be ${{\ell }}_{{\rm{c}}}\approx 10\;\mu {\rm{m}}$. The concentration of actin, defined as ${\rho }_{\mathrm{actin}}\equiv 2{{\ell }}_{{\rm{c}}}{N}_{\mathrm{poly}}/({s}_{\mathrm{actin}}{L}^{3}{N}_{{\rm{A}}})$ was taken to be $\approx 3\mu {\rm{M}}$, corresponding to typical experimental conditions [6, 29]. Here ${N}_{{\rm{A}}}$ is the Avogadro constant, and ${s}_{\mathrm{actin}}=5.46\;\mathrm{nm}$ is the typical axial separation between adjacent G-actin subunits along any strand of F-actin [32], and we chose ${{\ell }}_{\mathrm{cutoff}}=30{s}_{\mathrm{actin}}$, which is comparable to the size of filamin ($\approx 160\mathrm{nm}$), an F-actin crosslinker [33]. We generated 20 different networks, whose parameters are listed in table 1. In the rest of this article, unless otherwise stated, all energies are measured in units of temperature ${k}_{{\rm{B}}}T$, and all lengths in units of filament length ${{\ell }}_{{\rm{c}}}$.

Table 1.  A list of the parameters of the networks simulated for this article. ${N}_{\times }$ is the number of crosslinks; ${N}_{\mathrm{poly}}$, ${{\ell }}_{{\rm{p}}}$, ${{\ell }}_{{\rm{c}}}$ are the number, persistence length, and contour length of the filaments respectively; ${{\ell }}_{{ij}}$ is the average strand length; and L is the system-size. Networks $1-4$ yielded linear shear moduli $G\approx 0.2\;\mathrm{Pa}$ and bulk moduli $K\approx 0.3\;\mathrm{Pa}$. Networks 5–20 were derived by randomly removing only four-fold co-ordinated crosslinks from network 4 in order to keep ${N}_{\mathrm{poly}}$ and ${{\ell }}_{{\rm{c}}}$ fixed.

Network	${N}_{\times }$	${N}_{\mathrm{poly}}$	${{\ell }}_{{\rm{p}}}/{{\ell }}_{{ij}}$	${{\ell }}_{{\rm{c}}}/{{\ell }}_{{ij}}$	L
1	$1\;000$	300	15.7	6	${{\ell }}_{{\rm{c}}}$
2	$64\;000$	$19\;200$	15.7	6	$4{{\ell }}_{{\rm{c}}}$
3	$512\;000$	$153\;600$	15.7	6	$8{{\ell }}_{{\rm{c}}}$
4	$8000$	$2400$	15.7	6	$2{{\ell }}_{{\rm{c}}}$
5–20	$4600$–$7800$ in steps of 200	$2400$	—	—	$2{{\ell }}_{{\rm{c}}}$

Starting from a reference state with vanishing Cauchy stress we apply a shear deformation

$$\begin{eqnarray}{\boldsymbol{\Lambda }}=\left(\begin{array}{lll}1 & \gamma & 0\\ 0 & 1 & 0\\ 0 & 0 & 1\end{array}\right)\end{eqnarray} \tag{ 7 }$$

and investigate the nonlinear elastic properties of networks 1 and 2 (see table 1).

3.1.1. Shear modulus

The global elastic shear modulus $G=\frac{1}{V}\frac{{\rm{d}}{}^{2}{F}_{\phantom{{xx}}}}{{\rm{d}}{\Lambda }_{{xy}}^{2}}$ of network 2 is shown in figure 3 as a function of shear stress ${\sigma }_{{xy}}$, measured in units of ${\sigma }_{{\rm{c}}}$, which is the intersection between the low- and high-stress linear asymptotes of G. The linear modulus is approximately ${G}_{0}\sim 0.4$ Pa, which is on the lower end of values reported in experiments, but agrees well with, for example, Gardel et al [6]. As has been seen in previous work [31], a strain-stiffening behavior is observed (that is, G increases with increasing stress). Loss of numerical precision prevents us from increasing shear strains beyond 0.12.

Zoom In Zoom Out Reset image size

We also split the modulus in contributions due to the three-body (bending) and two-body (compression and stretching) interactions. Interestingly, both the linear and the nonlinear response are strongly dominated by the bending mode, while stretching only contributes very little.

A bending dominated linear response is by now well understood to be due to the excitation of non-affine 'floppy' modes [8]. For the nonlinear response, previous work [9] argues in favor of a transition from such a soft nonaffine bending mode to a much stiffer affine stretching mode. The stiffening then being due to the formation of percolating paths in the tensile direction of the shear deformation. Under large enough deformation the filament strands in these paths are brought to nearly full stretching, giving rise to a modulus that diverges with stress as $G\sim {\sigma }^{3/2}$ [6, 7]. In contrast, here we find that the bending response itself initially dominates nonlinear behavior. To further elucidate this apparently contradictory behavior we turn to a smaller system (network 1) with ${N}_{\times }=1000$ crosslinks but with otherwise the same parameters. In this system the shear strain can be pushed to much larger values and numerics is not an issue.

In figure 4 we plot the shear modulus $G({\sigma }_{{xy}})$ of this network as a function of shear stress ${\sigma }_{{xy}}$ together with the contributions from bending and stretching mode. As can be seen, $G({\sigma }_{{xy}})$ is dominated by bending at small and intermediate stresses, and by stretching only at the largest stresses. In this range, the characteristic scaling $G\sim {\sigma }_{{xy}}^{3/2}$ is observed, indicating the affine stretching of percolating paths (see figure 10 for an example of such a path). The bending-dominated nonlinear regime from figure 3 now shows up at intermediate stresses and displays a much stronger stress dependence than in the asymptotic ${\sigma }_{{xy}}^{3/2}$ regime.

Zoom In Zoom Out Reset image size

Between the low- and high-strain limits, non-monotonic changes of the two-body and three-body contributions may sometimes be also observed (as in figure 4, between the dashed vertical lines). In the case of the network sample of figure 4, this non-monotonic behavior is the response of the network to the buckling of a single short filament segment at that strain, whereupon the large local strain in the vicinity of the buckled filament causes surrounding filaments to further stretch and straighten leading to a larger increase in two-body stress and a correspondingly larger decrease in three-body stress (see the inset and the accompanying discussion in appendix A). Note, that the modulus itself hardly changes during this buckling event, only the relative weights for bending and stretching modes show rapid variation.

3.1.2. Normal stresses

Also of interest are the remaining components of the stress tensor ${\sigma }_{{xx}}$, ${\sigma }_{{yy}}$, ${\sigma }_{{zz}}$, ${\sigma }_{{yz}}$, and ${\sigma }_{{zx}}$, which are shown in figure 5 as a function of γ. We have restricted ourselves to small strains $\gamma \sim 1\%$. Hence the normal stress components (${\sigma }_{{xx}}$, ${\sigma }_{{yy}}$, ${\sigma }_{{zz}}$) are to a good approximation quadratic in the strain γ, furthermore ${\sigma }_{{yz}}={\sigma }_{{zx}}$ and all five components are small as compared to ${\sigma }_{{xy}}$—as expected.

Zoom In Zoom Out Reset image size

The normal stress components are observed to be positive, which with our sign convention, indicates that networks tend to contract when under shear deformation. This is in agreement with previous simulation results [8, 39, 40] and seems to be a general property of network-like materials [41]. The first normal stress difference ${N}_{1}\equiv {\sigma }_{{xx}}-{\sigma }_{{yy}}$ is the only other quantity (apart from the shear stress) readily accessible in cone-and-plate torsion experiments, and is directly proportional to the force exerted by the sample on the plates (for parallel-plate geometries this force is rather $-{\sigma }_{{yy}}$). N₁ has been reported as being negative for biopolymer networks at high strain [42]. Interestingly, in our system N₁ is positive and at small strains given by ${N}_{1}={\sigma }_{{xx}}-{\sigma }_{{yy}}=\gamma {\sigma }_{{xy}}$ (see inset figure 7), as expected for isotropic elastic solids [43]. Similar results were obtained in the theoretical considerations of [39], which also suggest a positive N₁. There, as well as in [42], it is argued that a negative N₁ in cone-and-plate geometries could be due to a relaxation of the hoop stresses (which correspond here to the ${\sigma }_{{xx}}$ component. A mechanism for this relaxation, however, is still elusive.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.1.3. Nonaffinity

Onck et al [9] have argued that stiffening is a network property rather than a single filament property. The important network rearrangements in 2D networks are idenfied as filament reorientations (rotations) which give rise to a maximum in the nonaffinity of the response.

To quantify the nonaffinity Γ, we compute the strand length due to a small increment in strain, ${{\ell }}_{{ij}}(\gamma +\delta \gamma )$ and compare it to the affine value ${{\ell }}_{{ij},\mathrm{aff}}(\gamma +\delta \gamma )$:

$$\begin{eqnarray}&&\Gamma \equiv \displaystyle \frac{1}{{N}_{\mathrm{bond}}}\displaystyle \sum _{i,j}^{{N}_{\mathrm{bond}}}\displaystyle \frac{{{\ell }}_{{ij}}^{2}(\gamma +\delta \gamma )-{{\ell }}_{{ij},\mathrm{aff}}^{2}(\gamma +\delta \gamma )}{{{\ell }}_{{ij}}^{2}(\gamma )}.\end{eqnarray} \tag{ 8 }$$

A plot of $| \Gamma |$ versus ${\sigma }_{{xy}}$ is shown in figure 6 and observed to increase monotonically with applied strain. However, restricting ourselves to tensile strands only, we do indeed observe a maximum roughly at the crossover (${\sigma }_{{xy}}\sim {\sigma }^{*}$) of the modulus to nonlinear behavior. Actually, ${\Gamma }_{\mathrm{tensile}}$ is negative, which indicates that strands under tension initially satisfy the demand for more strain chiefly by rotating until they can no longer do so and later resort chiefly to stretching at higher strains.

As discussed in the previous section, the computed shear modulus G displays the scaling $G\sim {\sigma }^{3/2}$, predicted from affine stretching of percolating paths, only at very large strains. At the onset of nonlinear behavior the modulus is still dominated by bending forces and for intermediate stresses we observe a much stronger increase of G with applied stress than expected. Both phenomena are not understood. In the following subsection, we show that generic nonlinear elasticity can account for the onset of nonlinear behavior. In the next subsection, we discuss a simple model for the orientation of fibers at high strain to explain the observed strong shear stiffening in the intermediate regime.

3.2.1. Isotropic nonlinear elasticity

When the deformation is no longer small, one has to take into account nonlinear contributions to the free energy. To analyze the onset of nonlinearities and the weakly nonlinear regime, we consider the strain energy density function of an isotropic medium, expanded in ${\boldsymbol{\eta }}$ up to cubic order:

$$\begin{eqnarray}f({\boldsymbol{\eta }}) & = & \displaystyle \frac{\lambda }{2}{(\mathrm{tr}{\boldsymbol{\eta }})}^{2}+\mu \mathrm{tr}{{\boldsymbol{\eta }}}^{2} & & +{c}_{1}\mathrm{tr}{{\boldsymbol{\eta }}}^{3}+{c}_{2}(\mathrm{tr}{\boldsymbol{\eta }})\mathrm{tr}{{\boldsymbol{\eta }}}^{2}+{c}_{3}{(\mathrm{tr}{\boldsymbol{\eta }})}^{3}.\end{eqnarray} \tag{ 9 }$$

λ and μ are the usual Láme constants from the linear theory. (Note: in this article, the symbols μ and G are often used interchangeably to denote the shear modulus.) At cubic order there are three invariants of ${\boldsymbol{\eta }}$ corresponding to three nonlinear elastic constants: ${c}_{1},{c}_{2}$ and c₃. The Cauchy stress tensor follows from equation (9) by differentiation:

$$\begin{eqnarray}&&{\sigma }_{\alpha \beta }=\displaystyle \frac{1}{\mathrm{det}\;{\boldsymbol{\Lambda }}}\ \displaystyle \frac{\partial {f}_{\phantom{\alpha \gamma }}}{\partial {\Lambda }_{\alpha \gamma }}{\Lambda }_{\gamma \beta }^{{\rm{T}}}.\end{eqnarray} \tag{ 10 }$$

We consider three simple cases: simple shear, pure shear and uniaxial extension, performing simulations in all cases with our network models, in order to determine the values of the model's elastic constants. For simplicity we always set the strain in the z-direction to zero.

Simple shear corresponds to a deformation gradient and corresponding strain tensor:

$$\begin{eqnarray}{\boldsymbol{\Lambda }}=\left(\begin{array}{lll}1 & \gamma & 0\\ 0 & 1 & 0\\ 0 & 0 & 1\end{array}\right)\quad \mathrm{and}\quad {\boldsymbol{\eta }}=\displaystyle \frac{1}{2}\left(\begin{array}{lll}0 & \gamma & 0\\ \gamma & {\gamma }^{2} & 0\\ 0 & 0 & 0\end{array}\right)\end{eqnarray} \tag{ 11 }$$

respectively. The resulting non-zero stresses are

$$\begin{eqnarray}&&{\sigma }_{{xy}}(\gamma )=\mu \gamma +{\mathcal{O}}({\gamma }^{3}),\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{\sigma }_{{xx}}(\gamma )=\displaystyle \frac{{\gamma }^{2}}{2}(\lambda +{\chi }_{1}+4\mu )+{\mathcal{O}}({\gamma }^{4}),\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{\sigma }_{{yy}}(\gamma )=\displaystyle \frac{{\gamma }^{2}}{2}(\lambda +{\chi }_{1}+2\mu )+{\mathcal{O}}({\gamma }^{4}),\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{\sigma }_{{zz}}(\gamma )=\displaystyle \frac{{\gamma }^{2}}{2}(\lambda +{c}_{2})+{\mathcal{O}}({\gamma }^{4}),\end{eqnarray} \tag{ 15 }$$

with ${\chi }_{1}={c}_{2}+3{c}_{1}/2$. Using these equations, the values of μ, $\lambda +{\chi }_{1}$, and $\lambda +{c}_{2}$ may be determined to fit the simulational data for small strains. Later, through uniaxial deformation, the the other two parameters of the model will be determined. Note from above that the relation $({\sigma }_{{xx}}^{\prime\prime }-{\sigma }_{{yy}}^{\prime\prime })/(2{\sigma }_{{xy}}^{\prime })=1$ holds, thus one may obtain a rough estimate of the range of γ within which this weakly nonlinear model is valid for the simulational model (see inset of figure 7). Resolving the stress tensor along the principal stretch directions, we obtain

$$\begin{eqnarray}&&{\sigma }_{\parallel }=\mu \gamma +\displaystyle \frac{{\gamma }^{2}}{2}(\lambda +{\chi }_{1}+3\mu )+{\mathcal{O}}({\gamma }^{3}),\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&{\sigma }_{\perp }=-\mu \gamma +\displaystyle \frac{{\gamma }^{2}}{2}(\lambda +{\chi }_{1}+3\mu )+{\mathcal{O}}({\gamma }^{3}),\end{eqnarray} \tag{ 17 }$$

for the normal stresses in the tensile and compressive principal stretch directions respectively. In figure 7 we compare data of the simulations for ${\sigma }_{\parallel }$ and ${\sigma }_{\perp }$ to the simple nonlinear model equation (9). The results are seen to agree at least qualitatively; in particular the onset of nonlinearity at ${\sigma }_{{xy}}(\gamma )={\sigma }^{*}$ is seen to coincide with a change in sign of ${\sigma }_{\perp }$, which is compressional only for small stresses and becomes tensile in the nonlinear region. Intuitively, the sign reversal of ${\sigma }_{\perp }$ suggests a tendency for the network to expand at small strains and contract at high strains for uniaxial extension, and this behavior is indeed observed in the next paragraph (see figure 8 and equation (20)). Even though demonstrating this ultimately involved the computation of an additional constant, namely c₃, we note here that the form of equation (20) demonstrates that all that is required is for c₃ to be positive, a property already apparent from the truncated strain energy density expansion in equation (9).

Zoom In Zoom Out Reset image size

For uniaxial strain we apply an extensional or compressional force in the x-direction, but allow the system to relax freely in the y-direction. The nonlinear strain is given by

$$\begin{eqnarray}{\boldsymbol{\eta }}=\left(\begin{array}{lll}\gamma & 0 & 0\\ 0 & -\nu \gamma & 0\\ 0 & 0 & 0\end{array}\right)\end{eqnarray} \tag{ 18 }$$

with $\nu =\nu (\gamma )$ determined from the condition of a free surface, that is, ${\sigma }_{{yy}}=0$. During such a deformation, the volume V of the elastic material sample changes according to

$$\begin{eqnarray}&&{(V/\mathring{V})}^{2}=(1+\gamma )(1-\nu (\gamma )\gamma ),\end{eqnarray} \tag{ 19 }$$

where $\mathring{V}$ is the initial volume in the stress-free state. Within linear elasticity, the volume change is determined by the 2D Poisson ratio $\nu (\gamma =0)={\nu }_{2{\rm{D}}}=\frac{\lambda }{\lambda +2\mu }$ (a constant). And from the simulations, we observe that ${\nu }_{2{\rm{D}}}\lt 1$ (see figure 8) indicating an increase in volume for small uniaxial extension, and a decrease in volume for small uniaxial compression. For higher strains, it turns out that the nonlinear terms tend to decrease the volume (see inset of figure 8). In particular, we find

$$\begin{eqnarray}&&{(V/\mathring{V})}^{2}=1+\gamma (1-{\nu }_{2{\rm{D}}})-{\gamma }^{2}(\rho +{\nu }_{2{\rm{D}}})+{\mathcal{O}}({\gamma }^{3}),\end{eqnarray} \tag{ 20 }$$

where $\rho =\frac{2{\mu }^{2}{\chi }_{2}+{\lambda }^{2}{\chi }_{1}}{{(\lambda +2\mu )}^{3}}$ and ${\chi }_{2}=3{c}_{3}+{c}_{2}$. Note that we have only kept the lowest order nonlinear terms in equation (20). Higher order terms are expected to limit the volume change at large strain.

Having already determined the value of μ through simple shear data, it is now possible, using the linear term of the equation above, to determine the value of λ from data on uniaxial extension. In this way we obtain the (3D) Poisson ratio $\nu =\frac{\lambda }{2(\lambda +\mu )}$ as $\nu \approx 0.28$,¹ indicating a ratio of bulk to shear modulus of $\approx 1.93$. Combining this with the other results from simple shear, one obtains the values of c₂ and ${\chi }_{1}$. Furthermore, using the quadratic term of the equation above, the value of ρ may be determined from the data. This in turn enables the determination of c₁ and c₃, thus completing the entire set of best-fit values for the elastic parameters of the model. As a check, we perform a pure shear deformation to observe how well the determined values fit the data.

Pure shear corresponds to a deformation gradient

$$\begin{eqnarray}{\boldsymbol{\Lambda }}=\left(\begin{array}{lll}1+\varepsilon & 0 & 0\\ 0 & {(1+\varepsilon )}^{-1} & 0\\ 0 & 0 & 1\end{array}\right)\end{eqnarray} \tag{ 21 }$$

and strain tensor

$$\begin{eqnarray}\eta =\left(\begin{array}{lll}\varepsilon +\displaystyle \frac{1}{2}{\varepsilon }^{2} & 0 & 0\\ 0 & -\varepsilon +\displaystyle \frac{3}{2}{\varepsilon }^{2} & 0\\ 0 & 0 & 0\end{array}\right),\end{eqnarray} \tag{ 22 }$$

which involves strains in the x- and y-directions such that $\mathrm{det}\;{\boldsymbol{\Lambda }}=1$ up to and including quadratic terms in . The resulting non-zero stresses are explicitly given by

$$\begin{eqnarray}&&{\sigma }_{{xx}}=2\mu \varepsilon +{\varepsilon }^{2}(2\lambda +2{\chi }_{1}+5\mu )+{\mathcal{O}}({\varepsilon }^{3}),\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&{\sigma }_{{yy}}=-2\mu \varepsilon +{\varepsilon }^{2}(2\lambda +2{\chi }_{1}+7\mu )+{\mathcal{O}}({\varepsilon }^{3}),\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&{\sigma }_{{zz}}=2{\varepsilon }^{2}(\lambda +{c}_{2})+{\mathcal{O}}({\varepsilon }^{3}).\end{eqnarray} \tag{ 25 }$$

The simulations show that the stress in the compressional (y-) direction of ${\boldsymbol{\eta }}$ is compressive (negative) only for small strains but then becomes tensile (positive) for large strains, which is in fair agreement with the weakly nonlinear model (see figure 9). A similar behavior is seen for the stress component in the compressional direction for simple shear (see figure 7).

Zoom In Zoom Out Reset image size

In summary, the weakly nonlinear model demonstrates that (1) it is possible to obtain a volume-change sign-reversal for large enough uniaxial strains, even for positive c₁, c₂, and c₃, and (2) that once the relevant parameters of the model have been determined through fitting with data, it is also able to predict changes of the sign of stress in other deformation types, such as the normal stress in pure shear, for example.

As the nonlinearities become more prominent, the above expansion ceases to be applicable. Besides the strength of the deformation which becomes too large to be treated perturbatively, it is the assumption of isotropy which breaks down at large strain. To illustrate this point, we highlight in figure 10 the force chains under the highest tension (in red) as well as the strands which sustain high compressional forces (in green). For such large shear deformations the stretching of single fibers gives rise to the known scaling, $G\sim {\sigma }^{3/2}$. In this section, we propose a simple toy model to account for the comparatively faster rate of stiffening observed before the asymptotic scaling is reached. The basic mechanism under consideration is the rotation of fibers or segments into the tensile direction. In the simplest version we model the strands as inextensible, but extensibility can be taken into account easily.

For simplicity, we will focus instead on the scaling dependence of the Young's modulus $E\sim {\sigma }_{\parallel }^{\alpha }$ on uniaxial stress ${\sigma }_{\parallel }$ at intermediate values of stress. Here, $\alpha$ is the scaling exponent. Figure 11 shows how in our simulations $\alpha ={\rm{d}}\;\mathrm{log}\;E/{\rm{d}}\;\mathrm{log}\;{\sigma }_{\parallel }$ changes with ${\sigma }_{\parallel }$. As expected, at low values of ${\sigma }_{\parallel }$, the scaling is linear, and $\alpha \;\sim \;0$. At intermediate values of ${\sigma }_{\parallel },\alpha$ may rise to a maximum which can be as high as five, before at last falling to a value of $3/2$ which indicates that a few strands (forming a force-chain percolating the network, as in figure 10) are being stretched affinely and are dominating the stress. The non-specific value for the maximum of $\alpha$ hints at an exponential regime at intermediate stresses.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In order for percolating tensile paths to be formed, filaments or parts of filaments (strands) have to align themselves with the tensile direction. Note that not all strands of the network will eventually be aligned with the tensile direction, otherwise the lengths of those network dimensions transverse to the tension would collapse to zero at full alignment, which is not what is observed (see figure 10). It turns out that we can understand the exceptionally strong stiffening in the intermediate nonlinear regime via this mechanism of filament alignment. To this end, consider an inextensible strand of length $b$ oriented at a force-dependent angle $\phi (f)$ to the tensile direction. Changing the pulling force $f\to f+{\rm{d}}$ f gives rise to a differential torque ${\rm{d}}M=-b\mathrm{sin}\phi \ {\rm{d}}f$, which tries to rotate the filament even further, ${\rm{d}}\phi ={\rm{d}}M/\tau$, where $\tau$ is an effective material-constant that, for simplicity, we take to be constant and independent of the pulling force. It plays the role of the (bending dominated) linear modulus $E(\sigma \to 0)$.

Thus, we obtain the following differential equation:

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\phi }{{\rm{d}}f}=-\frac{b}{\tau }\mathrm{sin}\phi .\end{eqnarray} \tag{ 26 }$$

which can be integrated to give

$$\begin{eqnarray}&&\mathrm{tan}\frac{\phi (f,\ {\phi }_{0})}{2}={{\rm{e}}}^{-{bf}/\tau }\mathrm{tan}\frac{{\phi }_{0}}{2},\end{eqnarray} \tag{ 27 }$$

where ${\phi }_{0}$ is the initial orientation of the strand. The crucial aspect is the exponential dependence on force $f$. The strain $\gamma$ can be computed from the average $\langle \mathrm{cos}\phi \rangle$, where the average is over the initial angles ${\phi }_{0}$ of the strands in the percolating path:

$$\begin{eqnarray}&&\gamma (f)\propto \frac{\langle \mathrm{cos}\phi \rangle -\langle \mathrm{cos}{\phi }_{0}\rangle }{\langle \mathrm{cos}{\phi }_{0}\rangle },\end{eqnarray} \tag{ 28 }$$

For strongly aligned strands, we approximate $\mathrm{tan}\frac{\phi }{2}\sim \frac{\phi }{2}\;{\rm{and}}\;\mathrm{cos}\phi \sim 1-\frac{{\phi }^{2}}{2}$ and compute

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\gamma }{{\rm{d}}f}\sim \frac{4b}{\tau }{{\rm{e}}}^{-2{bf}/\tau }\frac{\langle {\mathrm{tan}}^{2}\frac{{\phi }_{0}}{2}\rangle }{\langle \mathrm{cos}{\phi }_{0}\rangle }\end{eqnarray} \tag{ 29 }$$

resulting in an exponentially increasing modulus $E\sim (\tau /b){{\rm{e}}}^{+2{bf}/\tau }$. Note that the argument here is asymptotic, and one might expect this behavior to carry over to the shear modulus as well. But in simple shear, an additional complication arises, namely, compression is applied along one-dimension transverse to the maximum strain in order to keep $V$ constant. Extending the above argument to include this extra behavior is nontrivial and will be the subject of another article.

Similar exponential stiffening has been reported in [46] in the context of a network of rigid rods crosslinked by flexible linkers. The authors present an alternative explanation that is based on treating the flexible linkers as an effective elastic medium [47].

To check the robustness of the intermediate regime against parameter variations we varied the average strand length $b\equiv \langle {{\ell }}_{{ij}}\rangle$, i.e. the average distance between consecutive crosslinks along filaments. With the overall filament length fixed, a smaller strand length means more crosslinks per filament. This increases the modulus, $E\sim {b}^{-4}$, as can readily be seen in figure 12. Such a variation is in agreement with the theory of non-affine floppy bending modes [14, 22] of wavelength $\sim$ b and amplitude $\sim \gamma$ independent of $b$. The resulting modulus then is given by ${E}_{0}\sim (\kappa /{b}^{3})({{\ell }}_{{\rm{c}}}/b){N}_{\mathrm{fil}}/V\sim {b}^{-4}$, where the first factor is the spring constant of the bending mode, the second accounts for the number of strands per filament, and ${N}_{\mathrm{fil}}$ counts the number of filaments in the volume $V$. For the nonlinear response we observe that the onset of intermediate exponential stiffening is shifted to smaller stresses in the weakly crosslinked networks. This indicates that reorientation effects are stronger in weakly crosslinked networks, thus filaments are more easily rotated when not highly crosslinked to its neighbors.

Zoom In Zoom Out Reset image size

In appendix B we generalize to $n$-body interactions the often-quoted local stress tensor for pairwise interactions. The point-wise local stress for the networks of this article (which have only two- and three-body interactions) is therefore

$$\begin{eqnarray}&&{\boldsymbol{\sigma }}({\bf{r}})\equiv {{\boldsymbol{\sigma }}}^{(2)}({\bf{r}})+{{\boldsymbol{\sigma }}}^{(3)}({\bf{r}}),\end{eqnarray} \tag{ 30 }$$

where

$$\begin{eqnarray}&&{{\boldsymbol{\sigma }}}^{(2)}({\bf{r}})\equiv \displaystyle \sum _{\langle {ij}\rangle }^{N}{{\bf{f}}}_{i,j}^{(2)}\otimes {{\bf{r}}}_{{ij}}\;{\displaystyle \int }_{0}^{1}{\rm{d}}s\;{\phi }_{w}({\bf{r}}-[{{\bf{r}}}_{i}+s{{\bf{r}}}_{{ij}}])\end{eqnarray} \tag{ 31 }$$

and

$$\begin{eqnarray}&&{{\boldsymbol{\sigma }}}^{(3)}({\bf{r}})\equiv \frac{1}{3}\displaystyle \sum _{[{ijk}]}^{N}{{\bf{f}}}_{i,\{j,k\}}^{(3)}\otimes {{\bf{r}}}_{{ij}}\;{\displaystyle \int }_{0}^{1}{\rm{d}}s\;{\phi }_{w}({\bf{r}}-[{{\bf{r}}}_{i}+s{{\bf{r}}}_{{ij}}]).\end{eqnarray} \tag{ 32 }$$

The sum in equation (31) is over all filament segments $i\leftrightarrow$ j (that is, the edges), and the sum in equation (32) is over all bending angles i $\;\leftrightarrow \;j\leftrightarrow k$ but allows for all 3 permutations of the dummy indices $i,j,k{\text{}};{{\bf{f}}}_{i,j}^{(2)}\;{\rm{and}}\;{{\bf{f}}}_{i,\{j,k\}}^{(3)}$ are two-body and three-body forces defined in appendix B. The local linear strain is defined as

$$\begin{eqnarray}&&{\varepsilon }_{\alpha \beta }({\bf{r}})=\frac{1}{2}\left(\frac{\partial {u}_{\alpha }({\bf{r}})}{\partial {r}_{\beta }}+\frac{\partial {u}_{\beta }({\bf{r}})}{\partial {r}_{\alpha }}\right),\end{eqnarray} \tag{ 33 }$$

where the point-wise coarse-grained local linear displacement is

$$\begin{eqnarray}&&{u}_{\alpha }({\bf{r}})=\frac{{\displaystyle \sum }_{i}{u}_{i\alpha }\;{\phi }_{w}({\bf{r}}-{{\bf{r}}}_{i})}{{\displaystyle \sum }_{i}{\phi }_{w}({\bf{r}}-{{\bf{r}}}_{i})}.\end{eqnarray} \tag{ 34 }$$

Here ${{\bf{u}}}_{i}$ is the displacement of the $i$ th crosslink after a (global) deformation specified by, say $[{\Lambda }_{\alpha \beta }]$, has been applied to the system-boundary and the network subsequently minimized.

In figure 13 we show the distribution of local stress-component ${\sigma }_{{xx}}$ in a reference state without global prestress. These data demonstrate that on a local scale the network is strongly prestressed, giving rise to a broad distribution of local stresses. The width of the distribution shrinks as a function of growing coarse-graining length as one would expect.

Zoom In Zoom Out Reset image size

One final remark concerns the range of coarse-graining lengths accessible to our simulational data. The computation of local stresses and strains requires at least two mass points. Hence the coarse-graining scale should be larger than the distance between two 'particles', or crosslinks. This limits us to $w\geqslant \langle {{\ell }}_{{ij}}\rangle$.

We have seen in the previous section that the network supports considerable stresses on a local scale. In order to compute the local elastic constants, we thus need to consider a prestressed reference state. In the following, we will ignore crosslink positions, assuming that they are always at their equilibrium positions, and denote a reference state simply as $\mathring{{\boldsymbol{\Lambda }}}$ with a generally non-vanishing stress $\mathring{\sigma }$ and volume $\mathring{V}$. We will apply an additional deformation ${\boldsymbol{\Lambda }}$ to compute the elastic response of the reference system. (Note that an observable $a$ in the reference state is denoted by $\mathring{a}$.)

The stiffness of any state of the system is described most generally by a fourth-order modulus tensor ${C}_{\alpha \beta \gamma \delta }$ defined by [50]:

$$\begin{eqnarray}&&\frac{F-\mathring{F}}{\mathring{V}}={\mathring{\sigma }}_{\alpha \beta }\;{\eta }_{\alpha \beta }+\frac{1}{2}{C}_{\alpha \beta \gamma \delta }\;{\eta }_{\alpha \beta }\;{\eta }_{\gamma \delta }+O({[{\eta }_{\alpha \beta }]}^{3}),\quad {\rm{or}}\end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray}&&{S}_{\alpha \beta }={\mathring{\sigma }}_{\alpha \beta }+{C}_{\alpha \beta \gamma \delta }\;{\eta }_{\gamma \delta }+O({[{\eta }_{\alpha \beta }]}^{2}),\end{eqnarray} \tag{ 36 }$$

where ${S}_{\alpha \beta }\equiv {\mathring{V}}^{-1}\partial F/\partial {\eta }_{\alpha \beta }$ is the second Piola–Kirchoff stress tensor written as a function of ${\eta }_{\alpha \beta }$. Consequently, ${C}_{\alpha \beta \gamma \delta }$ possesses the following symmetries: ${C}_{\alpha \beta \gamma \delta }={C}_{\beta \alpha \gamma \delta }={C}_{\alpha \beta \delta \gamma }={C}_{\gamma \delta \alpha \beta }$ amounting to 21 independent components in general. In simulations of random polymer networks, it is common to assume elastic isotropy of the network and accordingly compute only 2 elastic constants: namely the bulk modulus $K$ and the shear modulus $G$ of the network. However, this assumption is true only in the thermodynamic limit for unstressed networks. Thus in our investigations of local moduli, we choose to maintain the most general symmetry possible in continuum elasticity, namely, the triclinic symmetry, and observe how this symmetry tends to isotropy (or some other symmetry) as the system-size or coarse-graining scale increases (while keeping the intensive properties of the system constant).

It has been shown [50] that the Cauchy stress in the new reference configuration ${\boldsymbol{\Lambda }}\mathring{{\boldsymbol{\Lambda }}}$ is to first-order given by

$$\begin{eqnarray*}&&{\sigma }_{\alpha \beta }={\mathring{\sigma }}_{\alpha \beta }+({C}_{\alpha \beta \gamma \delta }+{\mathring{\sigma }}_{\beta \delta }{\delta }_{\alpha \gamma }+{\mathring{\sigma }}_{\alpha \delta }{\delta }_{\beta \gamma }-{\mathring{\sigma }}_{\alpha \beta }{\delta }_{\gamma \delta }){u}_{\gamma \delta }.\end{eqnarray*}$$

Writing ${u}_{\alpha \beta }={\varepsilon }_{\alpha \beta }+{\omega }_{\alpha \beta }$, where ${\varepsilon }_{\alpha \beta }\;{\rm{and}}\;{\omega }_{\alpha \beta }$ are the symmetric and anti-symmetric parts of ${u}_{\alpha \beta }$ respectively, it follows that

$$\begin{eqnarray}&&{\boldsymbol{\sigma }}-\mathring{{\boldsymbol{\sigma }}}+\mathring{{\boldsymbol{\sigma }}}\;\mathrm{Tr}\;{\boldsymbol{\varepsilon }}-(\mathring{{\boldsymbol{\sigma }}}\cdot {\boldsymbol{\varepsilon }}+{\boldsymbol{\varepsilon }}\cdot \mathring{{\boldsymbol{\sigma }}})+(\mathring{{\boldsymbol{\sigma }}}\cdot {\boldsymbol{\omega }}-{\boldsymbol{\omega }}\cdot \mathring{{\boldsymbol{\sigma }}})=\underline{{\bf{C}}}:{\boldsymbol{\varepsilon }},\end{eqnarray} \tag{ 37 }$$

in tensor notation. (If there is little pre-stress [${\mathring{\sigma }}_{\alpha \beta }\;\approx \;0$], all the terms after the first term on the left-hand side of equation (37) may be taken to vanish. For this reason, these terms have often been neglected in simulations of glasses [12], where the so-called 'quenched stresses' are small.)

Our starting point is equation (37) which we apply locally at any observation point ${\bf{r}}$. The local coarse grained displacement field ${u}_{\alpha }({\bf{r}})$ is in general non-affine and given, to linear order, by equation (34). For small global strains, the local linear strain ${\varepsilon }_{\alpha \beta }({\bf{r}})$ due to ${u}_{\alpha }({\bf{r}})$ is expected to approach the global strain ${\eta }_{\alpha \beta }$ as the coarse-graining scale $w$ increases, approaching the system size and beyond (see figure 13). Indeed, since ${\sigma }_{\alpha \beta }({\bf{r}})$, ${\mathring{\sigma }}_{\alpha \beta }({\bf{r}}),{\omega }_{\alpha \beta }({\bf{r}})$, and ${\varepsilon }_{\alpha \beta }({\bf{r}})$ are all readily computable from the microscopic coordinates $\{{{\bf{r}}}_{i}\}$ via equations (31), (33) and (34), one may use equation (37) to solve for ${C}_{\alpha \beta \gamma \delta }({\bf{r}})$. The method proceeds as follows: for any reference state $(\{\mathring{{{\bf{r}}}_{i}}\},{\mathring{\Lambda }}_{\alpha \beta })$ of the network, we apply six small homogeneous deformations ${\Lambda }_{\alpha \beta }^{(d)}(d=1,\ldots ,6)$ along the six possible linearly independent directions of the strain tensor, and minimize with respect to the crosslink coordinates in order to obtain six different states $(\{{{\bf{r}}}_{i}^{(d)}\},{\Lambda }_{\alpha \gamma }^{(d)}{\mathring{\Lambda }}_{\gamma \beta })$. That is, for each deformation we set ${\Lambda }_{\alpha \beta }^{(d)}={\delta }_{\alpha \beta }+{e}_{\alpha \beta }$ where all the components ${e}_{\alpha \beta }$ are zero except for one which is very small in value (see table 2 for detailed information).

Table 2.  List of the six linearly independent deformations used in this article to compute the elastic modulus tensor. The other three components of ${\Lambda }_{\alpha \beta }^{(d)}$, that is, ${\Lambda }_{21}^{(d)},{\Lambda }_{32}^{(d)}$, and ${\Lambda }_{13}^{(d)}$ were set to 0 for all $d$.

Deformation d	${\Lambda }_{{xx}}^{(d)}$	${\Lambda }_{{yy}}^{(d)}$	${\Lambda }_{{zz}}^{(d)}$	${\Lambda }_{{yz}}^{(d)}$	${\Lambda }_{{zx}}^{(d)}$	${\Lambda }_{{xy}}^{(d)}$
1	1	1	1	0	0	${e}_{12}^{(1)}$
2	1	1	1	0	${e}_{31}^{(2)}$	0
3	1	1	1	${e}_{23}^{(3)}$	0	0
4	1	1	$1+{e}_{33}^{(4)}$	0	0	0
5	1	$1+{e}_{22}^{(5)}$	1	0	0	0
6	$1+{e}_{11}^{(6)}$	1	1	0	0	0

At any observation point ${\bf{r}}$, the local stress ${\mathring{\sigma }}_{\alpha \beta }({\bf{r}})$ before deformation and the local quantities ${\sigma }_{\alpha \beta }^{(d)}({\bf{r}})$, ${\omega }_{\alpha \beta }({\bf{r}})$, and ${\varepsilon }_{\alpha \beta }({\bf{r}})$ after each of the six deformations are computed. The local modulus at any observation point ${\bf{r}}$ may then be computed as follows: using equation (37) six times (one for each deformation), an overdetermined system of 36 linear equations in 21 unknowns (the components of the local modulus tensor) is obtained. This number of equations should be reduced, through careful selection (see, for example appendix C), to 21 linearly independent equations and solved. The ${C}_{\alpha \beta \gamma \delta }({\bf{r}})$ may be computed at observation points ${\bf{r}}$ located at the nodes of a regularly spaced grid of points within the system box. A similar procedure may be employed to determine the global modulus tensor defined in the previous subsection: in this case, it is perhaps simpler to directly use equation (36), instead of equation (37), with the replacements ${S}_{\alpha \beta }\to {S}_{\alpha \beta }^{(d)}\;{\rm{and}}\;{\eta }_{\alpha \beta }\to {\eta }_{\alpha \beta }^{(d)}=\displaystyle \frac{1}{2}\left({\Lambda }_{\alpha \gamma }^{(d){\rm{T}}}{\Lambda }_{\gamma \beta }^{(d)}-{\delta }_{\alpha \beta }\right)$, since it is known [51] that

$$\begin{eqnarray}&&{S}_{\alpha \beta }^{(d)}=\mathrm{det}[{\Lambda }_{\alpha \beta }^{(d)}]\;{\Lambda }_{\alpha \mu }^{(d)-1}\;{\sigma }_{\mu \nu }^{(d)}\;{\Lambda }_{\nu \beta }^{(d)-{\rm{T}}},\end{eqnarray} \tag{ 38 }$$

where the global Cauchy stresses ${\sigma }_{\mu \nu }^{(d)}\;\mathrm{and}\;{\mathring{\sigma }}_{\mu \nu }$ may be obtained using equation (6).

A particular linear combination of the 21 elastic constants gives the so-called 'best isotropic approximation' [52, 53] ${C}_{\alpha \beta \gamma \delta }^{\mathrm{iso}}=\lambda {\delta }_{\alpha \beta }{\delta }_{\gamma \delta }+2\mu ({\delta }_{\alpha \gamma }{\delta }_{\beta \delta }+{\delta }_{\alpha \delta }{\delta }_{\beta \gamma })$ to the modulus tensor ${C}_{\alpha \beta \gamma \delta }$, yielding an approximate shear modulus $\mu$ and bulk modulus $K=\lambda +2\mu /3$, whose expressions are given in appendix D. There, a quantity ${I}_{{\rm{A}}}$, called the 'index of anisotropy', measuring how far the system is from isotropy, is also defined [52]. ${I}_{{\rm{A}}}$ ranges from 0 to 1, and it vanishes for a completely isotropic medium. Note that in our computational models even the global moduli are anisotropic due to finite system size.

In figure 14 we show the average index of anisotropy as a function of coarse-graining length for two sizes of networks. The figure demonstrates that the network is isotropic in the thermodynamic limit, but anisotropic on a local scale and in small networks. For the larger network ($N=64\;000$), the index of anisotropy ${I}_{{\rm{A}}}\;\sim \;0.25$ for a coarse-graining scale comparable to the mean distance between crosslinks and ${I}_{{\rm{A}}}\sim$ 0.1 for a coarse-graining scale comparable to the filament length. These are the length-scales, where we expect to see fluctuations in the moduli, as discussed in the next paragraph.

Zoom In Zoom Out Reset image size

The distribution of local bulk modulus is shown in figure 15 for five values of coarse-graining length. As expected, the distributions become sharp as the coarse-graining length approaches system size. This limiting behavior is quantified by plotting the variance relative to the mean as a function of $w$, as done in the inset of figures 15.

Zoom In Zoom Out Reset image size

The fluctuations of the shear modulus are strongly correlated with the fluctuations in the density. The correlation factor is

$$\begin{eqnarray}&&C\equiv \frac{\langle \delta \rho ({r}_{0})\delta \mu ({r}_{0})\rangle }{{\Delta }_{\rho }{\Delta }_{\mu }}\sim 0.8\end{eqnarray} \tag{ 39 }$$

for all coarse-graining lengths. (Here ${\Delta }_{\rho }^{2}=\langle \delta {\rho }^{2}({r}_{0})\rangle \;{\rm{and}}\;{\Delta }_{\mu }^{2}=\langle \delta {\mu }^{2}({r}_{0})\rangle$.)

Next, we check the dependence of the elastic heterogeneity on degree of crosslinking. Reducing the number of crosslinks by roughly $40\%$ from the initial value (six crosslinks per filament), the network loses its stiffness towards shear deformations. The variance is similarly decreasing in this range of ${N}_{\times }$ so that the variance relative to the mean is approximately independent of crosslink concentration. On the other hand, the elastic heterogeneities increase with increasing shear stress as can be seen in figure 16, where we plot the variance relative to the mean a a function of stress

$$\begin{eqnarray}&&{\Delta }_{\tilde{\mu }}({\sigma }_{{xy}})=\sqrt{{\Delta }_{\mu }^{2}}/\mu \end{eqnarray} \tag{ 40 }$$

for several coarse graining lengths. Naturally the variance strongly decreases with coarse-graining length, as the local modulus approaches the global value.

Zoom In Zoom Out Reset image size

We expect the nonaffine deformations to be correlated with the fluctuations in the elastic constants, such that hard regions suppress the deformation below the affine value, whereas soft regions allow for deformations larger than the affine value. However, we only find a very small effect, when quantifying these correlations.

Defining hard (soft) regions as those whose local modulus exceeds (falls below) the mean, one may alternatively look at distributions of nonaffine displacements for the soft and hard regions separately as suggested in [12]. The normalized distributions are shown in figure 17. Clearly the large displacements are significantly more probable in soft regions.

Zoom In Zoom Out Reset image size