Stress-stabilized sub-isostatic fiber networks in a rope-like limit

TL;DR

This paper investigates how external stress stabilizes sub-isostatic fiber networks with rope-like interactions, revealing a power-law relationship between shear modulus and stress, and showing prestress effects on critical strain.

Contribution

It demonstrates that external stress can stabilize sub-isostatic networks with tensile interactions, extending understanding of their mechanical phase behavior without bending rigidity.

Findings

External stress stabilizes rope-like fiber networks.

Shear modulus scales with stress via a non-mean-field power law.

Prestress shifts the critical strain to lower values.

Abstract

The mechanics of disordered fibrous networks such as those that make up the extracellular matrix are strongly dependent on the local connectivity or coordination number. For biopolymer networks this coordination number is typically between three and four. Such networks are sub-isostatic and linearly unstable to deformation with only central force interactions, but exhibit a mechanical phase transition between floppy and rigid states under strain. Introducing weak bending interactions stabilizes these networks and suppresses the critical signatures of this transition. We show that applying external stress can also stabilize sub-isostatic networks with only tensile central force interactions, i.e., a rope-like potential. Moreover, we find that the linear shear modulus shows a power law scaling with the external normal stress, with a non-mean-field exponent. For networks with finite bending rigidity, we find that the critical stain shifts to lower values under prestress.

Equations

$$\mathcal{H}=\frac{\mu }{2}\sum _{⟨ij⟩}\frac{(l_{ij}-l_{ij,0}{)}^2}{l_{ij,0}}+\frac{\kappa }{2}\sum _{⟨ijk⟩}\frac{({\theta }_{ijk}-{\theta }_{ijk,0}{)}^2}{\frac{1}{2}(l_{ij,0}+l_{jk,0})}$$

$$\mathcal{H}=\frac{\mu }{2}\sum _{⟨ij⟩}\Theta (l_{ij}-l_{ij,0})\frac{(l_{ij}-l_{ij,0}{)}^2}{l_{ij,0}}$$

$${\sigma }_{\alpha \beta }=\frac{1}{2A}\sum _{⟨ij⟩}{f}_{ij,\alpha }{r}_{ij,\beta }$$

$$K\sim {\chi }_P{\sigma }_{\perp },\gamma \le {\gamma }_c$$

$$\delta \Gamma =\frac{⟨\parallel \delta {\mathbf{u}}^{\text{NA}}{\parallel }^2⟩}{l^2\delta {\gamma }^2}$$

Full text

I Introduction

Networks of biopolymers are ubiquitous in biological systems involved in structural and mechanical stability. Examples include cross-linked cortical actin in the cytoskeleton and branched collagenous networks in the extracellular matrix. The underlying local network geometry and the nature of interactions between the constituent fibers play a key role in determining the stability of these networks. Typically, these networks have average coordination or connectivity between three and four, corresponding to branched or cross-linked geometries, respectively. As shown by Maxwell, the isostatic threshold connectivity for linear stability of an interconnected mechanical structure of simple springs is twice the dimensionality for large number of elements, i.e., $z_{c}=2d$ Maxwell (1864). Based on this argument, biological networks are intrinsically sub-isostatic. Therefore, if the fibers interact only via central forces such as tension and compression, then these networks are unstable with respect to small deformations. Nevertheless, sub-isostatic networks can be rigidified through various stabilizing effects such as strain Alexander (1998), fiber bending interactions Kroy and Frey (1996); Satcher and Dewey (1996); Head et al. (2003); Wilhelm and Frey (2003); Das et al. (2007); Broedersz et al. (2011), active stresses Mizuno et al. (2007); Koenderink et al. (2009); Winer et al. (2009); Jansen et al. (2013); Sheinman et al. (2012a), or thermal fluctuations Plischke and Joós (1998); Dennison et al. (2013); Jaspers et al. (2014) giving rise to a stable linear elastic response.

Biopolymer networks also exhibit striking nonlinear elasticity: with barely a 10% increase in strain, the stiffness increases by almost two orders of magnitude. Such nonlinear mechanics are observed for intracellular cytoskeletal filaments, extracellular fibrin clots, and even whole tissues MacKintosh et al. (1995); Shah and Janmey (1997); Shadwick (1999); Gardel et al. (2004); Storm et al. (2005); Pritchard et al. (2014); Arevalo et al. (2015). The nonlinear mechanics of athermal networks have been described theoretically in terms of a crossover from bending-dominated to a stretching-dominated response Onck et al. (2005); Heussinger and Frey (2006), normal stresses Licup et al. (2015), or strain-controlled critical phenomena Sharma et al. (2016a). In these theoretical approaches, bending rigidity provides stability in the linear regime. It is also known experimentally that polymerization of hydrogels such as collagen and fibrin generally results in prestress van Oosten et al. (2016). In fact, some prestress is almost inevitable as crosslinks form between fibers Huisman et al. (2008). But, it is still not well understood how such prestress, either externally applied or due to internal constraints, affects network stability and nonlinear mechanics Vos et al. (2017). In the case of active stresses, such as by myosin motors in the cytoskeletal or platelet contraction in blood clots, such active prestress can give rise to shear moduli that can exceed the passive shear modulus of the underlying substrate Koenderink et al. (2009); Winer et al. (2009).

Here, in order to investigate the stabilization effect of prestress, we study sub-isostatic rope networks without bending interactions. The elastic response of a rope-like fiber is governed purely by central-force interactions under tension and has vanishing resistance under compression. This represents a minimal model for athermal fibers with zero bending rigidity, for which the Euler buckling threshold vanishes. The external normal stress is applied by either bulk or uniaxial expansion. We show that the linear shear modulus scales as power law with the imposed external normal stress, with a non-mean-field exponent which leads to a divergent susceptibility. This suggests that sub-isostatic rope networks become infinitely susceptible to any stress that invokes fiber stretching modes, including the self-generated normal stresses. We also show that the mechanics of stress-stabilized rope networks can be captured in terms of these normal stresses. Furthermore, by calculating the non-affine fluctuations for both prestressed rope and bend-stabilized networks, we show that prestressing removes criticality from sub-isostatic rope networks and shifts the critical point to lower values for a bend-stabilized network.

II Model

We use the phantom triangular lattice model Broedersz and MacKintosh (2011); Broedersz et al. (2012); Broedersz and MacKintosh (2014) to study fibrous networks in the rope limit. The networks are generated on a periodic 2D triangular lattice with lattice spacing $l_{0}=1$ and freely-hinging crosslinks at intersection points. The lattice occupies an area $A=aW^{2}$, where $W$ is the network size and $a$ is the area of a unit cell. The full triangular network has local connectivity of 6. In real biopolymer networks, the local connectivity can be either 3, corresponding to a branching point, or 4, corresponding to a crosslink between two fibers, so the maximum average connectivity cannot exceed 4. To satisfy this constraint, we randomly detach one of three fibers at every crosslink on a triangular lattice, reducing the average connectivity from 6 to 4. Moreover, we cut a single bond on each fiber at a random location to remove the unphysical effects of network-spanning fibers. Since there are $3\times W$ fibers on a triangular lattice and the average network connectivity $z$ is calculated as twice number of bonds divided by number of nodes, this fiber cutting step gives a connectivity $z=4-\frac{6}{W}$ which approaches 4 in the limit of large $W$. We then remove random bonds to obtain the desired connectivity. Dangling ends, which have no effect on the mechanical behavior of the network, are removed. Figure 1a shows a small section of a full phantom triangular model.

In order to compare our results with other geometries, we use three additional 2D network structures: (i) fully branched hexagonal (honeycomb) lattice Rens et al. (2016) (ii) Delaunay triangulation network Delaunay (1934); Gurtner and Durand (2014) and (iii) a Voronoi network Dirichlet (1848); Voronoi (1908); Heussinger and Frey (2007). We note that the notion of fibers in these structures is less well-defined than in the context of a triangular lattice. The honeycomb network is easily derived from a full triangular lattice by cutting specific bonds. The Delaunay networks are constructed by placing $N$ random points in a $W\times W$ box and triangulate them in a way that there is no point inside the circumcircle of any triangle (the unique circle passing through the three vertices of the triangle) which maximizes the smallest angle among all triangulations of the given point set Delaunay (1934). We use $N=W^{2}$ to obtain an average bond length close to $1$ similar to the triangular lattice model. An interesting aspect of a full Delaunay structure is that it has, by construction, a non-uniform local connectivity, in contrast to a uniform structure of a full triangular lattice. The average connectivity of a full Delaunay network, however, is 6 similar to a triangular lattice. To achieve a sub-isostatic Delaunay network ($z<4$) , we randomly remove bonds with no phantomization. Voronoi networks are derived from Delaunay networks of $N=W^{2}/2$ random points by connecting the centers of the circumcircles which gives an average bond length close to $1$ similar to the honeycomb lattice structure. Like the honeycomb lattice, Voronoi networks have a uniform local connectivity of 3. To remove edge effects, we impose periodic boundary conditions in both directions for all networks and utilize Lees-Edwards boundary conditions Lees and Edwards (1972) to apply shear strain. We use the network size of $W=100$ for phantom triangular, $W=120$ for honeycomb, and $W=70$ for Delaunay and Voronoi models. In order to obtain sufficient statistics, we use $50$ different realizations for every simulation. Moreover, to remove possible effects of underlying anisotropy, we average the quantities over both the positive and negative strain directions.

The energy of these networks has two main contributions: stretching of individual bonds and bending between nearest-neighbor bonds (collinear bonds in phantom triangular networks). Therefore, the total energy of the network with stretching stiffness $\mu$ and bending stiffness $\kappa$ is written as

[TABLE]

where $l_{ij,0}$ and $l_{ij}$ are the initial and current bond length between crosslinks $i$ and $j$, respectively, and $\theta_{ijk,0}$ and $\theta_{ijk}$ are the initial and current angle between neighboring bonds $ij$ and $jk$, respectively. For networks with finite stiffness, we set $\mu=1$ and vary the dimensionless bending stiffness $\tilde{\kappa}=\frac{\kappa}{\mu l_{0}^{2}}$. In the rope limit, however, the total energy depends only on extension of bonds, i.e., we remove bending and also compressive terms from the above energy expression which gives

[TABLE]

where $\Theta(x)$ is the Heaviside step function. This is indeed an extreme limit of an asymmetric Hookean spring which has spring constant $\mu$ in the extended state and no resistance under compression. Figure 2a shows the force-extension curve for a rope segment with insets describing stretching and compression of the segment with original length $l_{0}$ under extension $\Delta\ell=\ell-\ell_{0}$. After applying a deformation, we minimize the total energy of the network using FIRE algorithm Bitzek et al. (2006). The stress components are calculated using the microscopic definition of stresses in a polymeric system as discussed in Doi and Edwards (1988); Shivers et al. (2018)

[TABLE]

where $A$ is the area of the simulation box, $f_{ij,\alpha}$ is the $\alpha$-component of the force exerted on crosslink $i$ by crosslink $j$, and $r_{ij,\beta}$ is the $\beta$-component of the displacement vector connecting crosslinks $i$ and $j$. In order to investigate the effect of external stresses on a sub-isostatic rope network, we induce finite normal stresses by applying either bulk or uniaxial extension. After applying bulk or uniaxial strain to induce a finite external normal stress, we investigate the shear rheology of the network by applying incremental shear strains to the prestressed network. This procedure is schematically shown in Fig. 2b.

The phase diagram of sub-isostatic networks in the rope limit is shown in Fig. 2c and Fig. 2d by looking at non-affine fluctuations and differential shear modulus respectively. In the case of simple shear strain only, the network is floppy below an applied critical strain that is a function of network connectivity and geometry. Beyond this critical point, stretching modes rigidify the network (Arrow A in Fig. 2c). Although volume-preserving shear deformation of sub-isostatic networks has been extensively studied in prior work, this strain-induced rigidification occurs under any type of applied extensional strain; indeed, we see similar phase transition between unstable and stable states under both isotropic expansion and uniaxial extension (Arrow B in Fig. 2c). This strain-controlled transition occurs due to tension propagation between boundaries that generates a state of self-stress, therefore stabilizes the network Vermeulen et al. (2017). To test this in the rope limit, we deform networks by either isotropic expansion or uniaxial extension until the network develops a finite (normal) prestress $\sigma_{P}$ and apply step-wise shear strains in the direction of arrow A in Fig. 2c. The isotropic case mimics the uniform active stresses generated by motor proteins or cell contractility on a fibrous network substrate. The second case is motivated by axial expansion or compression experiments of biopolymer gels surrounded by buffer, where the solvent freely flows in or out of the sample thereby preserving the gel boundaries van Oosten et al. (2016). Since the uniaxial extension is applied in $y$-direction, throughout this paper, $\sigma_{P}=\sigma_{yy}$ refers to the normal prestress generated by uniaxial exension prior to a step-wise shear deformation, likewise $\sigma_{P}=\sigma_{B}=\frac{1}{2}(\sigma_{xx}+\sigma_{yy})$ in the case of bulk expansion and $\sigma_{\perp}=\sigma_{yy}$ denotes the generated normal stress during shear deformations which is equal to $\sigma_{P}$ at $\gamma=0$.

III Results

Figure 3a shows the differential shear modulus $K=\frac{d\sigma_{\parallel}}{d\gamma}$ where $\sigma_{\parallel}$ is the shear stress, versus shear strain $\gamma$ for a sub-isostatic rope network ($z=3.2$) for different amounts of external normal stress $\sigma_{B}$ caused by bulk expansion. In the absence of external stress ($\sigma_{B}=0$), the sub-isostatic rope network has no resistance under small shear deformation. Applying sufficient bulk expansion to induce finite $\sigma_{B}$ (crossing the phase boundary along the $\epsilon$ axis in Fig. 2c) stabilizes the network, resulting in a finite shear modulus in the linear regime, and further increasing $\sigma_{B}$ leads to an increase in $K$ in the linear regime. Similar behavior is observed in sub-isostatic fiber networks with finite bending stiffness $\kappa$ when $\kappa$ is increased Sharma et al. (2016a); Licup et al. (2015, 2016); Sharma et al. (2016b) in the absence of applied bulk strain (below the phase boundary in Fig. 2c). Despite this similarity between stress-stabilized and bend-stabilized networks, the microscopic picture of these two mechanisms is intrinsically different. Bend-stabilized networks resist deformations in the linear regime due to their bending stiffness $\kappa$, whereas stress-stabilized rope networks have already crossed the phase boundary ($\epsilon\geq\epsilon_{c}$ in Fig. 2c) and thus show resistance under small shear strains because of highly stretched segments. We find that in rope networks under either bulk or uniaxial expansion, the linear shear modulus $G=\lim_{\gamma\to 0}K$ increases with the external normal stress $\sigma_{P}$ as a power law, $G\sim\sigma_{P}^{\alpha}$, with a non-mean-field exponent $\alpha$ (see Fig. 3b). As expected, stress-stabilized networks with higher average connectivity $z$ show a larger linear shear modulus under equivalent $\sigma_{P}$. As shown in the inset of Fig. 3b, we find a very weak dependence of this scaling exponent on the network average connectivity $z$. Moreover, this scaling exponent appears to be independent of the prestressing method we used, i.e., bulk or uniaxial extension. As shown in Fig. 3b, the linear shear modulus of different network connectivity under large prestress deviates substantially from the power law scaling and has a converging trend. This is due to the fact that network segments are massively stretched under large expansion steps and hence the linear shear modulus is primarily governed by this stretching load rather than connectivity or density of networks. The significance of a sublinear scaling suggests the role of prestress in stabilizing a rope network. We define a susceptibility to the applied prestress, $\chi_{P}\propto\frac{dG}{d\sigma_{P}}\sim\sigma_{P}^{\alpha-1}$, which diverges as $\sigma_{P}\rightarrow 0$ for $\alpha<1$. In the small strain regime $\gamma<\gamma_{c}$, we find that the excess normal stress $|\sigma_{\perp}-\sigma_{P}|$ generated under shear remains negligible ($\lesssim 10^{-5}$) with fixed axial strain $\epsilon$. This is in contrast with fiber networks stabilized by finite bending modulus $\kappa$ or super-isostatic networks with only central force interactions, where in the linear regime, the normal stress in the shear deformation is quadratic , i.e., $\sigma_{\perp}\sim\gamma^{2}$ Janmey et al. (2007); Conti and MacKintosh (2009); Licup et al. (2016). Interestingly, normal stresses in the shear deformation of a super-isostatic rope network, which is stable due to the large number of constraints, are proportional to shear strain in the linear regime, i.e., $\sigma_{\perp}\sim\gamma$ (see Fig. 3d). This is due to the fact that the symmetry of the potential is broken in the rope limit. From these observations, we propose the following nonlinear stiffening relation in the stress-stabilized sub-isostatic rope networks Licup et al. (2015)

[TABLE]

in which $\gamma_{c}$ refers to the critical strain of a rope network in absence of any prestress. Since in the linear regime $\sigma_{\perp}\sim\sigma_{P}$, the above relation results in the sublinear scaling relation $G\sim\sigma_{P}^{\alpha}$. The dashed lines in Fig. 3a show predictions of the above stiffening relation.

The scaling relation $G\sim\sigma_{P}^{\alpha}$ with a non-mean-field exponent $\alpha<1.0$ signifies the influence of the non-mean-field characteristics of sub-isostatic phantom triangular networks, i.e their disordered structure and inhomogeneous local connectivities. Networks with more ordered geometries are expected to show mean-field behavior. To investigate this, we simulate full honeycomb, full Voronoi, and diluted Delaunay networks in the rope limit. As expected, networks with uniform local connectivity of 3, i.e., full honeycomb and Voronoi, exhibit mean-field $\alpha=1$ in contrast to the diluted Delaunay network, with inhomogeneous local connectivity, which yields a non-mean-field $\alpha\cong 0.9$ as the diluted phantomized triangular networks (see Fig. 3c).

The transition between floppy (bending-dominated) and rigid (stretching-dominated) states of sub-isostatic fiber networks has been studied in the absence (presence) of bending interactions Shivers et al. (2018); Sharma et al. (2016a); Sheinman et al. (2012b); Licup et al. (2015); Feng et al. (2016); Sharma et al. (2016b); van Oosten et al. (2016); Vahabi et al. (2016); Vermeulen et al. (2017); Jansen et al. (2018); Wyart et al. (2008); Rens et al. (2016, 2018); Merkel et al. (2018); Silverberg et al. (2014). Many of these prior studies have shown that this transition is accompanied by critical signatures, such as the divergence of the non-affine (inhomogeneous) fluctuations in the strain field. Although prior work has discussed possible discontinuity of the modulus at the transition Vermeulen et al. (2017); Merkel et al. (2018), it is important to note that such a discontinuity does not alter the critical nature of this strain-controlled transition. Since it involves a second derivative of the energy with respect to the control variable strain, the modulus can be thought of as analogous to the heat capacity Shivers et al. (2018), which can be discontinuous at a critical point.

Similar to previous studies Sharma et al. (2016a), we define the differential non-affinity as

[TABLE]

where $l$ is the typical bond length of the network and $\delta\mathbf{u}^{\text{NA}}=\mathbf{u}-\mathbf{u}^{\text{affine}}$ is the differential non-affine displacement of a crosslink caused by applying a small amount of shear strain $\delta\gamma$. We find the average of this quantity over all crosslinks in the network. Like in sub-isostatic spring networks, sub-isostatic rope networks in the absence of any external stresses show a mechanical phase transition under simple shear deformation between floppy and rigid states (Arrow A in Fig. 2c) coinciding with a peak in the non-affine fluctuations. Without any applied prestress, $\delta\Gamma$ of a sub-isostatic rope network shows a large peak at the critical strain $\gamma_{c}$ where the stretching energy becomes finite and stabilizes the network (see Fig. 4a). The critical strain has a strong dependence on the network connectivity $z$ as has been discussed previously Sharma et al. (2016a); Licup et al. (2016); Shivers et al. (2018). Applying any finite external stresses by either isotropic or uniaxial expansion, however, removes the criticality signatures and therefore $\delta\Gamma$ shows no peak as shown in Fig. 4a. This is due to the fact that by imposing large expansion prior to shear deformation (Arrow B in Fig. 2c), we indeed move the network out of floppy state and stabilize it. Moreover, applying small extensional strain $\epsilon<\epsilon_{c}$ leads to a decrease in $\gamma_{c}$ with $\gamma_{c}=0$ for $\epsilon\geq\epsilon_{c}$. This effect can be captured more easily in a bend-stabilized network in which finite $\sigma_{P}$ occurs for $\epsilon<\epsilon_{c}$. As shown in Fig. 4b, for bend-stabilized sub-isostatic fiber networks, the applied prestress by bulk expansion $\sigma_{B}$ clearly shifts the critical strain to lower values until a point where $\sigma_{B}$ is large enough, i.e., $\epsilon\geq\epsilon_{c}$ to transform the bending-dominated to stretching-dominated state and removes criticality from the network.

Interestingly, the onset of strain stiffening $\gamma_{0}$ which we define the strain where differential shear modulus is twice as large the linear shear modulus $K\simeq 2G$, increases with increasing external normal stresses for a sub-isostatic network in a rope-like potential (see Fig. 5). In the case of bend-stabilized networks, however, $\gamma_{0}$ decreases by applying extensional strain $\epsilon<\epsilon_{c}$ and shows a similar behavior to rope networks under extensional strains larger than $\epsilon_{c}$. In the experimental studies of biopolymer networks, one could apply external normal stress $\sigma_{P}$ and measure the differential shear modulus $K$ at this stressed state. By measuring shear modulus curves for every applied prestress, onset of strain stiffening $\gamma_{0}$ versus the prestress $\sigma_{P}$ can be obtained and could serve as an indicator that such networks have been strained past the critical strain $\epsilon_{c}$. So dependence of $\gamma_{0}$ on $\sigma_{P}$ can serve as an indicator of whether the network is bending-dominated or stretching-dominated in the reference state. Indeed the experimental studies on the reconstituted collagen networks are shown to develop normal stresses due to polymerization and boundary effects Licup et al. (2015) which can significantly affect their mechanical response. In real tissues, the external stresses exist due to the interactions between different tissues as well as embedded cells.

IV Conclusion

Floppy sub-isostatic fiber networks are stabilized via various mechanisms such as applying large strain, introducing fiber bending interactions, imposing active stresses, and thermal fluctuations. Here we have shown that external normal stresses can rigidify linearly unstable fiber networks. We considered the case in which connected network nodes interact only through rope-like tensile forces. The stability of rope-like structures has been studied previously in the context of disordered networks Delaney et al. (2005); Shen and Wolynes (2006); Wang et al. (2011); Han et al. (2018); Ronceray et al. (2019), although the role of prestress on the shear rheology of such systems was not examined. We found that the linear shear modulus of these stress-stabilized networks scales as a power law with the applied external normal stress. The scaling exponent exhibits a non-mean-field value for connectively disordered networks and a mean field value for connectively homogeneous structures. We also investigated the effect of prestress on the criticality of the networks. In order to stabilize a sub-isostatic rope network, a non-zero prestress which corresponds to an extension strain larger than $\epsilon_{c}$ needs to be applied. This indicates that the network becomes rigid and the criticality corresponding to the transition between two floppy and rigid states is removed. Indeed, the non-affine fluctuations of a stress-stabilized rope network show no peak and are clearly supressed. For a bend-stabilized network, however, small prestress corresponding to a small extension $\epsilon<\epsilon_{c}$ shifts the critical strain $\gamma_{c}$ to lower values which is clearly observed by calculating the non-affine fluctuations. Moreover, we find that the onset of strain stiffening $\gamma_{0}$ monotonically increases by applying prestress to a rope network. For a bend-stabilized network, however, this behavior is non-monotonic with a decreasing trend in the bending-dominated regime and increasing trend similar to the stress-stabilized rope networks in the stretching-dominated regime. The distinctive behavior of $\gamma_{0}$ versus external stress $\sigma_{P}$ can be used to determine whether a fiber network is in the bending-dominated or stretching-dominated regime.

Acknowledgments

This work was supported in part by the National Science Foundation Division of Materials Research (Grant DMR1826623) and the National Science Foundation Center for Theoretical Biological Physics (Grant PHY-1427654).