Moduli and modes in the Mikado model.

We determine how low frequency vibrational modes control the elastic shear modulus of Mikado networks, a minimal mechanical model for semi-flexible fiber networks. From prior work it is known that when the fiber bending modulus is sufficiently small, (i) the shear modulus of 2D Mikado networks scales as a power law in the fiber line density, G ∼ ρα+1, and (ii) the networks also possess an anomalous abundance of soft (low-frequency) vibrational modes with a characteristic frequency ωκ ∼ ρβ/2. While it has been suggested that α and β are identical, the preponderance of evidence indicates that α is larger than theoretical predictions for β. We resolve this inconsistency by measuring the vibrational density of states in Mikado networks for the first time. Supported by these results, we then demonstrate analytically that α = β + 1. In so doing, we uncover new insights into the coupling between soft modes and shear, as well as the origin of the crossover from bending- to stretching-dominated response.

Mikado networks[1] are a widely studied off-lattice model for the mechanics of semiflexible fiber networks.[2-9] Two-dimensional Mikado networks consist of slender rods with stretching modulus μ and bending modulus κ, which are randomly deposited in a box and crosslinked at each intersection point. At a phenomenological level, the model successfully reproduces important features of real fiber networks, including the affine to non-affine crossover in linear response (discussed below) and the onset of strain stiffening at a critical strain scale.[10-13] The model also serves as a bridge to other soft matter systems where floppiness and non-affinity play an important role, such as solid foams and confluent tissues,[14-16] soft sphere packings,[17-22] and random spring networks.[23-28]

The goal of this article is to determine how low frequency (“soft”) vibrational modes control the shear modulus G in 2D Mikado networks. It is natural to expect the displacements of crosslinks during shear deformation to project strongly on soft modes, which implies a relationship to the modulus. This relationship has the potential to explain the strong non-affine fluctuations that are known to emerge in fiber networks, by connecting them to the floppy (zero frequency) vibrational modes of freely bending (κ = 0) networks, which are inherently non-affine. However, precisely how soft modes set the modulus remains unclear, as there is an apparent contradiction between existing results in the literature. While both G[1,3,4,7] and the typical soft mode frequency ω[5,6] display anomalous power law scaling with the fiber line density ρ (defined below), the best estimates of their respective exponents do not appear to be compatible. In order to describe this discrepancy more precisely, we must first summarize what is known about moduli and modes separately.

Shear modulus. We begin with the shear modulus of bi-periodic networks of N fibers in an area A, as depicted in Fig. 1a. There are three independent, microscopic length scales that characterize such a network. The first is the fiber length f. The second is the mean length of a fiber segment between two crosslinkers. A Poissonian fiber deposition process results in a mean segment length s = π/(2ρ), where the fiber line density ρ ≡ Nf/A is the total length of fiber per unit area. The third length b ≡ (κ/μ)1/2 is mechanical; it represents the natural length over which a free fiber bends. When mechanical properties of networks are expressed in dimensionless form, they must be functions of two dimensionless ratios constructed from these three length scales. We choose b/f and ρf.

Fig. 1

A Mikado network with ρf = 40.0.

The most relevant features of G can be illustrated by varying the bending modulus κ while holding all other model parameters fixed (i.e., varying b/f at fixed ρf), as shown using our own data in Fig. 2a (details below). If we look first at larger values of κ, then eventually the energetic penalty for bending will be so large that the fibers essentially do not bend. Numerics indicate that when deformation is stretching-dominated, the displacements of crosslinks follow the global deformation affinely.[1,3] The shear modulus is therefore given by its affine value G ≃ Gaff, which can be calculated exactly[4] and scales as Gaff ∼ μρ/f in densely crosslinked networks. If instead we look at small κ, bending eventually becomes energetically favorable to stretching, the response becomes strongly non-affine, and G becomes much smaller than Gaff. Hence the shear modulus displays a crossover between stretching- and bending-dominated deformation, which correlate with affine and non-affine response, respectively. The bending-dominated regime is of particular interest, both because affine deformation is comparatively simple to model, and because non-affinity represents strong deviations from the mean response.

Fig. 2

(a) The shear modulus G of a Mikado network plotted versus (b/f)2 for varying ρf indicated in the legend (see text for definitions). The solid line denotes the affine shear modulus for ρf = 48.9. (b) The same data as in (a), now collapsed by plotting G/Gaff as a function of (b/f)2(ρf) with α = 7.2. The dashed line has a slope of 1, and the solid line denotes G = Gaff.

It is evident from Fig. 2a that the value (b/f)* where the bending- to stretching-dominated crossover occurs depends itself on ρf. If we assume power law scaling (b/f)* ∼ (ρf)− for some undetermined exponent α, it is then natural to expect that the dimensionless shear modulus ≡ G/Gaff can be written as a function of the variable z = (b/f)2(ρf). The scaling function has the formwhere z0 is a constant. The approach to ≃ 1 for large z ensures that G ≃ Gaff is recovered, while the form ∼ z for small z is necessary so that G is proportional to κ and independent of μ in the bending-dominated regime. The value of α has been estimated numerically by several groups.† Wilhelm and Frey found α ≈ 5.7.[1] Head et al.[3,4] estimated α ≈ 8 for networks with ρf as large as 120 (comparable to Wilhelm and Frey), but suggested that α will approach a value of 7 when larger values of ρf are accessed. Shahsavari and Picu[7] simulated ρf as large as 300 and found α ≈ 7. Head et al.[4] also presented a mean field scaling argument that predicts α = 7; however, it does not invoke modes and so does not clarify their relationship to the modulus. As shown in Fig. 2b and described in greater detail below, we obtain good data collapse with our own numerical results for α ≈ 7.2.‡

The scaling function (1) follows the discussion by Shahsavari and Picu.[7] Other authors rescale their data differently, but still determine exponents that can be compared to α directly. Namely, α is equal to μ − 1 from ref. 1, and 2 + 2/z from ref. 3,4.

Our simulations access a narrower range in ρf because, unlike prior work, we determine each network's eigenfrequencies, which is significantly more computationally expensive than calculating their shear moduli.

For completeness, we note that Mikado networks undergo a rigidity percolation transition at a critical fiber line density ρc ≈ 5.7.[1,2] Like the prior studies discussed above, we estimate α using data well above ρc. G crosses over to a qualitatively different scaling regime as ρ → ρc from above. This regime is not the focus of the present study.

Soft modes. Let us now turn to soft vibrational modes in Mikado networks. If fibers are allowed to bend freely (κ = 0), then Mikado networks always possess so-called floppy modes.[2,5,8] These are non-rigid body motions that can be performed without doing work, and therefore correspond to vibrational modes with zero frequency. A finite bending stiffness stabilizes floppy motions, lifting them to finite frequency. As floppy motions are intrinsically non-affine, the number and vibrational frequencies of soft modes must be fundamentally related to non-affine fluctuations in mechanical response.

In two seminal papers, Heussinger and Frey (HF) developed what has come to be known as floppy mode theory (FMT).[5,6] The theory (i) explicitly constructs the floppy modes in freely bending Mikado networks, and (ii) uses these modes to build a set of trial modes, and thereby predict the characteristic soft mode frequency ω at small but finite bending modulus using the variational method. Each trial mode is constructed by displacing one of the N fibers a distance δx along its axis. This motion is opposed by the rest of the network, and therefore carries an energetic cost W(δx) = (1/2)Kδx2. The average effective “spring constant” K sets the characteristic frequency, ω ∼ K1/2. By using the crosslink displacements from a floppy mode and the known distribution of segment sizes, HF foundwhere ω0 = (μ/f)1/2. Depending on the specific approximations that are made, FMT predicts β = 6 or β ≈ 5.75.

Rather than directly measuring the frequency ω, HF tested their theory indirectly by measuring the shear modulus. To predict G, they equated the energy of the shear deformation, GAγ2, with the energy NW(δx) stored by the soft modes under shear. However, it is not apparent what the typical axial displacement δx should be. In order to close this gap, they further assumed that the center of mass of each fiber displaces affinely. We shall refer to this as the affine fiber approximation (AFA). It implies that each fiber undergoes an axial displacement δx ∼ fγ relative to its environment. The shear modulus is then

Taken together, eqn (1)–(3) imply α = β. Based on this relation, Heussinger and Frey concluded that the numerical result α = 5.7 from Wilhelm and Frey[1] validates their theoretical predictions for β. However, as noted above, independent studies[3,4,7] conclude instead that α ≈ 7. In our opinion, the preponderance of evidence indicates that α is not equal to the FMT prediction for β.

Hence we have a puzzle. The discrepancy between α and the theoretically predicted value of β indicates that there is a gap in our understanding of non-affine response in the Mikado model – and, by extension, in fiber networks in general. Where is the missing physics? The persistence of this gap is due, at least in part, to the paucity of data for the vibrational density of states (and hence ω) in Mikado networks and related models.[29] Since ω has not been measured directly, one possibility is that β is larger than the HF theory predicts. We test this by presenting the first direct measurement of ω in the Mikado model, as shown in Fig. 3a (details below). We obtain the best collapse of the data for β = 6.1, with reasonable collapse in the range 5.8 ≤ β ≤ 6.5. This result for β is in reasonable agreement with FMT, and also too small to account for the discrepancy between α and β.

Fig. 3

(a) The characteristic frequency ω of soft vibrational modes in a Mikado network plotted versus (b/f)2(ρf) with β = 6.1. See Fig. 2 for legend. The dashed line has a slope of 1. (b) The ratio of the dimensionless shear modulus shear modulus G/Gaff to the characteristic frequency ω is not constant when plotted versus fiber line density, which implies α ≠ β.

The second possible explanation for the discrepancy, which we favor, is that the affine fiber approximation is incorrect, and hence that the relation α = β is unjustified. This view is supported by Fig. 3b, which plots the ratio (G/Gaff)/(ω2/ω02). eqn (3) predicts this ratio to be independent of the fiber line density, but instead we find a roughly linear trend. In the following sections, we will present evidence that in factThis relation is in reasonable agreement with our numerical results in Fig. 2 and 3. More importantly, we will derive (4) in two different ways, neither of which invokes the AFA. This resolves the apparent discrepancy between moduli and modes in the Mikado model.

Shear response in the Mikado model

We now define the Mikado model in greater detail, and demonstrate how to calculate the shear modulus from the crosslink displacements in response to an imposed shear stress.

Mikado networks consist of N filaments of equal length f, which are deposited randomly in a biperiodic square system of linear size L and area A = L2. Their positions and orientations are drawn from uniform probability distributions. If, during the deposition process, two filaments intersect, they form a crosslink around which the filaments can freely rotate. The crosslinks cannot be broken. In the initial condition, i.e. after deposition and prior to shear, the fibers are straight lines subdivided into co-linear segments by the crosslinks. We trim so-called “dangling ends,” i.e. segments at the end of a fiber that are attached via a single crosslink, because they neither stretch nor bend. We choose units of length and energy such that f = 1 and μ = 1 and fix A = 9f2. Nevertheless, we continue to write f and κ in scaling relations for transparency.

The energy of a network configuation is given by

Here 〈ij〉 denotes the segment connecting crosslinks i and j, and 〈ijk〉 labels two consecutive segments 〈ij〉 and 〈jk〉 on the same fiber. The quantity δ〈 is the change in the length of segment 〈ij〉 from its value in the initial condition. The angle θ〈 measures the deviation from co-linearity of adjacent segments 〈ij〉 and 〈jk〉, and 〈 is the average of 〈 and 〈.

By construction, the network is initially in a stress-free state. When subjected to a shear stress σ, the box undergoes a shear strain γ and the crosslinks displace from their initial positions. We consider a simple shear in which the lattice vectors of the initially square unit cell becomewhere F̂ is the deformation gradient. The crosslink displacements can be expressed as a sum of affine and non-affine terms. Denoting the position of crosslink i as X⃑, its affine displacement is A⃑ = F̂X⃑. Collecting the 2N× Cartesian components of the crosslink positions in a vector |X〉 (and similarly for other quantities), the crosslink positions after shearing are |X〉 + |A〉 + |U〉, where |U〉 represents the non-affine part of the displacements. While each element of |U〉 is an independent degree of freedom, specifying the shear strain selects the value of |A〉 on each crosslink. Hence a shear displacement of the network is described by 2N× + 1 degrees of freedom in total.

The response to shear can be determined by expanding the energy (5) in the components of |U〉 and γ. In the harmonic approximation, the result isHere Ĥ is the “standard” Hessian, i.e. the 2N× × 2N× matrix of second derivatives of the energy with respect to each Cartesian component U of the crosslink displacements |U〉,Similarly, the vector |Ξ〉 is defined viaIt represents the force imbalance on each crosslink that results from affinely deforming the entire system (unit cell plus crosslinks). Finally, the affine modulus isWe refer to the entire matrix on the lefthand side of eqn (7) as the extended Hessian Ĥext.

The shear modulus can be determined numerically by selecting a value of the shear stress σ and inverting eqn (7) numerically to find |U,γ〉. The modulus can then be read off from G = σ/γ. The results of this procedure, which were already discussed in the Introduction, are plotted in Fig. 2 after averaging over 100 samples per condition. To estimate the exponent α, we introduce a trial value αtrial and select the n data points (z,) satisfying z(αtrial) < 0.1. According to eqn (1), the variable g = /z should approach a constant value, independent of z. Deviations from this expectation are quantified by the cost function , which quantifies the relative amplitude of fluctuations about the mean ḡ. At αtrial = 7.2, the function obtains its minimum value Emin = 8.2 × 10−3, meaning data collapse to within less than 1% of the average. E remains within 50% of its minimum for 6.9 ≤ αtrial ≤ 7.6.

Density of states and soft modes

To make further progress, we now follow the approach developed in ref. 20,27,30, and consider the eigenmodes of Ĥext. These serve as a convenient eigenbasis in which to express the quasistatic displacements of a network's nodes. Eigenmodes of the Hessian are also eigenmodes of the dynamical matrix for the case where each crosslink is assigned unit mass, which motivates our choice to describe them as (undamped) “vibrational modes.” Nevertheless, eigenvalues of the Hessian have units of stiffness (force/length), rather than (frequency)2. The same Hessian can be used to model, e.g., overdamped relaxations in the presence of a solute, see ref. 20,27.

Each mode |U,Λ〉 satisfies the eigenvalue equationwhere ω is the eigenfrequency of mode n. We then express the shear response |U,γ〉 as a superposition of the modes,To better understand the quantity Λ, it is useful to note that the shear strain isFor two modes with equal coefficients c, the mode with the larger Λ will contribute more strongly to the strain. As the sign of Λ can fluctuate, its square can be understood as a measure of how strongly the mode couples to shear.

We now solve for G = σ/γ by inserting eqn (12) in (7) and exploiting the orthonormality of modes to findThe restriction to positive frequencies excludes trivial rigid body modes. We will evaluate (14) for G in Section 4. But first, we focus on the statistics of modes.

It is useful to rewrite eqn (14) in the continuum limit. To do this we assume that Λ2N2 is independent of N (which will be justified below) and define Λ(ω)/N2 to be the average of Λ2 in the interval [ω,ω + dω]). Similarly, we introduce the density of states D(ω), which is the probability density of ω in the same interval. We then have

The Mikado density of states D(ω) is plotted in Fig. 4 for one value of ρ and varying κ. It features two distinct bands. The lower band shifts with varying κ, while the upper band remains unchanged. We conclude that the lower band contains bending-dominated modes, and the upper band contains stretching-dominated modes.

Fig. 4

The Mikado density of states for ρf = 20.0 and κ = 10−14…10−8 in decade steps. The frequencies ω and ω are indicated with arrows for the lowest value of κ. The dashed line has slope − 3.

The bending-dominated band is characterized by its peak frequency ω, which we identify with the frequency scale studied by HF. We have estimated ω by fitting log D to a parabola in log ω in the vicinity of the peak. These data are plotted in Fig. 3a. As noted in the Introduction, we find that ω2 scales with (bf)2(ρf). Estimating β in the same way as α above, we find the best collapse for βtrial = 6.1, with Emin = 0.01. The cost function remains within 50% of its minimum for 5.8 ≤ βtrial ≤ 6.5.

The stretching-dominated band is characterized by a low frequency shoulder at a frequency ω. While ω is independent of κ, it shifts to lower frequencies when the fiber line density is increased, as shown in Fig. 5a. We discuss ω further in the following section. The large-ω tail of the stretching-dominated band is due to longitudinal vibrations ω‖ ∼ (μ/)1/2 of individual fiber segments with length . These frequencies are spread out by the distribution of segment lengths P(), which is exponential due to the Poissonian deposition process. The tail is then D ∼ P[(ω‖)]|d/dω‖| ∼ ω−3 (Fig. 4, dashed line).

Fig. 5

(a) Density of states for κ = 10−14 and varying fiber line density. (b) The same data, now with ω rescaled by the prediction of eqn (18) for ω.

Two types of soft modes

We will demonstrate that α = β + 1 in two separate ways, neither of which invoke the AFA. In the following section, we will directly calculate G in terms of the peak frequency ω in the bending-dominated regime. But first, we present a scaling analysis of the characteristic frequencies in the density of states. The main question will be determining when the bending- and stretching-dominated bands merge, i.e. when ω becomes comparable to ω. Therefore we must first determine how ω scales with fiber line density.

The stretching-dominated band is independent of the bending modulus κ, and therefore persists in freely-bending networks with κ = 0. The vibrational modes of central force networks have been studied extensively,[23-25,27,28,31-33] and we will now adapt these results to Mikado networks to determine the scaling of ω.

The mean coordination Z of a Mikado network obeys an upper bound Z ≤ Zc = 4, where Zc is known as the central force isostatic point.§ It represents the lowest coordination where there are enough springs to constrain all possible non-rigid body motions of the network. It is straightforward to show that ΔZ ≡ Zc − Z ∼ 1/(ρf) in Mikado networks. Therefore densely crosslinked, freely bending Mikado networks are central force networks close to, but below, their rigidity transition.[2,8,29]

The value Zc = 4 is the prediction of a mean field counting argument due to Maxwell. While not all random networks obey mean field counting, Mikado networks do.

Central force networks on either side of Zc display an abundance of soft modes whose spatial character strongly resembles floppy motions,[24,33]i.e. the relative motion of each pair of nodes i and j connected by a spring is predominantly transverse to the unit vector n̂ pointing from i to j. We shall refer to these as stretching-dominated soft modes, to distinguish them from the soft modes of FMT. The softest stretching-dominated soft modes have frequency ω. When Z < Zc, the density of states below ω is gapped, i.e. there are no modes between ω and the floppy modes at ω = 0.[24] The dependence of ω on ΔZ can be calculated with variational arguments reminiscent of FMT,[33] or with effective medium theory.[24] The result iswhere keff is an effective spring constantHere we have used the fact that the mean segment length scales as 1/ρ. The typical stretching-dominated soft mode is thereforeThis prediction is verified in Fig. 5b, where the low frequency shoulder in the density of states is collapsed for varying ρf by rescaling the frequency ω with the prediction of eqn (18).

Having determined the scaling of ω, we can now motivate the crossover between bending- and stretching-dominated shear response. We expect network response to project strongly on the softest modes. Therefore the response to shear will be dominated by fiber bending when the density of states has a well-defined bending-dominated band. This condition is satisfied so long as ω ≪ ω. The crossover from bending- to stretching-dominated response will occur when ω and ω become comparable,This defines the bending-to-stretching crossover in terms of the exponent β. In the Introduction, we showed that the crossover occurs when (b/f)2(ρf) = z0. As these two conditions must be equivalent, we conclude that α = β + 1. This is our first main result.

Directly relating the shear modulus to soft modes

We now derive the relation α = β + 1 by a different but compatible route. Namely, we calculate G in the bending-dominated regime directly from the sum rule (14).

We begin by rewriting the sum rule as a bound on G, by restricting the sum to the O(N) soft modes,We expect G to follow this bound closely, because the soft modes are the lowest non-trivial eigenfrequencies. Assuming each soft mode has the same frequency ω2 with a corresponding coupling strength Λ2 gives

In order to evaluate G, we must determine Λ2. Acting on eqn (11) from the left with 〈0,1|, solving for Λ2, and taking the thermodynamic limit givesAssuming that correlations between the force imbalances |Ξ〉 and the displacements |U〉 of a soft mode are negligible, the squared inner product isThe modes are normalized, so their mean squared displacement |U⃑,i|2 is O(1/N×).¶

Normalization guarantees |U⃑,i|2 ∼ O(1/N×), provided the mode is extended. As the floppy modes in freely bending Mikado networks are localized (motion is restricted to a central fiber and the secondary fibers with which it is crosslinked), one might worry that the soft modes are localized, as well. However, this is not the case. When the bending modulus is small but finite, the secondary fibers trigger motion in ternary fibers, and so forth. HF demonstrated that the axial displacements of the non-central fibers are of the same order of magnitude as the central fiber.

The quantity requires closer consideration. Recall that |Ξ〉 represents the force imbalance on each crosslink after an affine deformation of the network.|| Affine displacements do not bend fibers, so the forces they generate are only due to stretching. Let us consider a single fiber. Affine shear changes the length of each segment by an amount δ〈 = A(θ)γ, where A(θ) is a prefactor that depends only on the orientation of the fiber and is therefore the same for all segments on the same fiber. The force in segment α is f〈 = μδ〈/〈 = μA(θ)γ. Note that the force is independent of the segment length 〈. As a result, affine shear induces an identical force in each segment of a fiber. The segments are co-linear prior to shearing, and therefore the net force on a crosslink vanishes unless the crosslink is at one of the two ends of a fiber (or is attached to a dangling end, if they have not been trimmed). In other words, is zero on interior crosslinks and O(μ2) on terminal crosslinks. As there are O(N) terminal crosslinks, rather than O(N×), the sum in eqn (23) can be rewrittenReplacing the subscript n in eqn (22) with κ to indicate a typical soft mode, we findThe fact that Λ2 scales inversely with N2 guarantees that the shear modulus is intensive. Crucially, the coupling strength also depends on ρf. Inserting (25) in (21), we findfrom which it follows that α = β + 1. This is our second main result.

It may seem odd that the force imbalance due to an affine shear plays a role in determining G in the bending regime, where deformations are strongly non-affine. This can be understood by thinking of shear as a two-step process. In the first step, the system is deformed affinely, which, due to disorder, generates force imbalances on nodes. In the second step, the nodes displace non-affinely to relieve these force imbalances.

Discussion

We have demonstrated how soft modes determine the shear modulus of densely crosslinked Mikado networks in the bending-dominated regime. In so doing, we have resolved an apparent discrepancy in the literature between the predictions of floppy mode theory and independent measurements of the shear modulus. We find no reason to question the predictions of FMT regarding the scaling of the characteristic frequency ω in eqn (2). However, eqn (3) and the prediction that α = β, are not correct. As these relationships rely on the affine fiber approximation, the natural conclusion is that the AFA does not hold. Our results imply that the typical fiber displacement during shear is δx ∼ (ρf)(fγ = (ρf)1/2fγ, which is much larger than the value fγ predicted by the AFA. We have been unable to find an intuitive explanation for this scaling relation.

It is natural to ask why the early work of Wilhelm and Frey yielded an estimate α ≈ 5.7 that is significantly smaller than subsequent estimates. In our opinion, this is because they fitted data for ρ > 15 to the form G/Gaff ∼ Δρ, where Δρ = ρ − ρc (recall ρc ≈ 5.7). While the difference between using ρ and Δρ in the scaling relation becomes negligible sufficiently far above the percolation threshold, fitting to a power law in Δρ gives systematically lower estimates of α. We have shown that FMT predicts that G scales as a power of ρ in the bending-dominated regime.

Our work has generated two additional insights into Mikado network mechanics. First, bending-dominated response correlates with the existence of a well-separated band of bending-dominated soft modes. We have shown that Mikado networks also possess stretching-dominated soft modes. The bending-dominated regime ends, and the gap between bending- and stretching-dominated modes closes, when the characteristic frequency of stretching-dominated soft modes, ω, becomes comparable to the bending-dominated soft mode frequency ω. The frequency ω is set by the network's proximity to the central-force isostatic state. Our results therefore indicate that the isostatic point plays a governing role in the bending- to stretching-dominated crossover in 2D Mikado networks. This is distinct from fiber networks in 3D, which can undergo a similar crossover while remaining well below the central force isostatic point.[34]

A second insight is that bending-dominated soft modes couple to shear in a non-trivial way. This coupling is encoded in the system-size dependence of Λ2. Simple dimensional analysis reveals that Λ2 must be inversely proportional to the square of an extensive quantity. We have shown that this squared extensive quantity is an admixture of the number of fibers and the number of crosslinks, Λ2 ∼ 1/(NN×). This is a consequence of the co-linearity of crosslinks along a fiber. This prediction is corroborated by numerical studies in which the crosslink positions in the initial condition are gradually perturbed, which qualitatively changes the scaling of G with fiber line density.[6]

Our results suggest several directions for future work. In addition to elastic moduli, the density of states can be used to calculate linear response properties such as the storage and loss moduli,[20] sound transmission,[35] and thermal expansion coefficients.[26] It can also predict the onset of nonlinear properties such as strain stiffening and shocks.[27,36,37] Finally, a detailed understanding of structure–property relationships can facilitate the design of network architectures with tunable properties, such as moduli and failure onset.[13,38-42]

Conflicts of interest

There are no conflicts of interest to declare.