Fibre crosslinking drives the emergence of order in a three-dimensional dynamical network model

The extracellular-matrix (ECM) is a complex interconnected three-dimensional network that provides structural support for the cells and tissues and defines organ architecture as key for their healthy functioning. However, the intimate mechanisms by which ECM acquire their three-dimensional architecture are still largely unknown. In this paper, we study this question by means of a simple three-dimensional individual based model of interacting fibres able to spontaneously crosslink or unlink to each other and align at the crosslinks. We show that such systems are able to spontaneously generate different types of architectures. We provide a thorough analysis of the emerging structures by an exhaustive parametric analysis and the use of appropriate visualization tools and quantifiers in three dimensions. The most striking result is that the emergence of ordered structures can be fully explained by a single emerging variable: the number of links per fibre in the network. If validated on real tissues, this simple variable could become an important putative target to control and predict the structuring of biological tissues, to suggest possible new therapeutic strategies to restore tissue functions after disruption, and to help in the development of collagen-based scaffolds for tissue engineering. Moreover, the model reveals that the emergence of architecture is a spatially homogeneous process following a unique evolutionary path, and highlights the essential role of dynamical crosslinking in tissue structuring.


Introduction
The adequate architecture of any organ is mandatory for their efficient physiological function and any changes are associated with function impairment and putative developing dysfunctions and diseases [1][2][3].All biological tissues contain scaffolds of non-cellular components called extracellular matrices (ECMs) [4].Despite the great variability of proteins that make up the ECM (macromolecules such as collagen, glycoproteins etc.), it can be seen as a dynamic physical network of fibres interconnected by molecular bonds, i.e. crosslinks, generating a connected and elastic environment for the surrounding cells [2].
The network structure is in a state of constant remodelling, which is crucial to maintain tissue integrity and function.Crosslinks, however, can unbind spontaneously or under tension, which leads to viscoplastic material responses, such as softening and tension relaxation [5].Fibrosis and ageing are also characterized by an increase of enzymatic and non-enzymatic crosslinks [6,7] and this increase in crosslinking prevents ECM degradation by matrix metalloproteinases, both events leading to a decrease of ECM remodelling [8].Altogether, these events induce greater stiffness and the arrangement of the collagen fibres becomes less organized and more loose and fragmented, hence weakening tissue integrity and strength [9,10].An understanding of the basic organizing principles of ECM structure in three dimensions also helps in apprehending the complex dynamics of pathological tissues from degenerative diseases or tumour [8].
Because the global architecture of fibre networks seems to be fundamental for controlling tissue functions, modelling the process of ECM structure emergence will greatly improve our understanding of tissue biology and plasticity in physiological or pathological conditions.Numerous models of fibre networks can be found in the literature.Due to their simplicity and flexibility, the most widely used models are individual based models (IBM), which describe the behaviour of each agent (e.g. a fibre element) and its interactions with the surrounding agents over time [11,12].However, IBMs have a high computational cost which can become intractable when studying systems composed of too many agents, or systems at large scales, either spatial or temporal.In such cases, continuous or mean-field kinetic models may be preferred [13][14][15][16][17] since they are less costly, but at the expense of a loss of information at the individual level.Since it is well acknowledged that microstructure configurations modulate the macroscopic properties of crosslinked fibre networks [18], preserving the microscopic level description is of great importance to model tissue emergence.
Most of the computational models developed thus far for mimicking ECM networks are twodimensional [14,16,[19][20][21][22][23][24][25][26][27].Few studies have been conducted on three-dimensional models [28][29][30][31][32][33][34][35], although these are expected to yield different, more realistic results than two-dimensional ones since they better mimic biological structures themselves immersed in three-dimensional environments.One of the reasons for fewer three-dimensional models is the great increase in the number of agents needed to achieve a given spatial density and thus in the associated computational cost.Another reason is the lack of high quality data on ECM organization in three dimensions.However, the latter is becoming less and less of an issue with recent improvements in high resolution three-dimensional imaging and its availability.Among existing three-dimensional models, few of them feature dynamical crosslinking of ECM components.In [30,32,36], various models of three-dimensional fibrous networks composed of permanent or transient crosslinks (remodelling) are proposed.However, most of these models feature ECM remodelling in reaction to external factors (applied load [30,37], migrating cells [32], contractile cells [36]), and the literature so far provides few cues on the mechanisms underlying fibre self-organization.
In the present paper, we test the hypothesis that fibre macrostructures could spontaneously emerge without appealing to contact guidance or external mechanical challenges, as a result of simple mechanical interactions between the fibre elements composing the ECM network.We assess this hypothesis by means of a simple three-dimensional model where ECM fibres are discretized into unit fibre elements, consisting of non-stretching and nonflexible spherocylinders with the ability to spontaneously link to and unlink from their close neighbours.This dynamical crosslinking mechanism allows us to model both the overall temporal plasticity of the network and the complex physical properties of biological fibres such as elongation, bending, branching and growth, thus compensating our minimalistic description of the fibre units.We stress that in this paper, rather than developing a very complex model to reproduce the whole complexity of real tissues (at the cost of losing explicability), a large part of which corresponds to the redundancy of mechanisms to ensure the robustness of structures and regulations, we aim to keep a mathematical framework as simple as possible in order to break the complexity and shed light on some main and basic components at play royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.11: 231456 in the emergence of fibrous structures.The relevance of such an approach was previously validated in the frame of adipose tissue morphogenesis and regeneration in two dimensions [27,38].
Through computational simulations and exhaustive parametric analysis, we demonstrate that organized macrostructures can spontaneously emerge without external guidance.Overall, this study provides a comprehensive view on the role of ECM connectivity on tissue architecture emergence: -The model reveals that tissue architecture at equilibrium is simply controlled by the number of crosslinks per fibre in the network, an emerging variable not directly linked to the model parameters.If further validated on real tissues, this simple emerging variable could become an important putative target to control and predict the development of the architecture of biological tissues.Because of its simplicity, this variable is amenable to experimental measurements and could represent a major target for the development of therapeutic drugs to induce tissue recovery after injury, prevent tissue degradation during ageing, or help in the design of engineering collagen scaffolds for tissue regeneration.-A deep exploration of the model parameters reveals that this emerging variable, and therefore the global organization abilities of tissues, depend on a complex interplay between the model parameters related to the crosslinks, i.e their remodelling speed and their linked fibre fraction.These results rationalize how even subtle changes in fibre networks dynamical crosslinking can drive tissue reorganization and suggest that the development of biological crosslinkers to control ECM connectivity as a target for tissue reconstruction must carefully account for different parameters such as tissue remodelling activities.-Finally, a temporal analysis of the model simulations reveals that the different tissue architectures follow a simple and unique evolutionary path on timescales controlled by their remodelling characteristics, providing new insights into the temporal evolution of tissue structures as a function of the ECM remodelling properties.

Description of the model
The three-dimensional ECM is discretized into unit fibre elements consisting of line segments of fixed and uniform length, represented by their centres and directional unit vectors.We consider the following biological and mechanical features: (i) Fibre resistance to pressure: We suppose that fibre elements repel each other at short distances, which models size-exclusion effects.This is achieved via a repulsive force between close fibres based on Hertzian theory [39].This amounts to model fibres as spherocylinders of a given radius and length, that can interpenetrate each other.The intensity of the repulsion force α rep controls the amount of overlapping between fibres.(ii) Fibre elongation and breakage: In addition to carrying a unit of ECM fibre strength, fibre elements also carry a unit of fibre length.However, we provide a way to create longer fibres by allowing two nearby fibres to form a link.A crosslink is modelled as a linear spring with a given spring stiffness, connecting the two closest points of the fibre pair at the time of its creation.There is no prescription for the location of the crosslinks along the body of the fibres they connect.Several consecutively cross-linked fibre elements would model a long and flexible fibre having the ability to adopt complex geometries.Therefore, the cross-linking process models fibre elongation [40].The stiffness constant of the springs α rest controls the possible extension of the long fibres.Symmetrically, pairs of cross-linked fibres can spontaneously unlink, allowing for fibre breakage describing ECM remodelling processes [41].Linking and unlinking processes follow Poisson processes with frequencies ν link and ν unlink , respectively.As a result, the linked fibre ratio x link ¼ n link =ðn link þ n unlink Þ represents the equilibrium fraction of linked fibres among the pairs of neighbouring fibres.(iii) Crosslink fibre alignment: To model the ability of long fibres (those made of several cross-linked fibre units) to offer a certain resistance to bending, linked fibres are subjected to a potential torque at their junction.This torque vanishes when the fibres are aligned, and consequently acts as a linked-fibre alignment mechanism.This torque is characterized by a stiffness parameter α align playing the role of a flexural modulus.(iv) Large friction regime: As the Reynolds number in most biological tissues is very small [42], we suppose that inertial forces can be neglected and we consider an over-damped regime for fibre motion and rotation.Each of the mechanical interactions due to fibre-fibre repulsion (i), fibre-fibre attachment due to crosslinks (ii) and crosslinked fibre-fibre alignment (iii) generate elementary forces and torques between fibre pairs.The total force (resp.torque) acting on a fibre is then computed as the sum of all royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.11: 231456 the elementary forces (resp.torques) generated by the elements interacting with this fibre.The motion and rotation of the fibre is then deduced from Newton's equation of motion in an over-damped regime.More specifically, the N fib fibre elements are represented by straight lines of fixed length L fib represented by their centres X k ðtÞ [ V , R 3 and their non-oriented directional unit vectors v k ðtÞ [ S þ 2 .Moreover, from the fibre-fibre repulsion interaction, fibres may be seen as soft spherocylinders of radius R fib .We denote by ( p k,m (t)) k,m the fibre connectivity matrix, that is p k,m (t) is equal to 1 if fibres k and m are linked at time t and to 0 otherwise.
The motion and rotation of fibre k are then given by

Description of the simulation set-up and biological relevance of the model parameters
The spatial domain V is a cuboid of side lengths L x , L y and L z , respectively, in the x, y and z-dimension, centred on the origin For the sake of simplicity, we assume periodic boundary conditions: an agent exiting the domain by one side re-enters immediately from the opposite side, and interactions between agents are computed using the periodicized Euclidean distance.Fibres are initially randomly inseminated inside the domain according to a uniform law for both position and orientation.The differential system (2.1) is then numerically solved using a discrete upwind Euler scheme with adaptive time step, which has a very low computational cost.Details of the numerical implementation are given in appendix A.2.
The physical scaling of all the parameters of the model, as well as the values used in the simulations, are described in table 1.A few points may be noted: (a) the perception distance for link creation d max link and the link unloaded length d eq link are both equal to the diameter of a fibre 2 R fib .This means that the fibre units (spherocylinders) connect with the fibres they are in contact with or closer, and that the link tries to keep the bodies of the spherocylinders touching.In this regime, the presence of the links therefore participate in a non-overlapping configuration of the fibres.(b) The size of the domain is approximately four times the size of a fibre along its main axis (numerical checks were made by-hand to select a size of domain which optimizes between computation time and boundary effects) and (c) the fibre aspect-ratio L fib /2 R fib = 6 is quite small compared to the values used in other models of the ECM, which usually varies between 250 and 10 4 [15,33,34].This compensates for the fact that these models directly account for fibre bending and/ or fibre elongation, while our long fibres correspond to a sequence of crosslinked fibre units.On the same note, we stress the fact that our fibre units do not aim at modelling the individual collagen fibrils making up collagen fibres in ECM, but rather correspond to an intermediate scale where one fibre unit of our model is already a set of twined collagen fibrils that run in parallel to form a larger bundle [43].
We denote by ϕ fib the fibre density of the network, that is the ratio between the total volume of fibres (without overlapping) and the volume of the spatial domain: The quantity ϕ fib can be compared to the packing density, that is the maximal fraction of the domain that can be occupied by densely packed fibres.In the case of an ordered packing, the packing density of spherocylinders is ϕ order = 0.89, while for random or amorphous packing of spherocylinders with an aspect ratio of 6, the maximal random packing density ϕ random ≈ 0.4 [20].Thus, we may say that a system is 'sparse' if its fibre density is below ϕ random , 'dense' if it is between ϕ random and ϕ order , and 'hyperdense' if it is above ϕ order .In the following, we will study two types of systems: dense systems royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.11: 231456 containing N fib = 3000 fibres (ϕ fib = 0.58) and sparse systems with N fib = 1500 fibres (ϕ fib = 0.29, corresponding to measurements of extracellular volume fraction in muscle or myocardial fibrosis, see e.g.[44,45]).
For each of the three types of mechanical forces in the system, we define the 'characteristic interaction time' as the time needed for two isolated fibres interacting only via this force and initially positioned in the most unfavourable configuration to reach 99% of the equilibrium state.For repulsion, T rep is the time needed for two fully overlapped fibres (X 1 = X 2 and ω 1 = ω 2 ) to move apart by 99% of their equilibrium distance 2 R fib (i.e.‖X 1 − X 2 ‖ = 0.99 × 2 R fib ).Similarly, for the elastic spring T rest is the time needed for two fibres that are initially fully overlapping and crosslinked at their centre to move apart by 99% of their equilibrium distance d eq link .On the other hand, for nematic alignment T align is the time needed for two perpendicularly intersecting fibres (X 1 = X 2 and ω 1 ⊥ω 2 ) crosslinked at their centre to reach a relative angle arccosðv 1 Á v 2 Þ ¼ 0:9 .
Explicit computation leads to the following formula (numerical values are given for the parameters presented in table 1) T align ¼ 4:27 It may be noted that the alignment interaction is much slower than the repulsive and elastic restoring forces.In this regime, fibre elements are quite rigid and connected by strong springs (crosslinks), enabling us to prevent local accumulation of fibres and overstretching of long fibres (those made of several crosslinked fibre units).

Matrix crosslinking drives the local alignment of three-dimensional dynamical fibre networks
In figure 1a-c, we show various structures that can be obtained with our model by playing on the parameters in the ranges indicated in table 1.The fibres are represented by double arrows and coloured as a function of their local alignment with their neighbours.We refer the readers to appendix B.1 for more details on the computation of this quantifier, and just mention that the local alignment of fibre k, denoted Al k , is equal to 1 (fibre coloured in red) if all the neighbouring fibres display the exact same direction as fibre k, and to 0 (fibre coloured in blue) if the neighbouring fibres display uniformly distributed directional vectors.Moreover, we show in appendix B.1 that this quantifier is able to discriminate between fibres located in randomly oriented environments (corresponding to Al k < 0.5), fibres located in nearly planar environments (leading to Al k around 0.7), and fibres located in nearly uni-directional environments (leading to Al k above 0.8).
As one can observe, the fibre structures obtained at equilibrium range from highly aligned systems (mainly composed of red fibres, see figure 1a) to disordered systems with a low local alignment (mainly composed of fibres coloured in blue, see figure 1c).The model can also produce intermediate states composed of fibres with a median local alignment (figure 1b).
In order to assess the alignment states of our different fibre networks, we computed the mean of the local alignment indicator Al k over all the fibres of the system, denoted by Al sim .To account for stochastic variability (due to the random initial condition and the stochastic linking and unlinking processes), we computed the mean and standard deviation of Al sim over 10 simulations conducted with the same set of parameters, denoted by Al mean and Al STD .Similarly, we denote by N linkperfib = N links /N fib the number of royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.11: 231456 links per fibres in a network and by N mean linkperfib and N STD linkperfib its average and standard deviation over 10 simulations.We stress the fact that N linkperfib = 0.5 if all the fibres are connected to a neighbouring fibre.
By plotting the alignment quantifier Al mean as a function of the number of links per fibre N mean linkperfib (both computed on the systems at equilibrium), we discovered a striking and major correlation between these two quantities.This correlation is shown in figure 1b, with horizontal and vertical error-bars indicating the inter-simulation standard deviations N STD linkperfib and Al STD , respectively.The different markers indicate different fibre densities (dots for dense systems and triangles for sparse ones), the different colours refer to different networks dynamics ν link , and inside each colour series χ link is increasing with N mean linkperfib .Figure 1d reveals that the values of Al mean and N mean linkperfib at equilibrium are highly correlated.When N mean linkperfib is inferior to a critical threshold N critic ≈ 0.7 (indicated with a grey dashed line on figure 1d ), there is a logarithmic correlation between the number of links per fibre in the network and its mean alignment indicator (black dashed lines in figure 1d) with -α = 0.037, β = 1.006 and coefficient of determination r 2 = 0.87 for dynamical systems (non-blue markers); -α = 0.129, β = 0.651 and coefficient of determination r 2 = 0.96 for sparse non-dynamical networks (blue triangles); -α = 0.042, β = 0.433 and coefficient of determination r 2 = 0.985 for dense non-dynamical networks (blue dots).
Then, when N mean linkperfib .N critic we observe an abrupt drop of the equilibrium value of Al mean .Surprisingly and very interestingly, for dynamical systems (ν link > 0) there is no difference in alignment induced by the fibre density or the link characteristics ν link and χ link : the correlation observed is the same for all sets of points.
The second major observation from figure 1d is the difference between non-dynamical and dynamical networks at equilibrium.Indeed non-dynamical networks, composed of a fixed number of links, are systematically less aligned than dynamical ones (compare the values of Al mean between the blue markers and the other colours).Moreover, although we do recover the same type of correlation between the fibre local alignment and the number of links per fibre in the network, for nondynamical networks this correlation significantly depends on the fibre density.However, the critical number of links N critic allowing for larger alignment is the same for non-dynamical networks, either dense or sparse, and for dynamical networks.Therefore, N critic seems to be a general critical network connectivity value controlling the local alignment abilities of various networks.
Altogether, these results show that the emergence of organized networks (i) requires some remodelling abilities of the ECM matrix and (ii) is mainly controlled by the number of links per fibre.

ECM architecture emergence is driven by a complex interplay between remodelling speed and linked fibre fraction
The previous section took a particular focus on the local arrangement of the fibre units composing our three-dimensional fibre network, with little information on the global structures at the population scale.
In this section, we aimed to characterize quantitatively the macrostructures that emerge in our networks.
To this end, we used the stereographic projection of the fibre directional vectors.Disregarding the spatial position of a fibre, we represented its directional vector as a point on the surface of the unit half-sphere in three dimensions and then projected it onto the unit disk in two dimensions (see appendix B.2 for a detailed explanation).
As shown in figure 2, this representation enabled us to characterize the different global organizations of our fibre networks.Indeed, we observed three different types of stereographic projections in our simulations: fibres' directional vectors very concentrated around the centre of the disc, corresponding to a global alignment of the system (figure 1a, with stereographic projection shown as inset in figure 2a), fibres' directional vectors homogeneously distributed on the disc corresponding to a global disorder (figures 1c and 2e), and fibres' directional vectors distributed along a preferential axis, with royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.11: 231456 complete depletion in the direction perpendicular to this axis, corresponding to global curved/plane structures (figure 2b-d).
Together with the local alignment quantifier Al sim , we were now able to quantitatively characterize the different local and global fibre organizations inside our networks.We considered a system to be locally aligned if Al sim was above 0.7 (see appendix B.1 for justification of this value).At the same time, we considered that a system was globally aligned if its stereographic projection covariance ellipse had a semi-major axis smaller than 0.45 (implying that the point cloud covers less than 20% of the whole projection disk).
We therefore classified the simulations outcomes into three different states (unorganized, curved and aligned) using table 2. We ran a total of 1080 numerical simulations, exploring various values of the parameters ν link , χ link and N fib in the broad ranges indicated in table 1, and counted among their outcomes: -180 unorganized states (all occurring in non-dynamical systems, i.e. ν link = 0), -661 curved states, -239 aligned states (among which only 12 occurred in sparse systems).
Figure 2a shows the equilibrium values of quantifiers Al sim and A max for dense systems (see electronic supplementary material, appendix C.2 for the equivalent figure on sparse systems).The points are coloured according to the states defined previously (blue dots correspond to unorganized states, orange diamonds to curved states and red crosses to aligned states).The simulations already From figure 2a, we first observe that the unorganized states (blue dots) form a small, compact group of points with large semi-major axis length, while the aligned states (red crosses) make a long thin group with very high alignment indicator.On the other hand, the curved states (orange diamonds) form a scattered cloud of points with a broad range of values for both the semi-major axis length and the alignment indicator.Moreover, we observe that the transition between unorganized and curved states is very sharp (note the gap between the blue dots and orange diamonds in a).Indeed, no simulation displays an average alignment indicator at equilibrium between 0.65 and 0.77 (including sparse systems, see electronic supplementary material, appendix C.2), and there is a marked difference between the least organized of the curved states (illustrated in figure 1a) and the most organized of the unorganized states (illustrated in figure 2e).This confirms our choice of 0.7 for the threshold value between unorganized and curved states.
On the contrary, the transition from curved to aligned states is not a clear switch but a continuum of structures that can be illustrated by the two borderline cases in figure 2b,c.Thus, one must be aware that the partition between curved and aligned states is partly arbitrary and depends on the choice of the threshold.However, this classification into three states allowed us to distinguish between unorganized networks, globally aligned networks and networks locally aligned with twisting capacities at the population level, enabling us to go deeper into the model parameters controlling tissue architecture emergence at different scales.
We first found that the sharp transition between unorganized and curved states was fully controlled by the remodelling speed of the network ν link .Indeed, unorganized states were only and systematically observed for non-dynamical networks (ν link = 0), while dynamical networks (ν link > 0) never equilibrated in unorganized states but self-organized into either curved or aligned states, and this independently on the fibre density of the network (see electronic supplementary material, appendix C.2).By contrast, the transition between curved and aligned states is not controlled by a unique model parameter but is the interplay between several parameters.Indeed, figure 2f shows a heatmap of the percentage of simulations ending in an aligned state (versus a curved state), for dynamical dense networks (see electronic supplementary material, appendix C.2 for results on sparse networks), depending on the values of the network remodelling speed ν link (in ordinate) and the equilibrium linked fibre fraction χ link (in abscissa).As one can observe in figure 2f, there is a nonlinear relationship between the global alignment capacities of the networks and the parameters ν link and χ link .Indeed, analysis of the heatmap reveals that (i) reduced linked fibre fraction χ link can increase global alignment outputs because, for low-remodelling networks, the formation of crowded interconnected fibre structures inhibiting fibre motion is relieved by reduced link density.Moreover, (ii) the global alignment of networks with intermediate remodelling rates may undergo little change with reduced linked fibre-fraction and (iii) the global alignment ability of fast-remodelling networks will likely be impaired by reduced linked fibre-fraction.These results show that the different types of tissue architectures (aligned, curved or unorganized) depend on an interplay between parameters ν link and χ link .While ECM local alignment can be explained by the simple emerging variable that is the number of links per fibre in the network (as shown in §3.1), its direct relation with model parameters N fib , ν link and χ link is more complex.Indeed, figure 2g shows a heatmap of the number of links per fibre in the network N mean linkperfib as a function of ν link and χ link for dense dynamical networks (same simulations as f ).It demonstrates that N mean linkperfib is indeed an emerging variable, in the sense that it is not directly linked to the parameters ν link and χ link but rather is the result of a complex interplay between the two.Indeed, the number of links per fibre in the network increases along the diagonal, as ν link decreases and χ link increases (from top left to bottom right corner of the heatmap), crossing the critical threshold N critic doing so (the cells where N mean linkperfib % N critic are framed in green in g).These results explain why the proportion of aligned structures are maximal along the diagonal from the bottom left to the top right (i.e broadly perpendicular to the gradient of the emergent parameter).These results extend to the case of sparse networks (see electronic supplementary material, fig.8 in appendix C.2), confirming the strong correlation between ECM alignment abilities and the number of links per fibre they contain.
These results show that our networks can be seen as corresponding to different phases of physical materials depending on their remodelling abilities.If non-dynamical networks can be seen as solid structures unable to spontaneously reorganize, dynamical networks have properties reminiscent of fluid materials, the global architecture of which being controlled by an interplay between their remodelling speed and their linked fibre fraction.In the next section, we study the evolution in time of the structures, enabling us to give more insights into the role of these parameters in tissue structuring in time.

ECM architecture emergence follows a unique evolutionary path on timescales controlled by their remodelling characteristics
In this section, we study the temporal evolution of the spatial structures.Our very first observation is that, for all sets of parameters, the evolution in time of the quantifier Al mean follows a logarithmic growth (see electronic supplementary material, appendix C.4 for more details).We will use as a time reference the time-constant of this growth, denoted τ Al , which corresponds to the time needed for the quantifier to reach 63% of its asymptotic value (Al mean ≈ 0.7 in our case).
Movies displaying the full temporal evolution of a few simulations are available in supplementary data (see appendix C.1).In figure 3a-a"' and b-b"', we show the stereographic projection of a few well-chosen time frames (namely 0.5τ Al , τ Al , 3τ Al and T final ) for two of these simulations (respectively from Movie3.mp4 and Movie4.mp4).They correspond to dense systems with χ link = 0.8 and two different crosslink dynamics: fast remodelling network ν link = 0.1 (a-a"', Movie3.mp4) and slow remodelling network ν link = 0.001 (b-b"', Movie4.mp4).These screenshots enable us to answer the important question of how the network global structure emerges.It is not by accretion around a few structured areas that gradually merge together, but by an overall homogeneous structuring.Indeed, one can observe that the directional vectors gradually concentrate around a main direction without creating clustered points that merge together.This behaviour can be observed both for very aligned networks (A-A"') or curved states (B-B"'), and in fact in all our simulations, independently on the network density.Therefore, our model suggests that the emergence of tissue architecture occurs on a global scale.
We now turn towards the analysis of the time trajectories of the quantifiers of the structures.We show in figure 3c, the trajectory in the phase plane A max versus Al sim of simulations for low-dynamical dense networks ν link = 0.001 with various linked fibre fractions χ link (different colours, see electronic supplementary material, appendix C.4 for more dynamical networks).We observe that all the trajectories follow a common pattern.It begins with a sharp increase of the alignment indicator (from 0.15 to between 0.4 and 0.5) while maintaining a quasi-constant semi-major axis length: this corresponds to the partial depletion of one direction (denoted d 1 ) in the family of the fibres' directional vector, thus shifting from the initial uniform distribution to a mainly two-directional distribution (see appendix B.2 for more details on this interpretation).Non-dynamical networks do not go past that first stage (data not shown).
The trajectories then diversify: the alignment indicator keeps increasing while the semi-major axis length either decreases, stays constant or slightly increases.The first case is the most common and indicates that, while direction d 1 keeps depleting until near extinction, one of the two remaining directions starts to deplete as well.This diversification happens on the scale of the time-constant τ Al of the alignment indicator (marked on the trajectories of figure 3c with a black circle).
Lastly, simulations ending in an aligned state and part of those ending in a curved state display a stage of condensation of the fibres directional vectors around a main direction.This is marked by a shrinking of the covariance ellipse and a slow increase of the alignment indicator, which has already nearly reached its steady state (compare with the stabilization of Al mean in figure 10).This last point comes from the local quality of the quantifier Al sim (and by extension Al mean ): a system can be very royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.11: 231456 aligned locally, but not globally, if the main direction of the local structures varies smoothly across space.Thus, the transition between a curved and an aligned state is mostly characterized by a gradual shifting of multiple local structures towards the same direction.
Finally, we observe that the number of links per fibre (displayed in figure 3d) undergoes a transient increase followed by a two-stage exponential decay in time (appearing as a piece-wise linear decrease on the semi-logarithmic scale).For low dynamical networks, the initial accumulation of crosslinks is more pronounced, in the sense that the peak is higher and the subsequent decrease slower, when χ link is high.For the extreme case of large linked fibre fraction χ link = 0.9 ( pink curve in figure 3d), the phenomenon is so strong that only the first stage of exponential decay is observed during the time of the simulation.On the other hand, for small equilibrium linked fibre fraction (χ link = 0.1, blue curve), we do not observe any crosslinks accumulation or fast remodelling networks (see electronic supplementary material, appendix C.4 for more dynamical networks).
This behaviour can be explained by comparing the linking dynamics to the characteristic time of the repulsive interaction T rep = 4.32 U t .Parameter χ link describes the proportion of linked fibres among all linkable fibres at equilibrium, but this equilibrium takes time to establish (inversely proportional to ν link ).If the repulsion interaction operates faster than the links remodelling (i.e.T rep ≪ 1/ν link ), then the linkable configurations will change before the linking/unlinking processes could equilibrate on the current configuration: new links will appear between newly overlapping fibres while former overlapping fibres will still be linked even if not overlapping anymore, leading to an accumulation of links in the system.This happens all the more if the disparity between the frequencies ν link and ν unlink is more favourable to linking than unlinking (ν link > ν unlink , i.e. if χ link > 0.5).royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.11: 231456 The system thus exhibits a global, macroscopic relaxation phenomenon which emerges from its various local, microscopic properties.It can be seen that the characteristic time-scale of this relaxation is comparable to the time-constant of the alignment indicator τ Al (see position of the black circles on the curves in figure 3d, which indicates the value of τ Al for the corresponding set of parameters).
These results demonstrate a nonlinear dependence of the network properties on the type of links and the number of crosslinks per fibre.A high number of long-lasting crosslinks promotes crosslink accumulation resulting in medium/low alignment, while fast remodelling reduces the mechanical action of the individual links on the overall network, resulting in lowly connected networks being unable to align.Together with the results of §3.2, we showed that the network alignment abilities require a number of links adapted to their remodelling speed: fast remodelling networks need a high equilibrium linked fibre fraction to quickly reach a high alignment indicator, while slow remodelling networks need a medium/low equilibrium linked fibre fraction to prevent crosslink accumulation and avoid the formation of crowded interconnected fibre structures inhibiting fibre motion.

Discussion
In this work, we have implemented a three-dimensional model for fibre networks composed of fibre elements capable of dynamically crosslinking or unlinking from each other, to align with each other at the crosslinks and to repel their nearest neighbours to prevent cluttering.We showed that this model can spontaneously generate various types of macrostructures whose emergence can be finely described.The model reveals that the different macrostructures (i) can be easily explained by a single emerging intermediate variable, namely the number of links per fibre in the ECM network, (ii) are controlled by a nonlinear relationship between the linked fibre fraction and remodelling rate and (iii) follow the same unique evolutionary path for all structures and not multiple paths.
To our knowledge, this work is the first exhaustive study questioning the mechanisms of tissue architecture emergence via a simple mechanical model of dynamical fibre networks in three dimensions.The equilibrium structures obtained with our model can be classified into three types: (a) aligned states with a strong organization around one main direction, (b) curved states with a median, locally heterogeneous alignment indicator and a wide range of directional vectors living in a plane, named curved patterns and (c) unorganized states with very low alignment indicator and no preferential direction.These different types of macro architectures show that the model can cover a wide range of biological tissues, from highly aligned fibre structures reminiscent of muscular tissues [46] to disturbed alignment of collagen fibres observed in the first phase of wound healing [47].Unorganized states were exclusively obtained for non-dynamical networks composed of permanent crosslinks (ν link = 0), whose plasticity was very low due to their inability to rearrange their crosslinks.By contrast, dynamical networks exhibited a mixture of aligned and curved states.These results point to the essential role of matrix remodelling in ECM structuring, consistent with several results in the literature (see [48] and references therein).
This framework reveals that the different tissue architectures at equilibrium are directly controlled by a simple intermediary variable, the number of links per fibre (see §3.1).Our interpretation is that, when the number of links per fibre is inferior to the critical threshold N critic , the network is weakly constrained.In this configuration, an increase in the number of links per fibre improves the transmission of information in the network and thus enhances the alignment process.The logarithmic scaling indicates that the higher the number of links per fibre, the less prominent this feature becomes, until the gain (in terms of the equilibrium alignment indicator) becomes null.The system then shifts into a constricted regime where each new link adds to the constriction of the network and impedes its reorganization, leading to a decrease of the local alignment.
The fact that we observe the same correlation for all dynamical networks means that, as long as a network is slightly dynamical, its final alignment is mostly controlled by its number of links per fibre rather than by its remodelling dynamics or its density.On the other hand, non-dynamical networks are locked in mechanically constrained configurations, preventing the system from reorganizing efficiently compared to dynamical ones and leading to a much lower level of alignment.However, we showed that non-dynamical networks still contain some degrees of freedom allowing for spatial matrix reorganization, and that this organization is controlled again by the number of links per fibre in the network but also by the matrix density, which becomes an important factor.Our interpretation is that dense non-dynamical networks are more spatially constrained than sparse networks.Therefore, royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.11: 231456 adding new links to a sparse network can be more beneficial for the networks' overall alignment than to a dense network which has less degrees of freedom.
Altogether, this simple model suggests that different tissue architectures (different levels of fibre alignment) can already emerge as a result of simple interactions between dynamically linked fibres without the need for supplementary complex interactions involving external factors.Such simplified systems highlight the essential role of matrix remodelling on the tissue structuring.The existence of a simple emerging variable such as the number of links per fibre to control tissue structuring could have major therapeutic implications in systems where the architecture of the ECM is impacted (scarring, fibrosis, ageing), but could also prove very useful in the field of tissue engineering.Indeed, because of its simplicity, this variable is amenable to experimental measurements and represents a new putative target for the development of therapeutic drugs one could develop to restore the architecture of various biological tissues after external or internal alterations.It is noteworthy that this variable is not prescribed by model parameters but emerges from the initial simple rules as a combination of ECM remodelling dynamics, linked fibre fraction and fibre spatial organization.
The second major contribution lies in the analysis of the link between this emerging variable and the model parameters related to the crosslinks.Our model reveals that the number of links per fibre in the network, and therefore the global alignment abilities of dynamical fibre networks, results from a complex interplay between their linked fibre fraction and their remodelling speed.From such results, it is apparent that changes in linked fibre fraction will increase or decrease the global alignment abilities of the network, depending on the network remodelling rate.Thus, biological contexts in which fibre crosslinking activity undergoes changes may play an underappreciated role in driving tissue restructuring.Moreover, these results suggest that the development of biological crosslinkers controlling ECM crosslinking as a target for tissue reconstruction must be carefully accounting for ECM remodelling dynamics.
Finally, the third major contribution of the paper lies in the fine time evolution of the spatial structures.This documents the different temporal evolution of the structures as function of the ECM remodelling speeds and reveals an unique trajectory for all architectures combined with internal and transient temporal windows during which they self-organize.The model revealed that dynamical networks composed of long-lasting links exhibited a phase of crosslink accumulation followed by a long 'relaxation' phase (reduction of the number of links per fibre in the network) associated with a spatial reorganization of its fibres, while fast remodelling networks exhibited only the 'relaxation' phase.The long relaxation phase associated with slow realignment of the fibre units observed for slowly remodelling networks is reminiscent of the realignment phase observed on long time scales in later stages of wound healing [47].The crosslink accumulation phase has been observed in different ECM networks, for instance in ageing tissues [9].These new insights into the temporal evolution of the structures as function of the ECM remodelling properties could prove useful in the field of tissue engineering, where there is a need to design efficient biological crosslinkers [49,50].
In emerging systems, the characteristics of the final outcome cannot be predicted from the initial rules of the system and the paths from the initial interactions to the final equilibrium can be numerous and complex, corresponding to a stochastic evolution.This is not completely the case in our model because, if indeed the emerging macrostructures cannot be predicted from the initial rules and the emergence must be understood as a whole, the path is simple and unique and can be strongly predicted by an intermediate emerging variable (the number of links per fibre in the ECM).Altogether, our study suggests that the very aligned structures observed in fibrotic tissues could be mainly due to excess accumulation of crosslinks, consistent with the alterations of ECM structure observed as a consequence of increased crosslinking in lung fibrosis [51] or cancer [8], or again with previous studies on tissue-induced alignment of fibrous ECM [3,52].Such deciphering of the emergence would open numerous perspectives for future investigations.
In this study, several simplifications were made to break the complexity of real ECM systems.For instance, the dynamical remodelling of the fibre network (random linking/unlinking of fibres) can be seen as an indirect way to account for the presence of remodelling cells.This abstract way of looking at cellular activity on the ECM enables us to study independently the effect of matrix remodelling on its structure, instead of pre-imposing some cellular dynamics.By leaving ECM linking/unlinking as free and independent parameters of the model, we are then able to study the respective importance of these parameters and the importance of matrix crosslinking on its architecture.The model not only reveals that remodelling is essential in the production of aligned fibre structures, it also suggests that ECM architecture could be mainly driven by the number of links per fibre in the matrix.Of course, in royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.11: 231456 vivo experiments must be conducted to definitively validate this hypothesis and are out of the scope of this manuscript.On the modelling viewpoint, several perspectives can be considered.First, future works will be devoted to the study of the mechanical properties of these dynamical networks under tensile/ compressive stress and shear and study the viscoelastic properties of the different networks [53].Moreover, it is noteworthy that our model features networks composed of only one type of crosslink (permanent or transient with a given link-life).A natural perspective would be to study the selforganization abilities of networks composed of heterogeneous crosslinks, following the works of [36].Moreover, our network features active crosslinks, i.e crosslinks that generate an alignment of the fibres they are attached to.As a result, our fibre networks are not subject to any external mechanical stimuli.It would be interesting to add cells having the ability to generate locally biophysical cues such as tension, stiffness and fibre production/degradation [54] and study these effects on the structure and mechanical properties of the ECM networks.In the same direction, a complete cell/fibre threedimensional model could account for biomechanical feedback loops such as force-dependent fibre binding/unbinding [30,37].
Finally, we note that our fibre networks are reminiscent of nematic materials [55].A fundamental difference compared to previous studies is the active nematic alignment of the rod-like elements due to the alignment torque at the (dynamical) crosslinks.As a result, contrary to mixtures of passive particles interacting via volume exclusion, the alignment ability of our material depends crucially on the crosslinking dynamics and the network configurations.In this framework, crosslink remodelling can be seen as thermal fluctuations enabling network alignment.Therefore, this simple model could serve as a basis for studying the influence of the environment in the collective motion of isotropic or anisotropic cells (bacteria) [56][57][58].
The outline of this appendix is the following: in appendix A.1, we give the details of the mechanical interaction forces and torques acting on the fibres.In appendix A.2, we give details on the numerical implementation and appendix A.3 details the computation of the closest points of two finite segments used to compute the interactions and the crosslinks.

A.1. Model components
In this section, we give details of the computation of the forces and torques generated by the mechanical interactions described in §2.1 of the main text.
(i) Computation of the fibre-fibre repulsion forces.The force between two spherocylinders k and m is approximated by the force between two spheres of radius R fib , placed along the major axis of the fibre elements at such positions X k,m and X m,k that their distance is minimal (see appendix A.3 for the actual computation of this point).Denoting by F rep k,m the pairwise interaction force between the spherocylinders k and m and using Hertzian theory [39] where α rep is the maximal intensity of the fibre-fibre repulsion and R fib the threshold beyond which the force field vanishes (this can be regarded as the 'width' of the fibre).This force is applied at point X k,m , thus inducing a rotational torque on fibre k.
(ii) Computation of the fibre attachment forces due to crosslinks.Fibres closer than the threshold d max link can create a crosslink, modelled as a linear spring of stiffness α rest and unloaded length d eq link fixed to the two points of the crosslinked fibres that were closest at the time of its creation.Using Hooke's Law, the elastic restoring force sustained by fibre k due to its link with fibre m reads where X l k,m denotes the point of fibre k that was closest to fibre m at the time of the link creation.This force induces a rotational torque on fibre k To ensure coherence between the different features of the model, we require that 2R fib d eq link d max link .(iii) Computation of the linked fibre alignment forces.It is characterized by a stiffness parameter α align > 0 playing the role of a flexural modulus: the larger α align , the more rigid the fibre network.Given two linked fibres k and m, the torque sustained by the fibre k is such that, 8u [ R 3 , where e v m ¼ signðv k Á v m Þ Á v m so that there is no preferential orientation.(iv) Computation of fibre friction.We assume that the friction sustained by an infinitesimal element of a fibre follows a Stokes Law with friction coefficient μ fib [61].The total friction force sustained by a fibre k, computed by integrating this law on the whole length of the fibre, is equal to and the associated rotational torque is equal to For the sake of simplicity and because our model features only one type of element, we assumed a single friction coefficient μ fib for the fibre elements of our model.However, we note that several methods for computing the hydrodynamic properties of rigid macromolecules have been developed in the royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.11: 231456 literature.For instance in [62], the authors present a systematic method for computing the frictional tensor of rigid particles modelled as assemblies of beads.

A.2. Numerical implementation
The differential system (2.1) is numerically solved using a discrete upwind Euler scheme, with adaptive time step.The linking and unlinking Poisson processes are updated between each time step.We assume that a pair of fibres cannot change its linking state more than once in a single time step: this is reasonable if the length of the time-step dt is small enough compared to the mean occurrence time 1/ν of the Poisson process, so we prescribe dt ≤ dt link with dt link ¼ min 0:5 n link , 0:5 The probability for two fibres k and m to develop a crosslink between time t n and time t n+1 = t n + dt n is then given by while the probability for a crosslink to break is given by To ensure that agents do not swap position without even seeing each other, we also restrict the instantaneous translation of each fibre to half its radius R fib and its rotation to arctanð0:1Þ % 6 .This implies the following upper limits for the time step: Reduction of the computational cost is achieved by dividing the domain of simulation into cubes whose side-length is higher than the maximal range of the interactions: thus, interactions need only be computed for pairs of agents located in neighbouring cubes.The loops calculating the interactions are parallelized for further speeding up of the simulations.
One iteration of the Euler scheme proceeds as follows: -Parallel computation of all forces and torques sustained by the agents at time t n (right-hand part of equation (2.1)).-Computation of the adaptive time step (equations (A 8) and (A 11)) -Motion of the agents to their new position -Account for periodic boundary conditions.
-Attribution of each agent to a simulation box.

A.3. Closest points of two finite segments
Given two fibres k and m, we denote by X k,m = X k + l k,m ω k the point of fibre k closest to fibre m (figure 4).
The couple (l k,m , l m,k ) is the minimizer of the distance ‖X royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.11: 231456 2)].If ω k = ω m , there is an infinity of solutions of the form v = u + (X k − X m ) • ω k ; in this case, we arbitrarily chose the solution with the smallest |u| value.Otherwise, there exists a unique solution whose analytical expression is where C a denotes the cut-off function between −a and a.

Appendix B. Quantifiers and visualization tools for the fibre structures
The goal of this section is to define quantifiers allowing to quantitatively describe the local and global organization of the fibre structures obtained with our computational model.Figure 5a shows a typical simulation (almost) at equilibrium, in which fibres are represented as grey double arrows.As one can observe, this simulation shows two levels of organization: a high local alignment and globally twisting, curving patterns located near the centre of the domain.In order to quantitatively describe these states, we now define appropriate numerical quantifiers.

B.1. Local alignment indicator
Let R align denotes the sensing distance up to which a fibre may interact with its neighbours: in our model, it is equal to L fib + 2R fib .For any fibre k, we define its neighbourhood B k as the set of all fibres located at a distance less than R align and its local alignment indicator Al k as the fractional anisotropy of the fibres' directional vectors within B k .It is computed as follows.We denote by p m = ω m ⊗ ω m the projection matrix on the directional vector of fibre m.The mean of the projection matrices of the fibres inside B k is given by where jB k j denotes the number of fibres in B k .royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.11: 231456 The matrix P k is symmetric positive-definite, so its three eigenvalues λ 1 (P k ), λ 2 (P k ) and λ 3 (P k ) are real positive.The alignment indicator or fractional anisotropy in the neighbourhood B k is then equal to the mean of the eigenvalues.Figure 5b shows the same simulation as figure 5a, but here the fibres have been coloured as a function of their local alignment indicator, from blue (Al k = 0) to red (Al k = 1).As one can see, the curved patterns are much easier to distinguish.Thus, the local alignment quantifier also serves as a visualization tool by supporting the qualitative, visual observation of locally organized states.
Note that Al k = 1 if all the fibres in B k have the same directional vector.If the directional vectors are uniformly distributed then theoretically Al k = 0, but this is not always the case.Indeed, the actual sampling of a random distribution may not be fully isotropic, especially if the number of elements in the sample is small.Figure 6 displays the value of the alignment indicator obtained for various distribution of fibres and various sample sizes: it can be seen that a uniform distribution produces alignment indicator ranging from 0.1 (when the sample size is large) to as much as 0.55 (when the sample size is small), and that there is a large discrepancy between different samples.
In our simulations, the number of neighbours of a fibre is very stable: between 20 and 25 for dense systems and between 10 and 15 for sparse systems.Non-dynamical networks display mean alignment indicators between 0.3 and 0.45 for dense systems and between 0.4 and 0.65 for sparse systems: these values are comparable to those observed in our calibration tests for a uniform distribution with similar sample size.
It can be seen from figure 6 that these biases are much smaller for non-isotropic distributions: for mainly two-or one-directional distributions, the values computed are nearly the same regardless of the sample size and the discrepancy between different samples is small.For a two-directional royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.11: 231456 distribution (i.e. when the fibre directional vectors describe a disc), the eigenvalues on the mean projection matrix are theoretically λ 1 (P k ) = λ 2 (P k ) = 1/2 and λ 3 (P k ) = 0, leading to a theoretical alignment indicator of 1= ffiffi ffi 2 p % 0:707.This is very close to the value observed in our calibration tests (see yellow curve on figure 6).Nearly two-directional distributions, where the fibre directional vectors describe a 'band' or thick disc, give lower and lower alignment indicator as the prominence of the third direction (i.e. the band width) increases (see green curves on figure 6).Likewise, conical distributions, which are mainly one-directional, give an alignment indicator close to 1 which becomes lower and lower as the aperture angle of the cone increases (see red curves on figure 6).

B.2. Stereographic projection
The directional vectors of the fibres belong to the half unit sphere S þ 2 .This subset of R 3 can be projected onto the unit disc in two dimensions using a stereographic projection, as explained below.
We define the main direction of a system as the eigenvector associated with the largest eigenvalue of its total projection matrix If the system contains two or three equally represented directions (associated with equal eigenvalues), one of them is randomly selected.
We rotate the set of directional vectors so that this main direction lies on the z-axis or 'north-south axis'.Since the fibres orientation is not relevant in our model, the set of directional vectors can be restricted to the 'north hemisphere' of the sphere.A point ω = (x, y, z) on this hemisphere can then be projected onto the equatorial plane via the following transformation: The whole process is illustrated in figure 7.
Figure 5c shows the stereographic projection of the simulation displayed in figure 5a,b.As one can observe, the dots are not uniformly distributed but densely packed at the centre of the figure, indicating the existence of a main preferential direction in the system.However, not all fibres have a directional vector close to this main direction: a non negligible number of dots are distributed all around the circle, meaning that all possible directions are represented in the system.Furthermore, the presence of a 'circular branch' in the top-right part of the point cloud allows us to identify the locally twisting structure that can be observed in figure 5b: in this part of the system, nearby fibres have royalsocietypublishing.org/journal/rsosR. Soc.Open Sci.11: 231456 similar but gradually shifting directional vectors such that, on the scale of the whole structure, the fibres' directional vectors describe a circle (in the domain S þ 2 ).Thus, this representation enables us to quickly grasp the distribution of the fibres directional vectors around one or more poles.It must be noted that proximity on the stereographic projection indicates similar directional vectors, but not necessarily spatial proximity.Nonetheless, we can gain insights into the overall architecture of the network by drawing the covariance ellipse of the point-cloud (in red dashed line on figure 5c) and computing its semi-major axis length A max .As shown in the §3.2, this enables us to identify many type of 'states' or structures that can also hint on the spatial organization of the network.

8k [ ½½1, N fib , 9
ðtÞ and T rep k,m ðtÞ are the force and torque associated with the repulsion between fibres k and m, F rest k,m ðtÞ and T rest k,m ðtÞ are the force and torque due to the presence of a spring (crosslink) connecting fibres k and m, and T align k,m ðtÞ is the alignment torque generated by this crosslink.We refer to appendix A.1 for the detailed computations of these forces and torques.

1 ν link = 1 ν link = 10 Figure 1 .
Figure 1.(a-c) Illustration of the various structures that can be observed at equilibrium.Fibres are represented by double-headed arrows and coloured according to their local alignment with their neighbours (from blue: Al k = 0 to red: Al k = 1).The structures range from systems with uniformly high local alignment indicator (a) through systems with heterogeneous, intermediate local alignment indicator (b) to disordered systems with uniformly low local alignment indicator (c).(d ) Value of Al mean according to N mean linkperfib at equilibrium, with colour depending on the remodelling speed ν link and horizontal and vertical error-bars indicating the standard deviation N STD linkperfib and Al STD , respectively.The grey dashed-line indicates the critical value of N mean linkperfib and the black dashed lines the three logarithmic fits obtained for N mean linkperfib , N critic .

Figure 3 .
Figure 3. Temporal evolution of dense systems (N fib = 3000) with various linking dynamics.Panels a-a"': Stereographic projection of the system at times 0.5τ Al (a), τ Al (a'), 3τ Al (a") and T final (a"'), for one simulation with ν link = 0.1 and χ link = 0.8.Panels b-b"': Stereographic projection of the system at times 0.5τ Al (b), τ Al (B'), 3τ Al (b") and T final (b"'), for one simulation with ν link = 0.001 and χ link = 0.8.(c) Trajectory in the phase plane A max versus Al sim of individual simulations for slow-remodelling dense networks ν link = 0.001 and various linked fibre fractions χ link .The initial position is indicated with a black square, the final position with a black star and the time-constant τ Al with a black circle.The limits between each class of structures are drawn in dashed lines.(d) Evolution of N mean linkperfib for slow-remodelling dense networks ν link = 0.001 and various linked fibre fractions χ link , with shading indicating the inter-simulation standard deviation N STD linkperfib .The critical value N critic is indicated with a dashed line and the timeconstant τ Al with a black circle.

Figure 4 .
Figure 4. Scheme of two spherocylindrical fibres k and m indicating the position of the closest points X k,m and X m,k of their central segment (in a three-dimensional perspective) relative to their respective centre.

Figure 5 .
Figure 5. Illustration of the various way to visualize the state of a system, using as example the final state of a simulation.(a) Three-dimensional representation of each fibre as a grey double-headed arrow, with edges of the spatial domain V drawn in black.(b) Same representation, with fibres coloured according to their local alignment indicator (blue: Al k = 0, red: Al k = 1).See appendix B.1 for the actual computation.(c) Stereographic projection of the fibres directional vectors.See appendix B.2 for the actual computation.(d ) Stereographic projection of the fibres directional vectors, with the covariance ellipse drawn in red dashed line and its semi-major axis drawn in blue solid line.

Figure 6 .
Figure 6.Calibration of the alignment indicator quantifier Al on random sets of orientation vectors, for various distribution laws and sample sizes.The displayed values correspond to the average and standard deviation over 10 random draws with the same characteristic.

Figure 7 .
Figure 7. Illustration of the stereographic projection.The orientation axes are shown for reference.(a) Natural distribution of the fibres directional vectors on the unit sphere S 2 , with main direction indicated by a red line.(b) Rotation of the vectors set so that its main direction (in red) now lies along the z-axis.The definition-space of the vectors have been reduced to the 'North Hemisphere', that is to the subset S þ 2 in the new rotated coordinates system.The equatorial plane is shown in dark grey.(c) Projection of the vectors onto the equatorial plane, shown in three-dimensional perspective.royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.11: 231456

Table 2 .
Classification of the simulations outcomes into different states based on the local quantifier Al sim and the global quantifier A max .The case fAl sim , 0:7 & A max 0:45g never occurs in our simulations and is thus unnamed.