Structural hierarchy confers error tolerance in biological materials

Geometrical Characteristics of the Model System

We consider an extension of the well-studied dilute, triangular lattice in 2D. Nodes arranged in a triangular Bravais lattice are connected to nearest neighbors, and bonds are then removed at random such that some fraction, referred to as the bond portion, p, remains. The infinite triangular lattice has a connectivity of six bonds per node when p=1, while in 2D Maxwell counting dictates a minimum connectivity of four bonds per node; thus, the infinite, dilute triangular lattice should lose stiffness when p falls below 23. Lattices of finite size would require a slightly higher bond portion, due to the presence of underconstrained nodes at the boundary. This prediction has been thoroughly confirmed for ball-and-spring networks, via simulation and mean-field theoretic approaches (10, 13⇓–15).

We create hierarchical triangular lattices through an iterative process, in which the bonds of the lattice at one length scale are in turn crafted from smaller-scale triangular lattices. In principle, this process can be iterated ad infinitum; in practice, if the number of large bonds is held constant, the total number of nodes grows geometrically with the number of hierarchical levels. This places a practical limit on the number of levels that can be considered in simulations. We have numerically constructed and simulated lattices with one, two, and three levels of structural hierarchy (Fig. 1); bond portion is independently set on each hierarchical level.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 1.

(A) A dilute triangular lattice with one level of structure. Missing bonds are indicated with dashed lines. (B) A two-level triangular lattice, with each bond replaced by a smaller-scale, dilute triangular lattice. (C) An extension to three levels of structural hierarchy.

Hierarchical Model of Stiffness

We propose to model the stiffness of our networks by generalizing the scaling law proposed by Garboczi and Thorpe (10), for a single-scale, diluted triangular lattice. Garboczi and Thorpe (10) found that, for the ball-and-spring case, components of the elastic constant tensor should have the formK=kp−pc1−pc,p≥pc0,p<pc,[1]where pc is the minimum bond portion necessary for marginal stiffness, and k is the value of the modulus when the network is fully connected. We propose to describe large-scale bonds using an effective stiffness with the form of Eq. 1 and introduce plarge and psmall, the portion of bonds retained on the large and small scales. Because of the finite width of large-scale bonds, we do not assume that the minimum small- or large-scale bond portions needed for marginal stiffness are 23. The stiffness of a large-scale bond will be inherited from its small-scale structure, such that the overall stiffness scales asK=kplarge−pc,largepsmall−pc,small1−pc,large1−pc,small,[2]where K is tensile stiffness and k is the stiffness for a network fully connected on all scales. Now consider a general network with N distinct length scales. The stiffness ki of bonds at scale i will be inherited from the structure at scale i−1, so thatki∝pi−1−pc,i−1.[3]The overall stiffness for some general number N levels of structural hierarchy will then beK=k∏i=1Npi−pc,i1−pc,i.[4]

Simulation Procedure

Networks were simulated in 2D, with ball-and-spring interactions between pairs of connected nodes. To measure stiffness, nodes along the tops of networks were uniformly displaced along the vertical direction, after which the y coordinates of the top and bottom nodes were fixed. Next, the x coordinates of top and bottom nodes, as well as both coordinates of all other points, were relaxed using the fast inertial relaxation engine algorithm (28). Each bond is modeled as a fiber that resists stretching with a 1D stretching modulus, μ. Let each bond be a parametric curve r(s), where s ranges from 0 to the length of the bond, and let u→(s) be a field describing the displacement of a point on the bond in the strained state. The energy for a bond of length l is described by the functionalU=∫0lμ2du→ds2ds.[5]A uniform, 1D stretching modulus of unity was assigned to each bond, and nodes were relaxed until the rms residual force in the network was less than 1×10−10 in units of stretching modulus in the case of the one- and two-level networks and 1×10−9 in units of stretching modulus for three-level networks. We extract the tensile stiffness by fitting plots of elastic energy vs. strain to parabolas.

Simulation Results

We simulated one-, two-, and three-level lattices with a wide range of bond portions and measured their stiffness (Fig. 2). Large-scale bonds are 10 small-scale bonds long and contain three rows of small-scale bonds. Before comparing to our model, we must identify each critical bond portion for one-, two-, and three-level lattices.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

(A) Simulated stiffness plotted vs. bond portion for one length scale. (B) Heat map of simulated stiffness as a function of small and large bond portion for a two-level network. Stiffness is plotted in simulation units, as indicated in the key. (C) Heat map of simulated stiffness for a slice in bond portion space for three-level networks with full connectivity on the largest scale. Stiffness is plotted in simulation units, as indicated in the key. (D) Heat map of simulated stiffness for a slice in bond portion space for three-level networks with full connectivity on the smallest scale. Stiffness is plotted in simulation units, as indicated in the key. (E) Stiffness normalized by maximum stiffness for networks with one, two, and three levels of structure.

For the one-level lattice, we recover the expected linear relationship between stiffness and bond portion (Fig. 2A). We find the critical bond portion to be 0.670 (SI Appendix).

For the two-level lattice, we find that stiffness increases as either the small or the large bond portion is increased (Fig. 2B). To determine the threshold bond portion for the large (small) scale, we find those points in bond portion space for which the small (large) bond portion is equal to 1. We then fit a plot of large (small) bond portion vs. stiffness and fit these data to an equation of the form Eq. 1. We find that the critical bond portions are 0.83 and 0.58 for the small and large scales, respectively (SI Appendix). Initially, it may be surprising to find that the critical bond portion on the large scale is less than 0.67. However, as large-scale bonds are endowed with a finer-scale structure, they acquire an effective bending stiffness, rather than being governed strictly by harmonic, central force interactions. Notably, networks with bonds possessing bending stiffness have been demonstrated to be rigid even below the classic isostatic point (14, 15). Thus, it is crucial that pc be directly measured and not assumed from Maxwell counting.

For the three-level lattice, we find that stiffness increases as any bond portion is increased (Fig. 2 C and D and SI Appendix). Following an approach similar to that used for the two-level lattice, we find the critical bond portions are 0.83, 0.72, and 0.62 for the small, medium, and large scales, respectively.

To test our model, we compare the stiffness obtained from simulations to the stiffness predicted by our model. In Fig. 2E, we plot stiffness values from our model against stiffness values from simulations for one, two, and three levels. For ease of comparison, stiffness values for an N-level network are normalized to the maximum attainable stiffness for an N-level network. Each point in the plot compares a simulation result with its model prediction counterpart. There are no free parameters in our model; we normalize stiffnesses by the maximum values for one-, two-, or three-level lattices and use the identified critical bond portions. We find remarkable agreement; linear fits between simulated and predicted quantities have r2 values and slopes, respectively, of 0.983 and 1.003 for one level, 0.989 and 1.01 for two levels, 0.988 and 0.98 for three levels, and 0.993 and 1 for all data combined.

Experimental Prospects

While the attempt to manipulate actual tissues as we have manipulated our networks may prove elusive, we envision, as an accessible proof-of-principle experiment, cutting networks described in this work from elastic sheets. Such sheets would ideally have an out-of-plane thickness that is large in comparison with the in-plane width of bonds, to frustrate buckling, as described in ref. 29. Existing work has also established tapering of bonds at junctions as a useful technique for controlling the ratio of bending to stretching stiffness (30).

Model Applicability

We next performed a series of tests and analyses to determine how broadly our model can be applied. As a first step, we varied the size and aspect ratio of large-scale bonds. Up to this point we have shown results for networks in which large-scale bonds are 10 small-scale bonds long and contain three rows of small-scale bonds; however, we have also considered two-level networks whose large-scale bonds are 20 small-scale bonds long and contain five rows of small-scale bonds. In both cases, we find similarly strong agreement. More details are provided in SI Appendix.

While these two systems captured the same phenomena, we cannot test all possible configurations and must instead carefully consider the circumstances under which our model is applicable. Large-scale bonds can be too wide or too narrow. If the width of large-scale bonds is too great in comparison with their length, then the resulting structure will effectively be a single-scale lattice with a regular array of small perforations. If large-scale bonds contain two or fewer rows of small-scale bonds, the coordination number of the network with all bonds present will be four, making the network susceptible to immediate loss of stiffness upon dilution of bonds. Moreover, stretching of networks will lead to transverse contraction, producing a lateral load on horizontal bonds. If large-scale bonds buckle, our picture of a primarily stretching-stabilized network will be invalid. We here propose methods to identify the range of aspect ratios for which our model should be expected to apply.

We performed tensile simulations to determine the change in the displacement field of a single-scale lattice under tension due to the presence of a hole of width w. Note that a two-scale lattice with bonds of length l and width w, in units of small-scale bond length, may be viewed as a single-scale network with a periodic array of equilateral triangular holes with side length s=l−w3. We argue that our model, with a separation of scales, is appropriate when the change to the displacement field is of the order of the applied strain times the small-scale bond length over a length equal to or greater than the separation between neighboring holes. This analysis, explained in detail in SI Appendix, yieldsl/w≳2+3/w.[6]Large-scale bonds may also be too narrow and thus buckle. Using the theory of slender rods as discussed in ref. 31, we find buckling will become a concern when the ratio of length to width is of the order of the square root of the applied strain (ε); thus, to avoid buckling, the ratio of length to width must satisfylw≲ε−1/2.[7]In our case, large-scale bonds have length 10 and width 3, in units of small-scale bond length, for an aspect ratio of 5.77. This is above the limit separating hierarchical networks from perforated sheets (3.73). Further, as strains used in simulations do not exceed 5×10−3, this aspect ratio is below the limit at which buckling may be expected to occur (14.1). This places us comfortably within the constraints on the applicability of our model discussed above.

Thus, our proposed formulation for stiffness of a hierarchical lattice is expected to apply to a wide range of networks. Crucially, we observe that it accurately describes 1-, 2-, and 3-level lattices (Fig. 2E), suggesting that the smallest length scale can always be replaced by a network of even smaller bonds, and the stiffness will remain the product of all excess bond portions. This general formulation facilitates investigation of highly hierarchical (e.g., 10-level) lattices.

Structural Error Tolerance

Now that we have obtained a general relationship for the stiffness of a hierarchical structure, we are primed to consider the possibility of random errors in assembly, a likely complication in any real assembly process. We focus in particular on how stochastic deviation from a targeted set of connectivities (on all relevant length scales) results in a deviation in network stiffness. We consider two distinct regimes. In the first case, we consider a target point near the isosurface along which stiffness vanishes. In the second one, we consider a target point in bond portion space far from both the limiting case of full connectivity on any length scale and the contour of vanishing stiffness.

We first consider target points in bond portion space corresponding to marginally stiff structures. Such points are of interest, as highly compliant materials have critical biological roles (1) and are attracting increasing attention for applications (11). Such materials typically must not be susceptible to critical transitions in their elastic moduli as a result of small fluctuations in their fine-scale structure.

We use a numerical technique to estimate the distribution of stiffness arising from random errors. First, we choose a nominal point in bond portion space, such that a certain target stiffness is achieved, and the excess bond portion is the same for each length scale. We next add Gaussian random noise to the bond portion on each length scale. The stiffness of the resulting “noisy” point is then estimated by means of an interpolated function computed from simulation data for one- and two-level lattices and a fitted model for three-level lattices (Fig. 2); for lattices with more than three levels, we use Eq. 4. This process is carried out for 50,000 trials.

Strikingly, the variance in stiffness is greatly reduced with each additional level of hierarchical structure (Fig. 3A). Further, the stiffness distribution for the single-level network exhibits a large peak at zero, which is absent in the stiffness distribution of the two- and three-level networks. Thus, despite having opportunities for errors at three separate stages of assembly, the three-level network is more reliably constructed than the one-level network and can more reliably avoid stochastically generating a floppy network.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

(A) Stiffness probability distribution functions (PDFs) estimated from histogram data for networks with one, two, and three levels of structural hierarchy. Points in 1D, 2D, and 3D bond portion with the same minute nominal stiffness were chosen, and Gaussian random variables were added to each bond portion. Note the spike in the PDF for the one-level network at zero stiffness. In this case, we consider a nominal stiffness of 0.025 and Gaussian random noise for each bond portion with zero mean and a SD of 0.005. Other cases are addressed in SI Appendix. (B) Relative error in stiffness vs. levels of hierarchy is plotted for K¯=0.001, σ=0.0001. As additional levels of structural hierarchy are added, the relative error in the tensile stiffness decreases precipitously at first, and the effect saturates at a certain number of levels. Provided the assumptions leading to Eq. 9 hold, our analytical theory and numerical approach are in close agreement (r2≈1). Here, the product of excess bond portions is 0.001, and the noise has amplitude 0.0001 on each scale. (C) The number of all small-scale bonds in a network, divided by the area enclosed by the outer perimeter of the network, is shown vs. number of hierarchical levels for stiffness values of 0.02, 0.03, and 0.04, in units of stretching modulus over small-scale bond length. Networks were chosen to have the same bond portion on all three levels. In each case, increasing hierarchy leads to markedly lower density of small-scale bonds, attesting to the ability of hierarchy to confer both robustness and efficiency.

Next, we seek a general understanding of target points far from any boundary in bond portion space (see SI Appendix for detailed derivation). For N levels, there is a nominal excess bond portion, pe, for each level; we consider identically and independently distributed deviations from the nominal bond portions according to a normal distribution with zero mean and standard deviation (SD) σ. Referring to Eq. 4, we define the reduced stiffness:K¯=Kk∏i=1N1−pc,i.[8]With this definition we find the expected deviation in the stiffness of a network with N hierarchical levels to beΔK¯K¯≈NσK¯1/N,[9]where the approximation holds when σ≪pe. For a target K¯, this functional form predicts the optimal number of levels to beN*=⌊−2⁡lnK¯⌋,[10]where “⌊⌋” denotes the floor.

To test our analytical result, we again use the above numerical approach to calculate SD in stiffness for networks with 1–15 hierarchical levels. For lattices with more than 3 levels, we use Eq. 4 to estimate the stiffness of the resulting noisy point. We find very good agreement between the numerically generated data and our analytical prediction (r2≈1 for the case shown in Fig. 3B).

The above derivations were performed assuming identically and independently distributed random errors in bond portion; however, the results presented above do not strictly require such stringent conditions. Similar robustness against fluctuation in stiffness can occur for error rates that vary on different scales and for distributions of random errors in which errors in bond portion on different length scales are correlated (see SI Appendix for an in-depth treatment).

Interestingly, investigating networks with different error rate distributions allows us to identify a useful heuristic. Consider a two-level network with the same error rate in its large-scale bond portion as a one-level network. The two-level network will have a smaller variance in its stiffness than the one-level network as long as its small bond portion error rate is less than three times larger than its large bond portion error rate (note that this argument also works if the roles of the small and large scales are swapped; see SI Appendix for details).

These observations do not simply arise from the generic properties of random variables. Although a relation of the form Eq. 9 with N=1 would hold for a general product of random numbers, the dependence upon N of the divisor indicates a scaling of the relative fluctuation specific to a stiffness with the functional form of Eq. 4. To produce a supple solid of a particular stiffness, the portion of bonds that must be retained on each scale to achieve this stiffness increases with increasing hierarchy. This can ensure that the connectivity of the network on each scale exceeds the threshold connectivity by an amount large in comparison with the the typical fluctuations in connectivity.

We further find that hierarchical structures offer superior economy of material, in keeping with previous studies (16). To quantify this benefit, we define the number density, ns, of small-scale bonds as follows: Let ls be the length of a single small-scale bond, Ns denote the total number of small-scale bonds, and A denote the area enclosed by the perimeter of the network. Then, ns is given byns=NsA[11]and has units of ls−2. We regard small-scale bonds as the fundamental structural units of a network, such that ns provides a measure of the number of connections per unit area needed to achieve a target set of mechanical properties. To compare one-, two-, and three-level networks, we chose target stiffness values and identified the points in bond portion space with equal bond portions on each level corresponding to these stiffness values. We then estimated accompanying values of ns for networks with these bond portions by interpolation. For each stiffness chosen, there was a sharp decline in the number of bonds per unit area needed with increasing numbers of hierarchical levels. This result emphasizes that hierarchical materials’ greater error tolerance is a benefit inherent to their geometry, rather than deriving from redundancy. A nested, modular assembly process thus enables the creation of a structure that is both more reliable and more frugal.

Discussion and Conclusion

Contrary to expectation, the elastic moduli of hierarchical materials are more reliably controlled than the elastic moduli of materials with one relevant length scale. Thus, with respect to random errors, it is easier to make hierarchical structures than to make single–length-scale structures. This finding may have wide-ranging implications for evolutionary biology and materials design.

It may seem that evolving progressively larger, more complex bodies is accompanied, and impeded, by a growing assortment of mechanical challenges. To the contrary, this work provides evidence that adding hierarchical complexity can actually reduce stochastic variation in material properties. This effect decreases the need for coevolving error-correcting mechanisms, thus facilitating the evolution of new traits that are “good enough” for an organism to survive.

Upon successful assembly, structural hierarchy is known to endow materials with a host of desirable properties that are unattainable with single–length-scale structures. Our finding that hierarchy also reduces susceptibility to stochastic errors suggests that bottom–up production processes for synthetic hierarchical materials may produce finished products which perform as intended despite modest, but discernible errors.

While we have explicitly addressed 2D, lattice-based networks in this study, we anticipate that our findings will have broader applicability. Feng et al. (13) have demonstrated that an analogous model to Eq. 1 describes the dependence of the elastic moduli of a dilute face-centered cubic lattice in 3D on the portion of bonds retained. Moreover, connectivity has been established as a useful control parameter for predicting the mechanical properties of many disordered solids, such as jammed packings of hard (32, 33) and soft (34) particles and semiflexible biopolymer networks (21).

More generally, we have broadened the scope of Maxwell’s visionary means of characterizing frames to account for structural hierarchy, facilitating the understanding of cases in which it has proved difficult to relate materials’ emergent properties to their fine-scale structure. While recent computational advancements have enabled study of biological macromolecules over experimentally relevant time scales (35), comprehensive understanding at the level of an organism demands a coarse-graining procedure for which our model may offer a useful road map. A generalized counting heuristic may also offer a means of expediting feasible, yet cumbersome calculations in materials design. Hierarchical materials necessarily have many design attributes, but our accessible model may considerably narrow the search of parameter space needed to reach a goal.