## New Research In

### Physical Sciences

### Social Sciences

#### Featured Portals

#### Articles by Topic

### Biological Sciences

#### Featured Portals

#### Articles by Topic

- Agricultural Sciences
- Anthropology
- Applied Biological Sciences
- Biochemistry
- Biophysics and Computational Biology
- Cell Biology
- Developmental Biology
- Ecology
- Environmental Sciences
- Evolution
- Genetics
- Immunology and Inflammation
- Medical Sciences
- Microbiology
- Neuroscience
- Pharmacology
- Physiology
- Plant Biology
- Population Biology
- Psychological and Cognitive Sciences
- Sustainability Science
- Systems Biology

# A minimal-length approach unifies rigidity in underconstrained materials

Edited by Tom C. Lubensky, University of Pennsylvania, Philadelphia, PA, and approved February 20, 2019 (received for review September 6, 2018)

## Significance

What do a guitar string and a balloon have in common? They are both floppy unless rigidified by geometric incompatibility. The same kind of rigidity transition in underconstrained materials has more recently been discussed in the context of disordered biopolymer networks and models for biological tissues. Here, we propose a general approach to quantitatively describe such transitions. Based on a minimal length function, which scales linearly with intrinsic fluctuations in the system and quadratically with shear strain, we make concrete predictions about the elastic response of these materials, which we verify numerically and which are consistent with previous experiments. Finally, our approach may help develop methods that connect macroscopic elastic properties of disordered materials to their microscopic structure.

## Abstract

We present an approach to understand geometric-incompatibility–induced rigidity in underconstrained materials, including subisostatic 2D spring networks and 2D and 3D vertex models for dense biological tissues. We show that in all these models a geometric criterion, represented by a minimal length

Amaterial’s rigidity is intimately related to its geometry. In materials that crystallize, rigidity occurs when the constituent parts organize on a lattice. In contrast, granular systems can rigidify while remaining disordered, and arguments developed by Maxwell (1) accurately predict that the material rigidifies at an isostatic point where the number of constraints on particle motion equals the number of degrees of freedom.

Further work by Calladine (2) highlighted the important role of states of self-stress, demonstrating that an index theorem relates rigidity to the total number of constraints, degrees of freedom, and self-stresses. Recent work has extended these ideas in both ordered and disordered systems to design materials with geometries that permit topologically protected floppy modes (3⇓–5).

A third way to create rigidity is through geometric incompatibility, which we illustrate by a guitar string. Before it is tightened, the floppy string is underconstrained, with fewer constraints than degrees of freedom, and there are many ways to deform the string at no energetic cost. As the distance between the two ends is increased above the rest length of the string, this geometric incompatibility together with the accompanying creation of a self-stress rigidifies the system (3, 6). Any deformation will be associated with an energetic cost, leading to finite vibrational frequencies. This same mechanism has been proposed to be important for the elasticity of rubbers and gels (6) as well as biological cells (7).

In particular, it has been shown to rigidify underconstrained, disordered fiber networks under applied strain, with applications in biopolymer networks (8⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓–22). Just as with the guitar string, rigidity arises when the size and shape of the box introduce external constraints that are incompatible with the local segments of the network attaining their desired rest lengths. For example, when applying external shear, fiber networks strongly rigidify at some critical shear strain

Rigidity transitions have also been identified in dense biological tissues (29⇓⇓⇓–33). In particular, vertex or Voronoi models that describe tissues as a tessellation of space into polygons or polyhedra exhibit rigidity transitions (34⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓–49), which share similarities with both particle-based models, where the transition is driven by changes to connectivity (48), and fiber (or spring) networks, which can be rigidified by strain. Therefore, an open question is how both connectivity and strain can interact to rigidify materials (22).

Very recently, some of us showed that the 3D Voronoi model exhibits a rigidity transition driven by geometric incompatibility (46), similar to fiber networks. This has also been demonstrated for the 2D vertex model, using a continuum elasticity approach based on a local reference metric (42). For the case of the 3D Voronoi model, we found that there was a special relationship between properties of the network geometry and the location of the rigidity transition, largely independent of the realization of the disorder (46).

Here, we show that such a relationship between rigidity and geometric structure is generic to a broad class of underconstrained materials, including spring networks and vertex/Voronoi models in different dimensions (Table 1 and Fig. 1). We first demonstrate that all these models display the same generic behavior in response to isotropic dilation. Understanding key geometric structural properties of these systems allows us to predict the precise values of a discontinuity in the bulk modulus at the transition point. We then extend our approach to include shear deformations, which allows us to analytically predict a discontinuity in the shear modulus at the onset of rigidity. Moreover, we can make precise quantitative predictions of the values of critical shear strain

We also compare our predictions to previously published experimental data and highlight additional predictions, including a prefactor of 3 that we expect to find generically in a scaling collapse of the shear modulus, shear stress, and critical strain.

We achieve these results by connecting macroscopic mechanical network properties to underlying geometric properties. In the case of the guitar string, the string first becomes taut when the distance between the two ends attains a critical value

Here, we formulate a geometric compatibility criterion in terms of the constrained minimization of the average spring length

Just as with the guitar string, the description of the geometry given by

## Models

Here we focus on four classes of models, which include 2D subisostatic random spring networks without bending rigidity (9, 50⇓⇓⇓–54) and three models for biological tissues: the 2D vertex model (34, 37), the 2D Voronoi model (38, 44), and the 3D Voronoi model (46) (Table 1).

Two-dimensional spring networks consist of nodes that are connected by in total N springs, where the average number of springs connected to a node is the coordination number z. We create networks with a defined value for z by translating jammed configurations of bidisperse disks into spring networks and then randomly pruning springs until the desired coordination number z is reached (9, 27). We use harmonic springs, such that the total mechanical energy of the system is**1** in terms of a mean spring rest length *SI Appendix*, section IA). In simple constraint-counting arguments, each spring is treated as one constraint, and here we are interested in subisostatic (i.e., underconstrained, also called hypostatic) networks with

The tissue models describe biological tissues as polygonal (2D) or polyhedral (3D) tilings of space. For the Voronoi models, these tilings are Voronoi tessellations and the degrees of freedom are the Voronoi centers of the cells. In contrast, in the 2D vertex model, the degrees of freedom are the positions of the vertices (i.e., the polygon corners). Forces between the cells are described by an effective energy functional. For the 2D models, the (dimensionless) energy functional is

All four of these models are underconstrained based on simple constraint counting, as is apparent from the respective numbers of degrees of freedom and constraints listed in Table 1. We stress that Calladine’s constraint-counting derivation (2, 3) also applies to many-particle, non–central-force interactions.

Throughout this article, we often discuss all four models at once. Thus, when generally talking about “elements,” we refer to springs in the spring networks and cells in the tissue models. Similarly, when talking about “lengths ℓ” (of dimension d), we refer to spring lengths ℓ in the spring networks, cell perimeters p in the 2D tissue models, and cell surface areas s in the 3D tissue model (Table 1). Finally, when talking about “areas a” (of dimension D), we refer to cell areas a in the 2D tissue models as well as cell volumes v in the 3D tissue model.

Here we study the behavior of local energy minima of all four models under periodic boundary conditions with fixed dimensionless system size N; i.e., the model is nondimensionalized such that the average area per element is one (41, 44, 46). Under these conditions, a rigidity transition exists in all models even without area rigidity. In particular, for the 2D vertex and 3D Voronoi models, we discuss the special case

## Results

### Rigidity Is Created by Geometric Incompatibility Corresponding to a Minimal-Length Criterion.

We start by comparing the rigidity transitions in the four different models using Fig. 1, where we plot both the differential bulk modulus B and the differential shear modulus G vs. the preferred length

In all models, we find a rigid regime (*B*, *Inset*), as previously similarly discussed in ref. 10. Something similar has also been reported for a 2D vertex model (48).

For the cellular models, we find that the transition point for the case without area rigidity, *D* and *F* and Table 1). Moreover, our 2D vertex model transition point for *SI Appendix*, section IVC), and the location of the transition in vertex models depends somewhat on the energy minimization protocol (44), a feature that is shared with other models for disordered materials (55). Also, in Fig. 1 *D* and *F* the averaged shear modulus always becomes zero at a higher value than the respective average transition point listed in Table 1. This is due to the distribution of transition points having a finite width (see also finite width of *C* and *E*).

We find that in all these models, the mechanism creating the transition is the same: Rigidity is created by geometric incompatibility, which is indicated by the existence of prestresses. We have already shown this for the 3D Voronoi model (46) and the 2D Voronoi model with *SI Appendix*, section IIA).

We find something similar for the disordered 2D vertex model for *SI Appendix*, section IIA), to simplify our discussion here, we consider only configurations without such typically localized prestresses.

We observe that in all of these models, a geometric criterion, which we describe in terms of a minimal average length **2** into (*SI Appendix*, section IA)

According to Eq. **5**, energy minimization corresponds to a simultaneous minimization with respect to *SI Appendix*, section IIA). In contrast, in the rigid regime,

For the cellular models with **5** (46). Again, this is a geometric criterion, which also explains why the transition point *D* and *F*). Moreover, we can understand why the transition point is smaller for

### The Minimal Length Scales Linearly with Fluctuations.

We next study the scaling of the minimal length in the rigid vicinity of the transition. In the rigid regime, the system must compromise between minimizing *SI Appendix*, sections IC–IE).

In *SI Appendix*, section IB, we show analytically that in the absence of prestresses in the floppy regime, the minimal length

To check this prediction, we numerically simulate these models and observe indeed a linear scaling of the

For 2D spring networks, *A*, *Inset*). This scaling behavior of *SI Appendix*, section IF).

For cellular models where area plays a role, Eq. **6** is extended (Fig. 2 *B* and *C*):**6** and **7** are linear expansions of the function

### Prediction of the Bulk Modulus Discontinuity.

Knowing the behavior of the minimal-length function *SI Appendix*, section IC),*A*–*C* and *SI Appendix*, section IE):*SI Appendix*, section IE). A comparison of the predicted

### Nonlinear Elastic Behavior Under Shear.

As shown before (8⇓–10, 12, 14⇓–16, 18⇓⇓–21), underconstrained systems can also be rigidified by applying finite shear strain. We now incorporate shear strain γ into our formalism and test our predictions on the 2D spring networks. However, we expect our findings to equally apply to the cell-based models (*SI Appendix*, sections IC and ID). We also numerically verified that our analytical predictions also apply to 2D fiber networks without bending rigidity (*SI Appendix*, section IIC) (57).

To extend our approach, we take into account that the minimal-length function *SI Appendix*, sections ID and IV). While at the moment we have no formal proof that **10** comes from a numerical check (next paragraph), we hypothesize that for most systems

For a fixed value of γ, the interface between the floppy and the rigid regime is again given by *A*. Indeed, we also numerically find a quadratic scaling for the transition line, *C* and *SI Appendix*, section IIB). We find that for spring networks the coefficient b depends on *C*, *Inset*), which can be understood from properties of the density of states (*SI Appendix*, section IF). To optimize precision, values of b have been extracted from the relation *F*).

Knowing the functional form of *SI Appendix*, section IC):*B*).

When shearing the system starting in the floppy regime (i.e., for **12** predicts a discontinuous change in the shear modulus of *D*, *Inset* and the value of the scaling coefficient in *SI Appendix*, section IIB. Moreover, Eq. **12** also correctly predicts the behavior beyond *D*.

Eq. **12** also correctly predicts the shear modulus behavior for *E*), while for *F*; see *SI Appendix*, section ID for the cellular models), as reported before for many of the cellular models (37, 46, 56). In both cases, we verified that the respective coefficients coincide with their expected values based on the values of

In particular for

We also obtain explicit expressions for both shear stress *SI Appendix*, sections ID and IE). For the latter, we find a negative Poynting effect with coefficient *Inset*), where we use that close to the onset of rigidity,

## Discussion

In this article, we propose a unifying perspective on underconstrained materials that are stiffened by geometric incompatibility. This is relevant for a broad class of materials (6) and has more recently been discussed in the context of biopolymer gels (8, 12⇓–14, 21) and biological tissues (31, 37, 42, 46). Just as with a guitar string, we are able to predict many features of the mechanical response of these systems by quantifying geometric incompatibility—we develop a generic geometric rule

Our work is relevant for experimentalists and may explain the reproducibility of a number of generic mechanical features found in particular for biopolymer networks (12, 17, 21, 25). While we neglect here a fiber-bending rigidity that is included in many biopolymer network models (12⇓⇓–15, 21), future work that includes such a term will further refine our theoretical results (discussion below) and the following comparison with experiments. For shear deformations with

We can also apply these predictions to typical rheometer geometries (*SI Appendix*, section IG). We predict that an atypical tensile normal stress **13** and *SI Appendix*, section IG). This is precisely what has been found for many biopolymer gels like collagen, fibrin, or matrigel (12, 21, 25, 26). However, in contrast to ref. 21, our work suggests that the scaling factor between

Our work also highlights the importance of isotropic deformations when studying prestress-induced rigidity, as demonstrated experimentally in ref. 17. While previous work (8, 9, 12, 14, 15, 18, 20, 21) focused almost (10) entirely on shear deformations, we additionally study the effect of isotropic deformations represented by the control parameter *SI Appendix*, section IG). Second, we also correctly predict that the critical shear strain *A*). While we also predict an increase of the shear modulus G under extension, which was observed as well (17), additional effects arising from the superposition of pure shear and simple shear very likely play an important role in this case. While we consider this outside the scope of this article, it will be straightforward to extend our work by this aspect.

In summary, we have developed an approach to understand how many underconstrained disordered materials rigidify in a manner similar to that of a guitar string. While it is clear that the 1D string becomes rigid precisely when it is stretched past its rest length, we show that in 2D and 3D models, rigidity is governed by a geometrical minimal-length function

In addition, these predictions help unify or clarify several scaling collapses that have been identified previously in the literature. For 2D spring networks derived from jammed packings, we studied the dependence of our geometric coefficients on the coordination number z and find that approximately *SI Appendix*, Fig. S5*A*, *Inset*), we obtain that the critical shear strain

Moreover, we analytically predict and numerically confirm the existence and precise value of a shear modulus discontinuity

One major obstacle in determining elastic properties of disordered materials is the appearance of nonaffinities, which can lead to a breakdown of approaches like effective medium theory close to the transition (10). In our case, effects by nonaffinities are by construction fully included in the geometric coefficients

There are a number of possible future extensions of this work. First, we have focused here on transitions created by a minimal length, where the system is floppy for large

Second, while we studied here the vicinity of one local minimum of

Third, it will be important to study what determines the exact values of the geometric coefficients

Fourth, because we separated geometry from energetics, it is in principle possible to generalize our work to other interaction potentials, e.g., the correct expression for semiflexible filaments (58, 60), and to include the effect of active stresses (54, 65⇓–67). Note that our work directly generalizes to any analytic interaction potential with a local minimum at a finite length. Although in this more general case Eq. **5** would include higher-order cumulants of

Fifth, this work may also provide foundations to systematically connect macroscopic mechanical material properties to the underlying local geometric structure. For example for biopolymer networks, properties of the local geometric structure can be extracted using light scattering, scanning electron microscopy, or confocal reflectance microscopy (21, 68, 69). In particular, our simulations indicate that in models without an area term the *SI Appendix*, section IID), which suggests that local geometry may indeed be sufficient to characterize the large-scale mechanical properties of such systems. Remaining future challenges here include the development of an easy way to compute our geometric coefficients from simple properties characterizing local geometric structure without the need to simulate and to find ways to detect possible residual stresses that may have been built into the gel during polymerization.

Finally, our approach can likely be extended to also include isostatic and overconstrained materials. For example, it is generally assumed that the mechanics of biopolymer networks are dominated by a stretching rigidity of fibers that form a subisostatic network, but that an additional fiber-bending rigidity turns the network into an overconstrained system (12⇓⇓–15, 21, 22). The predictions we make here focus on the stretching-dominated limit where fiber-bending rigidity can be neglected, which is attained by a weak fiber-bending modulus and/or in the more rigid parts of the phase space. A generalization of our formalism toward overconstrained systems will allow us to extend our predictions beyond this regime and thus refine our comparison with experimental data.

## Materials and Methods

The 2D spring networks were initialized as packing-derived, randomly cut networks (9, 27). To improve the precision compared with the cellular models, we created our own implementation of the Polak–Ribière version of the conjugate gradient minimization method (70), where for the line searches we use a self-developed Newton method based only on energy derivatives. All states were minimized until the average force per degree of freedom was less than *A* and *B* and to find the *SI Appendix*, section IVB.

For the 2D vertex model simulations, we always started from Voronoi tessellations of random point patterns, generated using the Computational Geometry Algorithms Library (CGAL) (https://www.cgal.org/), and we used the BFGS2 implementation of the GNU Scientific Library (GSL) (https://gnu.org/software/gsl/) to minimize the energy. We enforced three-way vertices and the length cutoff for T1 transitions was set to

For the 2D Voronoi model simulations, we started from random point patterns and minimized the system energy using the BFGS2 routine of the GSL, each time using CGAL to compute the Voronoi tessellations. Due to limitations of CGAL, configurations were not shear stabilized.

For the 3D Voronoi model simulations, we used the shear-stabilized, energy-minimized states generated in ref. 46 using the BFGS2 multidimensional minimization routine of the GSL.

Details on the different simulation protocols (*SI Appendix*, section IV.

## Acknowledgments

We thank Daniel M. Sussman for fruitful discussions. M.M. and M.L.M. acknowledge funding from the Simons Foundation under Grant 446222, the Alfred P. Sloan Foundation, the Gordon and Betty Moore Foundation, and the Research Corporation for Scientific Advancement though the Cottrell Scholars program and computational support through NSF Grant ACI-1541396. M.L.M. also acknowledges support from the Simons Foundation under Grant 454947 and NSF-DMR-1352184 and NSF-PHY-1607416. K.B. and B.P.T. acknowledge funding from the Netherlands Organization for Scientific Research.

## Footnotes

- ↵
^{1}To whom correspondence should be addressed. Email: mmerkel{at}syr.edu.

Author contributions: M.M., B.P.T., and M.L.M. designed research; M.M. performed research; M.M. and K.B. analyzed data; and M.M., B.P.T., and M.L.M. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1815436116/-/DCSupplemental.

Published under the PNAS license.

## References

- ↵
- Maxwell JC

- ↵
- Calladine CR

- ↵
- Lubensky TC,
- Kane CL,
- Mao X,
- Souslov A,
- Sun K

- ↵
- Zhou D,
- Zhang L,
- Mao X

- ↵
- Mao X,
- Lubensky TC

- ↵
- ↵
- ↵
- ↵
- ↵
- Sheinman M,
- Broedersz CP,
- MacKintosh FC

- ↵
- ↵
- Licup AJ, et al.

- ↵
- Licup AJ,
- Sharma A,
- MacKintosh FC

- ↵
- ↵
- Sharma A, et al.

- ↵
- Feng J,
- Levine H,
- Mao X,
- Sander LM

- ↵
- ↵
- Vermeulen MFJ,
- Bose A,
- Storm C,
- Ellenbroek WG

- ↵
- Shivers JL,
- Feng J,
- Sharma A,
- MacKintosh FC

- ↵
- Shivers J,
- Arzash S,
- Sharma A,
- MacKintosh FC

- ↵
- Jansen KA, et al.

- ↵
- Rens R,
- Villarroel C,
- Düring G,
- Lerner E

- ↵
- ↵
- Rens R,
- Vahabi M,
- Licup AJ,
- MacKintosh FC,
- Sharma A

- ↵
- ↵
- de Cagny HCG, et al.

- ↵
- Baumgarten K,
- Tighe BP

- ↵
- ↵
- Angelini TE, et al.

- ↵
- ↵
- ↵
- Garcia S, et al.

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Matoz-Fernandez DA,
- Martens K,
- Sknepnek R,
- Barrat JL,
- Henkes S

- ↵
- Barton DL,
- Henkes S,
- Weijer CJ,
- Sknepnek R

- ↵
- Yang X, et al.

- ↵
- Moshe M,
- Bowick MJ,
- Marchetti MC

- ↵
- Giavazzi F, et al.

- ↵
- Sussman DM,
- Merkel M

- ↵
- Sussman DM,
- Paoluzzi M,
- Cristina Marchetti M,
- Lisa Manning M

- ↵
- Merkel M,
- Manning ML

- ↵
- Boromand A,
- Signoriello A,
- Ye F,
- O’Hern CS,
- Shattuk MD

- ↵
- Yan L,
- Bi D

- ↵
- Teomy E,
- Kessler DA,
- Levine H

- ↵
- ↵
- Yucht M,
- Sheinman M,
- Broedersz C

- ↵
- ↵
- Düring G,
- Lerner E,
- Wyart M

- ↵
- Woodhouse FG,
- Ronellenfitsch H,
- Dunkel J

- ↵
- ↵
- ↵
- Arzash S,
- Shivers J,
- Licup AJ,
- Sharma A,
- MacKintosh FC

- ↵
- ↵
- ↵
- Cioroianu AR,
- Storm C

- ↵
- Cui B,
- Zaccone A

- ↵
- Broedersz CP,
- Sheinman M,
- MacKintosh FC

- ↵
- Amuasi H,
- Fischer A,
- Zippelius A,
- Heussinger C

- ↵
- Kim S,
- Wang Y,
- Hilgenfeldt S

- ↵
- Ronceray P,
- Broedersz CP,
- Lenz M

- ↵
- Stam S, et al.

- ↵
- Fischer-Friedrich E

- ↵
- ↵
- Lindström SB,
- Vader DA,
- Kulachenko A,
- Weitz DA

- ↵
- Nash SG,
- Sofer A

## Citation Manager Formats

## Sign up for Article Alerts

## Article Classifications

- Physical Sciences
- Applied Physical Sciences

- Biological Sciences
- Biophysics and Computational Biology