Rigidity and fracture of fibrous double networks

Composite fiber networks are ubiquitous in biological systems and synthetic materials with tunable and robust mechanical properties. For example, the cytoskeleton, the scaffolding that gives eukaryotic cells mechanical integrity and shape, is a self-organized composite network of protein filaments, including actin and microtubules Pollard and D. (2017). The distinct rigidity of actin and microtubules enables cells to exhibit complex stress responses and architectures essential for a wide range of functions Pelletier et al. (2009); Shea N.Ricketts and M.Robertson-Anderson (2018); Ricketts et al. (2019); Farhadi et al. (2020). As another example, the load-bearing capability of musculoskeletal tissues such as articular cartilage arises from a network-like extracellular matrix made of collagen fibers and proteoglycans Mow et al. (1992); Fung (1993); O’Driscoll (1998); Lai and Levenston (2010). Finally, several synthetic double network hydrogels have recently emerged as extraordinarily robust materials with considerable toughness and fracture resistance compared to conventional single network hydrogels. For instance, the PAMPS-PAAm double network hydrogel, which consists of interacting networks of poly(2-acrylamide-2-methyl-propane sulfonic acid) and polyacrylamide, has a tearing energy $\sim 4400J/m^{2}$, which is several hundred to a thousand times that of single network PAAm and PAMPS hydrogelsHaque et al. (2012); Nonoyama and Gong (2015). The exceptional mechanical response of these double network systems derives from the synergistic interplay between two networks with very different single-filament and collective properties.

The rigidity of stiff networks made of a single type of fiber or biopolymer, henceforth called single networks (SN), have been vigorously investigated in the past two decades, uncovering mechanical phase transitions, distinct mechanical regimes, and novel non-linear mechanical properties Broedersz et al. (2011); Head et al. (2003); Wilhelm and Frey (2003); Gardel et al. (2004); Missel et al. (2010); Heussinger and Frey (2006); Das et al. (2007); Kang et al. (2009); Sheinman et al. (2012); Picu (2011). More recently, studies of the fracture mechanics of such networks have demonstrated that low network connectivity and system-wide distribution of damage can provide protective mechanisms against failure Zhang et al. (2017); Burla et al. (2020). The mechanics and fracture of composite networks, on the other hand, are only beginning to be explored, and theoretical studies to date have focused on composites made of rod-like inclusions in an SN Das et al. (2008); Das and MacKintosh (2010), bidispersity in single networks Burla et al. (2019), and continuum models of double network hydrogels Haque et al. (2012). The mechanical structure-function properties of double networks (DN) are less well understood, and there remain many open questions as to the mechanisms by which DNs achieve such remarkable mechanical performance. In particular, it is unknown how much the second network can affect the rigidity percolation threshold for the combined DN system, an important parameter for setting the stiffness. Nor is it known to what degree the second network can tune the strain necessary for network failure (extensibility), the maximum stress reached (strength), and the toughness under extension, all of which are important for determining the workable range of strains and stresses over which the system maintains its integrity. Addressing these questions will help guide the rational design of biomimetic soft materials with tunable mechanics and provide insights into the rigidity and fracture-resistance of load-bearing tissues such as cartilage Trivedi et al. (2008); Taylor et al. (2012).

Here, we address these questions by combining two structure-function frameworks, (i) a double network (DN) made of two interacting disordered networks with very distinct mechanics and (ii) rigidity percolation theory to construct a Rigidly Percolating Double Network model. Rigidity percolation theory models a biopolymer network as a disordered network of fibers consisting of flexible, sparsely connected regions and stiff, densely connected regions Silverberg et al. (2014); Feng et al. (1985); Das et al. (2007, 2012); Broedersz et al. (2011). When the network consists primarily of sparsely connected regions, it does not offer any resistance to shear deformations and has zero shear modulus. In contrast, when densely connected regions span the network, the network has a finite shear modulus. The system undergoes a mechanical phase transition from non-rigid to rigid at a certain fiber density known as the rigidity percolation threshold.

The rigidly percolating double network model is made of a stiff network interacting with a flexible network (Fig.1(a)). We study the shear response and crack propagation in this DN and show that the interplay of the mechanically distinct networks facilitates tunable mechanics and enhanced fracture resistance of the DN. Each of the two networks in the DN is modeled as a disordered kagome network and is constructed following the procedure described in Silverberg et al. (2014). The bonds in the two networks are randomly removed according to two different probabilities, $1-p_{1}$ for the stiff network, and $1-p_{2}$ for the flexible network, where $0<p_{1},p_{2}<1$ are the bond occupation probabilities. The energy cost of deforming this double network is given by:

where $E_{1}$ is the deformation energy of the stiff network, $E_{2}$ is the deformation energy of the flexible network, and $E_{3}$ is the deformation energy of the bonds connecting the two networks. In $E_{1}$, the first term corresponds to the energy cost of fiber stretching, and the second term to fiber bending Silverberg et al. (2014). In $E_{2}$, we have similar contribution for fiber stretching, but there is no energy cost of fiber bending. The stretching modulus of the fibers in the stiff and flexible networks are $\alpha_{1}$ and $\alpha_{2}$ respectively where $\alpha_{1}>\alpha_{2}$, the bending modulus of the fibers in the stiff network is $\kappa_{1}$. The networks interact via Hookean springs with spring constant $\alpha_{3}$ which connect the midpoints of bonds ($\textbf{x}_{1}$,$\textbf{x}_{2}$), and are only present when corresponding bonds are present in both networks. The indices $i,j,k$ refer to sites (nodes) in each lattice based network, such that $p_{ij}$ is 1 when a bond between those lattice sites is present, 0 if a bond is not present. The quantities $\textbf{r}_{ij}$ and $\textbf{s}_{ij}$ refer to the vector lengths between lattice sites $i$ and $j$ for the deformed stiff and flexible networks respectively, while $\textbf{r}_{ij0}$ and $\textbf{s}_{ij0}$ are the corresponding quantities for the initial undeformed networks. The angles $\Delta\theta_{ijk}$ correspond to the change in the angles between initially collinear bond pairs $ij$ and $jk$ for the deformed and undeformed network respectively. See Figure 1 (b) for illustration of the properties of the bonds in the networks. Unless otherwise noted, we have used the following biologically relevant parameters in the results presented: $\alpha_{2}/\alpha_{1}=0.1$, $\kappa/\alpha_{1}=.004$ Silverberg et al. (2014), and $\alpha_{3}=\alpha_{1}+\alpha_{2}$. ¹¹1The value of $\alpha_{3}$ was chosen to be the effective spring constant of two springs $\alpha_{1}$ and $\alpha_{2}$ in parallel. Simulations with a much small value of $\alpha_{3}/\alpha_{1}=0.01$ yielded qualitatively similar results and the percolation thresholds were unchanged.

For the shear response studies (Fig. 1(c)), we adopt a protocol where external deformations are applied along the top and bottom boundaries and periodic boundary conditions are used for the left and right sides of the network. Our simulations of the single network follow the same process, except the deformation energy consists of only $E_{1}$, since $p_{2}=0$. To obtain the linear mechanical response, we apply a shear strain of $5\%$ at the boundaries, minimize the deformation energy using a conjugate minimization scheme and calculate the shear modulus Silverberg et al. (2014).

We show the variation of the rigidity percolation threshold of a single network (SN) and four double networks (DN) by plotting the shear modulus versus bond occupation probability of the stiff network $p_{1}$ (Figure 2). The shear moduli $G$ are normalized by their respective values $G_{0}$ for fully occupied networks. The four DNs correspond to the four different values of the bond occupation probability $p_{2}=0.2,0.4,0.6,0.8$ of the flexible network. We find that the SN has a percolation threshold $p_{1,c}\sim 0.6$ in agreement with previous results Mao et al. (2013), while the DNs have a lower $p_{1,c}$, which decreases with increasing $p_{2}$, reaching $p_{1,c}\sim 0.35$ at $p_{2}=0.8$. This is a noteworthy result, because on their own a single stiff network based on a Kagome lattice has a percolation threshold $\sim 0.6$ Mao et al. (2013) and a flexible network based on such a lattice has a percolation threshold $\sim 1$, but when they form a double network, the resulting additional constraints due to their interaction lead to a lower, tunable percolation threshold. These constraints also allow the normalized shear rigidity of the DNs to be larger than that of the SN at the same value of $p_{1}$, and can be tuned by varying $p_{2}$. This result illustrates a mechanism for how the onset of rigidity for biological and synthetic double networks can be drastically modulated through very small changes in filament concentration in the secondary network.

For the crack propagation studies (Fig. 1(d)), we create a notch 20 times the bond rest length ($\sim 1/5$ times the system length) at the center of the top boundary following the protocol in Ref.Stok and Oloyede (2003), and study how the size of the notch increases as we apply larger and larger tensile strains along the left and right boundaries of the network.

The strains are applied quasi-statically in small increments of $1\%$ up to $70\%$, and after each application, the total energy is minimized to generate the new equilibrium configuration of the deformed DN. We have made the following assumptions regarding breaking and buckling of fibers: When a bond is stretched above a certain threshold, it will break, and when it is compressed above a certain threshold, it will buckle. Bonds in the stiff network break at $120\%$ of their rest length and buckle at $95\%$ of their rest length. Bonds in the flexible network break at $200\%$ of their rest length, but do not buckle. Broken or buckled bonds will no longer contribute to the deformation energy or rigidity of the network. Fiber breaking is an irreversible process, but buckled fibers in our model can “unbuckle” when the extra compression is removed.

To demonstrate, how the fracture mechanics change with the stiff network’s proximity to its rigidity percolation threshold, we present results for simulations of the DN close to ($p_{1}=0.62$) and away from ($p_{1}=0.80$) the rigidity threshold (Fig. 3). The values of $p_{1}$ were chosen, so that the DN has a finite rigidity, irrespective of $p_{2}$ ²²2We found, for example, that when $p_{1}$ was set to 0.55, the DN had zero shear rigidity and exhibited no stresses when $p_{2}$ was $0$, $0.2$, or $0.4$. We find that both the Young’s modulus (Fig.3(a) and (b)) as well as the stress (Fig. 3(c) and (d)) developed in the network initially increase with strain and reach a maximum as previously floppy regions become stretched and align to resist deformation. Once the fibers in the network start to experience strains larger than their stretching (or buckling) thresholds, however, they break (or buckle), causing the network to soften. Remarkably, we find that when the stiff network is close to its rigidity threshold, the normalized Young’s modulus, the maximum or peak stress, and the strain at failure can be shifted dramatically by the flexible network. This tunability arises because the sparsely populated stiff network allows the DN to undergo non-affine rearrangements Head et al. (2003); Das et al. (2007); Heussinger and Frey (2006), leading to large variations in rigidity.

We quantify these trends by comparing the peak stress $\sigma_{p}$, strain at maximum stress $\gamma(\sigma=\sigma_{p})$, and the network toughness $\tau$ versus $p_{2}$ for both DNs in figure 3. We find that the peak stress increases with $p_{2}$ for both DNs due to the additional constraints introduced by the secondary, flexible network. The strain at maximum stress decreases with $p_{2}$ when the stiff network is close to the percolation threshold and remains nearly constant when the stiff network is far from the percolation threshold. Thus, the additional constraints introduced by the secondary network plays a much greater role in restricting deformation when the stiff network is near the rigidity percolation threshold. Finally, we find that while the network toughness increases for both cases, the increase is greater for the network near the rigidity percolation threshold. This result is somewhat surprising since typically the toughness is proportional to the product of the peak stress and strain at failure which remains nearly constant for the $p_{1}=0.62$ data. Here, however, because the network fails gradually, the decrease in stress is far less abrupt than in typical materials and the network toughness is substantially increased.

These results show that the flexible network can modulate the mechanics of the DN far more effectively when the stiff network is just above its rigidity threshold. Additionally, they show how the DN can be modulated to either be extensible, breaking gradually, as is the case for low $p_{2}$ or be stronger, breaking in a more brittle fashion, as is the case for high $p_{2}$. The low $p_{2}$ limit is particularly important in biological tissues such as articular cartilage when it is undergoing osteoarthritis, where the proximity of the stiff (collagen) network to the rigidity percolation threshold varies as a function of tissue depth and the reinforcing flexible network is increasingly removed as the disease progressesGriffin et al. (2014); Jackson et al. .

To illustrate the above-mentioned trade-off visually, we present stills from simulations of crack propagation in DNs with $p_{1}=0.62$ and 0.80 as a function of the applied tensile strain $\gamma$ and $p_{2}$ (Fig.4). We find that when the stiff network is far above the rigidity threshold ($p_{1}=0.80$,Fig.4(a)), the DN ruptures abruptly at $\gamma\sim 0.2$ for all $p_{2}$ values, though the crack morphology is more uniform at higher $p_{2}$. In contrast, when the stiff network is close to the rigidity threshold ($p_{1}=0.62$, Fig.4(b)), we observed a wider range of responses. For $p_{2}=0$, 0.2, and 0.4 the networks are extensible, initially developing microcracks that are distributed throughout. With increasing strain these microcracks grow and the network decreases its rigidity while maintaining a percolated structure. For $p_{2}=0.6$ and 0.8 the networks are more brittle, rupturing less homogeneously and maintaining their rigidity up until the point of failure.

Importantly, this ability to tune the failure characteristics could have numerous important applications. For example, in biologically relevant scenarios the range of tensile strain that a tissue is exposed to may be limited. In this case, it would be more advantageous for the tissue to be extensible rather than brittle over this range of strains. As another example, it may be possible to construct DNs with varying compositions to guide the trajectory or even stall cracks propagating through the material. In fact, given the striking similarity of the crack morphology for $p_{2}=0.6,\gamma=0.35$ in Fig. 4b to the experimentally observed fracture of articular cartilage tissue (see for example Figure 4 in Stok and Oloyede (2003)) it is possible that cartilaginous tissues may already be employing such mechanisms. Finally, we note that the observed richness of behaviors presented here would be further enhanced by including additional tuning parameters such as structural hierarchy, network polarization, or bond polydispersity in either the stiff or flexible networks. This flexibility in resulting material properties and ease of implementation make double networks a very attractive platform for the fabrication of mechanically active artificial tissue constructs. The results presented here are an important step towards achieving this future.