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

Pauline Chassonnery, Jenny Paupert, Anne Lorsignol, Childérick Séverac,

Marielle Ousset, Pierre Degond, Louis Casteilla, Diane Peurichard

► To cite this version:

Pauline Chassonnery, Jenny Paupert, Anne Lorsignol, Childérick Séverac, Marielle Ousset, et al. Fibre crosslinking drives the emergence of order in a three-dimensional dynamical network model. Royal Society Open Science, 2024, 11 (1), pp.231456. 10.1098/rsos.231456. hal-04618698

HAL Id: hal-04618698

https://hal.science/hal-04618698

Submitted on 20 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
ROYAL SOCIETY
OPEN SCIENCE

royalsocietypublishing.org/journal/rsos

<figure><img src="image_1.png" /></figure>

Research

<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

Check for
updates

**Cite this article:** Chassonnery P, Paupert J, Lorsignol A, Séverac C, Ousset M, Degond P, Casteilla L, Peurichard D. 2024 Fibre crosslinking drives the emergence of order in a three-dimensional dynamical network model. *R. Soc. Open Sci.* **11**: 231456.

https://doi.org/10.1098/rsos.231456

Received: 7 October 2023

Accepted: 9 January 2024

**Subject Category:**

Mathematics

**Subject Areas:**

biophysics/biomechanics/computational biology

**Keywords:**

interaction networks, three-dimensional
mathematical modelling, self-organization,
extracellular matrix, dynamical crosslinking,
architecture emergence

**Author for correspondence:**

Diane Peurichard
e-mail: diane.a.peurichard@inria.fr

Electronic supplementary material is available
online at https://doi.org/10.6084/m9 figshare.c
7054386.

THE ROYAL SOCIETY
PUBLISHING

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

Pauline Chassonnery<fnref n="1" />,<fnref n="2" />, Jenny Paupert<fnref n="1" />,
Anne Lorsignol<fnref n="1" />, Childérick Séverac<fnref n="1" />, Marielle Ousset<fnref n="1" />,
Pierre Degond<fnref n="3" />, Louis Casteilla<fnref n="1" />,* and Diane Peurichard<fnref n="2" />

<figure><img src="image_4.png" /></figure>

PC, 0009-0002-5751-1533; DP, 0000-0002-0807-2266

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.

<footnote>*co-last author.</footnote>

<footnote n="1">¹RESTORE, Université de Toulouse, Inserm U1031, EFS, INP-ENVT, UPS, CNRS ERL5311, Toulouse, France</footnote>

<footnote n="2">²Inria Paris, team MAMBA, Sorbonne Université, CNRS, Université de Paris, Laboratoire Jacques-Louis Lions UMR7598, 75005 Paris, France</footnote>

<footnote n="3">³Institut de Mathématiques de Toulouse, UMR5219, Université de Toulouse, CNRS, UPS, 31062 Toulouse Cedex 9, France</footnote>
<figure><img src="image_2.png" /></figure>

Research

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

Check for
updates

**Cite this article:** Chassonnery P, Paupert J, Lorsignol A, Séverac C, Ousset M, Degond P, Casteilla L, Peurichard D. 2024 Fibre crosslinking drives the emergence of order in a three-dimensional dynamical network model. *R. Soc. Open Sci.* **11**: 231456.

https://doi.org/10.1098/rsos.231456

Received: 7 October 2023

Accepted: 9 January 2024

**Subject Category:**

Mathematics

**Subject Areas:**

biophysics/biomechanics/computational biology

**Keywords:**

interaction networks, three-dimensional
mathematical modelling, self-organization,
extracellular matrix, dynamical crosslinking,
architecture emergence

**Author for correspondence:**

Diane Peurichard
e-mail: diane.a.peurichard@inria.fr

Electronic supplementary material is available
online at https://doi.org/10.6084/m9 figshare.c
7054386.

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

Pauline Chassonnery<fnref n="1" />,<fnref n="2" />, Jenny Paupert<fnref n="1" />,
Anne Lorsignol<fnref n="1" />, Childérick Séverac<fnref n="1" />, Marielle Ousset<fnref n="1" />,
Pierre Degond<fnref n="3" />, Louis Casteilla<fnref n="1" />,* and Diane Peurichard<fnref n="2" />

PC, 0009-0002-5751-1533; DP, 0000-0002-0807-2266

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.

<footnote>*co-last author.</footnote>

<footnote n="1">¹RESTORE, Université de Toulouse, Inserm U1031, EFS, INP-ENVT, UPS, CNRS ERL5311, Toulouse, France</footnote>

<footnote n="2">²Inria Paris, team MAMBA, Sorbonne Université, CNRS, Université de Paris, Laboratoire Jacques-Louis Lions UMR7598, 75005 Paris, France</footnote>

<footnote n="3">³Institut de Mathématiques de Toulouse, UMR5219, Université de Toulouse, CNRS, UPS, 31062 Toulouse Cedex 9, France</footnote>
# 1. 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–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–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 two-dimensional [14,16,19–27]. Few studies have been conducted on three-dimensional models [28–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
# 1. 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–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–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 two-dimensional [14,16,19–27]. Few studies have been conducted on three-dimensional models [28–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
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.

## 2. Models and methods

### 2.1. 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 $\alpha_{\text{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 $\alpha_{\text{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 $v_{\text{link}}$ and $v_{\text{unlink}}$, respectively. As a result, the linked fibre ratio $\chi_{\text{link}} = v_{\text{link}}/(v_{\text{link}} + v_{\text{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 $\alpha_{\text{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
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.

## 2. Models and methods

### 2.1. 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 $\alpha_{\text{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 $\alpha_{\text{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 $v_{\text{link}}$ and $v_{\text{unlink}}$, respectively. As a result, the linked fibre ratio $\chi_{\text{link}} = v_{\text{link}}/(v_{\text{link}} + v_{\text{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 $\alpha_{\text{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
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 $\mathbf{X}_k(t) \in \Omega \subset \mathbb{R}^3$ and their non-oriented directional unit vectors $\omega_k(t) \in \mathbb{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

$$ 
\left. 
\begin{aligned}
\mu_{\text{fib}} L_{\text{fib}} \frac{d\mathbf{X}_k}{dt}(t) &= \sum_{m=1}^{N_{\text{fib}}} (\mathbf{F}_{\mathbf{k},m}^{\text{rep}}(t) + p_{k,m}(t) \mathbf{F}_{\mathbf{k},m}^{\text{rest}}(t)) \\
\mu_{\text{fib}} L_{\text{fib}}^3 \frac{d\omega_k}{dt}(t) &= \sum_{m=1}^{N_{\text{fib}}} \left( \mathbf{T}_{\mathbf{k},m}^{\text{rep}}(t) + p_{k,m}(t) (\mathbf{T}_{\mathbf{k},m}^{\text{rest}}(t) + \mathbf{T}_{\mathbf{k},m}^{\text{align}}(t)) \right) \wedge \omega_k(t)
\end{aligned}
\right\} \forall k \in [1, N_{\text{fib}}], \quad (2.1)
$$

where $\mathbf{F}_{\mathbf{k},m}^{\text{rep}}(t)$ and $\mathbf{T}_{\mathbf{k},m}^{\text{rep}}(t)$ are the force and torque associated with the repulsion between fibres $k$ and $m$, $\mathbf{F}_{\mathbf{k},m}^{\text{rest}}(t)$ and $\mathbf{T}_{\mathbf{k},m}^{\text{rest}}(t)$ are the force and torque due to the presence of a spring (crosslink) connecting fibres $k$ and $m$, and $\mathbf{T}_{\mathbf{k},m}^{\text{align}}(t)$ is the alignment torque generated by this crosslink. We refer to appendix A.1 for the detailed computations of these forces and torques.

## 2.2. Description of the simulation set-up and biological relevance of the model parameters

The spatial domain $\Omega$ 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

$$ \Omega = \left[ -\frac{L_x}{2}, \frac{L_x}{2} \right] \times \left[ -\frac{L_y}{2}, \frac{L_y}{2} \right] \times \left[ -\frac{L_z}{2}, \frac{L_z}{2} \right]. $$

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_{\text{link}}^{\text{max}}$ and the link unloaded length $d_{\text{link}}^{\text{eq}}$ are both equal to the diameter of a fibre $2R_{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}/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$ [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 $\phi_{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:

$$ \phi_{\text{fib}} = \frac{N_{\text{fib}} V_{\text{fib}}}{|\Omega|} = N_{\text{fib}} \frac{\pi R_{\text{fib}}^2 L_{\text{fib}} + (4/3) \pi R_{\text{fib}}^3}{L_x L_y L_z}. \quad (2.2) $$

The quantity $\phi_{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 $\phi_{order}=0.89$, while for random or amorphous packing of spherocylinders with an aspect ratio of 6, the maximal random packing density $\phi_{random} \approx 0.4$ [20]. Thus, we may say that a system is 'sparse' if its fibre density is below $\phi_{random}$, 'dense' if it is between $\phi_{random}$ and $\phi_{order}$ and 'hyperdense' if it is above $\phi_{order}$. In the following, we will study two types of systems: dense systems
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 $\mathbf{X}_k(t) \in \Omega \subset \mathbb{R}^3$ and their non-oriented directional unit vectors $\omega_k(t) \in \mathbb{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

$$ 
\left. 
\begin{aligned}
\mu_{\text{fib}} L_{\text{fib}} \frac{d\mathbf{X}_k}{dt}(t) &= \sum_{m=1}^{N_{\text{fib}}} (\mathbf{F}_{\mathbf{k},m}^{\text{rep}}(t) + p_{k,m}(t) \mathbf{F}_{\mathbf{k},m}^{\text{rest}}(t)) \\
\mu_{\text{fib}} L_{\text{fib}}^3 \frac{d\omega_k}{dt}(t) &= \sum_{m=1}^{N_{\text{fib}}} \left( \mathbf{T}_{\mathbf{k},m}^{\text{rep}}(t) + p_{k,m}(t) (\mathbf{T}_{\mathbf{k},m}^{\text{rest}}(t) + \mathbf{T}_{\mathbf{k},m}^{\text{align}}(t)) \right) \wedge \omega_k(t)
\end{aligned}
\right\} \forall k \in [1, N_{\text{fib}}], \quad (2.1)
$$

where $\mathbf{F}_{\mathbf{k},m}^{\text{rep}}(t)$ and $\mathbf{T}_{\mathbf{k},m}^{\text{rep}}(t)$ are the force and torque associated with the repulsion between fibres $k$ and $m$, $\mathbf{F}_{\mathbf{k},m}^{\text{rest}}(t)$ and $\mathbf{T}_{\mathbf{k},m}^{\text{rest}}(t)$ are the force and torque due to the presence of a spring (crosslink) connecting fibres $k$ and $m$, and $\mathbf{T}_{\mathbf{k},m}^{\text{align}}(t)$ is the alignment torque generated by this crosslink. We refer to appendix A.1 for the detailed computations of these forces and torques.

## 2.2. Description of the simulation set-up and biological relevance of the model parameters

The spatial domain $\Omega$ 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

$$ \Omega = \left[ -\frac{L_x}{2}, \frac{L_x}{2} \right] \times \left[ -\frac{L_y}{2}, \frac{L_y}{2} \right] \times \left[ -\frac{L_z}{2}, \frac{L_z}{2} \right]. $$

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_{link}^{\text{max}}$ and the link unloaded length $d_{link}^{\text{eq}}$ are both equal to the diameter of a fibre $2R_{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}/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$ [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 $\phi_{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:

$$ \phi_{\text{fib}} = \frac{N_{\text{fib}} V_{\text{fib}}}{|\Omega|} = N_{\text{fib}} \frac{\pi R_{\text{fib}}^2 L_{\text{fib}} + (4/3) \pi R_{\text{fib}}^3}{L_x L_y L_z}. \quad (2.2) $$

The quantity $\phi_{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 $\phi_{order}=0.89$, while for random or amorphous packing of spherocylinders with an aspect ratio of 6, the maximal random packing density $\phi_{random} \approx 0.4$ [20]. Thus, we may say that a system is 'sparse' if its fibre density is below $\phi_{random}$, 'dense' if it is between $\phi_{random}$ and $\phi_{order}$ and 'hyperdense' if it is above $\phi_{order}$. In the following, we will study two types of systems: dense systems
<table><caption>**Table 1.** Model parameters.</caption><thead><tr><th>name</th><th>value</th><th>units</th><th>description</th></tr></thead><tbody><tr><td colspan="4">agents</td></tr><tr><td>N<sub>fib</sub></td><td>[1500, 3000]</td><td>n.a.</td><td>number of fibres</td></tr><tr><td>L<sub>fib</sub></td><td>6</td><td>L</td><td>fibre length</td></tr><tr><td>R<sub>fib</sub></td><td>0.5</td><td>L</td><td>fibre radius</td></tr><tr><td colspan="4">mechanical interactions</td></tr><tr><td>α<sub>rep</sub></td><td>12.5</td><td>M · L<sup>-1</sup> · T<sup>-2</sup></td><td>magnitude of the repulsion force</td></tr><tr><td>α<sub>rest</sub></td><td>5.0</td><td>M · T<sup>-2</sup></td><td>magnitude of the elastic restoring force</td></tr><tr><td>α<sub>align</sub></td><td>2.0</td><td>M · L<sup>2</sup> · T<sup>-2</sup></td><td>magnitude of the alignment torque</td></tr><tr><td>d<sup>max</sup><sub>link</sub></td><td>1.0</td><td>L</td><td>perception distance for link creation</td></tr><tr><td>d<sup>eq</sup><sub>link</sub></td><td>1.0</td><td>L</td><td>link equilibrium length</td></tr><tr><td colspan="4">biological phenomena</td></tr><tr><td>ν<sub>link</sub></td><td>[0, 10]</td><td>T<sup>-1</sup></td><td>network remodelling speed</td></tr><tr><td>χ<sub>link</sub></td><td>[0.1, 0.9]</td><td>n.a.</td><td>equilibrium linked fibre fraction</td></tr><tr><td colspan="4">numerical parameters</td></tr><tr><td>L<sub>x</sub> = L<sub>y</sub> = L<sub>z</sub></td><td>30</td><td>L</td><td>side length of the cubic domain</td></tr><tr><td>T<sub>final</sub></td><td>5.10<sup>4</sup></td><td>T</td><td>total time of simulation</td></tr></tbody></table>

containing $N_{fib} = 3000$ fibres ($\phi_{fib} = 0.58$) and sparse systems with $N_{fib} = 1500$ fibres ($\phi_{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 $\omega_1 = \omega_2$) to move apart by 99% of their equilibrium distance $2 R_{fib}$ (i.e. $\|X_1 - X_2\| = 0.99 \times 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_{link}^{eq}$. On the other hand, for nematic alignment $T_{align}$ is the time needed for two perpendicularly intersecting fibres ($X_1 = X_2$ and $\omega_1 \perp \omega_2$) crosslinked at their centre to reach a relative angle $\arccos(\omega_1 \cdot \omega_2) = 0.9^\circ$.

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

$$ 
\left. 
\begin{aligned}
T_{\text{rep}} &= \frac{27 \mu_{\text{fib}} L_{\text{fib}}}{4\sqrt{2} R_{\text{fib}} \alpha_{\text{rep}}} = 4.32 U_t, \\
T_{\text{rest}} &= \ln(10) \frac{\mu_{\text{fib}} L_{\text{fib}}}{\alpha_{\text{rest}}} = 2.76 U_t, \\
T_{\text{align}} &= 4.27 \frac{\mu_{\text{fib}} L_{\text{fib}}^3}{\alpha_{\text{align}}} = 462 U_t.
\end{aligned}
\right\} \quad (2.3) 
$$

and

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).

## 3. Results

### 3.1. 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
<table><caption>**Table 1.** Model parameters.</caption><thead><tr><th>name</th><th>value</th><th>units</th><th>description</th></tr></thead><tbody><tr><td colspan="4">agents</td></tr><tr><td>N<sub>fib</sub></td><td>[1500, 3000]</td><td>n.a.</td><td>number of fibres</td></tr><tr><td>L<sub>fib</sub></td><td>6</td><td>L</td><td>fibre length</td></tr><tr><td>R<sub>fib</sub></td><td>0.5</td><td>L</td><td>fibre radius</td></tr><tr><td colspan="4">mechanical interactions</td></tr><tr><td>α<sub>rep</sub></td><td>12.5</td><td>M·L<sup>-1</sup>·T<sup>-2</sup></td><td>magnitude of the repulsion force</td></tr><tr><td>α<sub>rest</sub></td><td>5.0</td><td>M·T<sup>-2</sup></td><td>magnitude of the elastic restoring force</td></tr><tr><td>α<sub>align</sub></td><td>2.0</td><td>M·L<sup>2</sup>·T<sup>-2</sup></td><td>magnitude of the alignment torque</td></tr><tr><td>d<sup>max</sup><sub>link</sub></td><td>1.0</td><td>L</td><td>perception distance for link creation</td></tr><tr><td>d<sup>eq</sup><sub>link</sub></td><td>1.0</td><td>L</td><td>link equilibrium length</td></tr><tr><td colspan="4">biological phenomena</td></tr><tr><td>ν<sub>link</sub></td><td>[0, 10]</td><td>T<sup>-1</sup></td><td>network remodelling speed</td></tr><tr><td>χ<sub>link</sub></td><td>[0.1, 0.9]</td><td>n.a.</td><td>equilibrium linked fibre fraction</td></tr><tr><td colspan="4">numerical parameters</td></tr><tr><td>L<sub>x</sub> = L<sub>y</sub> = L<sub>z</sub></td><td>30</td><td>L</td><td>side length of the cubic domain</td></tr><tr><td>T<sub>final</sub></td><td>5.10<sup>4</sup></td><td>T</td><td>total time of simulation</td></tr></tbody></table>

containing $N_{fib} = 3000$ fibres ($\phi_{fib} = 0.58$) and sparse systems with $N_{fib} = 1500$ fibres ($\phi_{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 $\omega_1 = \omega_2$) to move apart by 99% of their equilibrium distance $2 R_{fib}$ (i.e. $\|X_1 - X_2\| = 0.99 \times 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_{link}^{eq}$. On the other hand, for nematic alignment $T_{align}$ is the time needed for two perpendicularly intersecting fibres ($X_1 = X_2$ and $\omega_1 \perp \omega_2$) crosslinked at their centre to reach a relative angle $\arccos(\omega_1 \cdot \omega_2) = 0.9^\circ$.

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

$$ 
\left. 
\begin{aligned}
T_{\text{rep}} &= \frac{27 \mu_{\text{fib}} L_{\text{fib}}}{4\sqrt{2} R_{\text{fib}} \alpha_{\text{rep}}} = 4.32 U_t, \\
T_{\text{rest}} &= \ln(10) \frac{\mu_{\text{fib}} L_{\text{fib}}}{\alpha_{\text{rest}}} = 2.76 U_t, \\
T_{\text{align}} &= 4.27 \frac{\mu_{\text{fib}} L_{\text{fib}}^3}{\alpha_{\text{align}}} = 462 U_t.
\end{aligned}
\right\} (2.3) $$

and

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).

# 3. Results

## 3.1. 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
<figure><img src="image_5.png" /><figcaption>**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_{linkperfib}^{mean}$ at equilibrium, with colour depending on the remodelling speed $v_{link}$ and horizontal and vertical error-bars indicating the standard deviation $N_{linkperfib}^{STD}$ and $Al_{STD}$, respectively. The grey dashed-line indicates the critical value of $N_{linkperfib}^{mean}$ and the black dashed lines the three logarithmic fits obtained for $N_{linkperfib}^{mean} < N_{critic}$.</figcaption></figure>

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
<figure><img src="image_5.png" /><figcaption>**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_{linkperfib}^{mean}$ at equilibrium, with colour depending on the remodelling speed $v_{link}$ and horizontal and vertical error-bars indicating the standard deviation $N_{linkperfib}^{STD}$ and $Al_{STD}$, respectively. The grey dashed-line indicates the critical value of $N_{linkperfib}^{mean}$ and the black dashed lines the three logarithmic fits obtained for $N_{linkperfib}^{mean} < N_{critic}$.</figcaption></figure>

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
links per fibres in a network and by $N_{\text{linkperfib}}^{\text{mean}}$ and $N_{\text{linkperfib}}^{\text{STD}}$ its average and standard deviation over 10 simulations. We stress the fact that $N_{\text{linkperfib}} = 0.5$ if all the fibres are connected to a neighbouring fibre.

By plotting the alignment quantifier $Al_{\text{mean}}$ as a function of the number of links per fibre $N_{\text{linkperfib}}^{\text{mean}}$ (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_{\text{linkperfib}}^{\text{STD}}$ and $Al_{\text{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 $v_{\text{link}}$, and inside each colour series $\chi_{\text{link}}$ is increasing with $N_{\text{linkperfib}}^{\text{mean}}$.

Figure 1d reveals that the values of $Al_{\text{mean}}$ and $N_{\text{linkperfib}}^{\text{mean}}$ at equilibrium are highly correlated. When $N_{\text{linkperfib}}^{\text{mean}}$ is inferior to a critical threshold $N_{\text{critic}} \approx 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)

$$Al_{\text{mean}} \approx \alpha \log(N_{\text{linkperfib}}^{\text{mean}}) + \beta, \quad (3.1)$$

with

- $\alpha = 0.037, \beta = 1.006$ and coefficient of determination $r^2 = 0.87$ for dynamical systems (non-blue markers);

- $\alpha = 0.129, \beta = 0.651$ and coefficient of determination $r^2 = 0.96$ for sparse non-dynamical networks (blue triangles);

- $\alpha = 0.042, \beta = 0.433$ and coefficient of determination $r^2 = 0.985$ for dense non-dynamical networks (blue dots).

Then, when $N_{\text{linkperfib}}^{\text{mean}} > N_{\text{critic}}$ we observe an abrupt drop of the equilibrium value of $Al_{\text{mean}}$. Surprisingly and very interestingly, for dynamical systems ($v_{\text{link}} > 0$) there is no difference in alignment induced by the fibre density or the link characteristics $v_{\text{link}}$ and $\chi_{\text{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_{\text{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 non-dynamical networks this correlation significantly depends on the fibre density. However, the critical number of links $N_{\text{critic}}$ allowing for larger alignment is the same for non-dynamical networks, either dense or sparse, and for dynamical networks. Therefore, $N_{\text{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.

## 3.2. 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
links per fibres in a network and by $N_{\text{linkperfib}}^{\text{mean}}$ and $N_{\text{linkperfib}}^{\text{STD}}$ its average and standard deviation over 10 simulations. We stress the fact that $N_{\text{linkperfib}} = 0.5$ if all the fibres are connected to a neighbouring fibre.

By plotting the alignment quantifier $Al_{\text{mean}}$ as a function of the number of links per fibre $N_{\text{linkperfib}}^{\text{mean}}$ (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_{\text{linkperfib}}^{\text{STD}}$ and $Al_{\text{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 $v_{\text{link}}$, and inside each colour series $\chi_{\text{link}}$ is increasing with $N_{\text{linkperfib}}^{\text{mean}}$.

Figure 1d reveals that the values of $Al_{\text{mean}}$ and $N_{\text{linkperfib}}^{\text{mean}}$ at equilibrium are highly correlated. When $N_{\text{linkperfib}}^{\text{mean}}$ is inferior to a critical threshold $N_{\text{critic}} \approx 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)

$$Al_{\text{mean}} \approx \alpha \log(N_{\text{linkperfib}}^{\text{mean}}) + \beta, \quad (3.1)$$

with

- $\alpha = 0.037$, $\beta = 1.006$ and coefficient of determination $r^2 = 0.87$ for dynamical systems (non-blue markers);

- $\alpha = 0.129$, $\beta = 0.651$ and coefficient of determination $r^2 = 0.96$ for sparse non-dynamical networks (blue triangles);

- $\alpha = 0.042$, $\beta = 0.433$ and coefficient of determination $r^2 = 0.985$ for dense non-dynamical networks (blue dots).

Then, when $N_{\text{linkperfib}}^{\text{mean}} > N_{\text{critic}}$ we observe an abrupt drop of the equilibrium value of $Al_{\text{mean}}$. Surprisingly and very interestingly, for dynamical systems ($v_{\text{link}} > 0$) there is no difference in alignment induced by the fibre density or the link characteristics $v_{\text{link}}$ and $\chi_{\text{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_{\text{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 non-dynamical networks this correlation significantly depends on the fibre density. However, the critical number of links $N_{\text{critic}}$ allowing for larger alignment is the same for non-dynamical networks, either dense or sparse, and for dynamical networks. Therefore, $N_{\text{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.

## 3.2. 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
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(c)</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(d)</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(e)</figcaption></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 2.** (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–e) displays the equilibrium state of a few other simulations, whose positions on the diagram are also indicated with a black star. (f) Heatmap of the percentage of simulations ending in an aligned state (versus a curved state), for dynamical dense networks as function of the values of the network remodelling speed $v_{link}$ (in ordinate) and the equilibrium linked fibre fraction $\chi_{link}$ (in abscissa). (g) Heatmap of the number of links per fibre $N_{linkperfib}^{mean}$ for dynamical dense networks as function of the values of the network remodelling speed $v_{link}$ (in ordinate) and the equilibrium linked fibre fraction $\chi_{link}$ (in abscissa). Zones where $N_{linkperfib}^{mean} \approx N_{critic}$ have been highlighted in green.</figcaption></figure>

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 $v_{link}$, $\chi_{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. $v_{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
<figure><img src="image_6.png" /><figcaption>**Figure 2.** (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–e) displays the equilibrium state of a few other simulations, whose positions on the diagram are also indicated with a black star. (f) Heatmap of the percentage of simulations ending in an aligned state (versus a curved state), for dynamical dense networks as function of the values of the network remodelling speed $v_{link}$ (in ordinate) and the equilibrium linked fibre fraction $\chi_{link}$ (in abscissa). (g) Heatmap of the number of links per fibre $N_{linkperfib}^{mean}$ for dynamical dense networks as function of the values of the network remodelling speed $v_{link}$ (in ordinate) and the equilibrium linked fibre fraction $\chi_{link}$ (in abscissa). Zones where $N_{linkperfib}^{mean} \approx N_{critic}$ have been highlighted in green.</figcaption></figure>

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 $v_{link}$, $\chi_{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. $v_{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
<table><caption>**Table 2.** Classification of the simulations outcomes into different states based on the local quantifier $A_{sim}$ and the global quantifier $A_{max}$. The case { $A_{sim} < 0.7$ & $A_{max} \le 0.45$ } never occurs in our simulations and is thus unnamed.</caption><thead><tr><th>A<sub>max</sub></th><th></th><th></th></tr></thead><tbody><tr><td></td><td>&le;0.45</td><td>&gt;0.45</td></tr><tr><td rowspan="2">A<sub>sim</sub></td><td>&ge;0.7</td><td>aligned state: alignment both local and global</td></tr><tr><td>&lt;0.7</td><td>(alignment global but not local)</td></tr></tbody></table>

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 three-dimensional view and stereographic projection in the *b* to *e*.

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 $v_{link}$. Indeed, unorganized states were only and systematically observed for non-dynamical networks ($v_{link} = 0$), while dynamical networks ($v_{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 $v_{link}$ (in ordinate) and the equilibrium linked fibre fraction $\chi_{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 $v_{link}$ and $\chi_{link}$. Indeed, analysis of the heatmap reveals that (i) reduced linked fibre fraction $\chi_{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 $v_{link}$ and $\chi_{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}$, $v_{link}$ and $\chi_{link}$ is more complex. Indeed, figure 2g shows a heatmap of the number of links per fibre in the network $N_{linkperfib}^{mean}$ as a function of $v_{link}$ and $\chi_{link}$ for dense dynamical networks (same
<table><caption>**Table 2.** Classification of the simulations outcomes into different states based on the local quantifier $A_{sim}$ and the global quantifier $A_{max}$. The case { $A_{sim} < 0.7$ & $A_{max} \le 0.45$ } never occurs in our simulations and is thus unnamed.</caption><thead><tr><th>A<sub>max</sub></th><th></th><th></th></tr></thead><tbody><tr><td></td><td>&le;0.45</td><td>&gt;0.45</td></tr><tr><td rowspan="2">A<sub>sim</sub></td><td>&ge;0.7</td><td>aligned state: alignment both local and global</td></tr><tr><td>&lt;0.7</td><td>(alignment global but not local)</td></tr></tbody></table>

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 three-dimensional view and stereographic projection in the *b* to *e*.

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 $v_{link}$. Indeed, unorganized states were only and systematically observed for non-dynamical networks ($v_{link} = 0$), while dynamical networks ($v_{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 $v_{link}$ (in ordinate) and the equilibrium linked fibre fraction $\chi_{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 $v_{link}$ and $\chi_{link}$. Indeed, analysis of the heatmap reveals that (i) reduced linked fibre fraction $\chi_{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 $v_{link}$ and $\chi_{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}$, $v_{link}$ and $\chi_{link}$ is more complex. Indeed, figure 2g shows a heatmap of the number of links per fibre in the network $N_{linkperfib}^{mean}$ as a function of $v_{link}$ and $\chi_{link}$ for dense dynamical networks (same
simulations as f). It demonstrates that $N_{\text{linkperfib}}^{\text{mean}}$ is indeed an emerging variable, in the sense that it is not directly linked to the parameters $v_{\text{link}}$ and $\chi_{\text{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 $v_{\text{link}}$ decreases and $\chi_{\text{link}}$ increases (from top left to bottom right corner of the heatmap), crossing the critical threshold $N_{\text{critic}}$ doing so (the cells where $N_{\text{linkperfib}}^{\text{mean}} \approx N_{\text{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.

### 3.3. 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 $A_{\text{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 $\tau_{\text{Al}}$, which corresponds to the time needed for the quantifier to reach 63% of its asymptotic value ($A_{\text{mean}} \approx 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\tau_{\text{Al}}$, $\tau_{\text{Al}}$, $3\tau_{\text{Al}}$ and $T_{\text{final}}$) for two of these simulations (respectively from Movie3.mp4 and Movie4.mp4). They correspond to dense systems with $\chi_{\text{link}} = 0.8$ and two different crosslink dynamics: fast remodelling network $v_{\text{link}} = 0.1$ (a–a''', Movie3.mp4) and slow remodelling network $v_{\text{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_{\text{max}}$ versus $A_{\text{sim}}$ of simulations for low-dynamical dense networks $v_{\text{link}} = 0.001$ with various linked fibre fractions $\chi_{\text{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 $\tau_{\text{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 $A_{\text{mean}}$ in figure 10). This last point comes from the local quality of the quantifier $A_{\text{sim}}$ (and by extension $A_{\text{mean}}$): a system can be very
simulations as f). It demonstrates that $N_{\text{linkperfib}}^{\text{mean}}$ is indeed an emerging variable, in the sense that it is not directly linked to the parameters $v_{\text{link}}$ and $\chi_{\text{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 $v_{\text{link}}$ decreases and $\chi_{\text{link}}$ increases (from top left to bottom right corner of the heatmap), crossing the critical threshold $N_{\text{critic}}$ doing so (the cells where $N_{\text{linkperfib}}^{\text{mean}} \approx N_{\text{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.

### 3.3. 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 $A_{\text{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 $\tau_{\text{Al}}$, which corresponds to the time needed for the quantifier to reach 63% of its asymptotic value ($A_{\text{mean}} \approx 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\tau_{\text{Al}}$, $\tau_{\text{Al}}$, $3\tau_{\text{Al}}$ and $T_{\text{final}}$) for two of these simulations (respectively from Movie3.mp4 and Movie4.mp4). They correspond to dense systems with $\chi_{\text{link}} = 0.8$ and two different crosslink dynamics: fast remodelling network $v_{\text{link}} = 0.1$ (a–a''', Movie3.mp4) and slow remodelling network $v_{\text{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_{\text{max}}$ versus $A_{\text{sim}}$ of simulations for low-dynamical dense networks $v_{\text{link}} = 0.001$ with various linked fibre fractions $\chi_{\text{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 $\tau_{\text{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 $A_{\text{mean}}$ in figure 10). This last point comes from the local quality of the quantifier $A_{\text{sim}}$ (and by extension $A_{\text{mean}}$): a system can be very
<figure><img src="image_13.png" /><figcaption>**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\tau_{Al}$ (a), $\tau_{Al}$ (a'), $3\tau_{Al}$ (a'') and $T_{final}$ (a'''), for one simulation with $v_{link} = 0.1$ and $\chi_{link} = 0.8$. Panels b–b''': Stereographic projection of the system at times $0.5\tau_{Al}$ (b), $\tau_{Al}$ (b'), $3\tau_{Al}$ (b'') and $T_{final}$ (b'''), for one simulation with $v_{link} = 0.001$ and $\chi_{link} = 0.8$. (c) Trajectory in the phase plane $A_{max}$ versus $A_{sim}$ of individual simulations for slow-remodelling dense networks $v_{link} = 0.001$ and various linked fibre fractions $\chi_{link}$. The initial position is indicated with a black square, the final position with a black star and the time-constant $\tau_{Al}$ with a black circle. The limits between each class of structures are drawn in dashed lines. (d) Evolution of $N_{linkperfib}^{mean}$ for slow-remodelling dense networks $v_{link} = 0.001$ and various linked fibre fractions $\chi_{link}$ with shading indicating the inter-simulation standard deviation $N_{linkperfib}^{STD}$. The critical value $N_{critic}$ is indicated with a dashed line and the time-constant $\tau_{Al}$ with a black circle.</figcaption></figure>

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 $\chi_{link}$ is high. For the extreme case of large linked fibre fraction $\chi_{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 ($\chi_{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 $\chi_{link}$ describes the proportion of linked fibres among all linkable fibres at equilibrium, but this equilibrium takes time to establish (inversely proportional to $v_{link}$). If the repulsion interaction operates faster than the links remodelling (i.e. $T_{rep} \ll 1/v_{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 $v_{link}$ and $v_{unlink}$ is more favourable to linking than unlinking ($v_{link} > v_{unlink}$, i.e. if $\chi_{link} > 0.5$).
<figure><img src="image_7.png" /><figcaption>**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\tau_{Al}$ (a), $\tau_{Al}$ (a'), $3\tau_{Al}$ (a'') and $T_{final}$ (a'''), for one simulation with $v_{link} = 0.1$ and $\chi_{link} = 0.8$. Panels b–b''': Stereographic projection of the system at times $0.5\tau_{Al}$ (b), $\tau_{Al}$ (b'), $3\tau_{Al}$ (b'') and $T_{final}$ (b'''), for one simulation with $v_{link} = 0.001$ and $\chi_{link} = 0.8$. (c) Trajectory in the phase plane $A_{max}$ versus $A_{sim}$ of individual simulations for slow-remodelling dense networks $v_{link} = 0.001$ and various linked fibre fractions $\chi_{link}$. The initial position is indicated with a black square, the final position with a black star and the time-constant $\tau_{Al}$ with a black circle. The limits between each class of structures are drawn in dashed lines. (d) Evolution of $N_{linkperfib}^{mean}$ for slow-remodelling dense networks $v_{link} = 0.001$ and various linked fibre fractions $\chi_{link}$ with shading indicating the inter-simulation standard deviation $N_{linkperfib}^{STD}$. The critical value $N_{critic}$ is indicated with a dashed line and the time-constant $\tau_{Al}$ with a black circle.</figcaption></figure>

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 $\chi_{link}$ is high. For the extreme case of large linked fibre fraction $\chi_{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 ($\chi_{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 $\chi_{link}$ describes the proportion of linked fibres among all linkable fibres at equilibrium, but this equilibrium takes time to establish (inversely proportional to $v_{link}$). If the repulsion interaction operates faster than the links remodelling (i.e. $T_{rep} \ll 1/v_{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 $v_{link}$ and $v_{unlink}$ is more favourable to linking than unlinking ($v_{link} > v_{unlink}$, i.e. if $\chi_{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 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.

## 4. 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 ($v_{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,
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.

## 4. 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 ($v_{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,
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
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
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 self-organization 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 three-dimensional 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–58].

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

**Data accessibility.** Data and relevant code for this research work are stored in GitHub: https://github.com/chassonnery/3D_DynamicalFiberNetwork and have been archived within the Zenodo repository: https://doi.org/10.5281/zenodo.8416498 [59].

Supplementary material is available online [60].

**Declaration of AI use.** We have not used AI-assisted technologies in creating this article.

**Authors' contributions.** P.C.: conceptualization, data curation, formal analysis, investigation, methodology, software, validation, visualization, writing—original draft, writing—review and editing; J.P.: conceptualization, formal analysis, investigation, methodology, visualization, writing—original draft, writing—review and editing; A.L.: conceptualization, formal analysis, investigation, methodology, visualization, writing—original draft, writing—review and editing; C.S.: conceptualization, formal analysis, investigation, methodology, visualization, writing—review and editing; M.O.: conceptualization, methodology; P.D.: conceptualization, methodology, supervision, visualization, writing—review and editing; L.C.: conceptualization, formal analysis, funding acquisition, investigation, methodology, project administration, resources, supervision, visualization, writing—original draft, writing—review and editing; D.P.: conceptualization, formal analysis, funding acquisition, investigation, methodology, project administration, resources, software, supervision, visualization, writing—original draft, writing—review and editing.

All authors gave final approval for publication and agreed to be held accountable for the work performed therein.

**Conflict of interest declaration.** We declare we have no competing interests.

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

**Acknowledgements.** P.D. holds a visiting professor association with the Department of Mathematics, Imperial College London, UK.

## Appendix A. Model

Here, we give details about the mathematical model presented in §2.1 of the main text. Let us recall the main features and introduce some notations. The $N_{fib}$ fibre elements are represented by their centres $\mathbf{X}_k(t) \in \mathbb{R}^3$ and their non-oriented directional unit vectors $\omega_k(t) \in \mathbb{S}_2^+$. The fibres repel their close neighbours by means of a soft repulsion mechanism modelling steric repulsion between spherocylinders of radius $R_{fib}$. Fibre elements have the ability to link to or unlink from each other to model fibre elongation or rupture. The linking and unlinking of fibres follow random (Poisson) processes in time. Fibres offer resistance to bending through an alignment torque acting between two linked fibre elements. Fibre motion and rotation is given by Newton's second law of motion, in an overdamped regime to model a medium with low Reynolds number.
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 self-organization 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 three-dimensional 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–58].

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

**Data accessibility.** Data and relevant code for this research work are stored in GitHub: https://github.com/chassonnery/3D_DynamicalFiberNetwork and have been archived within the Zenodo repository: https://doi.org/10.5281/zenodo.8416498 [59].

Supplementary material is available online [60].

**Declaration of AI use.** We have not used AI-assisted technologies in creating this article.

**Authors' contributions.** P.C.: conceptualization, data curation, formal analysis, investigation, methodology, software, validation, visualization, writing—original draft, writing—review and editing; J.P.: conceptualization, formal analysis, investigation, methodology, visualization, writing—original draft, writing—review and editing; A.L.: conceptualization, formal analysis, investigation, methodology, visualization, writing—original draft, writing—review and editing; C.S.: conceptualization, formal analysis, investigation, methodology, visualization, writing—review and editing; M.O.: conceptualization, methodology; P.D.: conceptualization, methodology, supervision, visualization, writing—review and editing; L.C.: conceptualization, formal analysis, funding acquisition, investigation, methodology, project administration, resources, supervision, visualization, writing—original draft, writing—review and editing; D.P.: conceptualization, formal analysis, funding acquisition, investigation, methodology, project administration, resources, software, supervision, visualization, writing—original draft, writing—review and editing.

All authors gave final approval for publication and agreed to be held accountable for the work performed therein.

**Conflict of interest declaration.** We declare we have no competing interests.

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

**Acknowledgements.** P.D. holds a visiting professor association with the Department of Mathematics, Imperial College London, UK.

## Appendix A. Model

Here, we give details about the mathematical model presented in §2.1 of the main text. Let us recall the main features and introduce some notations. The $N_{fib}$ fibre elements are represented by their centres $\mathbf{X}_k(t) \in \mathbb{R}^3$ and their non-oriented directional unit vectors $\omega_k(t) \in \mathbb{S}_2^+$. The fibres repel their close neighbours by means of a soft repulsion mechanism modelling steric repulsion between spherocylinders of radius $R_{fib}$. Fibre elements have the ability to link to or unlink from each other to model fibre elongation or rupture. The linking and unlinking of fibres follow random (Poisson) processes in time. Fibres offer resistance to bending through an alignment torque acting between two linked fibre elements. Fibre 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 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_{k,m}^{rep}$ the pairwise interaction force between the spherocylinders $k$ and $m$ and using Hertzian theory [39]

$$F_{k,m}^{rep} = \alpha_{\text{rep}}(2R_{\text{fib}} - \|\mathbf{X}_{k,m} - \mathbf{X}_{m,k}\|)^{3/2} \sqrt{2R_{\text{fib}}} \times \frac{\mathbf{X}_{k,m} - \mathbf{X}_{m,k}}{\|\mathbf{X}_{k,m} - \mathbf{X}_{m,k}\|}, \quad (A1)$$

where $\alpha_{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

$$T_{k,m}^{\text{rep}} = (X_{k,m} - X_k) \wedge F_{k,m}^{\text{rep}}, \quad (A2)$$

on fibre $k$.

(ii) **Computation of the fibre attachment forces due to crosslinks.** Fibres closer than the threshold $d_{link}^{max}$ can create a crosslink, modelled as a linear spring of stiffness $\alpha_{rest}$ and unloaded length $d_{link}^{eq}$ 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

$$F_{k,m}^{\text{rest}} = \alpha_{\text{rest}} \left( d_{\text{link}}^{\text{eq}} - \|\mathbf{X}_{k,m}^{l} - \mathbf{X}_{m,k}^{l}\| \right) \frac{\mathbf{X}_{k,m}^{l} - \mathbf{X}_{m,k}^{l}}{\|\mathbf{X}_{k,m}^{l} - \mathbf{X}_{m,k}^{l}\|}, \quad (A3)$$

where $\mathbf{X}_{k,m}^{l}$ 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$

$$T_{k,m}^{\text{rest}} = (X_{k,m}^{l} - X_k) \wedge F_{k,m}^{\text{rest}}, \quad (A4)$$

To ensure coherence between the different features of the model, we require that $2R_{fib} \le d_{link}^{eq} \le d_{link}^{max}$.

(iii) **Computation of the linked fibre alignment forces.** It is characterized by a stiffness parameter $\alpha_{align} > 0$ playing the role of a flexural modulus: the larger $\alpha_{align}$, the more rigid the fibre network. Given two linked fibres $k$ and $m$, the torque sustained by the fibre $k$ is such that, $\forall \mathbf{u} \in \mathbb{R}^3$,

$$T_{k,m}^{\text{align}} \wedge \mathbf{u} = \alpha_{\text{align}} \left( (\omega_k \wedge \tilde{\omega}_m) \wedge \mathbf{u} + \frac{1 - |\omega_k \cdot \omega_m|}{\|\omega_k \wedge \omega_m\|^2} (\omega_k \wedge \tilde{\omega}_m) \wedge ((\omega_k \wedge \tilde{\omega}_m) \wedge \mathbf{u}) \right), \quad (A5)$$

where $\tilde{\omega}_m = \text{sign}(\omega_k \cdot \omega_m) \cdot \omega_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 $\mu_{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

$$F_{k}^{\text{fric}} = -\mu_{\text{fib}} L_{\text{fib}}^{3} \omega_{k} \frac{d\mathbf{X}_{k}}{dt}, \quad (A6)$$

and the associated rotational torque is equal to

$$T_{k}^{\text{fric}} = -\mu_{\text{fib}} L_{\text{fib}}^{3} \omega_{k} \wedge \frac{d\omega_{k}}{dt}. \quad (A7)$$

For the sake of simplicity and because our model features only one type of element, we assumed a single friction coefficient $\mu_{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
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_{k,m}^{rep}$ the pairwise interaction force between the spherocylinders $k$ and $m$ and using Hertzian theory [39]

$$F_{k,m}^{rep} = \alpha_{\text{rep}}(2R_{\text{fib}} - \|\mathbf{X}_{k,m} - \mathbf{X}_{m,k}\|)^{3/2} \sqrt{2R_{\text{fib}}} \times \frac{\mathbf{X}_{k,m} - \mathbf{X}_{m,k}}{\|\mathbf{X}_{k,m} - \mathbf{X}_{m,k}\|}, \quad (A1)$$

where $\alpha_{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

$$T_{k,m}^{\text{rep}} = (X_{k,m} - X_k) \wedge F_{k,m}^{\text{rep}}, \quad (A2)$$

on fibre $k$.

(ii) **Computation of the fibre attachment forces due to crosslinks.** Fibres closer than the threshold $d_{link}^{max}$ can create a crosslink, modelled as a linear spring of stiffness $\alpha_{rest}$ and unloaded length $d_{link}^{eq}$ 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

$$F_{k,m}^{\text{rest}} = \alpha_{\text{rest}} \left( d_{\text{link}}^{\text{eq}} - \|\mathbf{X}_{k,m}^{l} - \mathbf{X}_{m,k}^{l}\| \right) \frac{\mathbf{X}_{k,m}^{l} - \mathbf{X}_{m,k}^{l}}{\|\mathbf{X}_{k,m}^{l} - \mathbf{X}_{m,k}^{l}\|}, \quad (A3)$$

where $\mathbf{X}_{k,m}^{l}$ 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$

$$T_{k,m}^{\text{rest}} = (X_{k,m}^{l} - X_k) \wedge F_{k,m}^{\text{rest}}, \quad (A4)$$

To ensure coherence between the different features of the model, we require that $2R_{fib} \le d_{link}^{eq} \le d_{link}^{max}$.

(iii) **Computation of the linked fibre alignment forces.** It is characterized by a stiffness parameter $\alpha_{align} > 0$ playing the role of a flexural modulus: the larger $\alpha_{align}$, the more rigid the fibre network. Given two linked fibres $k$ and $m$, the torque sustained by the fibre $k$ is such that, $\forall \mathbf{u} \in \mathbb{R}^3$,

$$T_{k,m}^{\text{align}} \wedge \mathbf{u} = \alpha_{\text{align}} \left( (\omega_k \wedge \tilde{\omega}_m) \wedge \mathbf{u} + \frac{1 - |\omega_k \cdot \omega_m|}{\|\omega_k \wedge \omega_m\|^2} (\omega_k \wedge \tilde{\omega}_m) \wedge ((\omega_k \wedge \tilde{\omega}_m) \wedge \mathbf{u}) \right), \quad (A5)$$

where $\tilde{\omega}_m = \text{sign}(\omega_k \cdot \omega_m) \cdot \omega_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 $\mu_{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

$$F_{k}^{\text{fric}} = -\mu_{\text{fib}}L_{\text{fib}}^{3}\omega_{k}\frac{d\mathbf{X}_{k}}{dt}, \quad (A6)$$

and the associated rotational torque is equal to

$$T_{k}^{\text{fric}} = -\mu_{\text{fib}}L_{\text{fib}}^{3}\omega_{k}\frac{d\omega_{k}}{dt}. \quad (A7)$$

For the sake of simplicity and because our model features only one type of element, we assumed a single friction coefficient $\mu_{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
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/v$ of the Poisson process, so we prescribe $dt \le dt_{\text{link}}$ with

$$dt_{\text{link}} = \min\left(\frac{0.5}{v_{\text{link}}}, \frac{0.5}{v_{\text{unlink}}}\right). \qquad (A 8)$$

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

$$\mathbb{P}(p_{k,m}(t_{n+1}) = 1 | p_{k,m}(t_n) = 0 \text{ and } \|\mathbf{X}_{k,m}(t_n) - \mathbf{X}_{m,k}(t_n)\| \le d_{\text{link}}^{\max} = 1 - e^{-v_{\text{link}}dt_n} \qquad (A 9)$$

while the probability for a crosslink to break is given by

$$\mathbb{P}(p_{k,m}(t_{n+1}) = 0 | p_{k,m}(t_n) = 1) = 1 - e^{-v_{\text{unlink}}dt_n}. \qquad (A 10)$$

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_{\text{fib}}$ and its rotation to $\arctan(0.1) \approx 6^\circ$. This implies the following upper limits for the time step:

$$\begin{gathered} \ddot{t}_{\text{trans}}(t_n) = \min_{1 \le k \le N_{\text{fib}}} \left( 0.5 \left| \frac{d\mathbf{X}_k}{dt}(t_n) \right| \right) \\ \ddot{t}_{\text{rot}}(t_n) = \min_{1 \le k \le N_{\text{fib}}} \left( \frac{0.1}{\left| \frac{d\omega_k}{dt}(t_n) \right|} \right). \end{gathered} \qquad (A 11)$$

and

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))

$$dt_n = \min(dt_{\text{trans}}(t_n), dt_{\text{rot}}(t_n), dt_{\text{link}}).$$

- Motion of the agents to their new position

$$\mathbf{X}_k(t_{n+1}) = \mathbf{X}_k(t_n) + dt_n \frac{d\mathbf{X}_k}{dt}(t_n)$$

and

$$\omega_k(t_{n+1}) = \omega_k(t_n) + dt_n \frac{d\omega_k}{dt}(t_n)$$

- Account for periodic boundary conditions.

- Attribution of each agent to a simulation box.

- Parallel update of linking configuration (equations (A 9) and (A 10)).

## A.3. Closest points of two finite segments

Given two fibres $k$ and $m$, we denote by $\mathbf{X}_{k,m} = \mathbf{X}_k + l_{k,m}\omega_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 $\|\mathbf{X}_k + u\omega_k - (\mathbf{X}_m + v\omega_m)\|$ for $(u, v) \in [-L_{\text{fib}}/2, L_{\text{fib}}]$
## 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/v$ of the Poisson process, so we prescribe $dt \le dt_{\text{link}}$ with

$$ dt_{\text{link}} = \min \left( \frac{0.5}{v_{\text{link}}}, \frac{0.5}{v_{\text{unlink}}} \right). \qquad (A 8) $$

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

$$ \mathbb{P}(p_{k,m}(t_{n+1}) = 1 | p_{k,m}(t_n) = 0 \text{ and } \|\mathbf{X}_{k,m}(t_n) - \mathbf{X}_{m,k}(t_n)\| \le d_{\text{link}}^{\max} = 1 - e^{-v_{\text{link}}dt_n} \qquad (A 9) $$

while the probability for a crosslink to break is given by

$$ \mathbb{P}(p_{k,m}(t_{n+1}) = 0 | p_{k,m}(t_n) = 1) = 1 - e^{-v_{\text{unlink}}dt_n}. \qquad (A 10) $$

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_{\text{fib}}$ and its rotation to $\arctan(0.1) \approx 6^\circ$. This implies the following upper limits for the time step:

$$ \mathbb{H}t_{\text{trans}}(t_n) = \min_{1 \le k \le N_{\text{fib}}} \left\{ 0.5 \left| \frac{d\mathbf{X}_k}{dt}(t_n) \right|^{R_{\text{fib}}} \right\} \qquad (A 11) $$

and

$$ dt_{\text{rot}}(t_n) = \min_{1 \le k \le N_{\text{fib}}} \left( \frac{0.1}{\left| \frac{d\omega_k}{dt}(t_n) \right|} \right). $$

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))

$$ dt_n = \min(dt_{\text{trans}}(t_n), dt_{\text{rot}}(t_n), dt_{\text{link}}). $$

- Motion of the agents to their new position

$$ \mathbf{X}_k(t_{n+1}) = \mathbf{X}_k(t_n) + dt_n \frac{d\mathbf{X}_k}{dt}(t_n) $$

and

$$ \omega_k(t_{n+1}) = \omega_k(t_n) + dt_n \frac{d\omega_k}{dt}(t_n) $$

- Account for periodic boundary conditions.

- Attribution of each agent to a simulation box.

- Parallel update of linking configuration (equations (A 9) and (A 10)).

## A.3. Closest points of two finite segments

Given two fibres $k$ and $m$, we denote by $\mathbf{X}_{k,m} = \mathbf{X}_k + l_{k,m}\omega_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 $\|\mathbf{X}_k + u\omega_k - (\mathbf{X}_m + v\omega_m)\|$ for $(u, v) \in [-L_{\text{fib}}/2, L_{\text{fib}}]$
<figure><img src="image_14.png" /><figcaption>**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.</figcaption></figure>

2)]. If $\omega_k = \omega_{m'}$, there is an infinity of solutions of the form $v = u + (X_k - X_m) \cdot \omega_{k'}$ in this case, we arbitrarily chose the solution with the smallest $|u|$ value. Otherwise, there exists a unique solution whose analytical expression is

$$ \left. \begin{aligned} l_{k,m} &= C_{L_{\text{fib}}/2} \left( \frac{\left( (\omega_k \cdot \omega_m) \omega_m \cdot (X_k - X_m) - \omega_k \cdot (X_k - X_m) \right)}{(1 - (\omega_k \cdot \omega_m)^2)} \right) \\ l_{m,k} &= C_{L_{\text{fib}}/2} \left( \frac{\left( (\omega_k \cdot \omega_m) \omega_k \cdot (X_m - X_k) - \omega_m \cdot (X_m - X_k) \right)}{(1 - (\omega_k \cdot \omega_m)^2)} \right), \end{aligned} \right\} \quad (A 12) $$

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 $A_l$ as the fractional anisotropy of the fibres' directional vectors within $B_k$.

It is computed as follows. We denote by $p_m = \omega_m \otimes \omega_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

$$ P_k = \frac{1}{|B_k|} \sum_{m_{s,t} \in B_k} p_m, \quad (B1) $$

where $|B_k|$ denotes the number of fibres in $B_k$.
<figure><img src="image_8.png" /><figcaption>**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.</figcaption></figure>

2)]. If $\omega_k = \omega_{m'}$, there is an infinity of solutions of the form $v = u + (X_k - X_m) \cdot \omega_{k'}$ in this case, we arbitrarily chose the solution with the smallest $|u|$ value. Otherwise, there exists a unique solution whose analytical expression is

$$ \left. \begin{aligned} l_{k,m} &= C_{L_{\text{fib}}/2} \left( \frac{\left( \left( \omega_k \cdot \omega_m \right) \omega_m \cdot (X_k - X_m) - \omega_k \cdot (X_k - X_m) \right)}{\left( 1 - (\omega_k \cdot \omega_m)^2 \right)} \right) \\ l_{m,k} &= C_{L_{\text{fib}}/2} \left( \frac{\left( \left( \omega_k \cdot \omega_m \right) \omega_k \cdot (X_m - X_k) - \omega_m \cdot (X_m - X_k) \right)}{\left( 1 - (\omega_k \cdot \omega_m)^2 \right)} \right), \end{aligned} \right\} \quad (A 12) $$

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 $A_l$ as the fractional anisotropy of the fibres' directional vectors within $B_k$.

It is computed as follows. We denote by $p_m = \omega_m \otimes \omega_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

$$ P_k = \frac{1}{|B_k|} \sum_{m_{s,t} \in B_k} p_m, \quad (B 1) $$

where $|B_k|$ denotes the number of fibres in $B_k$.
<figure><img src="image_15.png" /><figcaption>**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 Ω drawn in black. (b) Same representation, with fibres coloured according to their local alignment indicator (blue: Al<sub>k</sub> = 0, red: Al<sub>k</sub> = 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.</figcaption></figure>

The matrix $P_k$ is symmetric positive-definite, so its three eigenvalues $\lambda_1(P_k)$, $\lambda_2(P_k)$ and $\lambda_3(P_k)$ are real positive. The alignment indicator or fractional anisotropy in the neighbourhood $B_k$ is then equal to

$$ Al_k = \sqrt{\frac{3(\lambda_1(P_k) - \bar{\lambda})^2 + (\lambda_2(P_k) - \bar{\lambda})^2 + (\lambda_3(P_k) - \bar{\lambda})^2}{\lambda_1(P_k)^2 + \lambda_2(P_k)^2 + \lambda_3(P_k)^2}}, \quad (B2) $$

with $\bar{\lambda} = (\lambda_1(P_k) + \lambda_2(P_k) + \lambda_3(P_k))/3$ 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<sub>k</sub> = 0) to red (Al<sub>k</sub> = 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<sub>k</sub> = 1 if all the fibres in $B_k$ have the same directional vector. If the directional vectors are uniformly distributed then theoretically Al<sub>k</sub> = 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
<figure><img src="image_9.png" /><figcaption>**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 Ω drawn in black. (b) Same representation, with fibres coloured according to their local alignment indicator (blue: Al<sub>k</sub> = 0, red: Al<sub>k</sub> = 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.</figcaption></figure>

The matrix $P_k$ is symmetric positive-definite, so its three eigenvalues $\lambda_1(P_k)$, $\lambda_2(P_k)$ and $\lambda_3(P_k)$ are real positive. The alignment indicator or fractional anisotropy in the neighbourhood $B_k$ is then equal to

$$ Al_k = \sqrt{\frac{3(\lambda_1(P_k) - \bar{\lambda})^2 + (\lambda_2(P_k) - \bar{\lambda})^2 + (\lambda_3(P_k) - \bar{\lambda})^2}{\lambda_1(P_k)^2 + \lambda_2(P_k)^2 + \lambda_3(P_k)^2}}, \quad (B2) $$

with $\bar{\lambda} = (\lambda_1(P_k) + \lambda_2(P_k) + \lambda_3(P_k))/3$ 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<sub>k</sub> = 0) to red (Al<sub>k</sub> = 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<sub>k</sub> = 1 if all the fibres in $B_k$ have the same directional vector. If the directional vectors are uniformly distributed then theoretically Al<sub>k</sub> = 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
<figure><img src="image_16.png" /><figcaption>**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.</figcaption></figure>

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/\sqrt{2} \approx 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 $\mathbb{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

$$P_{\text{tot}} = \frac{1}{N_{\text{fib}}} \sum_{1 \le k \le N_{\text{fib}}} \omega_k \otimes \omega_k. \qquad (B3)$$

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 $\omega = (x, y, z)$ on this hemisphere can then be projected onto the equatorial plane via the following transformation:

$$p(\omega) = \left( \frac{x}{1+z}, \frac{y}{1+z} \right). \qquad (B4)$$

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
<figure><img src="image_10.png" /><figcaption>**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.</figcaption></figure>

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/\sqrt{2} \approx 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 $\mathbb{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

$$P_{\text{tot}} = \frac{1}{N_{\text{fib}}} \sum_{1 \le k \le N_{\text{fib}}} \omega_k \otimes \omega_k. \qquad (\text{B3})$$

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 $\omega = (x, y, z)$ on this hemisphere can then be projected onto the equatorial plane via the following transformation:

$$p(\omega) = \left( \frac{x}{1+z}, \frac{y}{1+z} \right). \qquad (\text{B4})$$

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
<figure><img src="image_17.png" /><figcaption>**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.</figcaption></figure>

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.

## References

1. Wohlgemuth R, Brashear S, Smith L. 2023 Alignment, cross-linking, and beyond: a collagen architect's guide to the skeletal muscle extracellular matrix. *Am. J. Physiol. Cell Physiol.* **325**, C1017–C1030. (doi:10.1152/ajpcell.00287. 2023)

2. Mays PK, McAnulty RJ, Campa JS, Laurent GJ. 1991 Age-related changes in collagen synthesis and degradation in rat tissues. Importance of degradation of newly synthesized collagen in regulating collagen production. *Biochem. J.* **276**, 307–313. (doi:10.1042/bj2760307)

3. Piotrowski-Daspit AS, Nerger BA, Wolf AE, Sundaresan S, Nelson CM. 2017 Dynamics of tissue-induced alignment of fibrous extracellular matrix. *Biophys. J.* **113**, 702–713. (doi:10.1016/j.bpj.2017.06.046)

4. Theocharis AD, Skandalis SS, Gialeli C, Karamanos NK. 2016 Extracellular matrix structure. *Adv. Drug Deliv. Rev.* **97**, 4–27. (doi:10.1016/j.addr.2015.11.001)

5. Malandrino A, Trepat X, Kamm RD, Mak M. 2019 Dynamic filopodial forces induce accumulation, damage, and plastic remodeling of 3D extracellular matrices. *PLoS Comput. Biol.* **15**, 1–26. (doi:10.1371/journal.pcbi.1006684)

6. Karamanos NK et al. 2021 A guide to the composition and functions of the extracellular matrix. *FEBS J.* **288**, 6850–6912. (doi:10.1111/ febs.15776)

7. Diller RB, Tabor AJ. 2022 The role of the extracellular matrix (ECM) in wound healing: a review. *Biomimetics (Basel, Switzerland)* **7**, 87. (doi:10.3390/biomimetics7030087)

8. Levental KR et al. 2009 Matrix crosslinking forces tumor progression by enhancing integrin signaling. *Cell* **139**, 891–906. (doi:10.1016/j. cell.2009.10.027)

9. Wilson SL, Guilbert M, Sulé-Suso J, Torbet J, Jeannesson P, Sockalingum GD, Yang Y. 2014 A microscopic and macroscopic study of aging collagen on its molecular structure, mechanical properties, and cellular response. *FASEB J.* **28**, 14–25. (doi:10.1096/fj.13- 227579)

10. Bailey AJ, Paul RG, Knott L. 1998 Mechanisms of maturation and ageing of collagen. *Mech. Ageing Dev.* **106**, 1–56. (doi:10.1016/S0047- 6374(98)00119-5)

11. Drasdo D. 2003 On selected individual-based approaches to the dynamics in multicellular systems. In *polymer and cell dynamics: multiscale modeling and numerical simulations* (eds W Alt, M Chaplain, M Griebel, J Lenz), pp. 169–203. Basel: Birkhäuser.

12. Hwang M, Garbey M, Berceli SA, Tran-Son-Tay R. 2009 Rule-based simulation of multi-cellular biological systems – a review of modeling techniques. *Cell. Mol. Bioeng.* **2**, 285–294. (doi:10.1007/s12195-009-0078-2)

13. Dallon JC, Sherratt JA, Maini PK. 1999 Mathematical modelling of extracellular matrix dynamics using discrete cells: fiber orientation and tissue regeneration. *J. Theor. Biol.* **199**, 449–471. (doi:10.1006/jtbi.1999. 0971)

14. Didonna BA, Levine AJ. 2006 Filamin cross-linked semiflexible networks: fragility under strain. *Phys. Rev. Lett.* **97**, 068104. (doi:10.1103/ PhysRevLett.97.068104)

15. Rizvi MS, Kumar P, Katti DS, Pal A. 2012 Mathematical model of mechanical behavior of micro/nanofibrous materials designed for extracellular matrix substitutes. *Acta Biomater.* **8**, 4111–4122. (doi:10.1016/j.actbio. 2012.07.025)

16. Unterberger MJ, Holzapfel GA. 2014 Advances in the mechanical modeling of filamentous actin and its cross-linked networks on multiple scales. *Biomech. Model. Mechanobiol.* **13**, 1155–1174. (doi:10.1007/s10237-014-0578-4)

17. Barre J, Degond P, Zatorska E. 2016 Kinetic theory of particle interactions mediated by dynamical networks. *Multiscale Model Simul.* **15**, 1294–1323. (doi:10.1137/ 16M1085310)

18. Lieleg O, Claessens MMAE, Bausch AR. 2010 Structure and dynamics of cross-linked actin networks. *Soft Matter* **6**, 218–225. (doi:10. 1039/B912163N)

19. Head DA, Levine AJ, MacKintosh FC. 2003 Distinct regimes of elastic response and deformation modes of cross-linked cytoskeletal and semiflexible polymer networks. *Phys. Rev. E* **68**, 061907. (doi:10.1103/PhysRevE. 68.061907)
<figure><img src="image_11.png" /><figcaption>**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.</figcaption></figure>

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.

## References

1. Wohlgemuth R, Brashear S, Smith L. 2023 Alignment, cross-linking, and beyond: a collagen architect's guide to the skeletal muscle extracellular matrix. *Am. J. Physiol. Cell Physiol.* **325**, C1017–C1030. (doi:10.1152/ajpcell.00287. 2023)

2. Mays PK, McAnulty RJ, Campa JS, Laurent GJ. 1991 Age-related changes in collagen synthesis and degradation in rat tissues. Importance of degradation of newly synthesized collagen in regulating collagen production. *Biochem. J.* **276**, 307–313. (doi:10.1042/bj2760307)

3. Piotrowski-Daspit AS, Nerger BA, Wolf AE, Sundaresan S, Nelson CM. 2017 Dynamics of tissue-induced alignment of fibrous extracellular matrix. *Biophys. J.* **113**, 702–713. (doi:10.1016/j.bpj.2017.06.046)

4. Theocharis AD, Skandalis SS, Gialeli C, Karamanos NK. 2016 Extracellular matrix structure. *Adv. Drug Deliv. Rev.* **97**, 4–27. (doi:10.1016/j.addr.2015.11.001)

5. Malandrino A, Trepat X, Kamm RD, Mak M. 2019 Dynamic filopodial forces induce accumulation, damage, and plastic remodeling of 3D extracellular matrices. *PLoS Comput. Biol.* **15**, 1–26. (doi:10.1371/journal.pcbi.1006684)

6. Karamanos NK et al. 2021 A guide to the composition and functions of the extracellular matrix. *FEBS J.* **288**, 6850–6912. (doi:10.1111/ febs.15776)

7. Diller RB, Tabor AJ. 2022 The role of the extracellular matrix (ECM) in wound healing: a review. *Biomimetics (Basel, Switzerland)* **7**, 87. (doi:10.3390/biomimetics7030087)

8. Levental KR et al. 2009 Matrix crosslinking forces tumor progression by enhancing integrin signaling. *Cell* **139**, 891–906. (doi:10.1016/j. cell.2009.10.027)

9. Wilson SL, Guilbert M, Sulé-Suso J, Torbet J, Jeannesson P, Sockalingum GD, Yang Y. 2014 A microscopic and macroscopic study of aging collagen on its molecular structure, mechanical properties, and cellular response. *FASEB J.* **28**, 14–25. (doi:10.1096/fj.13- 227579)

10. Bailey AJ, Paul RG, Knott L. 1998 Mechanisms of maturation and ageing of collagen. *Mech. Ageing Dev.* **106**, 1–56. (doi:10.1016/S0047- 6374(98)00119-5)

11. Drasdo D. 2003 On selected individual-based approaches to the dynamics in multicellular systems. In *polymer and cell dynamics: multiscale modeling and numerical simulations* (eds W Alt, M Chaplain, M Griebel, J Lenz), pp. 169–203. Basel: Birkhäuser.

12. Hwang M, Garbey M, Berceli SA, Tran-Son-Tay R. 2009 Rule-based simulation of multi-cellular biological systems – a review of modeling techniques. *Cell. Mol. Bioeng.* **2**, 285–294. (doi:10.1007/s12195-009-0078-2)

13. Dallon JC, Sherratt JA, Maini PK. 1999 Mathematical modelling of extracellular matrix dynamics using discrete cells: fiber orientation and tissue regeneration. *J. Theor. Biol.* **199**, 449–471. (doi:10.1006/jtbi.1999. 0971)

14. Didonna BA, Levine AJ. 2006 Filamin cross-linked semiflexible networks: fragility under strain. *Phys. Rev. Lett.* **97**, 068104. (doi:10.1103/ PhysRevLett.97.068104)

15. Rizvi MS, Kumar P, Katti DS, Pal A. 2012 Mathematical model of mechanical behavior of micro/nanofibrous materials designed for extracellular matrix substitutes. *Acta Biomater.* **8**, 4111–4122. (doi:10.1016/j.actbio. 2012.07.025)

16. Unterberger MJ, Holzapfel GA. 2014 Advances in the mechanical modeling of filamentous actin and its cross-linked networks on multiple scales. *Biomech. Model. Mechanobiol.* **13**, 1155–1174. (doi:10.1007/s10237-014-0578-4)

17. Barre J, Degond P, Zatorska E. 2016 Kinetic theory of particle interactions mediated by dynamical networks. *Multiscale Model Simul.* **15**, 1294–1323. (doi:10.1137/ 16M1085310)

18. Lieleg O, Claessens MMAE, Bausch AR. 2010 Structure and dynamics of cross-linked actin networks. *Soft Matter* **6**, 218–225. (doi:10. 1039/B912163N)

19. Head DA, Levine AJ, MacKintosh FC. 2003 Distinct regimes of elastic response and deformation modes of cross-linked cytoskeletal and semiflexible polymer networks. *Phys. Rev. E* **68**, 061907. (doi:10.1103/PhysRevE. 68.061907)
20. Wilhelm J, Frey E. 2003 Elasticity of stiff polymer networks. *Phys. Rev. Lett.* **91**, 108103. (doi:10.1103/PhysRevLett.91.108103)

21. Heussinger C, Frey E. 2006 Floppy modes and nonaffine deformations in random fiber networks. *Phys. Rev. Lett.* **97**, 105501. (doi:10.1103/PhysRevLett.97.105501)

22. Bai M, Missel AR, Levine AJ, Klug WS. 2011 On the role of the filament length distribution in the mechanics of semiflexible networks. *Acta Biomater.* **7**, 2109–2118. (doi:10.1016/j.actbio.2010.12.025)

23. Shahsavari A, Picu RC. 2012 Model selection for athermal cross-linked fiber networks. *Phys. Rev. E* **86**, 011923. (doi:10.1103/PhysRevE.86.011923)

24. Bertin E, Chaté H, Ginelli F, Mishra S, Peshkov A, Ramaswamy S. 2013 Mesoscopic theory for fluctuating active nematics. *New J. Phys.* **15**, 085032. (doi:10.1088/1367-2630/15/8/085032)

25. Boissard E, Degond P, Motsch S. 2013 Trail formation based on directed pheromone deposition. *J. Math. Biol.* **66**, 1267–1301. (doi:10.1007/s00285-012-0529-6)

26. Alonso R, Young J, Cheng Y. 2014 A particle interaction model for the simulation of biological, cross-linked fiber networks inspired from flocking theory. *Cell. Mol. Bioeng.* **7**, 58–72. (doi:10.1007/s12195-013-0308-5)

27. Peurichard D, Delebecque F, Lorsignol A, Barreau C, Rouquette J, Descombes X, Casteilla L, Degond P. 2017 Simple mechanical cues could explain adipose tissue morphology. *J. Theor. Biol.* **429**, 61–81. (doi:10.1016/j.jtbi.2017.06.030)

28. Stylianopoulos T, Barocas VH. 2007 Volume-averaging theory for the study of the mechanics of collagen networks. *Comput. Methods Appl. Mech. Eng.* **196**, 2981–2990. (doi:10.1016/j.cma.2006.06.019)

29. Huisman EM, van Dillen T, Onck PR, Van der Giessen E. 2007 Three-dimensional cross-linked F-actin networks: relation between network architecture and mechanical behavior. *Phys. Rev. Lett.* **99**, 208103. (doi:10.1103/PhysRevLett.99.208103)

30. Åström JA, Kumar PBS, Vattulainen I, Karttunen M. 2008 Strain hardening, avalanches, and strain softening in dense cross-linked actin networks. *Phys. Rev. E* **77**, 051913. (doi:10.1103/PhysRevE.77.051913)

31. Broedersz CP, Sheinman M, MacKintosh FC. 2012 Filament-length-controlled elasticity in 3D fiber networks. *Phys. Rev. Lett.* **108**, 078102. (doi:10.1103/PhysRevLett.108.078102)

32. Harjanto D, Zaman MH. 2013 Modeling extracellular matrix reorganization in 3D environments. *PLoS ONE* **8**, e52509. (doi:10.1371/journal.pone.0052509)

33. Lee B, Zhou X, Riching K, Eliceiri K, Keely P, Guelcher S, Weaver A, Jiang Y. 2014 A three-dimensional computational model of collagen network mechanics. *PLoS ONE* **9**, e111896. (doi:10.1371/journal.pone.0111896)

34. Lin S, Gu L. 2015 Influence of crosslink density and stiffness on mechanical properties of type I collagen gel. *Materials* **8**, 551–560. (doi:10.3390/ma8020551)

35. Picu RC, Sengab A. 2018 Structural evolution and stability of non-crosslinked fiber networks with inter-fiber adhesion. *Soft Matter* **14**, 2254–2266. (doi:10.1039/C7SM02555F)

36. Mak M. 2020 Impact of crosslink heterogeneity on extracellular matrix mechanics and remodeling. *Comput. Struct. Biotechnol.* **J. 18**, 3969–3976. (doi:10.1016/j.csbj.2020.11.038)

37. Cortes DB, Gordon M, Nédélec F, Maddox AS. 2020 Bond type and discretization of nonmuscle myosin II are critical for simulated contractile dynamics. *Biophys. J.* **118**, 2703–2717. (doi:10.1016/j.bpj.2020.03.033)

38. Peurichard D, Ousset M, Paupert J, Aymard B, Lorsignol A, Casteilla L, Degond P. 2019 Extracellular matrix rigidity may dictate the fate of injury outcome. *J. Theor. Biol.* **21**, 127–136. (doi:10.1016/j.jtbi.2019.02.017)

39. Shafer J, Dippel S, Wolf DE. 1996 Force schemes in simulations of granular materials. *J. Phys. I* **6**, 5–20. (doi:10.1051/jpl:1996129)

40. Ciarletta P, Ben Amar M. 2009 A finite dissipative theory of temporary interfibrillar bridges in the extracellular matrix of ligaments and tendons. *J. R. Soc. Interface* **6**, 909–924. (doi:10.1098/rsf.2008.0487)

41. Cox TR, Erler JT. 2011 Remodeling and homeostasis of the extracellular matrix: implications for fibrotic diseases and cancer. *Disease Models Mech.* **4**, 165–178. (doi:10.1242/dmm.004077)

42. Villa C, Chaplain M, Gerisch A, Lorenzi T. 2021 Mechanical models of pattern and form in biological tissues: the role of Stress-Strain constitutive equations. *Bull. Math. Biol.* **83**, 80. (doi:10.1007/s11538-020-00912-5)

43. Ushiki T. 2002 Collagen fibers, reticular fibers and elastic fibers. A comprehensive understanding from a morphological viewpoint. *Arch. Histol. Cytol.* **65**, 109–126. (doi:10.1679/aohc.65.109)

44. Yamada A, Kitagawa K, Nakamura S, Takafuji M, Goto Y, Okamoto R, Dohi K, Sakuma H. 2020 Quantification of extracellular volume fraction by cardiac computed tomography for noninvasive assessment of myocardial fibrosis in hemodialysis patients. *Sci. Rep.* **10**, 15367. (doi:10.1038/s41598-020-72417-5)

45. Dedieu V, Finat-Duclos F, Renou J, Joffre F, Vincensini D. 1999 *In vivo* tissue extracellular volume fraction measurement by dynamic spin-lattice MRI relaxometry: application to the characterization of muscle fiber types. *Invest. Radiol.* **34**, 185–9. (doi:10.1097/00004424-199903000-00005)

46. Purslow PP. 2020 The structure and role of intramuscular connective tissue in muscle function. *Front. Physiol.* **11**, 495. (doi:10.3389/fphys.2020.00495)

47. Xue M, Jackson CJ. 2015 Extracellular matrix reorganization during wound healing and its impact on abnormal scarring. *Adv. Wound Care* **4**, 119–136. (doi:10.1089/wound.2013.0485)

48. Bonnans C, Chou J, Werb Z. 2014 Remodelling the extracellular matrix in development and disease. *Nat. Rev. Mol. Cell Biol.* **15**, 786–801. (doi:10.1038/nrm3904)

49. Jiang YH, Lou YY, Li TH, Liu BZ, Chen K, Zhang D, Li T. 2022 Cross-linking methods of type I collagen-based scaffolds for cartilage tissue engineering. *Am. J. Transl. Res.* **14**, 1146–1159. (doi:10.1103/TranslRes.14.1146)

50. Xue X, Hu Y, Wang S, Chen X, Jiang Y, Su J. 2022 Fabrication of physical and chemical crosslinked hydrogels for bone tissue engineering. *Bioact. Mater.* **12**, 327–339. (doi:10.1165/rcmb.2016-03790C)

51. Philp CJ et al. 2018 Extracellular matrix cross-linking enhances fibroblast growth and protects against matrix proteolysis in lung fibrosis. *Am. J. Respir. Cell Mol. Biol.* **58**, 594–603. (doi:10.1165/rcmb.2016-03790C)

52. Garrison CM, Schwarzbauer JE. 2021 Fibronectin fibril alignment is established upon initiation of extracellular matrix assembly. *Mol. Biol. Cell* **32**, 739–752. (doi:10.1091/mbc.E20-08-0533)

53. Chaudhuri O, Cooper-White J, Janmey P, Mooney D, Shenoy V. 2020 Effects of extracellular matrix viscoelasticity on cellular behaviour. *Nature* **584**, 535–546. (doi:10.1038/s41586-020-2612-2)

54. Winkler J, Abisoye-Ogunlanyi A, Metcalf K, Werb Z. 2020 Concepts of extracellular matrix remodelling in tumour progression and metastasis. *Nat. Commun.* **11**, 5120. (doi:10.1038/s41467-020-18794-x)

55. Mušević I. 2017 Nematic liquid-crystal colloids. *Materials* **11**, 24.

56. Genkin MM, Sokolov A, Lavrentovich OD. Aranson IS. 2017 Topological defects in a living nematic ensnare swimming bacteria. *Phys. Rev. X* **7**, 011029. (doi:10.1103/PhysRevX.7.011029)

57. Berlyand L, Chi H, Potomkin M, Yip NK. 2023 Dynamics and steady state of squirrel motion in liquid crystal. *Eur. J. Appl. Math.* **2023**, 1–42. (doi:10.1017/S0956792523000177)

58. Mirza W, Corato MD, Pensalfini M, Vilanova G, Torres-Sánchez A, Arroyo M. 2023 Theory of active self-organization of dense nematic structures in actin gels. *arXiv - PHYS - Biological Physics* (preprint). (https://arxiv.org/abs/2306.15352)

59. Chassonnery P, Paupert J, Lorsignol A, Séverac C, Ousset M, Degond P, Casteilla L, Peurichard D. 2024 Fibre crosslinking drives the emergence of order in a three-dimensional dynamical network model. *Zenodo*. (doi:10.5281/zenodo.8416498)

60. Chassonnery P, Paupert J, Lorsignol A, Séverac C, Ousset M, Degond P, Casteilla L, Peurichard D. 2024 Fibre crosslinking drives the emergence of order in a three-dimensional dynamical network model. *Figshare*. (doi:10.6084/m9.figshare.c.7054386)

61. Leith D. 1987 Drag on nonspherical objects, aerosol science and technology. *Aerosol Sci. Technol.* **6**, 153–161. (doi:10.1080/02786828708959128)

62. Tirado MM. 2008 Translational friction coefficients of rigid, symmetric top macromolecules. Application to circular cylinders. *J. Chem. Phys.* **71**, 2581–2587. (doi:10.1063/1.4386131)
20. Wilhelm J, Frey E. 2003 Elasticity of stiff polymer networks. *Phys. Rev. Lett.* **91**, 108103. (doi:10.1103/PhysRevLett.91.108103)

21. Heussinger C, Frey E. 2006 Floppy modes and nonaffine deformations in random fiber networks. *Phys. Rev. Lett.* **97**, 105501. (doi:10.1103/PhysRevLett.97.105501)

22. Bai M, Missel AR, Levine AJ, Klug WS. 2011 On the role of the filament length distribution in the mechanics of semiflexible networks. *Acta Biomater.* **7**, 2109–2118. (doi:10.1016/j.actbio.2010.12.025)

23. Shahsavari A, Picu RC. 2012 Model selection for athermal cross-linked fiber networks. *Phys. Rev. E* **86**, 011923. (doi:10.1103/PhysRevE.86.011923)

24. Bertin E, Chaté H, Ginelli F, Mishra S, Peshkov A, Ramaswamy S. 2013 Mesoscopic theory for fluctuating active nematics. *New J. Phys.* **15**, 085032. (doi:10.1088/1367-2630/15/8/085032)

25. Boissard E, Degond P, Motsch S. 2013 Trail formation based on directed pheromone deposition. *J. Math. Biol.* **66**, 1267–1301. (doi:10.1007/s00285-012-0529-6)

26. Alonso R, Young J, Cheng Y. 2014 A particle interaction model for the simulation of biological, cross-linked fiber networks inspired from flocking theory. *Cell. Mol. Bioeng.* **7**, 58–72. (doi:10.1007/s12195-013-0308-5)

27. Peurichard D, Delebecque F, Lorsignol A, Barreau C, Rouquette J, Descombes X, Casteilla L, Degond P. 2017 Simple mechanical cues could explain adipose tissue morphology. *J. Theor. Biol.* **429**, 61–81. (doi:10.1016/j.jtbi.2017.06.030)

28. Stylianopoulos T, Barocas VH. 2007 Volume-averaging theory for the study of the mechanics of collagen networks. *Comput. Methods Appl. Mech. Eng.* **196**, 2981–2990. (doi:10.1016/j.cma.2006.06.019)

29. Huisman EM, van Dillen T, Onck PR, Van der Giessen E. 2007 Three-dimensional cross-linked F-actin networks: relation between network architecture and mechanical behavior. *Phys. Rev. Lett.* **99**, 208103. (doi:10.1103/PhysRevLett.99.208103)

30. Åström JA, Kumar PBS, Vattulainen I, Karttunen M. 2008 Strain hardening, avalanches, and strain softening in dense cross-linked actin networks. *Phys. Rev. E* **77**, 051913. (doi:10.1103/PhysRevE.77.051913)

31. Broedersz CP, Sheinman M, MacKintosh FC. 2012 Filament-length-controlled elasticity in 3D fiber networks. *Phys. Rev. Lett.* **108**, 078102. (doi:10.1103/PhysRevLett.108.078102)

32. Harjanto D, Zaman MH. 2013 Modeling extracellular matrix reorganization in 3D environments. *PLoS ONE* **8**, e52509. (doi:10.1371/journal.pone.0052509)

33. Lee B, Zhou X, Riching K, Eliceiri K, Keely P, Guelcher S, Weaver A, Jiang Y. 2014 A three-dimensional computational model of collagen network mechanics. *PLoS ONE* **9**, e111896. (doi:10.1371/journal.pone.0111896)

34. Lin S, Gu L. 2015 Influence of crosslink density and stiffness on mechanical properties of type I collagen gel. *Materials* **8**, 551–560. (doi:10.3390/ma8020551)

35. Picu RC, Sengab A. 2018 Structural evolution and stability of non-crosslinked fiber networks with inter-fiber adhesion. *Soft Matter* **14**, 2254–2266. (doi:10.1039/C7SM02555F)

36. Mak M. 2020 Impact of crosslink heterogeneity on extracellular matrix mechanics and remodeling. *Comput. Struct. Biotechnol.* **J. 18**, 3969–3976. (doi:10.1016/j.csbj.2020.11.038)

37. Cortes DB, Gordon M, Nédélec F, Maddox AS. 2020 Bond type and discretization of nonmuscle myosin II are critical for simulated contractile dynamics. *Biophys. J.* **118**, 2703–2717. (doi:10.1016/j.bpj.2020.03.033)

38. Peurichard D, Ousset M, Paupert J, Aymard B, Lorsignol A, Casteilla L, Degond P. 2019 Extracellular matrix rigidity may dictate the fate of injury outcome. *J. Theor. Biol.* **21**, 127–136. (doi:10.1016/j.jtbi.2019.02.017)

39. Shafer J, Dippel S, Wolf DE. 1996 Force schemes in simulations of granular materials. *J. Phys. I* **6**, 5–20. (doi:10.1051/jpl:1996129)

40. Ciarletta P, Ben Amar M. 2009 A finite dissipative theory of temporary interfibrillar bridges in the extracellular matrix of ligaments and tendons. *J. R. Soc. Interface* **6**, 909–924. (doi:10.1098/rsf.2008.0487)

41. Cox TR, Erler JT. 2011 Remodeling and homeostasis of the extracellular matrix: implications for fibrotic diseases and cancer. *Disease Models Mech.* **4**, 165–178. (doi:10.1242/dmm.004077)

42. Villa C, Chaplain M, Gerisch A, Lorenzi T. 2021 Mechanical models of pattern and form in biological tissues: the role of Stress-Strain constitutive equations. *Bull. Math. Biol.* **83**, 80. (doi:10.1007/s11538-020-00912-5)

43. Ushiki T. 2002 Collagen fibers, reticular fibers and elastic fibers. A comprehensive understanding from a morphological viewpoint. *Arch. Histol. Cytol.* **65**, 109–126. (doi:10.1679/aohc.65.109)

44. Yamada A, Kitagawa K, Nakamura S, Takafuji M, Goto Y, Okamoto R, Dohi K, Sakuma H. 2020 Quantification of extracellular volume fraction by cardiac computed tomography for noninvasive assessment of myocardial fibrosis in hemodialysis patients. *Sci. Rep.* **10**, 15367. (doi:10.1038/s41598-020-72417-5)

45. Dedieu V, Finat-Duclos F, Renou J, Joffre F, Vincensini D. 1999 *In vivo* tissue extracellular volume fraction measurement by dynamic spin-lattice MRI relaxometry: application to the characterization of muscle fiber types. *Invest. Radiol.* **34**, 185–9. (doi:10.1097/00004424-199903000-00005)

46. Purslow PP. 2020 The structure and role of intramuscular connective tissue in muscle function. *Front. Physiol.* **11**, 495. (doi:10.3389/fphys.2020.00495)

47. Xue M, Jackson CJ. 2015 Extracellular matrix reorganization during wound healing and its impact on abnormal scarring. *Adv. Wound Care* **4**, 119–136. (doi:10.1089/wound.2013.0485)

48. Bonnans C, Chou J, Werb Z. 2014 Remodelling the extracellular matrix in development and disease. *Nat. Rev. Mol. Cell Biol.* **15**, 786–801. (doi:10.1038/nrm3904)

49. Jiang YH, Lou YY, Li TH, Liu BZ, Chen K, Zhang D, Li T. 2022 Cross-linking methods of type I collagen-based scaffolds for cartilage tissue engineering. *Am. J. Transl. Res.* **14**, 1146–1159.

50. Xue X, Hu Y, Wang S, Chen X, Jiang Y, Su J. 2022 Fabrication of physical and chemical crosslinked hydrogels for bone tissue engineering. *Bioact. Mater.* **12**, 327–339.

51. Philp CJ et al. 2018 Extracellular matrix cross-linking enhances fibroblast growth and protects against matrix proteolysis in lung fibrosis. *Am. J. Respir. Cell Mol. Biol.* **58**, 594–603. (doi:10.1165/rcmb.2016-03790C)

52. Garrison CM, Schwarzbauer JE. 2021 Fibronectin fibril alignment is established upon initiation of extracellular matrix assembly. *Mol. Biol. Cell* **32**, 739–752. (doi:10.1091/mbc.E20-08-0533)

53. Chaudhuri O, Cooper-White J, Janmey P, Mooney D, Shenoy V. 2020 Effects of extracellular matrix viscoelasticity on cellular behaviour. *Nature* **584**, 535–546. (doi:10.1038/s41586-020-2612-2)

54. Winkler J, Abisoye-Ogunlanyi A, Metcalf K, Werb Z. 2020 Concepts of extracellular matrix remodelling in tumour progression and metastasis. *Nat. Commun.* **11**, 5120. (doi:10.1038/s41467-020-18794-x)

55. Mušević I. 2017 Nematic liquid-crystal colloids. *Materials* **11**, 24.

56. Genkin MM, Sokolov A, Lavrentovich OD. Aranson IS. 2017 Topological defects in a living nematic ensnare swimming bacteria. *Phys. Rev. X* **7**, 011029. (doi:10.1103/PhysRevX.7.011029)

57. Berlyand L, Chi H, Potomkin M, Yip NK. 2023 Dynamics and steady state of squirrel motion in liquid crystal. *Eur. J. Appl. Math.* **2023**, 1–42. (doi:10.1017/S0956792523000177)

58. Mirza W, Corato MD, Pensalfini M, Vilanova G, Torres-Sánchez A, Arroyo M. 2023 Theory of active self-organization of dense nematic structures in actin gels. *arXiv - PHYS - Biological Physics* (preprint). (https://arxiv.org/abs/2306.15352)

59. Chassonnery P, Paupert J, Lorsignol A, Séverac C, Ousset M, Degond P, Casteilla L, Peurichard D. 2024 Fibre crosslinking drives the emergence of order in a three-dimensional dynamical network model. *Zenodo*. (doi:10.5281/zenodo.8416498)

60. Chassonnery P, Paupert J, Lorsignol A, Séverac C, Ousset M, Degond P, Casteilla L, Peurichard D. 2024 Fibre crosslinking drives the emergence of order in a three-dimensional dynamical network model. *Figshare*. (doi:10.6084/m9.figshare.c.7054386)

61. Leith D. 1987 Drag on nonspherical objects, aerosol science and technology. *Aerosol Sci. Technol.* **6**, 153–161. (doi:10.1080/02786828708959128)

62. Tirado MM. 2008 Translational friction coefficients of rigid, symmetric top macromolecules. Application to circular cylinders. *J. Chem. Phys.* **71**, 2581–2587. (doi:10.1063/1.4386131)