The mechanics of phantom Mikado networks

The elasticity of networks comprised of semi-flexible polymers plays a vital role in regulating the mechanics of both intra- and extra-cellular matrices. The behaviour of this polymer scaffold will depend on the nature and density of cross-linking between constituent fibres. While modelling efforts have investigated the effects of cross-link density in biopolymer networks, this is often accompanied by changes in both the fibre density and the network structure. We investigate the elasticity of a quasi two-dimensional Mikado network of elastic rods, in which cross-link density is allowed to vary while polymer density is held constant. In particular, this model is extended by allowing constituent rods to cross without forming cross-links, while polymer density and network geometry are preserved. In doing so, the competing contributions to the shear modulus from cross-link density, mesh size, geometry and polymer density are decoupled. We find that previous scaling laws fail to capture the well-studied transition from bend- to stretch- dominated elasticity as cross-link density is varied. We identify a length scale which relates cross-link density to the transition between affine and nonaffine regimes, and which provides a collapse of simulation data curves for varying cross-link densities.


Introduction
Biopolymer networks are ubiquitous in living matter, from the collagen-dominant extracellular matrix (ECM) to the actin cortex of eukaryotic cells. The mechanics of such networks is vital to the proper functioning of many processes, including cellular motility [1,2], mechanical cell-cell communication [3,4] and stress management in articular cartilage [5]. Within the collagen scaffold of ECM, cell-mediated traction forces may be carried over many cell diameters [6,7], owing to the highly anisotropic fibre mechanics and the rearrangement and realignment of constituent fibres [8]. Within cells, dynamically shifting networks allow for structural stability in the presence of a fluctuating extracellular environment [9]. As stresses derived from the environment are transmitted to the cell membrane, rearrangements in stress fibres allow cells to adapt to various substrate mechanics and geometries [10]. The extracellular matrix itself can range from soft and homogeneous to stiff and aligned [11,12]. Clearly, biopolymer networks exist in many forms, and across many scales, the physics governing their mechanics being fundamental to their biological function. Despite the ubiquity of these polymer assemblies in vivo, a complete understanding of how network structural properties and polymer behaviours can produce the observed local and bulk mechanics is challenging.
A large family of observed biopolymer network systems, including many of those comprised from actin, collagen or intermediate filaments, are characterized as 'semi-flexible' [13]. Here, the tendency of a polymer to collapse into a random coil is slightly overcome by a small but significant bending stiffness, whence fibres can appear straight over long distances. This persistence of tangent-tangent correlations introduces hierarchical structure into the polymer assembly [14,15] and can result in more exotic elastic behaviours [16]. Furthermore, semi-flexible filaments often possess a profound mechanical anisotropy, being far more compliant in bending than in stretching. Given that the detailed structure of biopolymer networks is complex and intricate, with nontrivial organisation and long-range spatial correlations generated by their fibrous nature, the properties of individual fibres are not sufficient to describe the continuum-scale behaviour [17,18]. Instead, the macroscopic elastic moduli of such gels are also dependent on network properties, such as mesh spacing, network geometry, towards describing 3D network mechanics has been provided through the use of phantom face-centred cubic lattices [42,43]. In such models, while many fibres might pass through a lattice node, not all are considered cross-linked. Instead, only certain fibres form connections, such that the network structure and cross-link density can be somewhat decoupled. In reference to these lattice models, we term the model presented in this manuscript a 'phantom Mikado' model, in which cross-link and line density can be independently specified. This is achieved by relaxing the constraint that all fibre crossings represent cross-links. Instead, a proportion of cross-links with node coordination 4 from a standard Mikado geometry are removed, without otherwise adjusting the fibre constitutive properties or network geometry.
While lattice models have been used to investigate biopolymer network elasticity in many studies, their rigidly defined geometry and topology can limit their applicability in certain scenarios. In particular, regular lattice-based architectures exhibit preferential directions, whence unperturbed edges must adopt one of a small number of predefined orientations. Likewise, fibre segment lengths are constrained so as to only take certain values, rather than existing within a continuous distribution as observed in fully random architectures such as Voronoi or Mikado cases [15]. For these reasons, studies investigating the local micromechanics of biopolymer networks, for example the contractile behaviour of resident model cells [4], might find random architectures more suitable, as well as better representative of the physical biopolymer networks at this scale. There is therefore considerable utility in transferring ideas found within the lattice-based literature into models involving fully random networks, including the Mikado model.
Our objective is to investigate the mechanics and deformations of biopolymer networks in which cross-link and line density are decoupled. In particular, we aim to characterize the transition between the affine and nonaffine regimes for phantom Mikado networks. We determine the length scale where this transition occurs, and clarify how this differs from the standard Mikado model. Further, we aim to investigate whether energetic and displacement measures for affinity coincide within this model.

Model
We consider Mikado networks formed by randomly placing N monodisperse fibres of length L at uniformly distributed orientations within a square domain of area W 2 . We refer throughout to a non-dimensional polymer density ρ=NL 2 /W 2 , as used elsewhere in the literature [44,45]. Briefly, it provides a description of the number of cross-links along a single fibre as ρ;L/l c . On sufficiently large length-scales W, networks appear isotropic and homogeneous, such that the domain size does not affect results. In this study, domains of size W/ L=2.5-15 are employed, which ensured fluctuation errors remained small. When networks are generated, a permanent cross-link is assigned to the intersection of two fibres. Once the Mikado network has been formed, all dangling ends are removed, as these will not affect the static shear modulus of the network. In all the networks considered, the average coordination number was observed to be between 3 and 4, formed from nodes of coordination 2, 3 or 4.
To investigate the effects of cross-link density on network mechanics, further networks are produced from an initial Mikado geometry by removing a given proportion of cross-links. The resulting networks are described by a parameter p, which measures the ratio of initial to final cross-links in the sample as p≔n/n 0 , for cross-link number n, and subscript 0 denoting quantities in the standard Mikado case. In other models, the removal of cross-links usually corresponds to reductions in the overall network density and changes in the network geometry [26,46]. Mechanical stability can become an issue in these models, and it is often not feasible to remove more than fifty percent of model cross-links before network percolation is disrupted [46]. So as not to disrupt the network geometry, and to keep polymer density constant under cross-link removal, only nodes of coordination 4 are removed; if those with coordination 2 or 3 are removed, additional dangling ends are introduced, such that the effective network density has been lowered. Nodes of coordination 4 are removed uniformly at random according to the required extent. Following the removal of these cross-links, the crossing of two fibres is no longer necessarily identified with a cross-link, and many overlapping filaments do not constrain one another. An illustration of cross-link deletion, with cross-links highlighted in red, is shown in figure 1, demonstrating that polymer density is held constant as cross-link density varies. We note that networks of lower polymer density will have proportionally fewer nodes of coordination number 4, and thus proportionally fewer cross-links can be removed under our modelling assumptions.
Network fibres are modelled as shear-flexible beams. In particular, the energy contained within an individual fibre under deformation is given by: Here, E and G f are the fibre elastic and shear moduli respectively, while A and I give the area and second moment of area. The axial and transverse displacements along a fibre are given by u and v, for fibre arc length s. The rotation of the plane normal to the beam axis is given by ψ, while dv ds gives the rotation of the fibre crosssection (equal to ψ for Euler-Bernoulli beams), and du ds gives the axial strain. Finally, K is a shear correction factor, and the integral runs along the contour length of each fibre. The terms in the above energy equation can be partitioned as shown into bending, stretching and shearing modes, as indicated. Cross-links may be modelled as either freely rotating, where the translation of two fibres at the cross-link is fixed, yet rotation is unconstrained, or welded, where the relative orientation of crossing fibres is also fixed. We focus throughout on freely rotating cross-links, noting that differences arising from this modelling choice are slight [47]. As the origins of elasticity in semi-flexible biopolymers do not necessarily derive from classical linear elasticity, instead originating from different mechanisms for distinct species of biopolymer [48,49], we employ general beam sections to remove the notion of area and radius, and instead focus on ratios between the quantities μ ≔ EA and κ ≔ EI. Simulations investigating the effects of varying fibre shear modulus G f in the range , corresponding to a Poisson ratio of between 0 and 0.5, found that our conclusions are insensitive to this choice. We take the fibre shear modulus to be half the corresponding elastic modulus, and unit shear correction factor, giving a total network energy of: where summation is over network fibres. The total energy contained within the fibre assembly is used to produce the two-dimensional bulk shear modulus  through the relation W 2 total 2 2   g = . Each network is subjected to a small simple shear deformation as follows. Nodes on the top boundary are fixed both in rotation and vertical displacement and shifted horizontally to give 0.05% bulk shear strain γ. Nodes on the bottom boundary are fixed in all degrees of freedom. Periodic boundary conditions are applied on all degrees of freedom to nodes on vertical (left and right) boundaries. Computations are performed using the commercial, implicit finite element solver Abaqus/Standard [50]. Each fibre segment, that is the fibre extent between two model cross-links, is discretized with two linear Timoshenko beam elements, and a force balance is naturally maintained on all interior finite element nodes. The section stiffness for the beam elements is derived from equation (1) [51]. . (c) A partially cross-linked Mikado network is shown for the same sample as in (b). In particular, the proportion of remaining cross-links is given by p≔n/n 0 , for number of cross-links n. Longer fibre segment lengths l are observed than in the fully cross-linked case for equivalent polymer density, leading to a larger value of mean cross-link spacing l l c á ñ ≔ .
The relevant parameters governing the model networks are as follows. Alongside fibre number N, rod length L, cross-link number n and polymer density ρ=NL 2 /W 2 , there is mean segment length l l c = á ñ, bending length scale l b =(κ / μ) 1/2 , and cross-link proportion parameter p=n / n 0 .

Results
We now proceed to investigate results pertaining to the deformations of phantom Mikado networks. We begin by investigating the elasticity within the bending-dominated, nonaffine regime, before investigating the affinenonaffine transition.
3.1. Behaviour in the nonaffine regime As Mikado rods are placed within the domain, and the line density increases, a solution of rods becomes crosslinked into a network that acquires a static shear modulus, ceasing to be a liquid. This mechanical stability occurs with geometric network percolation when welded cross-links are employed, and at slightly higher fibre densities when freely-rotating links are used, as in this work [18]. Given the small bending length scales l b ) observed in semiflexible polymers, filaments in this regime are often far more compliant in bending than extension. Following this percolation, the mean cross-link spacing l c decreases with increasing polymer density [52]. As l c becomes small, the segments become relatively more difficult to bend, instead storing energy in stretching modes. Now, if cross-link density is lowered for networks of high polymer density, through the random removal of model cross-links, mean cross-link spacing increases (though the network never becomes a solution under the procedure described for this model). We expect that stretching modes will become unfavourable, and energy will be contained in the more compliant bending modes, whence the network returns to the nonaffine regime. As the above transition is predicted to occur for constant network geometry and constitutive choices, we predict that knowledge of the fibre elastic moduli and network fibre density will not be sufficient to determine the elastic regime in which a fibre network resides in the presence of variable crosslinking.
The nonaffine regime in fibre networks has been discussed previously [17,29], and is often characterized in Mikado networks through either a low polymer density, small bending length scale l b l ≔ ( ) [17], with other forms described in the literature [31]. We expect that, due to the dependence of fibre bending response on the mean segment length, cross-linking should be important to the network response in this regime. As the affine response of a fibre network is independent of fibre segment length, we do not expect a similar sensitivity to the cross-link density for networks within the affine regime.
We investigate nonaffine network behaviour by plotting the shear moduli for networks of two (lower) densities as cross-linking density is reduced, in figure 2(a). We find that decreasing cross-link density corresponds with a rapid reduction in the network shear modulus, as predicted. This effect becomes increasingly pronounced as the number of cross-links continues to decrease, such that the difference in shear modulus as p varies from p≈0.25 to p≈0.15 is far more pronounced than for p=1 to p≈0.9. In the former case, we find that for networks of polymer density ρ≈44.9 there is a reduction in the shear modulus of almost two orders of magnitude as cross-link parameter p changes by just 0.1. As such, relatively small changes to cross-link number can drastically alter the macroscopic elasticity of nonaffine samples, suggesting an important mechanism by which biopolymer network mechanics might be regulated. The changes for networks with almost maximal cross-linking are far slighter, such that these assemblies are relatively robust to small variations in cross-link density. For the lowest density networks considered 12.8 r » ( ) differences in the shear modulus remain slight except for those networks with the lowest extent of cross-linking, where again small changes in the nature of network cross-linking can dramatically alter the bulk shear modulus, as illustrated by figure 2(a), starred markers.
In figure 2(b) we normalize the shear modulus  by that found for maximal cross-linking, that is the shear modulus corresponding to the standard Mikado model 0  . In the lowest density case, ρ≈12.8, scaling with increasing cross-link removal is almost identical across an order of magnitude variation in the bending length l b . Small deviations are observed for the higher density network, with polymer density ρ≈44.9, though the scaling across both densities and all bending lengths considered remains similar, suggesting that this behaviour is conserved throughout the nonaffine regime. We also briefly compare to the behaviour in the affine regime, with large polymer density ρ≈109, in figure 2(b) (Inset). As predicted, networks within the affine regime are insensitive to decreasing cross-link density, displaying only minimal reductions in the bulk shear modulus as p decreases. However, at the very lowest values of p, corresponding to the removal of almost all constituent crosslinks, there is a sudden transition into the nonaffine regime, leading to a large reduction in network stiffness.
To better understand the elasticity of bending-dominated phantom Mikado networks, we briefly investigate a scaling argument, discussed by Head et al [31] and later Heussinger and Frey [53], for 'floppy' fibres within the nonaffine regime. In particular, we can consider networks formed from fibres with zero bending rigidity κ=0. As the networks are subisostatic, a geometric rearrangement of fibres can accommodate the imposed boundary shear without requiring the stretching of any fibres [13]. Now, introducing a small but non-zero fibre bending rigidity κ>0, networks acquire a non-vanishing static shear modulus  that can, by assumption, only depend on bending stiffness κ.
The bending behaviour of fibre segments of length l can be described through an effective transverse spring constant, given as k ⊥ ∼κ/l 3 by considering beam moments. Following an argument from Heussinger and Frey [53], we consider the filament deflections, of characteristic magnitude d, that arise in the translation of a single constituent fibre. These fibre translations occur on the scale of a whole fibre, such that we do not expect that their magnitude has a dependence on cross-link spacing l c . As such, the only remaining geometric length (for floppy fibres) is L, whence the expected energy in a floppy fibre with deflections of characteristic magnitude d∼L can be written: where L/l c gives the number of cross-links per fibre, l min is a minimum length, under which the fibre becomes too stiff to bend, and P(l) is the fibre segment length distribution. Now, as the number of cross-links n within the phantom Mikado network can be written n;NL/l c , compared with the initial number n NL l c , and cross-link removal is random, such that we expect the distribution of segment lengths to remain exponential. Combined, we can rewrite equation (4), using standard asymptotic techniques for small l l c min 0  , as: We can self-consistently determine the minimum length l min by considering the energy balance between a short fibre segment of this characteristic length under deflection, and the translation of an entire secondary filament, whose energy contribution is derived from the activation of the floppy modes there. In particular, ) . Significant reductions in  are observed for lower cross-link densities (decreasing p) in the higher density case ρ=44.9. In particular, the reduction in shear modulus accelerates as cross-link density decreases.
Substituting this minimum length expression back into equation (5), we arrive at the expected energy per fibre: Multiplying the above expression by the fibre number N produces the expected total network energy. In particular, we predict that the shear modulus should scale as: where we have neglected to include strain γ, and used the relationship NL W l 1 c 2 0 . We therefore expect that the shear modulus within the nonaffine regime should scale as p 6 ~. A comparison of the scaling of the shear modulus in the nonaffine regime with this prediction is shown on double logarithmic axes in figure 3. For values of p neither too close to unity nor the lower limit, we find that the shear modulus scales similar to our prediction p 6 . For the highest and lowest values of p, deviations from this scaling behaviour are observed. This is perhaps not surprising, given the assumptions made in the above arguments, which do not necessarily hold for the extreme values of p. In particular, we made the simplifying assumption that cross-link removal p occurred within each fibre. While this approximation is acceptable for the removal of a moderate number of cross-links, it is likely to be inaccurate when relatively few cross-links have been deleted, that is for p close to unity. Similarly, for the very lowest values of p = 1, the requirement that only cross-links of coordination 4 are removed might invalidate the assumption that the probability distribution of segment lengths remains exponential. Nevertheless, the good agreement within the intermediate regime indicated a strong scaling with p, as predicted by our above theory.

Transition between affine and nonaffine deformation
As filament density and bending stiffness increases, bending deformations of network fibres become increasingly unfavourable and progressively more energy is stored in stretching modes. As such, a fibre network will move from the bending-dominated regime towards an affine, stretch-dominated one, as discussed earlier in this manuscript. We briefly recap the scaling of the network shear modulus within the affine regime, before investigating whether previous characterizations of the transition continue to describe the scaling in the presence of varying cross-link density. When subjected to an affine deformation, a constituent fibre initially forming an angle θ with the x-axis will rotate and stretch according to a (shear) strain γ to incur an energy cost:  is plotted against cross-link proportion parameter. For moderate values of p, we find that the theoretical scaling prediction p 6 ~is a reasonable fit, suggesting that our treatment captures the general behaviour of the system well. However, for larger values of p close to unity, or small p=1, the scaling diverges.
for stretching stiffness μ. Averaging across all (uniformly distributed) angles θ, and summing over all network fibres gives the scaling of the shear modulus for high polymer-density networks within the affine regime:  l ≔ ( ) . Using this length scale, the ratios of network shear moduli with the corresponding affine prediction affine   collapse onto a common master curve across varying bending rigidities and polymer densities. We wish to investigate whether this length scale still governs the transition in the presence of varying cross-link density. The attempted summary of our data using this length scale, for varying κ, p and ρ is shown in figure 4(a). As cross-link density varies, the data no longer collapse onto a single curve. Instead, considerable deviations of several orders of magnitude arise between networks possessing different bending length scales and densities for p<1.
Given the dependence of the shear modulus in the nonaffine regime on p, we seek a data collapse under a corrected length scale c l of the form λ c =λ p x for scaling exponent x and cross-link proportion parameter p=n/n 0 . For large density networks, we find that the exponent is given empirically by x≈1/2, and for a given polymer density this exponent appears constant. However, for lower density networks, the best collapse if given by an exponent that deviates from 1/2 by the inverse of polymer density ρ, that is x≈1/2−α, where α∼ρ −1 . Figures 4(b)-(e) shows the attempted collapse for networks of each of the polymer densities considered here under a corrected length scale p c 1 2 l l = a -, where α=8/ρ. We find that plotting the shear modulus normalized by the affine prediction against L / λ c gives an excellent data collapse within networks of a given polymer density, for varying κ and cross-link density. We suggest that the deviations in the exponent observed for lower density networks are due to the larger proportion of nodes with coordination 2 or 3, resulting in significantly different network topologies as p is reduced. In particular, the number of constraints in a low density network for, say, p=0.5 will be far lower than for a high density network with equivalent cross-link removal parameter p, as a large proportion of cross-links remaining in the latter have coordination 4. As such, we expect that the correction to the length scale should be greater in high density cases. Note that the singularity in the exponent is not concerning, as ρ ceases to have meaning in these simulations below the geometric percolation threshold. As the quantity λ c returns to the affine length λ for p=1, we may view the previously reported affine length scale for fully cross-linked Mikado networks as a particular case of this more general result. We now proceed to investigate the origins of the affine-nonaffine transition as cross-link and polymer density, along with fibre bending stiffness, are varied. The length scale λ, discussed earlier for the fully crosslinked Mikado model, partitions networks into those in which the energetics are governed by bending deformations, and those where stretching dominates. We wish to see whether the new length scale λ c given above still describes the transition between bending and stretching as the dominant modes of energy storage. The ratio of total bending to stretching energy contained within a network is plotted against L/λ c in figures 5(a)-(b). We find that, under the new length scale, dominant energy storage switches at the transition from bending to stretching, so that the shear modulus will scale with μ in the affine regime and κ in the bending dominated regime. This confirms that the proposed length scale describes the same physics as in the fully cross-linked case. Beam shear energy does not change results dramatically, though leads to small-order corrections for the highest polymer density networks considered. This can be seen in the deviation in bend to stretch ratio for high density networks, in figure 5(b), where for stretch-dominated networks, shear and bending energy contributions can become comparable. Further confirmation that the transition is governed by a switch between energy storage modes is given in figures 5(c)-(d), where the energies stored in fibre shear, stretch and bending deformations are plotted independently. Here, the small perturbations due to shear, reported in other studies [47], are made explicit, and we can see that while there is some variation in beam shear and bending energies for networks in the affine regime, the sum of these energies remains small in comparison to that stored in fibre stretching. Interestingly, the changes in cross-link density considered here appear sufficient to induce the transition between regimes. Given that many physical biopolymer networks exist near the transition regime [31], this suggests a mechanism through which the stiffness of networks with constant polymer density can be extensively regulated to accommodate their biological function, while raising the possibility that pathological changes in cross-link density could have important mechanical implications in vivo.
Having discussed the affine-nonaffine transition, and the dominant energy storage modes within phantom Mikado networks, we now investigate the microscale distribution of fibre energies. In particular, we wish to understand how the network energy is partitioned between the many constituent fibres. We begin by  of each network is partitioned into the three deformation modes; stretch, shear and bending. In all cases, shear energy forms only a small part of total energy, suggesting that the mechanical network response is largely independent of fibre shearing in this study.
investigating the strain energy densities of filament segments. In this case, we can produce an analytic affine prediction for the cumulative strain energy density with the proportion of fibre segments. Consider an affine network deformation, producing strain energy densities  per fibre segment of: sin cos , 10 for initial fibre orientation θ with the horizontal x-axis. From this expression, fibres parallel to the x-or y-axis store no energy under bulk simple shear to leading order in small strain γ, and energy content is symmetric in θ=π/4. Further, energy density increases until fibres are oriented at an angle θ=π / 4. The cumulative strain energy density with fibre orientation F(θ) should therefore take the shape of:  ò q q q q q q~-( ) As networks are isotropic, such that all orientations are equally likely, this expression, plotted for θ in the range (0, π/4), and normalized to the maximum value F 4 q = p ( ) , is equivalent to plotting the affine cumulative strain energy density prediction against the proportion of fibres with those values. We compare this expression with the equivalent data from phantom Mikado networks in figures 6(a)-(c), with F(θ) plotted in dash. For the most affine networks considered (yellow in figure 6(a)), the simulation data are closer to the prediction, though energy is less equally distributed in the simulated networks than in the affine prediction. As the cross-link density is allowed to vary, we observe very little difference between the cumulative energy density curves in affine networks ( figure 6(b), markers) until p becomes very small, at which point a sudden switch to a nonaffine profile occurs. In nonaffine networks, a similar insensitivity is observed ( figure 6(c), markers), until almost all eligible cross-links are removed, whence there is a slight move towards a more affine profile. Taken together with our results for the network shear moduli, this suggests that the distribution of strains throughout the sample do not significantly Figure 6. (a)-(c) Cumulative fibre segment strain energy density is plotted, with comparison to the analytic affine prediction (dash line) given in equation (10). (a) Example results are displayed for affine (yellow) and highly nonaffine (blue) networks in the standard Mikado case, p=1. Affine networks more closely follow the theoretical prediction, while nonaffine networks are characterized by the presence of a small family of energetic fibre segments. The same result is plotted for (b) affine and (c) nonaffine networks as cross-links are removed. In affine networks, the strain energy density distribution is almost constant as cross-links are deleted, until the lowest values of p are reached, at which point a sudden move to the nonaffine profile is observed. The strain energy density distribution in nonaffine networks appears similarly insensitive to cross-link parameter p, though surprisingly a modest move towards the affine energy distribution is observed for the most sparsely cross-linked nonaffine networks. (d)-(f) Cumulative fibre energy content is plotted against proportion of fibres, with parameters inherited from (a)-(c), showing that similar results hold for fibre energy content as for segment energy density. Inset shows the energy fraction of the top 10% of fibres by energy content, highlighting the fact that a small family of fibres contributes a majority of total network energy. Markers are included to distinguish near identical data arising from networks within a given regime. change, even though the bulk stiffness can change by orders of magnitude. We expect that the change in profile for the initially affine networks as cross-link density becomes small is due to a cross-link-induced affinenonaffine transition, akin to that observed in figure 2(b) (Inset).
While we have quantified the fibre segment strain energy density distributions, it is not clear how the actual energy content in the constituent network fibres is distributed. In figure 6(d)-(f) we plot the cumulative fibre energy against the proportion of fibres, for the same network parameters used for the strain energy density distributions. That is, we order network fibres according to their energy content, and use this information to produce the cumulative energy, highlighting the energy share of different fibre families. A hallmark of nonaffine networks is the presence of an increasingly small family of fibres containing the majority of total network energy, while those networks that deform affinely, according to the length scale λ c discussed above, distribute energy across a greater number of fibres. Inset gives a closer look at the energy fraction of the top 10% of fibres by energy content, highlighting that the majority of network energy is contained in relatively few fibres. For initially affine networks (figure 6(e)), over half the total network energy is stored in less than 7.5% of fibres. This is not changed by large reductions in the network cross-link number until the cross-link parameter becomes small, p=0.07, whence 50% of total energy is contained in just ≈2.5% of fibres. The small move back towards a more affine profile as cross-link density becomes low in nonaffine networks is repeated here in figure 6(f). While approximately 1% of fibres contain the majority of network energy, the removal of over 80% of cross-links more equally distributed fibre energy throughout the network. We suggest that this change in energy profile is due to the changing proportion of coordination 2 and 3 nodes as cross-link density decreases.
In the above discussions, network deformations have been characterized through energy criteria. However, it is often more natural to characterize affinity in terms of the displacement and rotation of constituent fibres and cross-links. Such results may be harder to verify experimentally, requiring the analysis of tracer beads displaced within the biopolymer network, rather than measurements of the bulk shear modulus. However, it is useful to have such a description where microscale displacements are important, for example when investigating displacement fields derived from cell tractions. It is also important to question whether affinity, as characterized by displacement, coincides with energy-based measures such as those discussed above. Here, we investigate the displacement field in phantom Mikado networks by comparing the mean deviation of cross-link displacement from the corresponding affine predictions. In particular: where summation is over the number of cross-links n x , is the deformed position of a model cross-link, and u x X is the displacement of the i th cross-link, at initial position X i . We plot this measure, scaled by domain size W and strain γ, against the proposed length scale for affinity λ c in figure 7. While we do not find a data collapse for this measure, similar to other researchers [31], it is clear that a displacement-based affinenonaffine transition does occur. This transition is evidenced by the move between two plateau values for displacement measure D, large for nonaffine networks and small for affine cases. At a similar value of λ c observed for the energetic affine-nonaffine transition, there is a rapid switch between these plateau values. The removal of cross-links (lower values of p, indicated by colour) generally leads to larger nonaffine displacements throughout the sample, in all of the densities considered. We expect this to be relevant when considering the displacements experienced by cells resident in substrates with variable cross-linking [9].

Discussion
The mechanics of fibre networks has been studied under a variety of modelling assumptions, from lattice and disordered fibre network simulations, to mean-field theories. The Mikado model, on which this work is based, has provided key insights into the mechanics of biopolymer networks, where the fibre bending stiffness, network structure and polymer density leads to non-trivial elastic regimes. However, within this model all overlapping fibres are bound together at a permanent cross-link. Given recent improvements in electrospinning and photocross-linking technologies, in vitro platforms now have far more control over the biopolymer network structure, including the geometry and topology. Furthermore, many planar biopolymer networks involve fibres that overlap yet do not rigidly cross-link. To address the need for a more versatile approach to cross-link modelling, while recognising the importance and influence of the Mikado model within the discrete fibre network field, we extended the standard model to incorporate networks where not all fibre crossings result in the formation of a permanent cross-link.
We have shown in figure 2 that the scaling of the shear modulus with cross-link removal within the nonaffine regime, away from the affine transition, is qualitatively similar across different densities and bending stiffnesses. In particular, we found that the shear modulus of nonaffine, sparsely cross-linked phantom Mikado networks responded sensitively to the addition or removal of relatively few cross-links. For affine networks, we predicted and confirmed that, owing to the insensitivity of the affine model to the detailed nature or density of cross-links, the shear modulus should be unchanged as cross-links are removed. At the very lowest cross-link densities considered, even initially highly affine networks entered into the nonaffine regime, whence the shear modulus reduced sharply. We further produced a theoretical prediction, based on a fibre network floppy mode argument [53], that the scaling of the shear modulus should follow as p 6 ~within the nonaffine regime. This scaling found good agreement with simulation data for moderate values of p, yet diverged for the lowest and highest values considered. We expect that the cross-link removal criterion, that only cross-links of coordination 4 may be removed, leads to different fibre segment length distributions at the smallest values of p considered here. The differences observed for p close to unity likely arise due to a break down of the theoretical argument. In particular, when relatively very few cross-links have been removed, we do not expect each network fibre to experience a similar reduction of cross-links, such that the averaging argument employed no longer holds. We note that the strong sensitivity to p within the nonaffine regime is likely to have implications for 3D biopolymer networks, where the ratio of polymer to cross-link density is far higher than in planar networks.
It was not clear whether the bend-stretch transition in Mikado networks with varying cross-link density occurs under the same length scale as for the maximally cross-linked cases, namely l l l c c b 1 3 l = ( ) . However, we found that this length scale did not produce a data collapse for the shear moduli of phantom Mikado networks, as shown in figure 4(a). By considering an additional term relating to the cross-link parameter p, we proposed a corrected affine length scale p c 1 2 l l = a -, through which we found an excellent data collapse within a polymer density, as shown in figure 4(b)-(e). Within the simulations considered here, α depended inversely on density, α=8/ρ. For dense networks, this affine length scale becomes independent of ρ, though for low density assemblies, we observe sensitive differences according to polymer density. We suggested that this is likely due to the network differences that arise from different proportions of coordination 2 and 3 nodes, which are not removed under our model assumption. We note that for fully cross-linked networks p=1, this length scale reduces to the original quantity λ.
To ensure that the data collapse for phantom Mikado networks was still due to the dominant storage of energy within fibres transitioning from bending to stretching modes, we further investigated the energy partitions in network fibres. Under the new length scale, networks clearly displayed a transition similar to that observed in maximally cross-linked cases, where fibres became progressively less compliant in bending, and stored energy in stretching instead, as seen in figure 5. The ratio of shear energy to the dominant energy mode in all cases remained small. We continued the investigation into the constituent fibre behaviours by quantifying the cumulative segment strain energy density and fibre energy content in figure 6, compared with a theoretical affine prediction. In fully cross-linked Mikado cases, both segment energy density and fibre energy content were more equally distributed throughout the population than in nonaffine cases, where a small number of fibres contained the majority of total network energy. Generally, we found that the observed distributions were robust to changes in cross-link density, despite the large changes in network stiffness that we observed when cross-links were removed from, for example, nonaffine networks. The largest changes observed occurred when large numbers of cross-links were deleted in initially affine Mikado networks, inducing the aforementioned affine-nonaffine transition, and leading to a likewise change in fibre energy distributions. The fact that such a large proportion of the total network energy is stored in so few fibres in nonaffine cases suggests that such assemblies might be sensitive to damage, with the scission of a small family of fibres potentially altering the macroscale response. However, we note that there is potential for network redundancy in these situations, and that damage to energetic fibres could lead to deformations in nearby fibres instead. The fibre energy distribution could also be important when investigating the mechanical response of cells embedded in fibrous matrices, where the presence of relatively few energetic fibres could dictate different local responses from mechanically sensitive cells. Increasing the cross-link density could potentially protect against this effect, by more equally distributing energy across a larger group of fibres, suggesting another utility for varying cross-link density.
By measuring the mean deviations of network cross-links from the predicted affine displacement in figure 7, we investigated another measure for affinity that did not require any knowledge of fibre energetics. Displacement measures for affinity appear less well characterized in the literature than those derived from energy-based criteria [54]. However, our results for mean deviations from the affine deformation field suggest that affinity as measured using displacement follows the same proposed length scale λ c . In particular, there still exists a sharp transition between two regimes, where nonaffine displacement fluctuations are larger in energetically nonaffine networks. While we do not assign a definite quantitative interpretation to the values we observed, we do note that considerable differences in the displacement measure are observed as cross-link density varies, especially in those networks with high polymer densities. Here, larger nonaffine fluctuations arise as cross-links are removed, suggesting that displacement nonaffinity is sensitive to cross-link density.
We have not reported on results for welded networks, where cross-links prevent the rotation of incident fibres, as preliminary simulations suggested that only minor differences arise. In particular, while networks are consistently slightly stiffer owing to the increased number of constraints, the transition between regimes occurred under the same length scale, and the qualitative behaviour was similar. Quantitative studies investigating the relevance of cross-link type and cross-link density in comparison with experimental data are topics for future work.
The results presented here have considered only planar simulations, which have the advantage of being computationally cheaper and have been frequently utilized in the literature. As such, 3D Mikado networks have not been discussed. However, we note that the same procedure presented in this manuscript can also be extended to 3D Mikado models, where network geometry can be preserved through the removal of only those cross-links connecting nodes with coordination 3 (after dangling ends have been removed). This would allow for an investigation of the importance of 3D network geometry as decoupled from cross-link density. Our predictions from nonaffine phantom Mikado networks, namely that in this regime biopolymer assemblies become sensitive to small changes in cross-link density, are likely to be relevant to 3D network mechanics. Indeed, the phantom networks discussed here possess a far greater ratio of polymer to cross-link density, which might better represent 3D systems than the maximally cross-linked standard Mikado case. As 3D Mikado networks are considerably less constrained than planar cases, we also expect that many such networks are likely to reside closer to the nonaffine regime, such that cross-link density could play a vital role in regulating the macroscopic network response.
The Mikado model has been crucial in expanding our theoretical understanding of fibre assemblies, in part due to the natural representation of the biopolymer network. However, this work demonstrates that the condition of maximal cross-linking does not give a complete understanding of the network behaviour. Given the results presented here for Mikado geometries, it would seem important to study the effects of varying crosslinking in other geometries, to investigate whether a similar sensitivity is observed there, especially within the nonaffine regime. Our results provide further evidence that the mechanical response of fibrous networks is not insensitive to their detailed structure, even where isotropy is maintained, and suggest an important experimental target for photocross-linked, electrospun polymer networks. How variations in structure and cross-link density regulate the response of in vivo biopolymer networks remains a fascinating and important question.
In summary, we have investigated the mechanics of random networks of elastic rods in which cross-link density is varied. We observed that the scaling of the network shear modulus with cross-link removal in the nonaffine, bending dominated regime was the same across a range of bending rigidities and polymer densities, and provided theoretical predictions for the scaling of the shear modulus in this regime. We investigated previously characterized length scales for affinity in Mikado networks and demonstrated that they are not suitable for summarizing data from phantom Mikado architectures. We suggested an empirical length scale that characterizes the affine-nonaffine transition for a given polymer density as the extent of cross-linking is varied. On the scale of individual fibres, we confirmed that this transition occurred due to a bend-stretch change in energy storage, such that the corrected length scale described the same physics as in the maximally cross-linked cases. We also highlighted that in affinely deforming networks, a larger proportion of constituent fibres contribute to the total network energy, and that the distribution of fibre energies was largely robust to variations in network cross-linking. Finally, we confirmed that the affine-nonaffine transition occurred under the same length scale when characterized through cross-link displacements, with the magnitude of nonaffine fluctuations exhibiting a drastic increase as cross-link density was reduced.