A minimal model of elastic instabilities in biological filament bundles

We present a model of a system of elastic fibres which exhibits complex, coupled, nonlinear deformations via a connecting elastic spring network. This model can capture physically observed deformations such as global buckling, pinching and internal collapse. We explore the transitions between these deformation modes numerically, using an energy minimization approach, highlighting how supported environments, or stiff outer sheath structures, favour internal structural collapse over global deformation. We then derive a novel analytic buckling criterion for the internal collapse of the system, a mode of structural collapse pertinent in many biological filament bundles such as the optic nerve bundle and microtubule bundles involved in cell abscission.


Introduction
Filamentary bundles (fibres embedded in an elastic matrix) are a common biological structure. There is ample experimental evidence showing how their elastic response is critical for their functionality. For example, collagen bundles are composed of linked tropocollagen molecules with a range of mechanical properties from compliant to rigid [1][2][3]. The optic nerve bundle is composed of rigid bundles of myelinated neurons embedded in connective tissues, and damage to this bundle because of pressure-induced buckling is considered a likely cause of glaucoma [4]. Similar nerve bundle structures play a crucial role in spinal cord injuries [5]. Microtubules are known to moderate their flexibility and alter their load bearing capacity in muscles [6,7].
In this study, we develop a model that explicitly describes multiple filaments whose interactions are mediated through an elastic matrix, which we will show to be essential for capturing the surprising range of deformation modes exhibited by fibre bundles under compression. To our knowledge, there exists no fibre bundle model which has been demonstrated to be capable of describing these modes. In figure 1, we categorize these deformation modes into three types. Firstly, in figure 1a, the bundle may buckle globally, as observed for isolated actin and microtubule bundles [8,27]. Secondly, as shown in figure 1b, pinching (lateral contraction) may occur. This has been seen in the mouse embryonic fibroblast (MEF) [28], where the contraction of actin fibres was shown to be crucial in controlling nuclear shape. Finally, in figure 1c, buckling may instead occur internally, as is the case for microtubule bundles and the optic nerve bundle [4,29,31].
We are able to rationalize these deformation modes with three key parameters: the axial compression relative to the filament bending stiffness, the pinching load normalized against the interior filament bending stiffness and the non-dimensionalized ratio of the bending stiffness of the outer filaments to the inner filaments (representing, for example, the optic nerve bundle's stiff outer sheath).
The principal feature of our model is the explicit coupling between filaments of fixed length via both Hookean springs and an overlap penalty. This combination of interaction forces is essential for two reasons. The first is so that the model captures physical deformations. The second is so that the buckling of one filament can potentially trigger further buckling, and a nonlinear cascade of response for the whole bundle. This coherent buckling is critical, in particular, for explaining the internal collapse mode.
The internal buckling mode is of particular interest here, in part because there appears to be no relevant buckling criteria existing in the literature, but mainly because of their functional role in biological bundles. As our model will confirm, the localized pile-up type behaviour of these deformations leads to significant spikes in local pressure. For the optic nerve, this could be crucial and deleterious, as external pressure can inhibit neuronal signalling action [32]. By contrast, for cell abscission this pressure is a desirable property as it locally weakens the bundle, providing a point for cleavage [29]. The pinching deformation modes are shown to lead to a more uniformly distributed internal pressure which is of a lower order of magnitude, in some sense making it a 'safer' deformation. It is thus important to know what factors will favour one mode over another. Overall, a key outcome of our analyses is to show the mechanical forces relevant for mediating the transitions between different bundle morphologies.

The model 2.1. Model geometry
The basic model, illustrated in figure 2, comprises a set of m inextensible planar elastic rods r i , i ∈ 1, …m. Each rod, of the same arc length L, is represented as a discrete curve which is itself composed of n points (typically n = 200 here). Their Cartesian coordinates ðx ij , y ij Þ, j [ 1, . . .n, can be parametrized in terms of a set of angles θ ij The parameter W is the bundle's width. Thus u ij ¼ 0, 8i, j represents a set of straight filaments directed along thex-axis: the bundle's rest state. The m discrete rods have fixed positions at j = 0 for which x i0 = 0 and there is a vertical spacing y (i+1)0 − y i0 = W/(m − 1), as shown in figure 2a. The angles θ ij determine the bending of the curves, as shown in figure 2b. In this work, the filaments are confined to two dimensions, as already a range of deformation modes are manifest, with their biological relevance highlighted in figure 1. More complex threedimensional modes, such as torsional deformations, will be explored in future.

Energy functional
The set of filaments is then assigned the following energy functional: where the total energy E is the sum of the bending energy E b , elastic interaction potential E e , loading potential E l and geometric boundary constraint potentials E c . All parameters used to subsequently define the energy functional components are summarized in appendix A.

Bending energy E b
For each rod i, the bending energy takes the following form: This curvature measure is based on the radius of a circumscribed circle encompassing adjacent edges of the curve,  Figure 1. Three general classes of deformation modes: (a) a globally buckled state, here compared with the buckling of actin filaments [27] (adapted with permission, please be aware that this permission may not cover reuse under the OA agreement); (b) a pinched configuration, here compared with the basal organization of the actin filament network in an MEF fibroblast [28] (adapted with permission, please be aware that this permission may not cover reuse under the OA agreement); (c) internally buckled states, compared here with (left) microtubule bundle during cell abscission [29] (reproduced/adapted with permission which does not not cover reuse under the OA agreement) and (right) 'cupping' of the optic nerve in an eye that has suffered glaucoma [30] (adapted with permission, please be aware that this permission does not cover reuse under the OA agreement). and is a standard discrete form [33]. In the limit l → 0, it becomes a finite difference approximation for the curvature. The coefficients B ij represent the localized bending stiffnesses of the filament. Here, these values will be uniform along the filament length (B ij = B i ) and the only variation is between the outer bending stiffness B o = B 1 = B m and inner bending stiffness B in = B i , i = 2, 3…m − 1. Physically the ratio B r = B o /B in is used to model either surrounding support to the bundle or a stiff outer sheath, such as the epineurium and perineurium sheaths surrounding neuron fibre bundles [34].

Interaction energy E e
The interaction energy E e in (2.2) is the sum The component E 1 e is a set of Hookean spring potentials joining neighbouring filaments, dependent on the arc length paired distances d ij between neighbouring springs (i, i + 1) at location j along their length (figure 2b). Formally, E 1 e is the following Hookean potential: where d ij = |d ij | , with the vector d ij expressed as and the function Here d 0 = 1/(m − 1) is the value of d ij in the bundle's rest configuration (all θ ij = 0) and the constants K ij are the spring stiffnesses. The model can allow K ij to vary but in this study we choose a homogeneous connective tissue K ij = K. The factor l in (2.5) makes K spring energy a density per unit length. In addition, there are diagonal spring potentials E 1þ=À e dependent on distances d + ij between point pairs (i, j ) and (i, j ± 1), as indicated in figure 2b. The first set of potentials, corresponding to the distances d ij , model the stretching of connective tissue linking the filaments, while the potentials associated with the diagonal distances d + ij model the localized axial stretching of this tissue due to filament curvature. The diagonal spring potentials E 1+ e are also Hookean potentials ð2:7Þ is the value of d + ij in the bundle's rest configuration. Note there is no d À i0 or d þ in as these distances are formally set to zero.
The second interaction potential E 2 e prevents the overlap of the filaments by penalizing each area Aðd ij , d þ ij Þ and Aðd iðjþ1Þ , d þ ij Þ, illustrated in figure 2b, if an area gets close to zero. Aðd ij , d þ ij Þ is defined by the triangle between the two vectors d ij and d þ ij , such that ij Þ becoming zero, as for this to occur the length d ij or d þ ij is zero, or the angle between d ij and d þ ij is π. For Aðd iðjþ1Þ , d þ ij Þ, similar reasoning applies, with The term g(A) we choose to penalize area A getting close to zero is where A 0 = l/2(m − 1) is the size of each triangle in the rest configuration. The parameter C 1 is chosen to be C 1 = 10 6 in order to effectively prevent rod overlap. The parameter C 2 determines the extent of compression required for the potential to become large. In most of the examples in this study, we (i + 1, j royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 19: 20220287 use C 2 = 100, so that g(A) is only significant if the area of a given triangle is reduced to less than 10% of its original size. This means the elastic matrix between the filaments acts like a linear material until the filaments come close to contact. A lower C 2 value would mean the potential becomes larger for less significant compression (the value of the ratio A/A 0 ) making the material act more like one with pronounced strain hardening. In §3.6, we briefly explore the consequences of lowering C 2 and show this does not qualitatively alter our conclusions as to the nature of the system's phase space. The specific form we choose for E 2 e is We experimented with a similar overlap energy penalty, in which a repulsive potential was applied if d ij and d + ij approached zero length. However, this was not as effective as the area penalty, as points on the filament could fit between two points on a neighbouring filament without the spring distances ever becoming small. A final point to make is that this energy does not prevent non-local selfoverlap, occurring in globally buckled states when the load N c ≫ 1 and the bending coefficient ratio B r is of order unity. This can be remedied by adding an additional potential which penalizes non-local filament overlap, such as penalizing a point from approaching a line segment between two other points. However, such a term was found to significantly increase the computation time, and the range of forces used in this work mean that the local treatment was sufficient to prevent filament overlap. This allowed us to complete a significant range of parameter sweeps in a reasonable time frame.

Loading potentials E l
The loading energy E l in (2.2) is the sum of two terms: The potential E 1 l represents a compressive load ÀN cix applied to each rod i with the potential energy associated with the load taking the following form: In this work, we apply a uniform load N ci = N c for all i.
The potential E 2 l applies lateral compression force, perpendicular to the force N c (along the y-direction). This is defined as In this study, we prescribe a uniform compressive force on the outer filaments, so that and N p(ij) = 0 otherwise. We call this a pinching load N p in what follows, as it produces pinched geometries. The N p(ij) can also be readily chosen to model different systems. For example, pinching at only one end is relevant for optic nerve bundle shown in figure 1c where the scleral wall pinches the bundle. For other applications such as cell nuclei models, one could pinch at the middle to simulate the pinching of a nuclear cell.

Geometric boundary potentials E c
The geometric potentials impose boundary conditions such as pinned conditions (i.e. no end-to-end deflection), and clamped boundary conditions (i.e. no end-to-end deflection and fixed end tangents). The pinned case for which y i0 = y in is Hðy in À y i0 Þ, ð2:13Þ where Hð0Þ ¼ 0, HðsÞ . 0, 8s = 0. We use quadratic functions in the form H(s) = C 3 s 2 with C 3 a constant. Such a harmonic energy penalty is well established for imposing constraints for a wide range of systems [35,36]. In this study, a value of C 3 = 10 5 was used to ensure the conditions were strongly enforced. Clamped boundary conditions include the potential E 1 c and an additional potential E 2 c which constrains the filament ends to be parallel to the x-axis. It takes the form Hðu in Þ: ð2:14Þ The same quadratic form and C 3 value was used for this condition also. With regard to the biological examples shown in figure 2, the pinned boundary conditions applied alone would be appropriate for the actin filament buckling (a), as buckling is induced by compressing the filaments between two microscope slides. Indeed the single filament model in the cited study [27] used pinned boundary conditions. By contrast, the optic nerve bundle (c) begins at the lamina cribrosa, which is an elastically stiff porous material through which the neuron fibres are threaded [4]. This implies the curvature would be limited in this region, so clamped boundary conditions are appropriate. In addition, the pinching of the scleral wall restricts end movement of the filaments, so the addition of pinned conditions would also be appropriate. Similarly in the cell abscission example [29], there is an additional stiff elastic matrix structure (the mid-body structure) which deforms during abscission, which would also probably restrict curvature. We use both pinned and pinned and clamped conditions in what follows.

Non-dimensionalization
We non-dimensionlize the system parameters with a reference length corresponding to the width of the bundle W, and a reference bending stiffness B s = 0.01. This yields non-dimensionalized . The length scaling is chosen to ensure the parameter m determines the density of packing of the filaments through the rest separation d 0 = 1/(m + 1). Hereafter, we drop the primes and all quoted parameters are assumed non-dimensionalized unless otherwise stated.

Locating energy minima
To find stable states, E in equation (2.2) was minimized using the L-BFGS algorithm [37,38], a memory and computationally efficient minimization algorithm that we have shown previously to be suitable for high-dimensional buckling models [39,40]. The explicit form of the gradients ∂E/∂θ ij used to construct the gradient and approximate Hessian matrix are given in appendix B.

External and material parameters
The triplet (N c , N p , B r ) forms a set of externally acting parameters. This is clear for the loads N c , N p . The ratio B r is classed as external as it is used to model either a surrounding support to the bundle or a stiff outer sheath as discussed above. The parameters K (connective medium spring stiffness), number of rods m, filament length L and internal filament bending rigidity B in are the material parameters. In this work, the set (N c , N p , B r ) will parametrize the phase space of configurations and we find that they govern the transitions between the differing bundle geometries which minimize the system's energy. We will then see how the material parameter set (K, m, B in , L) alter the nature of this phase space, in particular how they determine the critical values of these parameters at which the transitions occur.

Deformation classification
Overall, four stable geometry classes are found, which are visualized in figure 3. These classes are the unbuckled system shown in figure 3a, and the three deformed systems (globally buckled, pinched and internally buckled), shown in figure 3b-d, respectively. We highlight the fact the globally and internally buckled configurations are 'symmetrically' buckled, in the sense they buckle with the same mode and direction. By contrast, in the pinched configurations, the upper and lower filaments buckle in opposing directions. The configurations shown in figure 3 are archetypal minima of the system, which occur either in the absence of pinching (globally and internally buckled configurations) or the compressive load (asymmetrically pinched case). In general, when both loads are present, we use these archetypal configurations to define three characteristics of each geometry, whose relative degree is used to categorize the configurations. They are: (i) B r 〈|κ|〉: the mean absolute moment, an average of the curvature |tan 2 ((θ ij − θ i( j−1) )/2)| along the filament length, then averaged over all filaments. We multiply by B r to account for the fact deformation is less for a given load in the outer filament. (ii) B r 〈|κ|〉 o and B r 〈|κ|〉 i , the quantity in (i) restricted, respectively, to the outer filaments and inner filaments only. (iii) P: the pinching parameter. The pinched configurations have filaments which buckle in the opposite directions either side of the centre, whilst for the internal/global buckling the filaments buckle in the same direction.
To measure the degree of pinching, we split the rods into the upper half i > (m + 1)/2 and the lower half i < (m + 1)/2 (or split around m/2 if m is odd). Then we calculate the weighted difference in angle between paired filaments from the upper and lower rods:  royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 19: 20220287 The mean absolute moment B r 〈|κ|〉 is simply used to discriminate unbuckled and buckled configurations. The globally buckled configurations have particularly large values of this parameter. If B r 〈|κ|〉 is greater than 0.00001 then the configuration is classed as buckled. We must then identify what type of buckling is present. A large value of the ratio B r 〈|κ|〉 i /B r 〈|κ|〉 o indicates structures with significant internal buckling characteristics (we use a ratio greater than 5). Otherwise, we must establish whether the structure is (globally) buckled or pinched, which is done using the pinching parameter P.
Generally we find the pinching parameter P is either very small or very close to 1, with small P indicative of systems with highly symmetric deformations such as global or internal buckling. We do find examples where there is a competition in classification due to the presence of both symmetric and asymmetric deformation modes. For example, globally buckled configurations with significant pinching forces are shown in figure 4a,b. The pinching load clearly affects the configuration, compared with figure 3b, but the dominant global characteristic is clear. Similarly in figure 4c, we see a configuration which is internally buckled but has some significant pinching. With a slight increase in pinching, figure 4c transforms to figure 4d and becomes classed as asymmetrically pinched. Note that significantly mixed pinched/internally buckled behaviour is only observed for mid-range B r values (we consider [1, 10 000]) in our study. For higher B r (say greater than 500) the symmetric and internally buckled configurations are significantly pronounced.

Surveying the phase space of deformation
We survey the available deformation modes over the parameter ranges N c ∈ [0. 5,10], N p ∈ [0, 0.4] and B r ∈ [1, 10 000] for the phase diagram in figure 3e. In figure 3e, we have also used L = 5, K = 2, m = 10, B = 1. This range incorporates biologically important example bundles such as the optic nerve, as we will demonstrate in §4. The spring constant K is non-dimensionalized with respect to the bending coefficient. This value has the same order of magnitude as the bending stiffness so that both effects are important to the system's behaviour (we vary this value shortly). We initialize the energy minimization in each of the four archetypal states, and then locate the local energy minima. When there are two local minima with energies within 1% of each other, these are marked with two colours.
We now explore the stability ranges and transitions in the phase diagram. We begin at (N c , N p , B r ) = (0, 0, 1) in the undeformed geometry. This is a bundle with no external protective support (B r = 1). As the pinching load N p is increased, the magnitude of the pinching deformation shown in figure 3c simply increases smoothly. There is no apparent instability. By contrast, if N c is increased instead, a sudden transition to the globally buckled state occurs, shown in figure 3b, governed by the classical Euler buckling criterion. Next, upon increasing the external support parameter B r , the globally buckled configuration is suppressed due to the extreme energetic penalty associated with deforming the stiff outer rods. At large B r > 100, as N c is increased, both the undeformed and pinched states suddenly collapse to an internally buckled state, such as that shown in figure 3d, at some critical value of N c , denoted by N Ã c . There is a tendency for domination of the pinching configurations when N p = 0 at lower N c . This is not a surprise, as the continuum version of elastic filament models does not have a straight equilibrium θ = 0 in the presence of a lateral force (see [41] eqn (99)). So the presence of a load N p should always lead to some deformation. As discussed above, we find the degree of pinching to grow steadily with N p (the example shown in figure 3c is one of the more exaggerated examples). For the high B r cases, the degree of pinching is severely restricted.

Additional examples
We focus on the same phase space, the parameter ranges  royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 19: 20220287 in figure 5a. This differs from the initial case which also has clamped constraints. For B r = 1, the globally buckled configuration is the only configuration present. This contrasts with the clamped case shown in figure 3e, where the Euler buckling criteria ensure that global buckling is suppressed until the Euler criterion N c > 4π 2 /25 is satisfied (as B = 1). For pinned boundary conditions, the critical criteria is π 2 /25 ≈ 0.395 which is always met within our N c range. A second contrasting feature to the clamped results is that, for B r = 1-33 (the first three layers of the plot), the globally buckled configuration dominates for lower N c . If we compare, for example, the results for B r = 33, there is little pinching above N c > 3, while pinching is observed up to loads of N c = 7 in the clamped case. At B r = 100 (fourth layer up) we also note that, in the pinned case, there are globally buckled configurations for high N c , whereas there are no such globally buckled minima in the clamped case. Hence the extra boundary constraint is important in this regime. By

Twenty rods
The minimum energy configurations for m = 20, L = 5, K = 2 and B = 1 with both pinned and clamped constraints are shown in figure 5b. These differ from the results shown in figure 3e only by the number of rods. There is significant similarity between the two results so there are only a couple of features to highlight. First, we note that for intermediate values of B r the pinching configurations disappear for lower N c , indicating a slight preference towards buckling when the outer support is not too strong. Second, for higher B r , while the N c values at which internal buckling becomes a minimal energy configuration are the same, there is a range of N c values over which both the pinching and global buckling configurations are effectively bistable. This is not present in the m = 10 case.

Varying matrix (spring) stiffness
The minimum energy configurations for m = 10, L = 5, K = 4 and B = 1 with both pinned and clamped constraints are shown in figure 5c. This has twice the spring stiffness of the results shown in figure 3e. The preferred configurations for low B r are nearly identical in both cases and the only significant difference manifests for high B r , where the critical N c value at which internal buckling becomes an energy minimum is significantly increased. The variation in the curvature of the filaments inherent to this mode will naturally lead to significant stretching of the spring matrix, which explains why the parameter K has a significant effect on the critical internal buckling condition. Further evidence of this is seen in figure 5d, which corresponds to the minimum energy configurations for m = 10, L = 5, K = 1 with both pinned and clamped constraints. The only significant difference is the lower N c value at which internal buckling occurs for high B r .

Forces in the system
One final aspect of these minima which we investigate is the internal pressure caused by the spring forces, which can be computed as the partial derivatives of the interaction energy E e with respect to the coordinates (x i , y i ). The magnitude of their distributions are shown in figure 6a-c. There are clear spikes where the filaments are pressed together. The peak pressure was found to be an order of magnitude higher for the buckled configurations compared with the pinched case, indicating such geometries would be particularly prone to localized damage for individual filaments. As discussed in the Introduction, this can be either a positive or negative feature of the bundle's mechanical response. The more uniform deformation of the pinch leads to more benign interactions.

Biological relevance
To further illustrate the potential biological relevance of these results, we argue that the material parameters we consider are relevant for the optic nerve. The optic nerve has a width and length of 3.55 and 28 mm, respectively [42], leading to the scaled length of L = 7.88 in our model. Young's modulus for the outer sheath has been estimated to be E o ≈ 44.6 MPa and the optic nerve itself E i ≈ 5.4 MPa, respectively [42]. These values can be translated to filament bending coefficients via πER 4 /2, assuming a circular cross section with radius R. For the outer sheath, R o = 0.37 mm [42]. The optic nerve thickness is typically around 3.55 mm, covering about 15 neuron sub-bundles, which leads to a nerve's width of R i = 3.55/30 = 0.12 mm. Thus, B r = B o /B in ≈ 10 4 . Finally, the typical inter-ocular pressures are 6.75 − 33.75 × 10 3 MPa, which gives a range of N ∈ [6,35] in our model. These parameters put the bundle in the regime where internal buckling dominates, suggesting optic nerve buckling in glaucoma can occur depending on the density of the collagen matrix which fills the space between the optic nerve bundles.

Elastic matrix model extensions
As discussed in §2.4, lowering the value of C 2 makes the inter-fibre matrix act more like a strain hardening material. To explore the influence of this, we compile energy minima for the same parameter set as in figure 3e, but with the lower C 2 value of C 2 = 40. This means the potential E 2 c becomes large when the matrix material is compressed below approximately 30% of its initial area. The results shown in figure 7 show only modest differences in comparison with figure 3e: at B r = 100 there is an increase in the dominance of pinching, and the force required for internal buckling increases slightly. The similarity between all phase diagrams presented throughout suggest that our qualitative conclusions generalize across different elastic matrix models. Furthermore, the energy approach to our model means that is straightforward to replace the Hookean potential in equation (2.5) with more complex models of biological elasticity, if specific systems are to be investigated in the future.

Internal buckling
As discussed in the Introduction, the internal buckling transition is of particular biological significance. An important royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 19: 20220287 example of this is shown in figure 1c, in which asymmetric cupping of the optic nerve is observed in glaucoma sufferers. This is not readily explained by linear elastic models [22,43]. The extreme localized pressures in the internally buckled state might explain the damage caused to single elements of the whole bundle, thus leading to localized sight loss, the defining property of glaucoma's insult.

The internal buckling criterion
As observed in the previous section, the model suggests the critical load at which this internal deformation occurs is consistent for a sufficiently strong external support parameter B r > 100, independent of the pinching force N p , and varies with the model parameters. This suggests a clear buckling criterion could be developed for such systems which depends on the material parameters K, B in , L, m. To characterize this criterion, we chart the critical load N Ã c at which the straight configuration collapses into the internally buckled state in figure 8a, for a variety of spring stiffnesses K, filament lengths L, and bending stiffnesses B. We fix m = 10, N p = 0 and B r = 10 000. The solid lines represent the following analytic prediction for the critical load N Ã c : where n e is minimum Euler buckling mode, n e = 1 for pinned boundary conditions and n e = 2 for clamped conditions, and α = 10.61 (we will shortly discuss the relevance of this value). Given the complexity of the model and the deformations involved, this is a remarkably simple formula and, as indicated in figure 8a, it is also highly accurate. We now discuss its derivation in some detail as it provides significant insight into the nature of the instability and the factors which determine it.

A discrete model
We make an energy estimation of the internally buckled state by assuming a triangular shape for the buckled filaments, as indicated in figure 8b, with the amplitude of buckling decaying laterally to zero, as in the numerically obtained minima shown in figure 8c. The indicated angle θ t dictates both the buckling frequency and, through the fixed length criteria, the amplitude of buckling D. Both are functions of the distance Δ through which the force N c does work. This allows us to explicitly calculate the energy in equation (2.2) and minimize it via a two-step process: first we find the load at which the system becomes unstable when Δ = 0, then we minimize the energy at this load for the wavemode n at which this occurs. We begin by only assuming m = 3 so that only the middle filament actually deforms. The general m formula is then a straightforward extension of this.
With reference to the undeformed state, the energy of the simple internally buckled state is ð4:2Þ Here E l is the work done by the load N e bending the middle filament. All parameters used to define the components of E int are summarized in appendix A. Note that in this simplified model we are not including the boundary constraints. It

The interaction energy E e
With three rods, d 0 = 1/(m − 1) = 1/2, and the total spring energy is: Here, the integral is the elastic deformation energy associated with the shaded triangle shown in figure 3a, which takes into account both the energy of compressing the springs below the triangle, and extending the springs above. The factor of 1/L is used to ensure that as we vary L, there are the same number of springs per unit length to match our simulation set-up. The prefactor of n accounts for the total deformation energy of the 2n copies of the shaded triangle and the fact that there are six spring connections for each point in our model (four diagonal, two direct, so that in our earlier notation this is the sum E 1 e þ E 1þ e þ E 1À e ). We deduce Finally we note that, as we seek to derive an instability criterion, we can assume the overlap potential would be negligible and hence equation (4.5) represents the total interaction energy of the system.

The energy E b
To derive E b , we take a similar approach as in the derivation of the discretized simulation bending energy. The acute angle θ t of the shaded triangle is We approximate the curvature energy as the difference between the two angles of the triangle as marked in figure 8a. Using the arctan formula of our model, we define the curvature κ t of each triangle as Thus the curvature energy, averaged over the 2n triangles is ð4:8Þ

Instability as a balance between bending and matrix stretching
We take the derivative of equation (4.2) with respect to Δ (using equations (4.5) and (4.8)) and find the critical excess force N e when Δ = 0 to be The first term represents the energy due to spring stretching. It decreases with n as the fixed length of the filaments means higher mode deformations have lower amplitude and hence less stretching of the elastic matrix joining them. The second term represents the bending energy which increases with n. Thus, we see the minimum of this function will be intermediate in n. An optimal n, denoted by n*, is found by minimizing N e with respect to n. Substitution of this lowest-load mode n* into equation (4.1) returns the critical buckling load royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 19: 20220287 stiffness. Reflecting on this, the coefficient 8 is entirely determined by our approximation of the shape of the buckled configuration as a series of triangular kinks, and so cannot be expected to accurately represent the true minimum energy buckled configuration. It is reasonable therefore to treat the pre-factor as a fitting parameter α. As mentioned above, we add on the critical Euler load which accounts for the boundary conditions to give equation (4.1).

Accounting for neighbouring deformations
To account for the decay in buckling amplitude from the central of the bundle when m > 3 (seen in figure 8c), we assume neighbouring rods are also buckled, as in the red curves shown in figure 8b. This leads to a lower effective stretching d* which will be less than the value d used in the three-rod case described above. We should expect the variation in neighbouring rod curvature to vary as a function of the undeformed separation 1/(m − 1), the equilibrium separation of the rods. Thus, we propose that the stretching d* takes the form for some constant α. Since d 0 = 1/(m − 1), the factor of 1/(m − 1) cancels in the ratio d*/d 0 used in equation (4.3). Hence, working through the same steps as in the above section, we obtain the same critical buckling criterion in equation (4.1), but now we recognize the parameter α is additionally accounting for the variation in curvature across the filaments in the bundle.

The constant α
For the results in figure 8b, a value of α = 10.61 provides an excellent fit to the data across a wide variety of parameter sets (B in , K, L) for a bundle with m = 10. Our simple assumption in equation (4.1) is that α is independent of all the system parameters. We now further test the variation in the value of α with the bundle number m. To do so we locate the critical buckling N Ã c numerically by steadily increasing N c until a perturbation about the undistorted bundle state is found to develop internal buckling. We then use the predicted internal buckling formula in equation (4.1) to calculate the value of α for each set (K, B in , L, m). The results for a range of parameters are shown in figure 9. For lower length filaments the parameter is approximately constant with m, whereas for the longer filaments there is some weak decay with m.
This implies that the simple first-order model misses some of the more subtle aspects of the general nonlinear model. For instance in the general model, the deformation shape is more complex relative to the triangular approximation, the deformation magnitude is suppressed at either end of the filament by the boundary conditions, the end displacement at the unpinned end is not uniform, and the outer rods may deform slightly. However, all values of α are of the same order of magnitude as the value of 8 in equation (4.10), indicating the model is predicting the correct scaling, resulting from balancing the elastic and bending energies.

A note on confinement and wavemode
It is well known that filaments in embedded media buckle at higher wavemodes than the free filament model (see [7][8][9]41]), with the mode determined by the stiffness of the embedding/ constraining medium. In this case, the stiff outer filaments/ sheath act to confine the filaments. The formula for the critical mode n* in equation (4.10) indicates that the spring stiffness K and bundle length to width ratio L promote higher mode buckling. The confinement means that the peak buckling amplitude experienced by the central filaments must decay in the filaments to the outer edge of the bundle and this differential deformation leads to stretching of the spring medium, which, as discussed above, is lower for high wavemodes. Increasing K means the benefit of adopting high wavemodes and decreasing this differential stretching is increased, thus promoting higher modes. One can also see from equation (4.3) that increasing the length L leads to an increase in the elastic energy for a given mode n, thus the energy can be lowered by increasing the mode n. Once again this energy is only important because of the confinement leading to differential deformation of the filaments and hence spring stretching. By contrast, it is unsurprising that increasing the stiffness of the filaments B in promotes lower wavemodes, a factor not directly related to the confinement.

Conclusion
We have developed a constitutively simple (in terms of its material parameters) model for interactive elastic filament bundles via an energy functional, equation (2.2), which accounts for: the filament bending, the filament interaction and non-overlapping constraints, bundle loading and boundary conditions. This minimal model exhibits a diverse range of free energy minima, which we categorized into four distinct classes: unbuckled, globally buckled, pinched and internally buckled. Importantly, these follow many biologically observed buckling phenomena, such as those shown in figure 1. These numerical experiments were shown to mimic biologically observed buckling responses.
The model is detailed in §2. A key feature we highlighted, to ensure numerical and physical reliability, was the non-local overlapping potential. This was critical to mediate the complex nonlinear interaction inherent to the model. We believe this insight could be of significant benefit to researchers developing similar theoretical models.
A common feature of many biological filament bundles is the embedding of the system in a stiff matrix or sheath. Interestingly, for sufficiently stiff outer sheaths, the buckling mode was observed to change dramatically from a Euler-like global buckling to an internal buckling mode. We derived both the royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 19: 20220287 scaling and analytical formula to predict the critical compressive load leading to this internal buckling transition. As internal buckling leads to large, highly localized pressures (relative to a simple pinching mode), we predict internal collapse may describe catastrophic mechanical failure in biological systems.
Furthermore, the derived critical compressive load, given in equation (4.1), is remarkably straightforward. We highlight that the relevant parameters of this prediction should be relatively uncomplicated to estimate. The bending stiffness B in can be derived from the filament's Young's modulus and Poisson ratio, while the spring/matrix stiffness K is a measure of the resistance to deformation of the material/ bonds connecting the filaments. Thus, we believe it has significant predictive potential for applications of biological filament bundles in supported environments, where internal filament collapse plays a key role.
Possible future investigations using this model might also address what happens if the filaments are of different types (they have variable bending stiffness), or have structural weakness (bending stiffness varying with arclength). One can also consider variations on the pinching force such as pinching at multiple points or the pinching/pulling of individual filaments ( plucking). In each case, one would have access to a richer variety of geometries than presented here, but we believe the principle of competition between elastic matrix stretching and filament bending highlighted here would still play a key role in determining these geometries.
The model presented here is two-dimensional, but has relevance to three-dimensional bundles under the condition that there are no torsional stresses. Under such conditions, the pinching modes presented here could represent a crosssection of an axially symmetric bundle pinching with uniformly applied pinching pressure. This would be relevant for the optic nerve where the hoop stress in the scleral wall applies a uniform pressure to the lamina cribosa [4]. The global buckling mode represents a planar buckling of a threedimensional bundle, this behaviour has been observed in related global buckling models [19]. Some caution can be found in results in a related model of a filament connected to a fixed surface found in [41]. In this case, there is some indication that the spring connection can lead to a competition between non-planar and planar buckling modes for the lowest-load collapse, although both have similar wavelengths. Similar considerations would apply to the internal buckling modes. However, it remains the case that the internal buckling wavelength would still be determined by the balance between the elastic stiffness penalizing higher wavemodes and the interaction force promoting them. This is a critical principle underlying internal buckling modes in such systems.