Pauline Chassonnery

Jenny Paupert

Anne Lorsignol

Childérick Séverac

Marielle Ousset

Pierre Degond

Louis Casteilla

Diane Peurichard
email: diane.a.peurichard@inria.fr

Fiber crosslinking drives the

Keywords: AMS Subject Classification 92-10, 92C10, 82C22, 93A16 Interaction networks, Three dimensional mathematical modelling, Self-organization, Extracellular matrix, Dynamical crosslinking, Architecture emergence. Subjects Biophysics, Biomechanics, Computational Biology

Introduction

The adequate architecture of any organ is mandatory for their efficient physiological function and any changes is associated with function impairment and putative developing dysfunctions and diseases [START_REF] Wohlgemuth | Alignment, cross-linking, and beyond: collagen architect's guide to the skeletal muscle extracellular matrix[END_REF][START_REF] Mays | Age-related changes in collagen synthesis and degradation in rat tissues. importance of degradation of newly synthesized collagen in regulating collagen production[END_REF][START_REF] Piotrowski-Daspit | Dynamics of tissue-induced alignment of fibrous extracellular matrix[END_REF]. The tissue architecture depends mainly on the mechanical forces exerted by the ExtraCellular Matrix (ECM) [START_REF] Theocharis | Extracellular matrix structure[END_REF]. 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 fibers interconnected by molecular bonds, i.e. crosslinks, generating a connected and elastic environment for the surrounding cells [START_REF] Mays | Age-related changes in collagen synthesis and degradation in rat tissues. importance of degradation of newly synthesized collagen in regulating collagen production[END_REF].

The network structure is in constant remodeling, 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 [START_REF] Malandrino | Dynamic filopodial forces induce accumulation, damage, and plastic remodeling of 3d extracellular matrices[END_REF]. Fibrosis and aging are also characterized by an increase of enzymatic and non-enzymatic crosslinks [START_REF] Karamanos | A guide to the composition and functions of the extracellular matrix[END_REF][START_REF] Diller | The role of the extracellular matrix (ecm) in wound healing: A review[END_REF] and this increase in crosslinking prevents ECM degradation by matrix metalloproteinases, both events leading to a decrease of ECM remodeling [START_REF] Levental | Matrix crosslinking forces tumor progression by enhancing integrin signaling[END_REF]. Altogether, these events induce greater stiffness and the arrangement of the collagen fibers becomes less organized and more loose and fragmented, hence weakening tissue integrity and strength [START_REF] Wilson | A microscopic and macroscopic study of aging collagen on its molecular structure, mechanical properties, and cellular response[END_REF][START_REF] Bailey | Mechanisms of maturation and ageing of collagen[END_REF]. 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 tumor [START_REF] Levental | Matrix crosslinking forces tumor progression by enhancing integrin signaling[END_REF].

Because the global architecture of fiber networks seems to be fundamental for controlling tissue functions, modeling the process of ECM structure emergence will greatly improve our understanding of tissue biology and plasticity in physiological or pathological conditions. Numerous models of fiber 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 behavior of each agent (e.g. a fiber element) and its interactions with the surrounding agents over time [START_REF] Drasdo | On selected individual-based approaches to the dynamics in multicellular systems[END_REF][START_REF] Hwang | Rule-based simulation of multi-cellular biological systems -a review of modeling techniques[END_REF]. However, IBM 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 [START_REF] Dallon | Mathematical Modelling of Extracellular Matrix Dynamics using Discrete Cells: Fiber Orientation and Tissue Regeneration[END_REF][START_REF] Didonna | Filamin cross-linked semiflexible networks: Fragility under strain[END_REF][START_REF] Rizvi | Mathematical model of mechanical behavior of micro/nanofibrous materials designed for extracellular matrix substitutes[END_REF][START_REF] Unterberger | Advances in the mechanical modeling of filamentous actin and its cross-linked networks on multiple scales[END_REF][START_REF] Barre | Kinetic theory of particle interactions mediated by dynamical networks[END_REF] 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 fiber networks [START_REF] Lieleg | Structure and dynamics of cross-linked actin networks[END_REF], 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 two-dimensional [START_REF] Head | Distinct regimes of elastic response and deformation modes of cross-linked cytoskeletal and semiflexible polymer networks[END_REF][START_REF] Wilhelm | Elasticity of stiff polymer networks[END_REF][START_REF] Didonna | Filamin cross-linked semiflexible networks: Fragility under strain[END_REF][START_REF] Heussinger | Floppy Modes and Nonaffine Deformations in Random Fiber Networks[END_REF][START_REF] Bai | On the role of the filament length distribution in the mechanics of semiflexible networks[END_REF][START_REF] Shahsavari | Model selection for athermal cross-linked fiber networks[END_REF][START_REF] Bertin | Mesoscopic theory for fluctuating active nematics[END_REF][START_REF] Boissard | Trail formation based on directed pheromone deposition[END_REF][START_REF] Alonso | A Particle Interaction Model for the Simulation of Biological, Cross-Linked Fiber Networks Inspired From flocking Theory[END_REF][START_REF] Unterberger | Advances in the mechanical modeling of filamentous actin and its cross-linked networks on multiple scales[END_REF][START_REF] Peurichard | Simple mechanical cues could explain adipose tissue morphology[END_REF]. Few studies have been conducted on 3D models [START_REF] Stylianopoulos | Volume-averaging theory for the study of the mechanics of collagen networks[END_REF][START_REF] Huisman | Three-Dimensional Cross-Linked F-Actin Networks: Relation between Network Architecture and Mechanical Behavior[END_REF][START_REF] Åström | Strain hardening, avalanches, and strain softening in dense cross-linked actin networks[END_REF][START_REF] Broedersz | Filament-Length-Controlled Elasticity in 3D Fiber Networks[END_REF][START_REF] Harjanto | Modeling Extracellular Matrix Reorganization in 3D Environments[END_REF][START_REF] Lee | A three-dimensional computational model of collagen network mechanics[END_REF][START_REF] Lin | Influence of crosslink density and stiffness on mechanical properties of type i collagen gel[END_REF][START_REF] Picu | Structural evolution and stability of non-crosslinked fiber networks with inter-fiber adhesion[END_REF], although these are expected to yield different, more realistic results than 2D ones since they better mimic biological structures themselves immersed in 3D environments. One of the reasons for fewer 3D 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 3D. However, the latter is becoming less and less of an issue with recent improvements in high resolution 3D imaging and its availability. Among existing 3D models, few of them feature dynamical crosslinking of ECM components. In [START_REF] Åström | Strain hardening, avalanches, and strain softening in dense cross-linked actin networks[END_REF][START_REF] Harjanto | Modeling Extracellular Matrix Reorganization in 3D Environments[END_REF][START_REF] Mak | Impact of crosslink heterogeneity on extracellular matrix mechanics and remodeling[END_REF], various models of 3D fibrous networks composed of permanent or transient crosslinks (remodeling) are proposed. However, most of these models feature ECM remodeling in reaction to external factors (applied load [START_REF] Åström | Strain hardening, avalanches, and strain softening in dense cross-linked actin networks[END_REF][START_REF] Daniel | Bond type and discretization of nonmuscle myosin ii are critical for simulated contractile dynamics[END_REF], migrating cells [START_REF] Harjanto | Modeling Extracellular Matrix Reorganization in 3D Environments[END_REF], contractile cells [START_REF] Mak | Impact of crosslink heterogeneity on extracellular matrix mechanics and remodeling[END_REF]), and the literature so far provides little cues on the mechanisms underlying fiber self-organization.

In the present paper, we test the hypothesis that fiber macrostructures could spontaneously emerge without appealing to contact guidance or external mechanical challenges, as a result of simple mechanical interactions between the fiber elements composing the ECM network. We assess this hypothesis by means of a three dimensional model where ECM fibers are discretized into unit fiber elements, consisting of non-stretching and nonflexible sphero-cylinders with the ability to spontaneously link to and unlink from their close neighbors. This dynamical crosslinking mechanism allows us to model both the overall temporal plasticity of the network and the complex physical properties of biological fibers such as elongation, bending, branching and growth, thus compensating our minimalistic description of the fiber units. The relevance of such discretization was previously validated in the frame of adipose tissue morphogenesis and regeneration in 2D [START_REF] Peurichard | Simple mechanical cues could explain adipose tissue morphology[END_REF][START_REF] Peurichard | Extracellular matrix rigidity may dictate the fate of injury outcome[END_REF].

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 fiber in the network, an emerging variable not directly linked to the model parameters. This simple emerging variable therefore becomes 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 enabling 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 fiber fraction. These results rationalize how even subtle changes in ECM 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 function of the ECM remodelling properties.

2 Models and methods

Description of the model

The 3D ECM is discretized into unit fiber elements consisting of line segments of fixed and uniform length, represented by their centers and directional unit vectors. We consider the following biological and mechanical features: (i) Fiber resistance to pressure: We suppose that fiber elements repel each other at short distances, which models size-exclusion effects. This is achieved via a repulsive force between close fibers based on Hertzian theory [START_REF]Force schemes in simulations of granular materials[END_REF]. This amounts to model fibers 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 fibers. (ii) Fiber elongation and breakage: In addition to carrying a unit of ECM fiber strength, fiber elements also carry a unit of fiber length. However, we provide a way to create longer fibers by allowing two nearby fibers to form a link. A crosslink is modelled as a linear spring with a given spring stiffness, connecting the two closest points of the fiber pair at the time of its creation. There is no prescription of the location of the crosslinks along the body of the fibers they connect. Several consecutively cross-linked fiber elements would model a long and flexible fiber having the ability to adopt complex geometries. Therefore, the cross-linking process models fiber elongation [START_REF] Ciarletta | A finite dissipative theory of temporary interfibrillar bridges in the extracellular matrix of ligaments and tendons[END_REF]. The stiffness constant of the springs α rest controls the possible extension of the long fibers. Symmetrically, pairs of cross-linked fibers can spontaneously unlink, allowing for fiber breakage describing ECM remodelling processes [START_REF] Cox | Remodeling and homeostasis of the extracellular matrix: implications for fibrotic diseases and cancer[END_REF]. Linking and unlinking processes follow Poisson processes with frequencies ν link and ν unlink respectively. As a result, the linked fiber ratio χ link = ν link ν link + ν unlink represents the equilibrium fraction of linked fibers among the pairs of neighbouring fibers. (iii) Crosslink fiber alignment: To model the ability of long fibers (those made of several cross-linked fiber units) to offer a certain resistance to bending, linked fibers are subjected to a potential torque at their junction. This torque vanishes when the fibers are aligned, and consequently acts as a linked-fiber 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 of biological tissues is very small [47], we suppose that inertial forces can be neglected and we consider an over-damped regime for fiber motion and rotation.

Each of the mechanical interactions due to fiber-fiber repulsion (i), fiber-fiber attachment due to crosslinks (ii) and crosslinked fiber-fiber alignment (iii) generate elementary forces and torques between fiber pairs. The total force (resp torque) acting on a fiber is then computed as the sum of all the elementary forces (resp. torques) generated by the elements interacting with this fiber. The motion and rotation of the fiber is then deduced from Newton's equation of motion in an overdamped regime. More specifically, the N fib fiber elements are represented by straight lines of fixed length L fib represented by their centers X k (t) ∈ Ω ⊂ R 3 and their non-oriented directional unit vectors ω k (t) ∈ S + 2 . Moreover, from the fiber-fiber repulsion interaction, fibers may be seen as soft spherocylinders of radius R fib . We denote by p k,m (t) k,m the fiber connectivity matrix, that is p k,m (t) is equal to 1 if fibers k and m are linked at time t and to 0 otherwise.

The motion and rotation of fiber k are then given by:

           µ fib L fib dX k dt (t) = N fib m=1 F rep k,m (t) + p k,m (t)F rest k,m (t) µ fib L 3 fib dω k dt (t) = N fib m=1 T rep k,m (t) + p k,m (t) T rest k,m (t) + T align k,m (t) ∧ ω k (t) ∀k ∈ 1, N fib , (1)
where

F rep k,m (t) and T rep k,m (
t) are the force and torque associated with fiber k and fiber m repulsion, F rest k,m (t) andT rest k,m (t) are the force and torque due to the presence of a spring (crosslink) connecting fiber k and fiber m, and T align k,m (t) is the alignment torque generated by this crosslink. We refer to apendix A.1 for the detailed computations of these forces and torques.

Description of the experiments and biological relevance of the model parameters

The spatial domain Ω is a cuboid of side length L x , L y and L z respectively in the x, y and z-dimension, centered on the origin :

Ω = - L x 2 , L x 2 × - L y 2 , L y 2 × - L z 2 , L z 2 .
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. Fibers are initially randomly inseminated inside the domain according to a uniform law for both position and orientation. The differential system (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 annex 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 fiber 2R fib . This means that the fiber units (spherocylinders) connect with the fibers 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 fibers. (b) The size of the domain is approximately 4 times the size of a fiber 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 fiber aspect-ratio

L fib
2R fib = 6 is quite small compared to the values used in other models of the ECM, which usually varies between 250 and 10 4 [START_REF] Rizvi | Mathematical model of mechanical behavior of micro/nanofibrous materials designed for extracellular matrix substitutes[END_REF][START_REF] Lee | A three-dimensional computational model of collagen network mechanics[END_REF][START_REF] Lin | Influence of crosslink density and stiffness on mechanical properties of type i collagen gel[END_REF]. This compensate the fact that these models directly account for fiber bending and/or fiber elongation, while our long fibers correspond to a sequence of crosslinked fiber units. On the same note, we stress the fact that our fiber units do not aim at modelling the individual collagen fibrils making up collagen fibers in ECM, but rather correspond to an intermediate scale where one fiber unit of our model is already a set of twined collagen fibrils that run in parallel to form a larger bundle [START_REF] Tatsuo | Collagen fibers, reticular fibers and elastic fibers. a comprehensive understanding from a morphological viewpoint[END_REF].

We denote by ϕ fib the fiber density of the network, that is the ratio between the total volume of fibers (without overlapping) and the volume of the spatial domain :

ϕ fib = N fib V fib |Ω| = N fib πR 2 fib L fib + (4/3)πR 3 fib L x L y L z . (2)
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 fibers. In the case of an ordered packing, the packing density of spherocylinders is ϕ order = 0.89, while for random or amorphous packing of sphero-cylinders with an aspect-ratio of 6, the maximal random packing density ϕ random ≈ 0.4 [START_REF] Wilhelm | Elasticity of stiff polymer networks[END_REF]. Thus, we may say that a system is "sparse" if its fiber 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 containing N fib = 3000 fibers (ϕ fib = 0.58) and sparse systems with N fib = 1500 fibers (ϕ fib = 0.29, [START_REF] Yamada | Quantification of extracellular volume fraction by cardiac computed tomography for noninvasive assessment of myocardial fibrosis in hemodialysis patients[END_REF][START_REF] Dedieu | In vivo tissue extracellular volume fraction measurement by dynamic spin-lattice mri relaxometry: application to the characterization of muscle fiber types[END_REF]).

For each of the three types of mechanical forces in the system, we define the "characteristic interaction time" the time needed for two isolated fibers 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 fibers (X 1 = X 2 and ω 1 = ω 2 ) to move apart by 99% of their equilibrium distance 2R fib (i.e. ||X 1 -X 2 || = 0.99 × 2R fib ). Similarly, for the elastic spring T rest is the time needed for two fibers that are initially fully overlapping and crosslinked at their center 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 fibers (X 1 = X 2 and ω 1 ⊥ ω 2 ) crosslinked at their center to reach a relative angle arccos(ω 1 • ω 2 ) = 0.9 • .

Explicit computation leads to the following formula (numerical values are given for the parameters presented in Table 1) :

               T rep = 9µ fib L fib R fib α rep = 8.64 U t , T rest = ln(100) µ fib L fib α rest = 5.53 U t , T align = 4.8 µ fib L 3 fib α align = 523 U t . (3)
It may be noted that the alignment interaction is much slower than the repulsive and elastic restoring forces. In this regime, fiber elements are quite rigid and connected by strong springs (crosslinks), enabling to prevent local accumulation of fibers and overstretching of long fibers (those made of several crosslinked fiber units). 3 Results

Matrix crosslinking drives the local alignment of 3D dynamical fiber networks.

In Figure 1.(A-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 fibers are represented by double arrows, colored as function of their local alignment with their neighbors. We refer the readers to annex B.1 for more details on the computation of this quantifier, and just mention that the local alignment of fiber k, denoted Al k , is equal to 1 (fiber colored in red) if all the neighbouring fibers display the exact same direction as fiber k, and to 0 (fiber colored in blue) if the neighbouring fibers display uniformly distributed directional vectors. Moreover in appendix B.1, we show that this quantifier is able to discriminate between fibers located in randomly oriented environments (corresponding to Al k < 0.5), fibers located in nearly planar environments (leading to Al k around 0.7), and fibers located in nearly uni-directional environments (leading to Al k above 0.8).

As one can observe, the fiber structures obtained at equilibrium range from highly aligned systems (mainly composed of red fibers, see Figure 1.A) to disordered systems with a low local alignment (mainly composed of fibers colored in blue, see Figure 1.C). The model can also produce intermediate states composed of fibers with a median local alignment (see Figure 1.B).

In order to assess the alignment states of our different fiber networks, we computed the mean of the local alignment indicator Al k over all the fibers 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 links per fibers in a network and by N mean linkperfib and N STD linkperfib its average and standard deviation over 10 simulations.

By plotting the alignment quantifier Al mean as a function of the number of links per fiber 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 1.D, with horizontal and vertical error-bars indicating the inter-simulation standard deviations N STD linkperfib and Al STD respectively. The different markers indicate different fiber densities (dots for dense systems and triangles for sparse ones), the different colors refer to different networks dynamics ν link , and inside each color series χ link is increasing with N linkperfib .

Figure 1.D 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 1.D), there is a logarithmic correlation between the number of links per fiber in the network and its mean alignment indicator (black dashed lines in Figure 1.D) :

Al mean ≈ α log(N mean linkperfib ) + β, (4) with
• α = 0.037, β = 1.006 and coefficient of determination r 2 = 0.87 for dynamical systems (nonblue 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 fiber 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 . 1.D 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 colors). Moreover, although we do recover the same type of correlation between the fiber local alignment and the number of links per fiber in the network, for non-dynamical networks this correlation significantly depends on the fiber density. However, the critical number of links N critic allowing for larger alignment is the same for non-dynamical networks, either dense of 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 fiber.

ECM architecture emergence is driven by a complex interplay between remodelling speed and linked fiber fraction

The previous section took a particular focus on the local arrangement of the fiber units composing our 3D fiber network, with little information on the global structures at the population scale. In this section, we aimed to characterize quantitatively the macro structures that emerge in our networks. To this aim, we used the stereographic projection of the fiber directional vectors. Disregarding the spatial position of a fiber, we represented its directional vector as a point on the surface of the unit half-sphere in 3D and then projected it onto the unit disk in 2D (see annex B.2 for a detailed explanation). As shown in Figure . 2, this representation enabled us to characterize the different global organizations of our fiber networks. Indeed, we observed three different types of stereographic projections in our simulations : fibers directional vectors very concentrated around the center of the disk, corresponding to a global alignment of the system (Figure 1.A and Figure 2.A), fibers directional vectors homogeneously distributed on the disk corresponding to a global disorder (Figure 1.C and Figure 2.E), and fibers directional vectors distributed along a preferential axis, with complete depletion in the direction perpendicular to this axis, corresponding to global curved/plane structures (Figure 2

.(B-D)).

Together with the local alignment quantifier Al sim , we were now able to quantitatively characterize the different local and global fiber 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).

A max ⩽ 0. 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 2 A shows the equilibrium values of quantifiers Al sim and A max for dense systems (see annex C.2 for the equivalent figure on sparse systems). The points are colored 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 displayed in Figure 1 are indicated with a black star and their stereographic projection shown as inset. Four other simulation outcomes are singled out with black stars on the phase diagram and illustrated with a 3D view and stereographic projection in the panels B to E bellow.

From Fig. 2.A, 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 (notice the gap between the blue dots and orange diamonds in panel A). Indeed, no simulation displays an average alignment indicator at equilibrium between 0.65 and 0.77 (including sparse systems, see annex C.2), and there is a marked difference between the least organized of the curved states (illustrated in Figure 1.A) and the most organized of the unorganized states (illustrated in Figure 2.E). 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 panels B and C of Figure 2. 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 enabled to distinguish between unorganized networks, globally aligned networks and networks locally aligned with twisting capacities at the population level, enabling 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 fiber density of the network (see Appendix C.2). In 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 2.F shows a heatmap of the percentage of simulations ending in an aligned state (vs a curved state), for dynamical dense networks (see Appendix C.2 for results on sparse networks), depending on the values of the network remodeling speed ν link (in ordinate) and the equilibrium linked fiber fraction χ link (in abscissa). As one can observe in Fig. 2.F, 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 fiber fraction χ link can increase global alignment outputs because, for low-remodelling networks, the formation of crowded interconnected fiber structures inhibiting fiber motion is relieved by reduced link density. Moreover, (ii) the global alignment of networks with modest remodelling rates may undergo little change with reduced linked fiber-fraction and (iii) the global alignment ability of fast-remodelling networks will likely be impaired by reduced linked fiber-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 fiber in the network (as shown in section 3.1), its direct relation with model parameters N fib , ν link and χ link is more complex. Indeed, Figure 2.G shows a heatmap of the number of links per fiber in the network N mean fib as function of ν link and χ link for dense dynamical networks (same simulations as panel F). It demonstrates that N mean fib 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 fiber in the network increases along the diagonals, 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 linkperfib ≈ N critic are framed in green in panel 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 Fig. C.2 in Appendix C.2), confirming the strong correlation between ECM alignment abilities and the number of links per fiber 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 fiber fraction. In the next section, we study the evolution in time of the structures, enabling to give more insights into the role of these parameters in tissue structuring in time.

ECM architecture emergence follows a unique evolutionnary 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 an inverted exponential growth (see annex 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 annex C.1). In Figure 3.A-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 remodeling network ν link = 0.1 (A-A"', Movie3.mp4) and slow remodeling 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 behavior 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 3.C the trajectory in the phase plane A max vs Al sim of simulations for low-dynamical dense networks ν link = 0.001 with various linked fiber fractions χ link (different colors, see 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 fibers directional vector, thus shifting from the initial uniform distribution to a mainly two-directional distribution (see annex 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 3.C 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 fibers 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 stabilisation 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 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 fiber (displayed in Figure 3.D) 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 fiber fraction χ link = 0.9 (Figure 3.D), 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 fiber fraction (χ link = 0.1, blue curves) we do not observe any crosslinks accumulation or fast remodeling networks (see 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 = 10. Parameter χ link describes the proportion of linked fibers among all linkable fibers at equilibrium, but this equilibrium takes time to establish (inversely proportional to ν link ). If the repulsion interaction operates faster than the links remodeling (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 fibers while former overlapping fibers 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).

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 3.D, 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 fiber. A high number of long lasting crosslinks promotes crosslink accumulation resulting in medium/low alignment, while fast remodeling 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 section 3.2, we showed that the network alignment abilities require a number of links adapted to their remodeling speed : fast remodeling networks need a high equilibrium linked fiber fraction to quickly reach a high alignment indicator, while slow remodeling networks need a medium/low equilibrium linked fiber fraction to prevent crosslink accumulation and avoid the formation of crowded interconnected fiber structures inhibiting fiber motion.

Discussion

In this work, we have implemented a 3D model for fiber networks composed of fiber elements capable to dynamically crosslink or unlink each others, to align with each others at the crosslinks and to repel their nearest neighbors to prevent fibers from 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 fiber in the ECM network, (ii) are controlled by a nonlinear relationship between the linked fiber 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 fiber networks in 3D. 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 (highly aligned fiber structures reminiscent of muscular tissues [START_REF] Peter | The structure and role of intramuscular connective tissue in muscle function[END_REF] to disturbed alignment of collagen fibers observed in the first phase of wound healing [START_REF] Xue | Extracellular matrix reorganization during wound healing and its impact on abnormal scarring[END_REF]). 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. In contrast, dynamical networks exhibited a mixture of aligned and curved states. These results point to the essential role of matrix remodeling in ECM structuring, consistent with several results in the literature (see [START_REF] Bonnans | Remodelling the extracellular matrix in development and disease[END_REF] 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 fiber (see section 3.1). Our interpretation is that, when the number of links per fiber is inferior to the critical threshold N critic , the network is weakly constrained. In this configuration, an increase in the number of links per fiber 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 fiber, 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 fiber 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 fiber 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, 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.

The existence of a simple emerging variable such as the number of links per fiber 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 fiber fraction and fiber spatial organization, independently of supplementary complex interactions involving external factors such as migrating cells, contractile forces etc.

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 fiber in the network, and therefore the global alignment abilities of dynamical fiber networks, result from a complex interplay between their linked fiber fraction and their remodelling speed. From such results, it is apparent that changes in linked fiber fraction will increase or decrease the global alignment abilities of the network, depending on the network remodelling rate. Thus, biological contexts in which fiber crosslinking activity undergoes changes may play an underappreciated role in driving tissue restructuring. Moreover, these results suggest that the development of biological crosslinkers controling 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 remodeling speeds and reveals an unique trajectory all architecture 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 fiber in the network) associated with a spatial reorganization of its fibers, while fast remodeling networks exhibited only the "relaxation" phase. The long relaxation phase associated with slow realignment of the fiber units observed for slowly remodelling networks is reminiscent of the realignment phase observed on long time scales in later stages of wound healing [START_REF] Xue | Extracellular matrix reorganization during wound healing and its impact on abnormal scarring[END_REF]. The crosslink accumulation phase has been observed in different ECM networks, for instance in ageing tissues [START_REF] Wilson | A microscopic and macroscopic study of aging collagen on its molecular structure, mechanical properties, and cellular response[END_REF]. 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 [START_REF] Jiang | Cross-linking methods of type i collagen-based scaffolds for cartilage tissue engineering[END_REF][START_REF] Xue | Fabrication of physical and chemical crosslinked hydrogels for bone tissue engineering[END_REF].

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 fiber 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 [START_REF] Philp | Extracellular matrix cross-linking enhances fibroblast growth and protects against matrix proteolysis in lung fibrosis[END_REF] or cancer [START_REF] Levental | Matrix crosslinking forces tumor progression by enhancing integrin signaling[END_REF], or again with previous studies on tissue-induced alignment of fibrous ECM [START_REF] Piotrowski-Daspit | Dynamics of tissue-induced alignment of fibrous extracellular matrix[END_REF][START_REF] Garrison | Fibronectin fibril alignment is established upon initiation of extracellular matrix assembly[END_REF]. Such deciphering of the emergence would open numerous perspectives for future investigations. In vivo experiments must be conducted to definitively validate this hypothesis and are out of the scope of this manuscript.

In this study, we demonstrated the ability of fiber networks to spontaneously self-organize as function of the kinetics of their crosslinks. It is noteworthy that our model features networks composed of only one type of crosslinks (permanent or transient with a given link-life). A natural perspective would be to study the self-organization abilities of networks composed of heterogeneous crosslinks, following the works of [START_REF] Mak | Impact of crosslink heterogeneity on extracellular matrix mechanics and remodeling[END_REF]. Moreover, our network features active crosslinks, i.e crosslinks that generate an alignment of the fibers they are attached to. As a result, our fiber networks are not subject to any external mechanical stimuli. Future works will be devoted to the study of the mechanical properties of these dynamical networks under tensile/compressive stress, shear, etc. Another interesting perspective would be to add cells having the ability to generate locally biophysical cues such as tension, stiffness and fiber production/degradation [START_REF] Winkler | Concepts of extracellular matrix remodelling in tumour progression and metastasis[END_REF] and study these effects on the structure and mechanical properties of the ECM networks.

A Model

Here, we give details about the mathematical model presented in Section 2.1 of the main text. Let us recall the main features and introduce some notations. The N fib fiber elements are represented by their centers X k (t) ∈ R 3 and their non-oriented directional unit vectors ω k (t) ∈ S + 2 . The fibers repel their close neighbors by means of a soft repulsion mechanism modeling steric repulsion between spherocylinders of radius R fib . Fiber elements have the ability to link to or unlink from each other to model fiber elongation or rupture. The linking and unlinking of fibers follow random (Poisson) processes in time. Fibers offer resistance to bending through an alignment torque acting between two linked fiber elements. Fiber motion and rotation is given by Newton's second law of motion, in an overdamped regime to model a medium with low Reynolds number.

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 fibers. 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 section 2.1 of the main text.

Computation of the fiber-fiber repulsion forces (i) 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 fiber elements at such positions X k,m and X m,k that their distance is minimal (see annex A.3 for the actual computation of this point). Denoting by F rep k,m the pairwise interaction force between the spherocylinders f and m and using Hertzian theory [START_REF]Force schemes in simulations of granular materials[END_REF] :

F rep k,m = α rep (2R fib -||X k,m -X m,k ||) 3/2 2R fib × X k,m -X m,k ||X k,m -X m,k || , ( 5
)
where α rep the maximal intensity of fiber-fiber repulsion and R fib the threshold beyond which the force field vanishes (this can be regarded as the "width" of the fiber). This force is applied at point X k,m , thus inducing a rotational torque :

T rep k,m = (X k,m -X k ) ∧ F rep k,m , (6)
on fiber k.

Computation of the fiber attachment forces due to crosslinks (ii) Fibers closer than the threshold d max link can create a crosslink, modeled as a linear spring of stiffness α rest and unloaded length d eq link fixed to the two points of the crosslinked fibers that were closest at the time of its creation. Using Hooke's law, the elastic restoring force sustained by fiber k due to its link with fiber m reads

F rest k,m = α rest d eq link -||X l k,m -X l m,k || X l k,m -X l m,k ||X l k,m -X l m,k || , (7)
where X l k,m denotes the point of fiber k that was closest to fiber m at the time of the link creation. This force induces a rotational torque on fiber k :

T rest k,m = X l k,m -X k ∧ F rest k,m . (8)
To ensure coherence between the different features of the model, we require that 2R fib ≤ d eq link ≤ d max link . Computation of the linked fiber alignment forces (iii). It is characterized by a stiffness parameter α align > 0 playing the role of a flexural modulus : the larger α align , the more rigid the fiber network. Given two linked fibers k and m, the torque sustained by the fiber k is such that, ∀u ∈ R 3 ,

T align k,m ∧ u = α align (ω k ∧ ω m ) ∧ u + 1 -|ω k • ω m | ||ω k ∧ ω m || 2 (ω k ∧ ω m ) ∧ ((ω k ∧ ω m ) ∧ u) , (9)
where

ω m = sign(ω k • ω m )
• ω m so that there is no preferential orientation. Computation of fiber friction (iv) We assume that the friction sustained by an infinitesimal element of a fiber follows a Stokes law with friction coefficient µ fib [START_REF] Leith | Drag on nonspherical objects[END_REF]. The total friction force sustained by a fiber k, computed by integrating this law on the whole length of the fiber, is equal to :

F fric k = -µ fib L fib dX k dt (10)
and the associated rotational torque is equal to :

T fric k = -µ fib L 3 fib ω k ∧ dω k dt . (11)
A.2 Numerical implementation

The differential system 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 fibers 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 ν link , 0.5 ν unlink . (12)
The probability for two fibers k and m to develop a crosslink between time t n and time t n+1 = t n + dt n is then given by :

P p k,m (t n+1 ) = 1 p k,m (t n ) = 0 and ||X k,m (t n ) -X m,k (t n )|| ≤ d max link = 1 -e -ν link dtn (13)
while the probability for a crosslink to break is given by :

P p k,m (t n+1 ) = 0 p k,m (t n ) = 1 = 1 -e -ν unlink dtn (14)
To ensure that agents do not swap position without even seeing each other, we also restrict the instantaneous translation of each fiber to half its radius R fib and its rotation to arctan(0.1) ≈ 6 • . This implies the following upper limits for the time-step :

               dt trans (t n ) = min 1≤k≤N fib   0.5 R fib dX k dt (t n )   , dt rot (t n ) = min 1≤k≤N fib   0.1 dω k dt (t n )   . (15)
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 follow :

• Parallel computation of all forces and torques sustained by the agents at time t n (right-hand part of equation ( 1)).

• Computation of the adaptive time-step (equations ( 12) and ( 15))

dt n = min(dt trans (t n ), dt rot (t n ), dt link ).
• Motion of the agents to their new position :

X k (t n+1 ) = X k (t n ) + dt n dX k dt (t n ) ω k (t n+1 ) = ω k (t n ) + dt n dω k dt (t n )
• Account for periodic boundary conditions.

• Attribution of each agent to a simulation box.

• Parallel update of linking configuration (equations ( 13) and ( 14)).

A.3 Closest points of two finite segments

Given two fibers k and m, we denote by X k,m = X k + l k,m ω k the point of fiber k closest to fiber m (see Figure 4). The couple (l k,m , l m,k ) is the minimizer of the distance

||X k + uω k -(X m + vω m )|| for (u, v) ∈ -L fib 2 , L fib 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 :

     l k,m = C L fib 2 ((ω k • ω m )ω m • (X k -X m ) -ω k • (X k -X m )) 1 -(ω k • ω m ) 2 , l m,k = C L fib 2 ((ω k • ω m )ω k • (X m -X k ) -ω m • (X m -X k )) 1 -(ω k • ω m ) 2 , (16)
where C a denotes the cut-off function between -a and a.

B Quantifiers and visualization tools for the fiber structures

The goal of this section is to define quantifiers allowing to quantitatively describe the local and global organization of the fiber structures obtained with our computational model. Figure 5.A shows a typical simulation (almost) at equilibrium, in which fibers are represented as gray 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 center 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 fiber may interact with its neighbours : in our model, it is equal to L fib + 2R fib . For any fiber k, we define its neighbourhood B k as the set of all fibers located at a distance less than R align and its local alignment indicator Al k as the fractional anisotropy of the fibers 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 fiber m. The mean of the projection matrices of the fibers inside B k is given by

P k = 1 |B k | m s.t. Xm∈B k p m , (17)
where |B k | denotes the number of fibers in B k .

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 :

Al k = 3 2 (λ 1 (P k ) -λ) 2 + (λ 2 (P k ) -λ) 2 + (λ 3 (P k ) -λ) 2 λ 1 (P k ) 2 + λ 2 (P k ) 2 + λ 3 (P k ) 2 (18)
with λ = (λ 1 (P k ) + λ 2 (P k ) + λ 3 (P k ))/3 the mean of the eigenvalues. Figure 5.B shows the same simulation as Figure 5.A, but here the fibers have been colored 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 fibers 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 fibers 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 fiber 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 distribution (i.e. when the fiber directional vectors describe a disk), 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/ √ 2 ≈ 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 fiber directional vectors describe a "band" or thick disk, 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 increase (see red curves on Figure 6).

B.2 Stereographic projection

The directional vectors of the fibers belong to the half unit sphere S + 2 . This subset of R 3 can be projected onto the unit disk in 2D using a stereographic projection, as explained below.

We define the main direction of a system as the eigenvector associated to the largest eigenvalue of its total projection matrix

P tot = 1 N fib 1≤k≤N fib ω k ⊗ ω k . (19)
If the system contains two or three equally represented directions (associated to 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 fibers 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 :

p(ω) = x 1 + z , y 1 + z . ( 20
)
The whole process is illustrated in Figure 7.

Figure 5.C shows the stereographic projection of the simulation displayed in Figure 5.A and B. As one can observe, the dots are not uniformly distributed but densely packed at the center of the figure, indicating the existence of a main preferential direction in the system. However, not all fibers 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 to identify the locally twisting structure that can be observed in Figure 5.B : in this part of the system, nearby fibers have similar but gradually shifting directional vectors such that, on the scale of the whole structure, the fibers directional vectors describe a circle (in the domain S + 2 ). Thus, this representation enables us to quickly grasp the distribution of the fibers 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 5.C) and computing its semi-major axis length A max . As shown in the section 3.2, this enables us to identify many type of "states" or structures that can also hint on the spatial organization of the network.

C Supplementary data C.1 Temporal evolution of the spatial structures : Movies

We describe here the movies showing the results of some of our simulations, available online at this address.

Each video is divided in three panels. On the left is a 3D representation of the system with fibers colored according to their alignment indicator (see colorbar on the right) and edges of the spatial domain drawn in black, in the middle the stereographic projection of the fibers directional vectors and on the right the trajectory of the simulation in the plane A max -Al sim . The current time (in U t ) is displayed at the top.

(Movie1) Simulation of a dense system (L fib = 3000) with fast remodeling dynamics (ν link = 0.1) and low equilibrium linked fiber fraction (χ link = 0.2).

This video shows a system quickly organizing : at t = 1000U t , the system has already transitioned from its initial unorganized state to a curved state. At t = 3000U t , the main direction can be seen emerging in the form of a large cluster of points in the stereographic projection. At t = 10 000U t , the stereographic projection displays a planar distribution of the directional vectors, with extra accretion of points in the main direction and total depletion in the perpendicular direction. The system has already nearly reached its maximal alignment indicator and, from that point onward, it will mainly undergo small local adjustments of the fibers position and orientation (see the 3D representation on the left panel). The alignment indicators of individual fibers harmonise, the mean alignment indicator increases slightly and the point cloud of the stereographic projection condensates in a clear straight band. During the entire simulation, the semi-major axis length A max of the stereographic projection covariance ellipse stays nearly constant.

(Movie2) Simulation of a dense system (L fib = 3000) with slow remodeling dynamics (ν link = 0.001) and low equilibrium linked fiber fraction (χ link = 0.2).

This video shows a system organizing more slowly than the previous one (approximately twice slower) but achieving a more aligned final state. The system reaches a curved state at t = 1900U t . The main direction can be seen emerging on the stereographic projection around time t ≈ 5000U t . The point cloud of the stereographic projection then begins to condensate around this main direction in a nearly symmetric manner while the various local structures rotate to align together (see left panel), reaching an aligned state at t = 23 000U t and continuing to align.

(Movie3) Simulation of a dense system (L fib = 3000) with fast remodeling dynamics (ν link = 0.1) and high equilibrium linked fiber fraction (χ link = 0.8).

This video shows a system organizing very quickly, with a stereographic projection adopting as early as t = 4000U t a band-like pattern which quickly gets thinner. At t = 6000U t , the 3D representation shows a clear wavy pattern with very uniform local alignment indicators. At that time the mean align-ment indicator is already high (> 0.9). The stereographic projection then begins to contract while the wavy pattern flatten, and the simulation ends in an aligned state.

(Movie4) Simulation of a dense system (L fib = 3000) with slow remodeling dynamics (ν link = 0.001) and high equilibrium linked fiber fraction (χ link = 0.8).

This video shows a system evolving very slowly. The mean alignment indicator reaches the 0.5 threshold around t = 11 000U t . At that time, the local alignment indicator of individual fibers displays wide discrepancies and the stereographic projection point cloud has not visibly changed. A main direction can be seen emerging at approximately t = 20 000U t , but the central cluster of points is very large and does not contract over time, as can be seen by the fact that the quantifier A max does nearly not decrease. The system ends in a curved state with heterogeneous local structures.

(Movie5) Comparison of two simulations with different fiber densities, both ending in an aligned state. The top row shows a dense system (L fib = 3000) with intermediate remodeling dynamic (ν link = 0.01) and moderate equilibrium linked fiber fraction (χ link = 0.3). The bottom row displays a sparse system (L fib = 1500) with intermediate remodeling dynamic (ν link = 0.01) and very high equilibrium linked fiber fraction (χ link = 0.9).

It is noteworthy that the two systems display a very similar temporal evolution. This comes from the fact that they have the same remodeling speed ν link and a comparable number of links per fiber N linkperfib . The latter is achieved by giving the sparse system a higher equilibrium linked fiber fraction χ link , which compensate for its lesser number of linkable configurations (i.e. overlapping fiber pairs).

C.2 Snapshots of sparse systems

In this section, we take a closer look at the spatial organization of sparse systems. Figure 8.A compares the values of quantifiers Al sim and A max when the simulation has reached equilibrium, with color depending on the type of state reached (blue dots correspond to unorganized states, orange diamonds to curved states and red crosses to aligned states). A few of simulations corresponding to either typical or borderline cases are singled out with black stars and illustrated with a 3D view and stereographic projection in the panels B to I.

We first observe that the group of unorganized states (blue dots) is less compact than it was for dense systems and reaches greater values of Al sim . The groups of curved states (orange diamonds) and aligned states (red crosses) have the same characteristics in term of Al sim and A max than before, but the first one is much more populated and the second much less (it only contains 10 simulations).

The most aligned state observed in sparse systems (panel B) is less straight than the typical aligned state for dense systems. Typical curved states (panels E and F) and unorganized states (panel I) however are very comparable to what was observed in dense systems. The transition between aligned and curved states is still continuous (compare panels C and D) and the transition between curved and unorganized states sharp (compare panels G and H), though the gap (in term of Al sim ) and the visual difference are lesser.

C.3 Correlation between the links life-expectancy and the ECM architecture

Here, we explore whether the network organization abilities could be controlled by the life expectancy of a link, which depends of both ν link and χ link via the following relation :

T link-life = 1 ν unlink = χ link (1 -χ link )ν link . (21)
Figure 9 displays the value of Al mean at equilibrium as a function of T link-life , with a vertical errorbar indicating the inter-simulation standard deviation Al STD . The value of ν link is indicated in color and, inside each color series, χ link is increasing with T link-life . The characteristic time of the alignment interaction T align (see section 2.2) is indicated for the sake of comparison. We observe that, in the case of dense systems (left panel), the Al mean a flat maximum for T link-life ∈ [10, 500] U t , while for sparse system (right panel) it reaches its highest value at T link-life ≈ 500 U t .

This can be explained by the fact that, when the average life expectancy of a link T link-life is very small compared to the characteristic time of the alignment force T align = 523 U t , the links do not persist long enough to fully exert their aligning influence and the equilibrium alignment indicator is lesser. This is especially true for sparse systems, which display a clear drop for T link-life < 500 U t . For dense systems the drop is slower and less pronounced.

On the other hand, when T link-life is large compared to T align , on average the links last longer than necessary to wield their full effect and lock the system in non-optimal configurations by obstructing the action of other links. Though these locally locked structures will disappear over time, others will appear -or, to put it another way, the transmission of information (i.e. the fiber direction) in the network is too slow for all the agents to synchronize and the system will not be able to reach an extremely aligned equilibrium state. Our main observation is that, for all parameters, the Al mean follows an inverted exponential growth, that is a quick initial growth followed by a slow convergence towards an asymptotic value. We computed the time-constant τ Al of this growth, that is the time needed to reach 63% of the asymptotic value, and plotted it on the corresponding curve with a black circle. It can be seen that, for a given value of ν link , the shorter the time-constant, the higher the equilibrium value of the alignment indicator (compare the curves inside each panel). By comparing the panels from left to right, we see that the faster the remodeling of the network, the faster the convergence of the system towards its equilibrium value. Moreover, by comparing the extreme cases χ link = 0.1 (blue curve) with χ link = 0.9 (pink curves) of panels A and C, we see that the dependence of the reorganization time τ Al on the equilibrium linked fiber fraction is not trivial. Indeed, fast remodeling networks (panel C) seem to reorganize faster when the equilibrium linked fiber fraction is large (pink curve) than low (blue curve), while the reverse is observed for slow remodeling networks (panel A). Altogether, these results suggest that for each network dynamics, there exists a most efficient range of equilibrium linked fiber fraction allowing for quicker convergence to equilibrium.

C.4 Temporal evolution of the quantifiers

To explore in more details the dependence between the convergence speed and the parameters of We can first see on the left panel of Figure 11 that τ Al decreases when ν link increases according to a non-linear relationship which saturates for ν link ≥ 0.1 (compare the different color points). These results show that fast remodeling networks relax faster to their steady-states than slow-dynamical networks. Moreover, sparse systems organize quicker than dense systems at low linking dynamics (ν link ≤ 0.01, compare the dot and triangle markers for the green and yellow populations), while there is no difference between dense and sparse systems for fast remodeling networks (ν link ≥ 0.1 where dot and triangle markers are superimposed).

For each value of ν link , there is a most efficient range of equilibrium linked fiber fraction χ link allowing for a lower value of τ Al and so a quicker convergence to equilibrium. For slow remodeling networks (ν link = 0.001, green markers) this range lays between χ link = 0.2 and χ link = 0.3, because systems with too much crosslinks will undergo stiffening and take longer to relax, but systems with too few crosslinks will have difficulty to align themselves. As one can observe, the range of χ link allowing the fastest convergence to equilibrium shifts towards 1 as the network remodeling speed ν link increases. As the network remodeling increases, a greater number of crosslinks will then promote a quicker alignment.

When looking only at parameter sets which, on average, lead to aligned equilibrium states (right panel of Figure 11), we can see that these parameter sets cover all remodeling dynamics and correspond to the range of equilibrium linked fiber fraction leading to fastest convergence for each remodeling speed. We conclude that the most efficient systems (which organize the fastest) are also those that align most.

Figure 1 :

 1 Figure 1: Panels A-C: Illustration of the various structures that can be observed at equilibrium. Fibers are represented by double-headed arrows and colored 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 (panel A) through systems with heterogeneous, intermediate local alignment indicator (panel B) to disordered systems with uniformly low local alignment indicator (panel C). Panel D: Value of Al mean according to N mean linkperfib at equilibrium, with color depending on the remodelling speed ν link and horizontal and vertical error-bars indicating the standard deviation N STD linkperfib and Al STD respectively. The gray 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 2 :

 2 Figure 2: Panel A: Alignment indicator Al sim versus semi-major axis length of the covariance ellipse of the stereographic projection A max , for each simulation of dense systems (N fib = 3000). Red crosses correspond to systems in an aligned state, orange diamonds to curved states and blue dots to unorganized states. The simulations previously displayed in Figure 1 are indicated with a black star and their stereographic projection given as inset. Panel B to E display the equilibrium state of a few other simulations, whose position on the diagram are also indicated with a black star. Panel F: Heatmap of the percentage of simulations ending in an aligned state (vs a curved state), for dynamical dense networks as function of the values of the network remodeling speed ν link (in ordinate) and the equilibrium linked fiber fraction χ link (in abscissa).

Figure 3 :

 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. Panel C : Trajectory in the phase plane A max vs Al sim of individual simulations for slowremodelling dense networks ν link = 0.001 and various linked fiber 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. Panel D : Evolution of N mean linkperfib for slow-remodelling dense networks ν link = 0.001 and various linked fiber fractions χ link , for 10 simulations with shading indicating the inter-simulation standard deviation N STD linkperfib . The critical value N critic is indicated with a dashed line and the time-constant τ Al with a black circle.

Figure 4 :

 4 Figure 4: Scheme of two sphero-cylindrical fibers k and m indicating the position of the closest points X k,m and X m,k of their central segment (in a 3D perspective) relative to their respective center.

Figure 5 :

 5 Figure 5: Illustration of the various way to visualize the state of a system, using as example the final state of a simulation. Panel A : 3D representation of each fiber as a gray double-headed arrow, with edges of the spatial domain Ω drawn in black. Panel B : Same representation, with fibers colored according to their local alignment indicator (blue : Al k = 0, red : Al k = 1). See annex B for the actual computation. Panel C : Stereographic projection of the fibers directional vectors. See annex B.2 for the actual computation. Panel D : Stereographic projection of the fibers directional vectors, with the covariance ellipse drawn in red dashed line and its semi-major axis drawn in blue solid line.

Figure 6 :

 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 :

 7 Figure 7: Illustration of the stereographic projection. The orientation axis are shown for reference. Panel A : Natural distribution of the fibers directional vectors on the unit sphere S 2 , with main direction indicated by a red line. Panel 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. Panel C : Projection of the vectors onto the equatorial plane, shown in 3D perspective.

Figure 8 :

 8 Figure 8: Panel A: Alignment indicator Al sim versus semi-major axis length of the covariance ellipse of the stereographic projection A max , for each simulation of a sparse system. Red crosses correspond to systems in an aligned state, orange diamonds to curved states and blue dots to unorganized states. Panels B to I display the equilibrium state of a few simulations (with 3D view and stereographic projection) to illustrate typical or borderline cases. Their position on the diagram are indicated with a black star.

Figure 9 :

 9 Figure 9: Value of Al mean at equilibrium according to the value of T link-life , with color depending on the remodeling speed ν link and vertical error-bars indicating the inter-simulation standard deviation Al STD . The characteristic time of the alignment interaction T align is indicated with gray dashed-lines for the sake of comparison.

Figure 10 :

 10 Figure 10: Temporal evolution the quantifier Al mean for dense systems (N fib = 3000) with various linking dynamics. Shading indicate the inter-simulation standard deviation Al STD . The time-constant of this growth is indicated with a black circle and the limit between unorganized and curved or aligned states is drawn with a dashed line.

Figure 10

 10 Figure10displays the temporal evolution of the quantifier Al mean for dense systems with various values of ν link and χ link .Our main observation is that, for all parameters, the Al mean follows an inverted exponential growth, that is a quick initial growth followed by a slow convergence towards an asymptotic value. We computed the time-constant τ Al of this growth, that is the time needed to reach 63% of the asymptotic value, and plotted it on the corresponding curve with a black circle. It can be seen that, for a given value of ν link , the shorter the time-constant, the higher the equilibrium value of the alignment indicator (compare the curves inside each panel). By comparing the panels from left to right, we see that the faster the remodeling of the network, the faster the convergence of the system towards its equilibrium value. Moreover, by comparing the extreme cases χ link = 0.1 (blue curve) with χ link = 0.9 (pink curves) of panels A and C, we see that the dependence of the reorganization time τ Al on the equilibrium linked fiber fraction is not trivial. Indeed, fast remodeling networks (panel C) seem to reorganize faster when the equilibrium linked fiber fraction is large (pink curve) than low (blue curve), while the reverse is observed for slow remodeling networks (panel A). Altogether, these results suggest that for each network dynamics, there exists a most efficient range of equilibrium linked fiber fraction allowing for quicker convergence to equilibrium.To explore in more details the dependence between the convergence speed and the parameters of

Figure 11 :

 11 Figure 11: Time-constant τ Al of the average alignment indicator (Al mean ) according to the equilibrium linked fiber fraction χ link , with color depending on the remodeling speed ν link and marker depending on the fiber density. Left : Results for all tested sets of parameters, whatever the outcomes of the simulations. Right : Results for all tested sets of parameters which lead, on average, to an aligned state.

Figure 12 :

 12 Figure 12: Temporal evolution of dense systems (N fib = 3000) with medium or high linking dynamics. Panels C-C' : Trajectory in the phase plane A max vs Al sim of individual simulations for mediumremodelling (ν link = 0.01) or quick-remodelling (ν link = 0.1) dense networks and various linked fiber 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. Panels D-D' : Evolution of N mean linkperfib for slow-remodelling dense networks ν link = 0.001 and various linked fiber fractions χ link , for 10 simulations with shading indicating the inter-simulation standard deviation N STD linkperfib . The critical value N critic is indicated with a dashed line and the time-constant τ Al with a black circle.

Table 1 :

 1 Model parameters.

	Name	Value	Units	Description
				Agents
	N fib	[1500, 3000]	N/A	Number of fibers
	L fib	6	L	Fiber length
	R fib	0.5	L	Fiber radius
			Mechanical interactions
	α rep	12.5	M.L -1 .T -2 Magnitude of the repulsion force
	α rest	5.0	M.T -2	Magnitude of the elastic restoring force
	α align	2.0	M.L 2 .T -2 Magnitude of the alignment torque
	d max link d eq link	1.0 1.0	L L	Perception distance for link creation Link equilibrium length
			Biological phenomena
	ν link	[0, 10]	T -1	Network remodeling speed
	χ link	[0.1, 0.9]	N/A	Equilibrium linked fiber fraction
			Numerical parameters
	L x = L y = L z	30	L	Side length of the cubic domain
	T final	5.10 4	T	Total time of simulation

Table 2 :

 2 Classification of the simulations outcomes into different states based on the local quantifier Al sim and the global quantifier A max . The case {Al sim < 0.7 & A max ⩽ 0.45} never occurs in our simulations and is thus unnamed.

		45	> 0.45
	⩾ 0.7	Aligned state : alignment both local and global.	Curved state : alignment local but not global.
	Al sim

< 0.7 (alignment global but not local) Unorganized state : no alignment, either local or global.

Acknowledgment : P. Degond holds a visiting professor association with the Department of Mathematics, Imperial College London, UK.

Ethics : This article does not present research with ethical considerations.

Data Access :

The code used to perform numerical simulations of our model can be found on GitHub. Supplementary data are also available online.

Competing interests :

The authors declare no competing interests.

Funding : This study has been partially supported through the grant EUR CARe N°ANR-18-EURE-0003 in the framework of the Programme des Investissements d'Avenir, by Sorbonne Alliance University with an Emergence project MATHREGEN, grant number S29-05Z101 and by Agence Nationale de la Recherche (ANR) under the project grant number ANR-22-CE45-0024-01.