Reentrant Rigidity Percolation in Structurally Correlated Filamentous Networks

Many biological tissues feature a heterogeneous network of fibers whose tensile and bending rigidity contribute substantially to these tissues' elastic properties. Rigidity percolation has emerged as a important paradigm for relating these filamentous tissues' mechanics to the concentrations of their constituents. Past studies have generally considered tuning of networks by spatially homogeneous variation in concentration, while ignoring structural correlation. We here introduce a model in which dilute fiber networks are built in a correlated manner that produces alternating sparse and dense regions. Our simulations indicate that structural correlation consistently allows tissues to attain rigidity with less material. We further find that the percolation threshold varies non-monotonically with the degree of correlation, such that it decreases with moderate correlation and once more increases for high correlation. We explain the eventual reentrance in the dependence of the rigidity percolation threshold on correlation as the consequence of large, stiff clusters that are too poorly coupled to transmit forces across the network. Our study offers deeper understanding of how spatial heterogeneity may enable tissues to robustly adapt to different mechanical contexts.


INTRODUCTION
Networks and network-like structures are ubiquitous in biological cells and tissues, and provide the basis for their mechanical properties and functions.Biopolymer networks are largely responsible for the mechanical response of the cytoskeleton of cells [1][2][3][4][5][6] and the extracellular matrix of tissues [7][8][9][10][11]; more recently, rigidly percolating connected networks of cells have been shown to account for the viscoelasticity of developing embryos [12].These networks are generally highly disordered and spatially inhomogeneous as a result of how they are assembled and disassembled.For example, cytoskeletal networks are highly dynamic and have a complex and heterogeneous spatial organization allowing for context-dependent cell remodeling and response [1].As a second example, the collagen II scaffold in articular cartilage is densest in the vicinity of chondrocytes, cells which secrete extracellular matrix material to construct and sustain collagen networks [13].In fact, previous work has established spatial heterogeneity as a crucial consideration in developing faithful cartilage replacements and scaffolds for tissue regeneration [14,15].
In the last two decades, there has been extensive study of disordered biopolymer networks through invitro experiments and simulations, which have provided a wealth of information about these networks' responses to mechanical stimuli [16][17][18][19][20][21][22][23][24].To date, however, almost all computational studies of biopolymer networks have focused on spatially homogeneous systems and ignored the presence of structural correlations, which can have significant consequences for the collective properties of the network.Here, we address this gap and present a novel investigation of the percolation of rigidity in structurally correlated fiber networks, which are found in many cells and tissues, using a lattice-based framework.
Lattice-based fiber networks are a prominent paradigm for modeling biopolymer scaffolds [19,21,22,24,25].These networks are constructed by laying down infinitely long fibers in a regular pattern, such that whenever two fibers cross, there is a crosslink which allows free rotation of the fibers but does not allow them to slide.The fibers can stretch and bend, but pay energy penalties for these deformations.Each fiber can further be thought of as a collinear series of connected bonds, such that random removal of bonds yields a broad distribution of fiber lengths.
The mechanical response of such a disordered fiber network can be mapped to the fraction of bonds present.Starting with a network in which all the bonds are present, one can progressively decrease network rigidity by removing bonds.Once the network reaches a certain threshold of bond occupation, its elastic moduli undergo a dramatic, many-decade decrease, dropping to negligibly small values.This mechanical phase transition is known as rigidity percolation [19,25].Dilute fiber networks have shown great promise as micromechanical models for invitro cytoskeletal networks [16][17][18][19], and more recently for extracellular matrix networks in tissues [8,11,26].
Previous computational studies of fiber networks have examined spatially homogeneous disordered networks, in which bonds are excluded or retained purely at random, and have not considered the possibility of correlations in the inclusion of bonds [16,17,[19][20][21].In this paper we introduce a model in which bonds are added in a structurally correlated manner to account for the pronounced heterogeneity in the distribution of material observed in cells and tissues.In this model, the likelihood of adding a bond is contingent upon the number of adjacent bonds already present.This protocol gives rise to networks in which already dense regions become further enriched with material, while sparse regions remain comparatively dilute.

II. MODEL A. Network Construction
As we are interested in biopolymer networks in which vertices correspond to crossings of adjacent filaments, we choose as our model network the kagome lattice [8], with a maximum coordination number of 4. We adjust the elastic moduli of networks by randomly retaining a subset of the bonds, such that some portion, p, is included.In the absence of structural correlation, p corresponds the probability that each bond is retained, such that there is an independently and identically distributed probability of keeping each bond.We instead use the term bond portion to reflect the fact that, for correlated dilute networks, once the network is seeded with an initial set of bonds, other bonds do not have an identical likelihood of inclusion.
We employ an iterative process, introduced in [27], in which, at each step, a candidate bond is chosen at random from those bonds that have not yet been included in the network.We then count the number of bonds adjacent to the candidate bond, n a that have already been retained, where adjacency is defined by the condition that two bonds share a common vertex.For the kagome lattice, the maximum possible number of adjacent bonds, n a,max , is 6.Given a correlation strength, c, where 0 ≤ c < 1, the candidate bond is added with a probability where c = 0 corresponds a purely random dilution, and c ≈ 1 corresponds to maximally correlated dilution, in which a bond is rejected unless all adjacent bonds have been retained.This process is repeated until the desired portion of bonds has been retained.In figure 1, we show representative samples of networks with varying degrees of dilution and correlation, including no correlation.

B. Mechanical Model
The bonds of the network resist stretch and compression with a Hookean spring stiffness, α, and resist bending with a bending rigidity, κ.Bending resistance is implemented by regarding adjacent, collinear bonds as consecutive segments of a fiber.If the bond connecting vertices i and j is collinear with and adjacent to the bond joining vertices j and k, then we penalize a change in the angle ∠ijk by an amount proportional to the square of the angular deflection.
We focus on the mechanical response in the linear response regime, so that the deformation energy consists of terms quadratic in the strain.Following [8], we truncate the deformation energy to leading order in the displacement of vertices from the reference configuration of the network, and model the energy as Here, < ij > denotes a sum over pairs of vertices sharing a bond, < ijk > denotes a sum over vertices of adjacent, colinear bonds, and p ij is defined to be 1 if the bond between bonds i and j is retained, and 0 otherwise.Further, u ij denotes the difference between the displacement vectors for vertices i and j, and rij denotes the direction vector of the bond between vertices i and j in the reference state.

C. Structural Relaxation Procedure
We simulate the shear mechanics of our model networks by imposing a series of small, simple shear displacements of the vertices at the top.
The displacements of the vertices at the bottom of the network are constrained to be zero, while along the sides we impose periodic boundary conditions.Given these constraints, we minimize the energy given in (2) for several different shear strains, s , and compute the shear modulus at a target shear strain, * s by numerically estimating the second derivative of the strain energy with respect to shear strain: Full details of this calculation may be found in Appendix A.

A. Shear Mechanics
We first examined how the rigidity percolation threshold can be tuned by varying the degree of c p . 4 .6.80 .6 .8

Candidate Neighbors
FIG. 1. Sample networks are shown in which correlation strength is varied from 0 to .8, and the bond portion is varied from .4 to .8.As the correlation strength is increased, networks transition from arrangements of homogeneously dispersed bonds to networks with dense clusters interspersed with sparse voids.While individual clusters may be rigid, they may be tenuously connected to neigboring clusters, such that, beyond some critical correlation strength, networks will become less effective in transmitting force over macroscopic distances.In the magnified inset at the left of the figure, we illustrate the correlated construction process.A candidate,marked in magenta, is considered for inclusion, and has five neighbors, marked in orange.
correlation.We considered nine distinct values of c, ranging from 0 to 0.8, in steps of .1, and 61 distinct values of p, ranging from 0.4 to 1, in steps of 0.05.We carried out structural relaxation for 10 realizations for each combination of c and p, and identified G for each combination as the geometric mean of the 10 values.The shear modulus G, normalized by its universal maximum, G 0 , is shown vs. p for several values of c in Figure 2a, accompanied by a full phase diagram in Fig. 2b.Figures 2a and 2b provide the first indications of an intriguing variation in the rigidity percolation threshold with the degree of correlation c of the disordered network.For each value of c, we find a qualitatively similar scaling of the shear modulus with the bond fraction.Interestingly, however, the rigidity percolation threshold i.e., the critical bond fraction, p c , at which the shear modulus first differs appreciably from zero, shifts markedly and non-monotonically as c is varied: While introducing a moderate correlation strength initially diminishes p c , this effect saturates at about c = .6,and p c increases for still larger values of c.
To quantitatively identify the rigidity percolation threshold for each value of c, we considered pairs of bond fraction and shear modulus for which the shear modulus ranged from 10 −9 to ∼ 10 −12 .This ensured that the shear modulus was greater than machine or algorithmic error, but still small in comparison with its maximal value, G 0 , at p = 1.We used the method of least squares b. a.
FIG. 2. In panel a, we show the scaling of shear modulus with bond portion for several structural correlation strengths.While the dependence of the shear modulus on the bond portion is qualitatively similar in each case, the point of rigidity percolation shifts initially to the left, then back to the right with increasing structural correlation.In panel b, we show the dependence of G on c and p for the full range of parameter space considered.The reentrance of the dependence of pc on structural correlation strength is clearly discernible in a contour of marginal stiffness on the left side of the heat map.We attribute this reentrance to two competing effects: the need for rigid clusters, and the need for strong coupling between adjacent clusters.
to fit each set of bond fraction-shear modulus pairs to a power law of the form For each value of c, we found a good fit to equation 4 over over at least seven decades of dynamic range in the shear modulus.In each case, the correlation coefficient, R 2 , between log 10 (G) and log 10 [k(p − p c ) β ] is ≥ .96.Further details are provided in Appendix B.
As shown in Fig. 3a, in which p c is plotted vs. c, our power law fits affirm the trend in p c previously identified by inspection in Fig 2 .Surprisingly, the scaling exponent, β, on the other hand, exhibits the opposite trend, increasing with c until about c = .6,and decreasing thereafter.This unexpected scaling of β can be explained as an earlier onset of rigidity percolation being associated with a more abrupt rise in the shear modulus at the point of percolation.
As shown in figure 3c, we find β to decrease linearly with p c (R 2 = .92).We attribute the decrease in β with p c to two competing factors determining percolation: the presence of large, rigid clusters, and sound mechanical coupling of adjacent clusters.As the correlation strength is increased beyond its optimum value, the network segregates into large, dense regions that are too poorly connected to enable optimal transmission of stress.The optimal correlation strength of .6 strikes an ideal balance, enabling the greatest gain in stiffness per unit material.To substantiate this conjecture, we turn to detailed analysis of network displacement fields.

B. Analysis of Non-affine Deformation
To gain further insight into the micromechanical mechanisms underlying the reentrance in the rigidity percolation threshold, we quantified the degree of non-affinity in each network.Non-affinity quantifies the departure of a displacement field of a strained material from the displacement expected in a simple, homogeneous elastic continuum.In an affinely deforming network subjected to simple shear, a vertex with an initial location r 0 = (x 0 , y 0 ) will be mapped to the final location r A = (x 0 + s y 0 , y 0 ), where s is the shear strain.The non-affine displacement field, u N A is defined as r − r N A , where r is the true displacement field.The non-affine parameter, Γ, is then defined as where l 0 is the length of an undeformed bond, and N is the number of vertices in the network [28].Non-affinity is a well established means of characterizing the difference in behavior of a purely entropic rubber and a gel of semiflexible polymers [24].We computed Γ for all ten network realizations for each combination of p and c, and calculated the final non-affine parameter as the arithmetic mean over all realizations.As reported in previous studies, we find a pronounced peak in the non-affine parameter near the rigidity percolation threshold, as shown in Figure 4a; however, for increasing structural correlation, this peak becomes progressively lower and broader.We posit that this peak broadening is associated with the formation of local rigid regions, joined by weakly connected interstitial regions, such that stress is distributed in a non-uniform manner over a larger range of bond fractions for networks with high structural correlation.We investigate this idea further by considering the dispersion in the distribution a.

FIG. 3. Scaling behavior near the rigidity threshold:
In a, we show the results of fitting shear modulus -bond portion pairs to Eq. 4. In case, we find a sound fit spanning 6 to 7 decades in the shear modulus.In b, we show fitting results for the critical bond portion, pc, and the scaling exponent, β, of the shear modulus with the excess bond portion.The scaling exponent exhibits the opposite trend from the critical bond portion, indicating that a low percolation threshold is accompanied by a more abrupt increase in the shear modulus.In c, we demonstrate that the critical bond portion is a reliable predictor of the critical exponent relating growth in the shear modulus to the excess bond portion.We find that the value c = .6yields an optimal trade-off between a need for large, rigid clusters, and a need for sound mechanical coupling of adjacent clusters to enable coordinated, system-spanning force propagation.For excessive correlation, dense clusters amount to stiff inclusions in an otherwise under-coordinated network.
of bond strains, and the spatial correlations in the nonaffine displacement field.We first consider a radial non-affine correlation function, g(r).We take the inner product of non-affine displacements for all pairs of points within some cutoff distance, r cut , of one another, and bin displacement vectors between pairs of points into annular sectors of thickness ∆r.We then define g(r) as where the first average runs over all points, and the second average runs over all distinct pairs of points i, j such that the positions of vertices i and j in the undeformed lattice are separated by a distance in the range [r, r + ∆r).This normalizes g(r) to be equal to 1 when r = 0.
For each combination of p and c, we find g(r) to be well fit (R 2 ≥ .98)by the form We show a representative set of curves in Fig. 4c, for c = .6and varying bond fractions.Initially, the floor of g(r) decreases, reaching its lowest point near the rigidity percolation threshold, then steadily increases for larger values of p.We attribute the early decrease in the floor to incipient rigid clusters that deform differently from the surrounding soft regions, such that there is no coordinated, long-range force transmission.The subsequent rise in the minimum of g(r) is associated with the emergence of system-wide spanning force chains beyond the rigidity percolation threshold.
From the decay distance, λ, in equation ( 7), we infer an effective mechanical length scale.Different correlation strengths yield qualitatively similar dependence of λ upon bond portion, but λ at a fixed bond portion steadily increases with increasing correlation.This further affirms the idea that the size of a coherently deforming cluster becomes progressively larger with growing correlation.Results are shown in Figure 4, with trend lines computed from cubic basis splines.
We finally seek to account for the reentrant scaling of p c with structural correlation strength by identifying the critical mechanical length scale, λ c , at the onset of rigidity percolation for each structural correlation strength, c.We estimate λ c by interpolation using the previously mentioned cubic basis splines.We find that p c varies non-monotonically with λ c , with an initial decrease until λ c exceeds about 5 bond lengths, after which p c once more increases.The optimal value of λ c is that obtained for a structural correlation strength of .6, in concert with our previous findings.We further observe that λ c appears to diverge according to a power law as c approaches 1.In the limit c = 1, either all bonds can be present or no bonds can be present, so that the only percolating network would be a fully connected network, in which vertex displacements are correlated over arbitrarily large distances.We thus find that, while small, rigid islands must nucleate to enable the most efficient percolation, excessively large rigid clusters leave too little material elsewhere to enable the formation of system-wide force chains.We show the non-affine parameter, as defined in Eq. 5, vs bond portion for a number of correlation strengths.In each case, the non-affine parameter exhibits a pronounced peak at the rigidity percolation transition.Larger correlation strengths correspond to non-affine re-arrangements occurring over a wider range of bond-fractions and are reflected in the peaks becoming progressively broader with increasing c. b We display azimuthally averaged non-affine correlations gNA(r) as a function of distance r between vertices.In each case, we find the decay to be exponential, allowing us to extract an emergent mechanical length scale, λ.This length scale monotonically increases with bond portion p, suggesting longer range correlations for networks with more fiber content.The non-affine correlations in the large r limit decrease with increasing bond portion below the rigidity percolation threshold; however, above this threshold they increase with increasing bond portion.While we choose a correlation strength of .6 here, networks with differing correlation strength exhibited qualitatively similar behavior.In c, we show emergent mechanical length scales λ for different combinations of p and c.At low bond portions, the length scale λ steadily increases with correlation strength for a given p, whereas at high bond portions, values of λ for different values of c converge and approach the system size.In d, we identify a critical mechanical length scale, λc, as the decay length of non-affine correlations at the onset of rigidity percolation, which exhibits a power law divergence as the structural correlation strength nears 1 (see inset).We show that the critical bond portion pc varies non-monotonically with λc in a manner reminiscent of the scaling of pc with with structural correlation strength c in Fig. 3(b).

IV. CONCLUSION
We have introduced and investigated a model of rigidity percolation in spatially correlated networks.Our study of the scaling of the shear modulus near percolation, coupled with our analysis of networks' strain fields, offers a straightforward physical picture accounting for the reentrant scaling of p c with c.While the length scale over which a network's displacement field is well coordinated grows monotonically with correlation strength, eventually neighboring rigid clusters become poorly coupled.Weak tethers between dense islands of bonds lead to strain being highly concentrated, rather than the load being distributed evenly throughout the network.
This work broadens the already successful rigidity percolation framework to better account for the mechanical response of structurally correlated, heterogeneous networks found in cells and tissues.We anticipate this work will usher in further studies exploring the role of anisotropy [29] observed in many extracelluler matrices.Our findings indicate that, rather using just an averaged, system-wide characterization of network topology, local spatial patterns should be considered to fully understand tissues' responses to We neglect values of the shear modulus less than 10 −9 , in units of the bond stretching modulus, to avoid numerical noise.We also neglect the final 20 points in each series, corresponding to bond portions greater than .8,as the shear modulus depends upon the bond portion in a qualitatively different manner as it reaches its upper plateau.For each correlation strength, we obtain a sound fit of the series of (G, p) pairs to Eq. B3 by the method of least squares.In each case, the correlation coefficient R 2 for log 10 (G) vs. log 10 [k(p − p c ) α ] is ≥ .96.

Appendix C: Quantifying Anisotropy in Non-affine Correlations
We further look for evidence of orientational order in non-affine correlations by averaging not over an annular sector, but rather over all pairs of points whose relative displacement has magnitude r, and makes an angle ϕ with the positive x axis.We define a measure of correlation ψ(r, ϕ): , (C1) where summations are over all vertices, δ is the Kronecker delta, and θ 1,2 is the angle between the displacement r 2 − r 2 and the positive x axis.Results are shown in Fig. 5.The color scale for each combination of correlation strength and bond portion is mapped to the range spanning the minimum and maximum values of ψ for that combination.
While for low bond portions, the correlation between non-affine displacement is highly isotropic, correlations decay far more gradually along the direction of applied strain at large bond portion.Well above the rigidity percolation threshold, networks deform in a nearly affine manner, as shown in Fig. 4 a.In this regime, the discrete rotational symmetry of the Kagome network introduces elastic anisotropy.This anisotropy is more pronounced for networks with less structural correlation, as highly correlated networks have greater fluctuation in local stiffness, a trait known to increase non-affinity [28].The decay in non-affine displacement coefficient becomes more gradual with growing bond portion for all correlation strength, but growth in decay length becomes markedly more rapid for highly correlated networks.We also note that, with growing bond portion, non-affine correlations exhibit increasing anisotropy, with decay in correlations becoming much more gradual along the direction of applied shear.

6 FIG. 4 .
FIG.4.a.We show the non-affine parameter, as defined in Eq. 5, vs bond portion for a number of correlation strengths.In each case, the non-affine parameter exhibits a pronounced peak at the rigidity percolation transition.Larger correlation strengths correspond to non-affine re-arrangements occurring over a wider range of bond-fractions and are reflected in the peaks becoming progressively broader with increasing c. b We display azimuthally averaged non-affine correlations gNA(r) as a function of distance r between vertices.In each case, we find the decay to be exponential, allowing us to extract an emergent mechanical length scale, λ.This length scale monotonically increases with bond portion p, suggesting longer range correlations for networks with more fiber content.The non-affine correlations in the large r limit decrease with increasing bond portion below the rigidity percolation threshold; however, above this threshold they increase with increasing bond portion.While we choose a correlation strength of .6 here, networks with differing correlation strength exhibited qualitatively similar behavior.In c, we show emergent mechanical length scales λ for different combinations of p and c.At low bond portions, the length scale λ steadily increases with correlation strength for a given p, whereas at high bond portions, values of λ for different values of c converge and approach the system size.In d, we identify a critical mechanical length scale, λc, as the decay length of non-affine correlations at the onset of rigidity percolation, which exhibits a power law divergence as the structural correlation strength nears 1 (see inset).We show that the critical bond portion pc varies non-monotonically with λc in a manner reminiscent of the scaling of pc with with structural correlation strength c in Fig.3(b).

8 FIG. 5 .
FIG.5.We show the correlation in non-affine parameter, as a function of the magnitude and orientation of separation.The color scale for each panel is normalized to the maximum value for a given combination of bond portion and correlation strength.The decay in non-affine displacement coefficient becomes more gradual with growing bond portion for all correlation strength, but growth in decay length becomes markedly more rapid for highly correlated networks.We also note that, with growing bond portion, non-affine correlations exhibit increasing anisotropy, with decay in correlations becoming much more gradual along the direction of applied shear.