Micromechanical theory of strain stiffening of biopolymer networks

Filamentous biomaterials such as fibrin or collagen networks exhibit an enormous stiffening of their elastic moduli upon large deformations. This pronounced nonlinear behavior stems from a significant separation between the stiffnesses scales associated with bending versus stretching the material's constituent elements. Here we study a simple model of such materials, floppy networks of hinged rigid bars embedded in an elastic matrix, in which the effective ratio of bending to stretching stiffnesses vanishes identically. We introduce a theoretical framework and build upon it to construct a numerical method with which the model's micro- and macromechanics can be carefully studied. Our model, numerical method and theoretical framework allow us to robustly observe and fully understand the critical properties of the athermal strain-stiffening transition that underlies the nonlinear mechanical response of a broad class of biomaterials.

Figure 1

The strain-stiffening transition is conventionally studied by deforming simple model networks in which the characteristic energies of bending and stretching interactions are well separated. We show the shear modulus $G$ of such a model (see text for details) as a function of shear strain, for different values of the ratio of stretching to bending stiffness $\tilde{\mu }$ (see Appendix pp1). The curve with circle symbols corresponds to the shear modulus of the model introduced in this work, which corresponds to the limit $\tilde{\mu }\to \infty$. By studying this limit, our model allows us to reveal the critical properties of the strain-stiffening transition that underlies the nonlinear mechanics of biopolymer networks.

Figure 2

Examples of the floppy frames employed for our study: (a) a packing-derived network obtained by diluting the edges of a packing's contact network, and (b) a diluted off-lattice honeycomb network with disorder of the node positions ($d_{\text{max}}=0.5$). See Appendix pp4 for a detailed description of the protocol used to construct these frames.

Figure 3

For undeformed networks ($N=12800$) of different coordinations $z$ we plot (a) the shear modulus $G$ (black triangles) and bulk modulus $K$ (gray circles), and (b) the characteristic scale of nonaffine velocities squared ${\dot{x}}^2\equiv \langle \dot{\mathit{x}}|\dot{\mathit{x}}\rangle /N$ for both shear and dilation.

Figure 4

For a honeycomb lattice ($N=3600$, $z=3.0$, $d_{\text{max}}=0.5$), (a) the approximation of the edge forces as given by Eq. (32) and denoted here by $\tilde{\tau }$ is tested by plotting the relative magnitude squared of their difference with the actual forces $\tau$ [see Eq. (14)], as a function of the distance to the critical strain $\delta \gamma$. (b) Validation of relation (34); for the same deformed network we plot $\langle {\varphi }_0\left|{\left(\mathcal{S}{\mathcal{A}}^{-1}{\mathcal{S}}^T\right)}^{-1}\right|{\varphi }_0\rangle$ against $\tau /{\omega }_0^2$, to find a linear relation.

Figure 5

Elastic moduli of deformed, elastically embedded floppy frames of rigid edges, as a function of the distance to the critical strain $\delta \gamma \equiv {\gamma }_c-\gamma$. Shown are the decomposition of the moduli's constituent terms, grouped by their scaling behavior with respect to strain, as expressed by Eqs. (40, 41, 42), and predicted by our theory. (a) The shear modulus $G$ of a sheared disordered honeycomb lattice ($N=3600$, $z=3.0$, $d_{\text{max}}=0.5$) with a bond-bending interactions as the embedding elastic energy; (b) the bulk modulus $K$ of a packing-derived disordered frame with $z=3.0$ embedded in an Hookean-spring elastic network.

Figure 6

(a) The distribution of eigenfrequencies $D\left(\omega \right)$ of the operator $\mathcal{S}{\mathcal{S}}^T$, calculated for sheared networks of $N=4096$ nodes with coordinations $z=3.6$ and $z=3.8$, at different stages of the deformation. Each curve shows the distribution of eigenfrequencies binned over 100 independent realizations. The distributions feature a plateau above an onset frequency ${\omega }^{*}\sim \delta z\equiv z_c-z$, with $z_c$ the Maxwell threshold. Upon deformation, a single, lowest frequency mode ${\omega }_0$ per realization escapes from the plateau and vanishes as ${\omega }_0\sim \sqrt{\delta \gamma }$ upon approaching the critical strain, as shown in panels (b) and (c). The mean (over realizations) frequency of the lowest mode ${\omega }_0$ is indicated with a vertical line, and its standard deviation is indicated by horizontal lines. (b) The lowest eigenvalue ${\omega }_0^2$ of the operator $\mathcal{S}{\mathcal{S}}^T$ as function of the distance to the critical point ${\gamma }_c$ for packing-derived networks ($N=1600$) of different coordination number $\delta z$. Each data point represents the median over 20 realization. (c) The same data as presented in panel (b), recast into the scaling form given by Eq. (43).

Figure 7

(a) The shear modulus $G$ as function of the distance to the critical point ${\gamma }_c$ for packing-derived networks of different coordination ($N=1600$) number $\delta z$. Each curve is the result of the median of 20 realisation. (b) The shear modulus $G$ rescaled by $\delta z$ to obtain a collapse of the same data. Panels (c) and (d) are the same as panels (a) and (b) but for the bulk modulus $K$ of networks under expansive deformation.

Figure 8

(a) Decaying correlations $C\left(r\right)$ as function distance between pairs for a single network realization $(N=40000)$ under increasing amount of deformation reveals the increasing length scale in the system. (b) The collapse of the correlation functions by the length scale $l_r\sim 1/\sqrt{\delta \gamma }$.

Figure 9

Example of a plastic instability in a packing-derived network ($N=1600$, $z=3.8$). (a) The shear modulus shows a characteristic dip where the network softens over a small strain interval. (b) The shear stress shows a sign of the instability as well.