Simple mechanical cues could explain adipose tissue morphology

TL;DR

This paper proposes that simple mechanical interactions between adipocytes and collagen fibers can explain the formation of adipose tissue structures, challenging the idea that vasculature is essential for this process.

Contribution

It introduces a two-dimensional individual-based model demonstrating that mechanical cues alone can lead to adipose tissue morphology without requiring vasculature.

Findings

Structures resemble experimental observations

Cell clustering emerges spontaneously from mechanical interactions

Vasculature not necessary for tissue organization

Abstract

The mechanisms by which organs acquire their functional structure and realize its maintenance (or homeostasis) over time are still largely unknown. In this paper, we investigate this question on adipose tissue. Adipose tissue can represent 20 to 50% of the body weight. Its investigation is key to overcome a large array of metabolic disorders that heavily strike populations worldwide. Adipose tissue consists of lobular clusters of adipocytes surrounded by an organized collagen fiber network. By supplying substrates needed for adipogenesis, vasculature was believed to induce the regroupment of adipocytes near capillary extremities. This paper shows that the emergence of these structures could be explained by simple mechanical interactions between the adipocytes and the collagen fibers. Our assumption is that the fiber network resists the pressure induced by the growing adipocytes and forces them to regroup into clusters. Reciprocally, cell clusters force the fibers to merge into a well-organized network. We validate this hypothesis by means of a two-dimensional Individual Based Model (IBM) of interacting adipocytes and extra-cellular-matrix fiber elements. The model produces structures that compare quantitatively well to the experimental observations. Our model seems to indicate that cell clusters could spontaneously emerge as a result of simple mechanical interactions between cells and fibers and surprisingly, vasculature is not directly needed for these structures to emerge.

Equations

$$\mathcal{W}(C,F)=W_{pot}(C,F)+W_{al}(F),$$

$$W_{pot}(C,F)=\sum _{1\le i\le N_A}\sum _{1\le k\le N_f}{W}_{i,k}(X_i,Y_k,{\theta }_k)$$

$$W_{al}(F)=c_1\sum _{(k,m)\in [1,N_f]}{p}_{km}^t{\sin }^2({\theta }_k-{\theta }_m).$$

$$(C,F)=\underset{\Phi (\tilde{C})\le 0,\Psi (\tilde{F})=0}{\mathrm{argmin}}\mathcal{W}(\tilde{C},\tilde{F}),$$

$$\begin{split}\Phi(C)&=\big{(}\Phi_{ij}(X_{i},X_{j})\big{)}_{(i,j)\in[1,N_{A}]^{2}},\\ \Psi(F)&=\big{(}\vec{\Psi}_{km}(Y_{k},Y_{m},\theta_{k},\theta_{m})\big{)}_{(k,m)\in\mathcal{N}_{f}}.\end{split}$$

$$E=\frac{{\sum }_{c=1}^{N_C}\text{Card}(\mathcal{R}\cap {\mathcal{C}}_c)}{{\sum }_{c=1}^{N_C}{n}_c},$$

$$A=\frac{1}{N_f}\sum _{1\le c_f\le N_{T_f}}{\lambda }_{c_f}^{+}{n}_{c_f},$$

$$\mathcal{W}(C,F)=W_{pot}(C,F)+W_{al}(F),$$

$$W_{pot}(C,F)=\sum _{1\le i\le N_A}\sum _{1\le k\le N_f}{W}_{i,k}(X_i,Y_k,{\theta }_k)$$

$$W_{al}(F)=c_1\sum _{(k,m)\in [1,N_f]}{p}_{km}^t{\sin }^2({\theta }_k-{\theta }_m),$$

$$W_{i,k}=\{\begin{array}{l}\frac{\tilde{W}({\lambda }_k^{+})}{d_{0,i}}(d_{0,i}-d(X,Y,\theta ))\text{ if d(X,Y,θ)≤d0,i}\\ 0\text{ otherwise}\end{array}$$

$$d(X,Y,\theta )=|X-Y+\frac{L_f}{2}\omega (\theta )|+|X-Y-\frac{L_f}{2}\omega (\theta )|-L_f,$$

$$\begin{array}{r}{d}_{0,i}=-L_f+2\sqrt{{\Delta }_i},{\Delta }_i=(\frac{L_f}{2}{)}^2+(\tau R_i{)}^2.\end{array}$$

$$P_k=\frac{1}{n_k}\sum _{m|Y_m\in B(Y_k,R_{al})}{\omega }_m\otimes {\omega }_m,$$

$${\lambda }_k^{+}=\frac{1+\sqrt{\Delta }}{2}$$

$$\begin{split}\Delta&=1+\frac{4}{n_{k}^{2}}\bigg{[}\bigg{(}\sum_{m|Y_{m}\in B(Y_{k},R_{al})}\cos\theta_{m}\sin\theta_{m}\bigg{)}^{2}\\ &-\sum_{m|Y_{m}\in B(Y_{k},R_{al})}(\cos\theta_{m})^{2}\sum_{m|Y_{m}\in B(Y_{k},R_{al})}(\sin\theta_{m})^{2}\bigg{]}\end{split}$$

$$\tilde{W}({\lambda }_k^{+})=(W_1-W_0){\lambda }_k^{+}+W_0,$$

$${\Phi }_{ij}(X_i,X_j)=(R_i+R_j{)}^2-|X_i-X_j{|}^2.$$

$${\vec{\Psi }}_{km}(Y_k,Y_m,{\theta }_k,{\theta }_m)=Y_k+{\ell }_{km}\omega ({\theta }_k)-Y_m-{\ell }_{mk}\omega ({\theta }_m),$$

$${\ell }_{km}=\frac{(x_m^0-x_k^0)\sin {\theta }_m^0-(y_m^0-y_k^0)\cos {\theta }_m^0}{\sin ({\theta }_m^0-{\theta }_k^0)}$$

$$\chi (X,R)=\frac{\sum _{i|X_i\in B(X,R)}\pi R_i^2}{\pi R_{\max }^2{N}_R},$$

$$\mathcal{P}(X,R)={\chi }^{\alpha },$$

$$t_{ins}=\frac{1}{{\nu }_{ins}},t_e=\frac{N_A}{{\nu }_{ins}}.$$

$$R_i^3(t+\Delta t)=R_i^3(t)+K_g(1+\eta {\rho }_g)$$

$$t_g=\frac{R_{\max }^3\Delta t}{K_g}.$$

$$\begin{array}{rl}\mathbb{P}(p_{km}^{t+\Delta t}=1|p_{km}^t=0)& =1-e^{-{\nu }_{\ell }\Delta t}\text{ if max(ℓkm,ℓmk)≤2Lf }\\ & =0\text{ otherwise}\\ \mathbb{P}(p_{km}^{t+\Delta t}=0|p_{km}^t=1)& =1-e^{-{\nu }_d^{km}\Delta t},\end{array}$$

$$t_{\ell }=\frac{1}{{\nu }_{\ell }},t_d=\frac{1}{{\nu }_d}.$$

$${\chi }_{\ell }=\frac{{\nu }_f}{{\nu }_f+{\nu }_d}$$

$$\begin{array}{rl}\mathcal{L}(C,F,\lambda ,\vec{\mu })& =\mathcal{W}(C,F)+{\Phi }_{\lambda }(C)+{\Psi }_{\vec{\mu }}(F),\\ {\Phi }_{\lambda }(C)& =\sum _{1\le i,j\le N_A}{\lambda }_{ij}{\Phi }_{ij}(C)\\ {\Psi }_{\vec{\mu }}(F)& =\sum _{{}_{(k,m)\in {\mathcal{N}}_f}}{\vec{\mu }}_{km}.{\vec{\Psi }}_{km}(F),\end{array}$$

$$X_i^{p+1}=X_i^p-{\alpha }_a^i{\nabla }_{X_i}\mathcal{L}(C^p,F^p,{\lambda }^p,{\vec{\mu }}^p)$$

$$Y_k^{p+1}=Y_k^p-{\alpha }_f^k{\nabla }_{Y_k}\mathcal{L}(C^p,F^p,{\lambda }^p,{\vec{\mu }}^p)$$

$${\theta }_k^{p+1}={\theta }_k^p-{\alpha }_{\theta }^k{\partial }_{{\theta }_k}\mathcal{L}({\mathcal{C}}^p,{\mathcal{F}}^p,{\lambda }^p,{\vec{\mu }}^p)$$

Full text

Author Summary

Because of the key role of adipose tissue in energy homeostasis and associated diseases, there is a great deal of interest in understanding the biology of this tissue. Very little is known about the key to understanding its structuration as lobules. We postulate that lobule emergence is the result of a self-organization process driven by bidirectional mechanical interactions between adipocytes and fibers. We test this hypothesis by means of a 2D individual based model of interacting adipocytes and fiber elements. Indeed, our model produces structures that compare quantitatively well to the experimental observations. This clearly shows that cell clusters of adipose tissue could spontaneously emerge as a result of simple mechanical interactions, with no direct involvement of vasculature.

1 Introduction

White adipose tissue (WAT) is the main energy store of the organism. It is interconnected with all physiological functions via its endocrine functions. It plays a key role in the energy homeostasis and weight of the organism. It is a highly plastic tissue composed of differentiated adipocytes that are able to store and release fatty acids as well as to secrete numerous cytokines and hormones [32]. Mature adipocytes represent only 40 to 60% of the whole cell population. The other cells form a heterogeneous population named the stroma-vascular fraction (SVF). Adipocyte progenitors are present in the SVF throughout adult life [44]. They can proliferate and/or be recruited according to physiological or pathological situations, participate in the turnover of adipocytes and are also believed to be supporting cells. Because of their important role and due to the explosive worldwide development of obesity, the molecular pathways driving adipocyte differentiation are now well investigated and described [10]. In contrast, the global organization at the tissue scale is poorly understood. Since Wassermann’s work in 1960 [47], very few investigations have been performed at this scale. These seminal investigations revealed that adipose tissue is constituted of distinct lobules containing clusters of adipocytes. Moreover, observing its development, Wassermann described the emergence of mature WAT from primitive structures constituted of an unstructured fiber network containing endothelial cells and fibroblast-like cells. The latter are believed to be preadipocytes. In adult adipose tissue, lobules housing adipocytes are separated from each other by well-structured separations (or septa) composed of extracellular matrix (ECM) [31]. Thereafter the number of lobular units seems to remain approximately constant. In excessive development of adipose tissue occurring during obesity, increased fibrosis (formation of excess fibrous tissue) is observed and many reports associate these changes with adipocyte dysfunctions [13, 45]. This suggests that a proper maintenance of adipose tissue architecture is critical for its normal functionality.

Because the global architecture of adipose tissue and its organization into lobules are robust throughout adult life and seem to be fundamental elements of adipose tissue homeostasis, modeling the process of lobule emergence will greatly improve our understanding of adipose tissue biology and plasticity in physiological or pathological conditions. Numerous models of tissue morphogenesis can be found in the literature, describing the emergence of self-organization of cells and fibers. Due to their simplicity and flexibility, the most widely used models are Individual Based Models (IBM) (see [14], [25] and references therein). They describe the behavior of each agent (e.g. a cell or a fiber element) and its interactions with the surrounding agents over time. Due to the high computational cost of IBM, mean-field kinetic or continuous models, which are more efficient to describe the large scales, are often preferred. All these models include one or several of the following interactions: (i) cell/cell, (ii) cell-fiber and (iii) fiber-fiber interactions. Models of interacting cells moving in ECM-free media such as [15] focus on interactions of type (i). Based on the mechanisms reviewed in [16], a wide variety of models incorporating interactions of type (ii) have been proposed such as: (a) mechanical models [30], (b) chemotaxis-type models ([29, 3] and references therein), where cell motion is driven by chemical gradients or (c) models of contact guidance [20] (see [12, 22, 23]) where the ECM gives directional information for cell motion. However, how the processes are coordinated to produce directed motion is not well understood. Fiber-fiber interactions have been explored in [2], where a model of a fibrous network composed of cross-linked fiber elements is proposed. Other authors treat the fibrous network as a continuum, such as a porous medium [39] or an active gel [26] for instance. However, the literature so far provides little clues on the mechanisms underlying contact guidance or fiber self-organization. In the present paper, we demonstrate that directionally organized cell and fiber structures can emerge without appealing to contact guidance or fiber directional interactions, as a result of simple mechanical interactions between the cells and the fiber network. To test this hypothesis qualitatively, it is sufficient in a first instance to consider a two-dimensional model, which we do for reasons of simplicity and computational efficiency. Our model is of more microscopic nature than previously proposed mechanical models [30, 37] and aims at describing the emergence of the lobular structures observed in adipose tissue. Our goal is to test the scenario that, due to the fibers resisting the pressure induced by the growing adipocytes, the latter are forced to regroup into clusters. Adipocyte clusters in turn force the fibers to merge into a well-organized network. To validate this scenario, we developed a two-dimensional IBM modelling adipocytes interacting with ECM fiber elements. The model and experiment data showed strikingly similar lobule-like structures and revealed that vasculature was not needed for lobular self-organization to emerge, although vasculature is a proven key factor of adipogenesis [7].

2 Material and methods

2.1 Experiments and image processing

Fig. 1 shows part of a sub-cutaneous adipose tissue from an adult mouse. The adipocytes were visualized by immuno-staining of perilipin, a protein that surrounds the unilocular lipid droplets. In Fig. 1A, the ECM fiber network appears in black as a background. This picture clearly reveals the organization of the adipose tissue in lobules. A 3D image of one isolated lobule is shown in supplementary information. We implemented classical image processing methods to extract the centers and radii of the cells (Fig. 1B) and the different lobules (Fig. 1C) (see D for details). The quality of the cell and lobule segmentation methods were carefully checked.

2.2 Description of the model

We postulated that, in WAT, the agents contributing the most to mechanical balance were the ECM fibers and the adipocytes. The ECM was discretized into unit fiber elements consisting of line segments of fixed and uniform lengths represented by their centers and their directional unit vectors. We supposed that two fiber elements that crossed each-other could form a link, thereby creating a longer fiber. The cells were described as 2D spheres represented by their centers and radii. At any given time, the two sets of agents were supposed to realize the minimum of the mechanical energy of the system (described below). We incorporated the following biological features : (i) Pre-adipocyte differentiation: Immature cells are much smaller than adipocytes. So, we supposed that they had negligible impact on the mechanical equilibrium and we did not incorporate them in the model. The transformation of an immature cell into an adipocyte was modelled as the creation (or “insemination”) of a new adipocyte. All new adipocytes were inseminated with the same small radius. New adipocytes were inseminated at random times following a Poisson process. The location of the insemination was also random with either uniform probability in the domain (below referred to as “random insemination”) or with a bias resulting in a higher insemination probability at locations where existing adipocytes were already present (below referred to as “biased insemination”). In this biased insemination case, the existing cell density in a disk of radius $R$ around the randomly chosen insemination point $X$ was computed and normalized by the maximal possible density (corresponding to adipocytes in contact with each other), resulting in a dimensionless parameter $\chi$ comprised between [math] and $1$. Then, the insemination probability at $X$ was taken proportional to $\chi^{\alpha}$, with biasing parameter $\alpha>0$. In this scenario, a pre-adipocyte sensed the adipocyte density $\chi$ up to a sensing distance $R$ and made a decision whether to differentiate into an adipocyte according to the value of the parameter $\chi^{\alpha}$ (the larger $\alpha$, the larger the local adipocyte density needed to be to trigger differentiation). By correlating the new adipocytes location to existing ones, biased insemination indirectly accounted for a vascular network bringing blood supplies in given locations of the tissue to trigger the appearance of new adipocytes (see Section ’Discussion’ below). (ii) Adipocyte growth: The ability to store and release energy according to the needs of the organism was modelled through the regular growth of the cells. Therefore, thanks to (i) and (ii), we incorporated both hyperplasia (cell number increase) and hypertrophy (cell size increase). As the turnover of adipocytes is small and not related to adipose tissue morphology [4], we neglected the apoptosis of adipose cells. We assumed that the volume of each adipocyte reached a maximal value beyond which it stayed constant. (iii) Adipocyte incompressibility and non-overlapping: Adipocytes are reservoirs of fat, which is an incompressible liquid, and they cannot overlap. Therefore, we assumed that the radius of each disk was unaffected by whatever mechanical efforts were exerted onto it and that two neighboring disks could not overlap. (iv) Fiber resistance to adipocyte pressure: We viewed the ECM as a soft medium composed of fibers having locally preferred directions. Each fiber element represented the elementary mechanical action that the ECM exerted on neighbouring cells. It carried a directional information corresponding to the preferred direction of the collagen fibers. The largest mechanical action of the fiber was exerted normally to this direction. The fibers repelled the cells by means of a soft potential allowing some penetration of the cells inside the ECM. A fiber element produced a unit of potential strength. Larger mechanical actions were achieved by having several of these unit fiber elements in a close neighbourhood. Obviously, there exists a limit to ECM strength, but we assumed that our system operated below that limit and we did not implement any upper bound on the number of fiber elements per unit area. In practice, we assumed that the fiber-cell repulsion potential iso-lines were ellipses with focii at the two ends of the fiber segment and that the potential vanished beyond a certain distance from the fiber. The choice of elliptic isolines corresponded to the simplest anisotropic shape in 2D. (v) Fiber growth, elongation and ability to bend: In addition to carrying a unit of ECM fiber strength, fiber elements also carried a unit of fiber length. However, we provided a way to create longer fibers by allowing two crossing fibers to create a cross-link. Any displacement of a cross-linked fiber pair would maintain the position of the cross-link relative to the center of each fiber. Several consecutively cross-linked fiber elements would model a long fiber having the ability to bend or even take possible tortuous geometries. As the number of fibers cross-linked to a given fiber was not limited, we could also account for fiber branching. Therefore, the cross-linking process would model fiber elongation [6] and symmetrically, pairs of cross-linked fibers could also spontaneously unlink, allowing for fiber breakage describing ECM remodelling processes. Linking and unlinking processes followed Poisson processes with frequencies $\nu_{\ell}$ and $\nu_{d}$ respectively, and the parameter $\chi_{\ell}=\frac{\nu_{\ell}}{\nu_{\ell}+\nu_{d}}$, where $\chi_{\ell}\in[0,1]$ represented a measure of the fraction of linked fibers among the pairs of intersecting fibers. (vi) Fiber alignment and resistance to bending: To model the ability of long fibers (those made of several cross-linked fiber units) to offer a certain resistance to bending, linked fibers were subjected to a potential torque at their junction. This torque vanished when the fibers were aligned, and consequently acted as a linked-fiber alignment mechanism. Any force exerted by a cell in the vicinity of a cross-link would result in the cross-linked fibers making an obtuse angle with respect to one another until the exerted torque at the cross-linked balanced the effect of the force, thereby providing a discrete representation of the bending of a continuous beam. This torque was characterized by a stiffness parameter $c_{1}>0$ playing the role of a flexural modulus. The larger the $c_{1}$ the more rigid the fiber assembly was.

The mechanical energy of the system included the cell-fiber interaction potentials (iv) and the linked fiber-fiber alignment potential (vi). At each time step, a minimum of this mechanical energy subject to the nonoverlapping constraint between cells (iii) and to the linkeage constraint between linked fibers (v) was sought. At the beginning of the next time step, new adipocytes were inseminated (i), adipocyte radii were increased (ii) and pairs of fibers were linked/unlinked (v). These phenomena disrupted the mechanical equilibrium and a minimum of this new energy was again sought, and so on.

More specifically, the $N_{A}$ adipocytes were modeled as 2D growing spheres of center $X_{i}$ and radius $R_{i}$ for $i$ in $[1,N_{A}]$, and the $N_{f}$ ECM fiber elements were represented by straight lines of fixed length $L_{f}$, of center $Y_{k}$ and orientational angle $\theta_{k}$ for $k\in[1,N_{f}]$. Given a configuration at a fixed time, we denoted by $C$ the set of cell center positions and radii: $C=\{(X_{i},R_{i})\;,\;i\in[1,N_{A}]\}$ and by $F$ the set of fiber center positions and fiber directional angles: $F=\{(Y_{k},\theta_{k})\;,\;k\in[1,N_{f}]\}$. The global mechanical energy of the system was written:

[TABLE]

where $W_{pot}$ and $W_{al}$ were the fiber-cell repulsion potential (iv) and the linked fiber-fiber alignment potential (vi) respectively:

[TABLE]

The time-dependent coefficients $p^{t}_{km}$ were set to 1 if fibers $k$ and $m$ were linked at time $t$ and to 0 otherwise (see B.3). The alignment potential $W_{al}$ of intensity $c_{1}$ consisted of the sum of elementary alignment potentials between fibers of a linked pair. The repulsion potential $W_{pot}$ was supposed to be the sum of two-particle potential elements, $W_{i,k}$, modeling the mechanical interaction between cell $i$ and fiber $k$. The iso-lines of these potential elements consisted of ellipses with focii located at the two ends of the fiber segment (see Fig. 2). We supposed that fiber $k$ could repel cell $i$ up to a distance $\tau R_{i}$ in its orthogonal direction, where $R_{i}$ is the radius of the cell and $\tau$ a model parameter. We stress the fact that fiber-cell repulsion aimed to model repulsion of the fiber by the border of the cell, while a cell was modeled by its center. Therefore, the repulsion distance $\tau R$ of the model corresponded to a repulsion distance $(\tau-1)R$ from the border of the cell.

At each time step, the minimization of the global mechanical energy under the constraints (iii) and (v) was written:

[TABLE]

where $\Phi(\tilde{C})\leq 0$ expressed cell-cell nonoverlapping inequality constraints (iii) and $\Psi(\tilde{F})=0$ expressed fiber-fiber linkeage equality constraints (v) (see B.2), where:

[TABLE]

Here, $\mathcal{N}_{f}$ denoted the set of linked fiber pairs: $\mathcal{N}_{f}=\{(k,m)\in[1,N_{f}]^{2},\,k<m,\,p_{km}^{t}=1\}$. We send the reader to B.4 for more details on the numerical algorithm used for solving the minimization problem. The model was implemented on a 2D square domain and boundary conditions were assumed periodic (i.e. each agent was assumed periodically repeated beyond the boundary of the square domain). The numerical simulations were initialized with a randomly distributed fiber network (according to a uniform distribution over all possible direction angles or over a sub-interval of directions angles centered about a given angle $\theta_{1}$ and of width $2\theta_{2}$). New cells were inseminated at a constant rate until reaching a cell volume fraction of 50%, a number consistent with the experimental observations. We refer the reader to B for a more complete description of the model ingredients.

2.3 Biological relevance of the model parameters

The reference time was chosen to be 10 times the insemination time $t_{ins}$ (inverse of insemination rate), i.e $t_{ref}=10t_{ins}$, and the time step was chosen to be $\Delta t=t_{ins}=0.1t_{ref}$. Simulations were stopped at dimensionless time $t=100t_{ref}$. The time to reach the maximal number of cells corresponded to $t_{e}=18t_{ref}$ and that for cells to reach their maximal sizes to $t_{g}=18t_{ref}$. Simulations were run five times the time needed for all cells to reach their maximal sizes to ensure that an equilibrium was attained in the end. The fiber unlinking time $t_{d}=1/\nu_{d}$ varied between $t_{d}=10t_{ref}$ and $t_{d}=10^{4}t_{ref}$, and the linking time $t_{\ell}=1/\nu_{\ell}$ ranged from about twice to ten times $t_{d}$ according to the value of $\chi_{\ell}$. Note that these are linking/unlinking times per fiber but the actual frequencies of linking/unlinking events must be multiplied by the number of fibers $N_{f}$ ($N_{f}=800$ in our simulations, see B.3 for details on the time scales). It is difficult to relate these time scales to actual biological time scales because adipocyte growth rate is highly variable according to the subject, its age, its diet, its energy consumption, etc. Nonetheless, adipocyte turnover is low [4] and so we can estimate that it takes about 100 days for a nascent adipocyte to grow to its maximum size. In our model, a nascent adipocyte needs about 20 time units to grow, so, it is consistent to assume that one time unit of our model is about five days. In our model, we can estimate that there is significant ECM remodeling when at least $10\%$ of the fibers (i.e about 100 fibers) have been linked or unlinked. Therefore, the ECM remodeling rate can be related to $100\,\nu_{d}$ and consequently, the ECM remodeling time $t_{d}/100$ ranges between $0.5$ and $510^{2}$ days. In [8], it is estimated that ECM remodeling takes about 15 days which falls between these two bounds.

By estimating the adipocyte radius to $30\mu m$ from the experiments (see Fig. 1), the fiber elements considered here were $60\mu m$ large and $140\mu m$ long. In [46], it is shown that the collagen fibers are organised into bundles of collagen fibrils. The bundle width varies from $1$ to $20\mu m$ depending on tissues and organs, and in loose connective tissues such as adipose tissue these bundles can run parallel to each other to be twined into larger bundles. As for the length of the fiber elements, we note that actual collagen fibers are able to bend around the cells. As in our model fiber elements are rigid, the bending ability of the fibers is restored by connecting fiber elements together without forcing them to be aligned (the finite value of the flexural modulus $c_{1}$ allows for an angle between two connected fiber elements to appear). However, if the fiber elements are too long, the curvature of the resulting fiber may be too restricted. In order to describe fibers able to bend around the cells, the fiber element length should be about a cell diameter, i.e 60 $\mu$m. In F, we report on simulations using the values of fiber element length of 54 $\mu$m and width 10$\mu$m consistent with these considerations. However, these values lead to computationally intensive simulations and for the purpose of building a phase diagram of the system, we instead have used fiber element length 140$\mu$m and width 60 $\mu$m that lead to tractable simulations. If not otherwise stated, the model parameters used for the simulations are given by table 3, the values of the numerical parameters are listed in table 11 and more information on the values of the parameters for the biological phenomena can be found in B.

The relevance of these parameters is assessed inF where we observe that the results with the more biologically relevant parameter values are similar to those obtained with the larger values, thus justifying the use of the latter for an extensive exploration of the parameter space.

3 Results

3.1 Influence of the model parameters

Fig. 4 (I) shows a phase diagram representing the various morphologies obtained using random adipocyte insemination and varying the linking/unlinking frequency $\nu_{d}$ (horizontal axis of the diagram) as well as the linked fiber fraction $\chi_{\ell}$ (vertical axis of the phase diagram). It illustrates that the linking-unlinking dynamics strongly influenced the morphology of the final structures. Indeed, varying $\nu_{d}$ and $\chi_{\ell}$, we observed that the cell cluster morphology changed from compactly shaped (Fig. 4 (I 1) or (I 4)) to elongatedly shaped (Fig. 4 (I 3) or (I 6)), while the fiber cluster morphology changed from disordered (Fig. 4 (I 1)) to aligned (Fig. 4 (I 2)). Finally the fibers self-organized into rigid long fiber threads (Fig. 4 (I 3)). For a slow linking-unlinking process (Fig. 4 (I 1)), the fiber network topology was almost frozen and this generated rigidly connected fiber structures that were too stiff to self-organize (see H.1 and H.2). A faster linking-unlinking process (Fig. 4 (I 2)) allowed for the remodeling of the network topology and for the local alignment of the fibers thanks to the alignment torque generated at the links (see H.3).

However, for a very fast linking-unlinking dynamics (Fig. 4 (I 3)), the aligned fiber threads were reinforced by the fast creation of links, increasing the rigidity of the network (see H.4). A preferred fiber direction locally emerged and favored the growth of cell clusters in that particular direction, thereby generating elongated cell clusters. For a smaller linked fiber fraction $\chi_{\ell}=0.1$ (Fig. 4 (I 5) or (I 6)), the fiber network was less rigid compared to the previous case. Fiber clusters were consequently smaller, and failed to surround the cell structures, generating bigger cell clusters.

In Fig. 5(A) and (B), we show simulation results obtained with the larger fraction of linked fibers $\chi_{\ell}=0.35$ and for $\nu_{d}=10^{-3}$ (A) and $\nu_{d}=210^{-2}$ (B). Similarly, Fig. 5(C) and (D) show simulation results using biased insemination with biasing parameter $\alpha=0.1$. Again, the larger linked fiber fraction $\chi_{\ell}=0.35$ is used, and two values of $\nu_{d}$ are investigated: $\nu_{d}=10^{-3}$ (C) and $\nu_{d}=210^{-2}$ (D). It shows that biased insemination did not have a significant influence on the final cell and fiber structures. For a well calibrated fiber linking-unlinking dynamics, the model was able to produce lobule-like structures without the need of biased insemination (see E.3 for more details).

3.2 Quantitative analysis

In order to quantify cell and ECM fiber structures, we defined a set of statistical descriptors (SQ). A cell cluster was defined as a group of at least $5$ adipocytes in contact with each other, and $N_{C}$ was the number of cell clusters per 100 adipocytes. Parameter $E$ measured the average elongation of cell clusters. It was comprised between [math] (for disk-like cell clusters) and $1$ (for cord-like clusters) and had the expression:

[TABLE]

where $\mathcal{R}\cap\mathcal{C}_{c}$ is the set of cells with less than five neighbors in cluster $c$, and $n_{c}$ its total number of cells (see C for more details). We verified that our conclusions did not depend on the chosen minimal size (here $5$) of the cell clusters. Due to computational constraints, we characterized cell clusters by their shape rather than by their number of cells. Indeed, comparing the number of cells with experimental data would have required to simulate a larger number of cells, which was computationally too intensive. The SQ $\Theta$ measured the standard deviation of the shape anisotropy direction of cell clusters. For this purpose, each cell cluster was best-matched to an ellipse (again because ellipses are the most simple anisotropic shapes in 2D) and the cluster shape anisotropy direction was defined as the angle of the ellipse semi-major axis with a reference direction. Small values of $\Theta$ indicated a preferred shape anisotropy direction of cell-clusters. Moreover, a fiber cluster was defined as a set of neighboring quasi-aligned fiber elements, and parameter $A$ measured the mean alignment of the fibers of the clusters. Parameter $A$ was comprised between 0 and 1: small values of $A$ indicated an isotropy in the fiber orientation angle distribution (’disorganized fiber network’) while large values of $A$ indicated a global alignment of the fibers in the clusters (’organized fiber network’). The expression of $A$ was given by:

[TABLE]

where $N_{f}$ is the total number of fibers, $\lambda_{c_{f}}^{+}$ is the mean alignment of fibers in cluster $c_{f}$ and $n_{c_{f}}$ its number of fibers. We refer the reader to C for details on the computation of the SQ.

3.3 Identification of different morphologies

For each set of model parameters, we computed the SQ on the obtained final structures and averaged them over $10$ realizations. In Fig. 4 (II), we plotted the mean alignment of fiber clusters $A$ (in orange) and the mean cell cluster elongation $E$ (in black) as functions of the fiber unlinking frequency $\nu_{d}$ for two values of the linked fiber fraction $\chi_{\ell}=0.1$ (Fig. 4 (II A)) and $\chi_{\ell}=0.35$ (Fig. 4 (II B)) (see E.1 for a complete diagram scanning the values of $\chi_{\ell}\in[0.1,0.35]$ and for the SQ $\Lambda$ and $N_{C}$). These plots revealed an increase of the cell cluster elongation $E$ as the unlinking frequency $\nu_{d}$ increases, whatever value the linked fiber fraction $\chi_{\ell}$ took. We identified two plateaus of values of $E$: $E\approx 0.55$ for $\nu_{d}\leq 10^{-2}$ and $E\approx 0.75$ for $\nu_{d}\geq 10^{-1}$. The corresponding cell structures were compactly shaped (Fig. 4 (I A) (I B), (I D) or (I E)) or elongatedly shaped (Fig. 4 (I C)) respectively. Fig. 4 (II) reveals that the mean alignment $A$ of fiber clusters has a plateau value for $\nu_{d}\leq 10^{-2}$. For a slow fiber linking-unlinking process ($\nu_{d}\leq 10^{-2}$), the fibers in the clusters are poorly aligned (low value of $A$), and the mean alignment of fiber clusters increases with $\nu_{d}$ from the value $A\approx 0.6$ to $A\approx 0.8$ as $\nu_{d}$ increases in $[10^{-2},1]$. This tends to show that the fiber linking-unlinking dynamics strongly influences the final structures. For a slow fiber linking-unlinking dynamics, due to fiber interconnections the fiber network is extremely rigid. In this case, fibers fail to align because of the high connectivity of the network which prevents any configurational change. As $\nu_{d}$ increases, the lifetime of each link decreases, allowing for the remodeling of the fiber network. Consequently, fiber structures are more flexible and more aligned. If the linking-unlinking process is too fast, the fiber structures easily align, forming fiber threads. These fiber threads are then reinforced by the fast creation of links, increasing the rigidity of the fiber patterns. Preferred directions locally emerge in the fiber network and favor the growth of cell clusters in these directions. The cell clusters are consequently elongated.

We therefore identified three different morphologies -which will be referred to as ’phases’- according to the values of the parameters: (A) Lobule-like cell clusters surrounded by a disorganized (unaligned) fiber network, (B) lobule-like cell structures in an aligned fiber network, and (C) elongated cell structures in a network composed of long and rigid fiber threads. Each type of structure is obtained in a specific range of the parameters $\nu_{d}$ and $\chi_{\ell}$ of the fiber linking-unlinking dynamics. We chose the mean fiber cluster alignment $A$ and the mean cell cluster elongation $E$ to quantify the passage from one morphology to another, and we identify the threshold values $A^{*}=0.68$ for $A$ and $E^{*}=0.65$ for $E$. Structures of type (A) correspond to $A<A^{*}$ and $E<E^{*}$. Type (B) is described by $A>A^{*}$ and $E<E^{*}$ and finally type (C) by $A>A^{*}$ and $E>E^{*}$. Fig. 4 (I) shows a phase diagram representing the various obtained morphologies in the $(E,A)$ plane. Each point in this phase diagram has been obtained by averaging the SQ over 10 simulations. The separatrix between phases (A) and (B) (i.e. the line $A=A^{*}$) is shown in red and the separatrix between phases (B) and (C) (i.e. the line $E=E^{*}$) is shown in blue. For each region, we show simulation results in the corresponding range of the parameters $\nu_{d}$ and $\chi_{\ell}$ for the sake of illustration.

The most biologically relevant structures, composed of well-separated lobule-like cell clusters in an organized fiber network, were characterized by a low value of the mean elongation $E$ and a moderate value of the fiber mean alignment $A$. From Fig. 4 I and II, we realized that the best combination of parameters was a linked fiber fraction of $\chi_{\ell}=0.35$ and an unlinking frequency $\nu_{d}\in[10^{-2},10^{-1}]$. For these parameters, the model produced biologically relevant cell and fiber structures (Fig. 4 (I 2)), without the need for biased insemination. We refer to E.3 for a discussion of the influence of the biasing parameter $\alpha$ and to E.2 for the role of the alignment force $c_{1}$.

In F, we have explored another range of parameters corresponding to those discussed at the end of Section 2.2 above. We have shown that they give rise to similar structures and phase transitions as in the case described here. Specifically, we have used shorter fiber elements, of the order of a cell diameter. The fiber elements were made thinner in the same proportion to keep the overall shape of their interaction potential similar. Simultaneously, a decreased interaction distance between cells and fibers was used, merely reducing to a contact interaction. Finally, a fiber-fiber interaction potential, of the same strength and interaction distance as the cell-fiber interaction potential was added. To mimic fibers initially made of several fiber elements, the latter were initially placed by groups of interlinked four elements. Within this range of parameter values, we obtained very similar structures as those generated with larger fibers, larger cell-fiber interaction distance and no fiber-fiber repulsion. These results support the conclusion that the observed phenomena are generic and not tightly connected with specific parameter values. Due to computational constraints, in the simulations reported in the present section, larger values of the fiber length and width were used allowing for a smaller number of fiber elements to be simulated. This led to computationally tractable simulations which could be exploited to explore the parameter space, generate the phase diagram and calibrate the parameters for comparisons with experimental data. Finally, the convergence time to equilibrium and the influence of the final adipocyte number have been explored. Studies on the convergence time showed that the SQ remained stable after an initial transient of less than 100 times the insemination time, which showed that the analyzed structures had correctly reached an equilibrium. Few quantitative differences were recorded when the final number of adipocytes was changed but more importantly, the important tendencies such as variations of the SQ’s as functions of the model parameters remained the same whatever final adipocyte number was chosen.

3.4 Comparisons with experimental data

The image processing enabled the computation of the SQ $N_{C}$, $E$ and $\Theta$ on the biological images, thus allowing for a quantitative comparison between the model and the experimental results. The data were compared with simulation results at equilibrium only. Indeed, to date, the in-vivo registration and tracking of cell positions and ECM location during the development of adipose tissue is not feasible. To compare the biological and numerical SQ, we non-dimensionalized the mean cluster number $N_{C}$ and the mean elongation of cell clusters $E$ by reference values referred to as $N_{ref}$ and $E_{ref}$ respectively. The numerical (respectively biological) reference values $N_{ref}$, $E_{ref}$ were defined as the mean value of $N_{C}$ or $E$ over all the numerical (respectively biological) experiments. As biological images suggest that parts of adipose tissue exhibit a preferred direction, we ran simulations for different initial fiber configurations such that the initial fiber angles $\theta_{init}$ were uniformly chosen in the range $[\theta_{1}\pm\theta_{2}]$. Parameter $\theta_{1}$ measured the mean initial fiber direction and $\theta_{2}$ was related to its standard deviation. The larger $\theta_{2}$, the more disordered the initial network was (with $\theta_{2}=\pi$ as the extreme case where the fiber initial distribution was fully isotropic). We also considered the case where the mean fiber direction could depend on the position in the form $\theta_{1}(x_{1},x_{2})=\theta_{1}^{+}$ if $x_{1}>0$ and $\theta_{1}(x_{1},x_{2})=\theta_{1}^{-}$ if $x_{1}<0$ where $(x_{1},x_{2})$ are the coordinates of a point in the computational domain and $x_{1}=0$ corresponds to the vertical middle line. This corresponded to the case where the initial fiber mean direction changed abruptly from $\theta_{1}^{-}$ to $\theta_{1}^{+}$ across the vertical middle line. We refer to E.4 for a thorough analysis of the influence of the initial network anisotropy on the shape of the final structures. Here for the purpose of comparing with actual data, we realized a database containing, for each set of model parameters $(\nu_{d},\theta_{1},\theta_{2})$ (or in the case of position-dependent initial mean fiber direction, $(\nu_{d},\theta_{1}^{-},\theta_{1}^{+},\theta_{2})$), the SQ $N_{C}$, $E$ and $\Theta$ for $10$ different simulations. The simulations were performed with random insemination and the other parameters were chosen to the value $\chi_{\ell}=0.35$ and $c_{1}=1$, according to the previous analysis. For each biological image, we first searched the database of numerical simulations to find the combination of parameters which minimized the quadratic difference between the experimental SQ and the mean of the model SQ. Then, within the $10$ simulations generated for this set of parameters, we selected the one minimizing the quadratic difference between the experimental and model SQ.

Fig. 6 (A) to (C) show three biological images before (I) and after (II) lobule detection, as well as the corresponding best simulation (III) applying the previously detailed method and the associated set of parameters $\nu_{d}$, $\theta_{1}$ (or $\theta_{1}^{\pm}$) and $\theta_{2}$. Table (D) provides the SQ corresponding to images (A) to (C) for both the experimental and numerical images. Finally, in Fig. 6 (E), a larger portion of the adipose tissue is shown (the white scale bar at the bottom right indicates $500\,\mu m$ for (E) and $100\,\mu m$ for (A) to (C)).

The lobule segmentation (column (II)) revealed that the elongation of the lobular structures increased from (A) to (C). We notice that the value of the unlinking frequency $\nu_{d}$ corresponding to the best numerical fit increases as well, confirming the analysis made above. In these three cases, the values of the SQ given by the model are very close to the experimental ones (see Table (D)). These results show that the model is able to reproduce the data in a fairly wide range of situations by simply modifying the model parameters. In real tissues the organizational level varies from high to low according to the position in the tissue. In Fig 6 (E), a large portion of adipose tissue is shown. It reveals well-organized lobular structures in its upper part and more disordered structures in its lower part. Simulation of the entire tissue using the model would be possible (although computationally intensive) by simply varying the parameters of the model to match the variation of the organizational level of the tissue.

Discussion

To our knowledge, this work is the first attempt to understand the emergence of the lobular structure of the adipose tissue by interfacing a mathematical model and experimental results. The originality of our work lies in the assumption that adipose tissue architecture results from a self-organization process principally driven by mechanical interactions between adipocytes and the ECM. This corresponds to a co-organization where the cell clusters and the fiber structures evolve simultaneously. Our mathematical model is able to reproduce the clustering of adipocytes into lobular units surrounded by the ECM fiber network. Simply varying a few parameters allowed us to span a large variety of morphologies. Our results suggest that adipose tissue organization could be principally driven by mechanical cues, in addition to a limited number of biologically-controlled phenomena such as fiber-fiber chemical linking. They highlight the importance of the ECM and the mechanical forces induced by it. They are consistent with recent papers demonstrating the impact of the ECM and of its dysregulation (i.e. fibrosis) on adipose tissue function [28, 13, 33].

Furthermore, the model can reproduce many different types of tissue morphologies by varying the parameters controlling the mechanical response of the ECM to the growth of adipocyte number and size. This gives us good confidence that although validated on mice, the model will be able to reproduce human adipose tissue. Indeed, there is increasing evidence that the differences between adipose tissue structures observed in different species or even at different locations within the same species can be related to mechanical pressure, adipocyte sizes and ECM characteristics. For instance, [38] proposes to define different types of subcutaneous WAT in humans based on their observed structural and ultrastructural features. They describe lobular subtypes for fibrous white adipose tissue that can be found in regions where mechanical constraints are large. Additionally they distinguish lobular and non lobular adipose tissue on the basis of the nature and richness of the ECM. All these differences in the mechanical and structural properties of the adipocytes and the ECM can be fed into the model to produce a wide array of different morphologies.

The structures that emerged from the mathematical model can be classified into three types: (a) middle-sized compact cell clusters surrounded by a disorganized fiber network, (b) middle-sized compact cell clusters surrounded by a well-organized network of thick fiber threads, or (c) elongated clusters surrounded by a network of thin and rigid fiber threads. Each type of structure corresponded to a range of model parameters of the fiber linking-unlinking dynamics. Structures of type (a) were obtained for a slow fiber linking-unlinking dynamics. The rigidly connected fiber structures could not self-organize, leading to a disordered fiber network. However, the system was able to produce middle sized cell clusters of lobular shape as a result of cell-fiber repulsion. This reflected the ability for a connected fiber network to exert a pressure on the cell structures and confine them into separated zones. For a faster fiber linking-unlinking dynamics, the remodeling of fiber structures was enabled, and fibers could arrange more easily into organized patterns, thanks to the torque acting on linked fibers. For a well chosen range of the unlinking frequency $\nu_{d}$ and of the flexural modulus $c_{1}$ (see discussion of $c_{1}$ in E.2), the model produced biologically relevant cell and fiber structures of type (b). Finally, structures of type (c) were observed for fast fiber linking-unlinking dynamics. In this case, fibers easily aligned with each other and fiber stiffness was reinforced by the links and associated torque acting on linked fibers. This imposed local directional constraints to cell cluster growth, favoring cell cluster elongation. Moreover, due to increased fiber rigidity, the fibers failed to surround the cell structures.

The present results, which show a segregation of cells and fibers and the emergence of alignment (aka nematic order) among the fibers is consistent with experimental and theoretical work about mixtures of isotropic and nematic particles. For instance, [24] examines a mixture of F actin and inert polymer Polyethylene glycol (PEG) and observes that actin filaments bundle together when a threshold concentration of PEG is reached. In [18], vibrofluidized steel rods in the presence of high concentrations of hard spheres self assemble into linear structures under similar conditions. These observations are consistent with Monte-Carlo simulations of mixtures of hard sphero-cylinders and spheres peformed in [27]. In these references, the coalescence and nematic ordering of the rods is analyzed in terms of the depletion force stemming from volume exclusion effects. The phenomena observed in our simulations bear evident analogies with these observations, even though the fibers and cells do not interact via volume exclusion but through a soft repulsive potential. The difference is important in that the interaction remains active even in the absence of noise (which was the case in our simulations) while the depletion force vanishes in the same circumstance. In particular, this explains why we observe segregation and nematic ordering in the absence of noise. Furthermore, the aforementioned references focus on the nematic ordering of the fibers, the spherical particles being introduced only for the purpose of reducing the free volume. In the present paper, we have provided evidence that this interaction could also affect the shapes of the clusters of spherical particles (the cells) and have provided quantitative analysis of this effect. The interesting finding is that the fiber nematic order reflects itself into the elongated shape of the cell clusters, an effect which we have not seen noticed elsewhere.

By correlating the appearance of new adipocytes to a larger concentration of existing adipocytes, biased insemination equipped our model with an indirect way to investigate the influence of vasculature. Indeed, blood supplies the substrates required for adipogenesis and favors the appearance of new adipocytes at the extremities of capillaries where existing adipocytes are already present, leading to a correlation between new adipocyte locations and those of existing ones. However, fully random insemination alone reproduced the global tissue organization qualitatively well, and little improvement was observed by turning biased insemination on. This result suggests that no preferred location for differentiation of immature cells into adipocytes is required. It looks inconsistent with Wasserman’s analysis [47] and other works demonstrating that vasculature plays a key role in adipogenesis. Indeed, it has been demonstrated that adipose progenitors reside in the adipose vasculature and that pericytes which enwrap capillaries could represent a reservoir of adipocyte progenitors [9, 40]. Thus it is classically considered that the global adipose tissue architecture is primarily driven by vasculature [47, 7]. Our findings open the possibility that the role of vasculature in the emergence of global adipose tissue architecture could be more complex than previously thought.

Other phenomena were explored, such as local fiber alignment (to model ECM reorganization by stem cells) or the suppression of isolated cells (to model isolated cell apoptosis). These additional phenomena led to a broad range of tissue structures. Although not necessarily relevant for adipose tissue, these structures could account for other organs (such as muscles, liver, etc.) or pathological adipose tissues (such as fibrotic ones). More quantitative work is needed and these questions will be developed in future works.

Further improvements of the model could be made. If the choice of modelling the cells as 2D spheres was motivated by our observations on the biological images (see Fig. 1), more complex cell shapes could be taken into account without profound modifications of the model. Indeed, cells could be modelled as ellipses, or as sets of several overlapping disks connected to each other. These modelling choices should have to be accompanied by appropriate expressions of the interactions/ constraints in the minimization algorithm, but the general methodology would be preserved. Moreover, vasculature formation could be explicitly introduced in the model and its role as a niche for pericytes and pre-adipocytes could be explored. Incorporating immature cells could help investigating their role in the reconstruction of the ECM [17]. Similarly, coupling cell apoptosis with spontaneous ECM reconstruction [19] would improve the treatment of the latter. The adipocyte growth law could be made dependent on the local stress as in [5, 36]. Macroscopically, introducing a disruption of the equilibrium by suppressing a part of the tissue would open up new applications of the model, such as wound healing. From a mathematical viewpoint, the development of a macroscopic model from the present IBM would allow us to perform simulations on larger domains and address the question of the organization of the whole tissue. A first step has been made in [11] where a macroscopic model for interacting fibers has been derived. This macrocopic model has been further analysed and its derivation has been numerically validated in [34]. From a biological viewpoint, experimental data on ECM protein deficient mice strains (such as collagen type VI) [28] will be helpful to further validate our model. Designing new experiments allowing for the in-vivo registration and tracking of cell positions and ECM location during the development of adipose tissue would be of great value as giving unprecedented access to the time dynamics of the morphogenesis process.

The model could also be extended to 3D without profound methodological modifications. In Fig. 7, we show a three-dimensional image of an adipose lobule. This image reveals that while the cross-section of the lobule in a plane parallel to the skin is vaguely reminiscent of an ellipse, its shape in any plane perpendicular to the skin is much more elongated and convoluted. The results of our 2D model correctly reproduce the lobular organization in a plane parallel to the skin. A 3D extension of the model is needed to correctly account for its complex 3D organization.

To extend the model to 3D, several options are possible. Fibers could still be considered as line segments and their connection condition would allow for close but non-coplanar segments to connect with each other. Another possibility to account for the fact that fiber septa are more akin to two-dimensional surfaces would be to introduce planar ECM elements in the form of disks. Interconnecting disks would have the possibility to connect along their intersecting segments. On the other hand, cells could be modeled as spheres in 3D like in 2D, which would not require any change from the present 2D model. As mechanisms for network generation exhibit some differences in 2D and 3D, an increase of dimension in our model might affect the form of cell and fiber clusters. Testing different modelling choices for individual fibers or introducing new phenomena in the model could be necessary to reproduce the 3D structures observed in vivo, as for example to distinguish between sheet-like or tubular fiber structures according to real experiments. To this aim, more information will be needed on the biological viewpoint. When passing the model to 3D, we expect to recover that the shape of the 3D lobules is mainly controlled by the mechanical properties of the extra-cellular matrix, expectations which need to be confirmed by numerical experiments.

Data availability

Data supporting this work are available on figshare: $https://figshare.com/projects/Simple\_mechanical$

$\_cues\_could\_explain\_adipose\_tissue\_morphology/16419$

Acknowledgments

DP gratefully acknowledges the hospitality of Imperial College London, where part of this research was conducted. DP wishes to thank J. Fehrenbach (IMT, Toulouse) for enlighting discussions on image processing, and C. Guissard (StromaLab, Toulouse) for providing 3D imaging of adipose lobule.

Appendix A Immunohistochemistry and confocal microscopy.

Whole mouse inguinal adipose tissues (AT) were fixed, embeded in agarose and cut into 300$\mu$m slices. Slices were permeabilized in PBS/2% Normal Horse serum 0.2% triton 4h at, room temperature (RT), and then incubated in anti-perilipin antibody (1/250, P1873 Sigma, 24h RT). After washing, slices were incubated in alexa488-conjugated goat anti-rabbit IgG (1/250, A11008 Molecular Probes). Imaging was performed using a Confocal Laser Scanning microscope (LSM510 NLO - Carl Zeiss, Jena, Germany) with an objective lens LCI “Plan- Neofluar” 25x/0,8 and excited using a 488 nm argon laser. Images were obtained by the stitching of 45 acquired images with 10% overlapping with Fiji defined by image metadata Grid/collection stitching plugins [35]. Automatic image segmentation techniques were developed for (i) cell segmentation and (ii) lobule segmentation. A fully-automatic method for (i) based on a marker-controlled watershed technique [41] was implemented (see D for details).

Appendix B Mathematical model

Here, we give details about the mathematical model presented in section 2.2 of the main text. Let us recall the main features and introduce some notations. The two-dimensional Individual Based Model (IBM) consists of $N_{A}$ adipocytes, which are modeled as 2D growing spheres of center $X_{i}$ and radius $R_{i}$ for $i$ in $[1,N_{A}]$, and $N_{f}$ extra-cellular-matrix fiber elements which are represented by straight lines of fixed length $L_{f}$, of center $Y_{k}$ and orientational angle $\theta_{k}$ for $k\in[1,N_{f}]$. Adipocytes are prevented to overlap to model volume exclusion between cells. Fiber elements have the ability to link to or unlink from each other to model fiber elongation or rupture. The resistance of fibers to growing adipocytes is modeled by a repulsion potential $W_{pot}$ between cells and fibers. Additionally, fibers offer resistance to bending through an alignment potential $W_{al}$ acting between two linked fiber elements. Cells and fibers seek to minimize their mechanical interaction energy resulting from these two potentials, subject to the non-overlapping constraint between cells and to the linkeage constraint between linked fibers. In the course of the simulation, adipocytes grow, new adipocytes appear as a result of stem-cell differentiation, new fiber links can appear and existing links can disappear. These phenomena disrupt the mechanical equilibrium and force the cells and fibers to move in order to restore the equilibrium.

The outline of this appendix is the following: in B.1 we give the details of the mechanical interaction potentials, and B.2 shows the functionals defined for the cell and fiber constraints. In B.3, we give details on the modelling of biological phenomena and detail the time scales used for the simulations. The numerical algorithm for solving the constrained problem and numerical considerations are given in B.4-B.6.

B.1 Mechanical interaction potential

We recall that given a configuration at a fixed time, $C$ denotes the set of cell center positions and radii: $C=\{(X_{i},R_{i})\;,\;i\in[1,N_{A}]\}$ and $F$ the set of fiber center positions and fiber directional angles: $F=\{(Y_{k},\theta_{k})\;,\;k\in[1,N_{f}]\}$. The global mechanical energy of the system reads (1), which we remind here for the reader’s convenience:

[TABLE]

where $W_{pot}$ and $W_{al}$ read:

[TABLE]

where $p^{t}_{km}$ are time-dependent coefficients that are equal to 1 if fibers $k$ and $m$ are linked and are equal to 0 otherwise. Their time-evolution is given in B.3. The repulsion potential $W_{pot}$ is supposed to be the sum of two-particle potential elements, $W_{i,k}$, modeling the mechanical interaction between cell $i$ and fiber $k$. For two given vectors $X$ and $Y$ of $\mathbb{R}^{2}$ and an angle $\theta\in[-\pi,\pi]$, $W_{i,k}=W_{i,k}(X,Y,\theta)$ reads:

[TABLE]

where:

[TABLE]

and $\omega(\theta)=\begin{pmatrix}\cos\theta\\ \sin\theta\end{pmatrix}$ is the unit vector associated to the fiber angle $\theta$. The fiber-cell repulsion potential iso-lines are ellipses with focii located at the two ends of the fiber segment. The potential vanishes beyond distance $d_{0,i}$ to the fiber center (see Fig. 2). Parameter $d_{0,i}$ is set such that the length of the ellipse semi-minor axis is $\tau R_{i}$ (see Fig. 2), with $\tau$ a parameter set equal to $3$ in the simulations. In this case, a fiber repels cell $i$ up to a distance $\tau R_{i}$ in its orthogonal direction. A direct computation gives:

[TABLE]

Finally, the factor $\tilde{W}(\lambda^{+}_{k})$ in (5) measures the strength of each repulsion potential element. In order to model the fact that a fiber network is stiffer when the fibers are aligned [16], $\tilde{W}(\lambda^{+}_{k})$ is assumed to be a linear increasing function of the fiber local alignment $\lambda^{+}_{k}$. The local alignment of fibers around fiber $k$ is computed in a neighborhood $B(Y_{k},R_{al})$, where $R_{al}$ is the sensing distance up to which fiber $k$ senses the direction of its neighbors. Let $P_{k}$ denote the mean of the projection matrices on the direction vectors of the fibers in $B(Y_{k},R_{al})$:

[TABLE]

where $n_{k}$ denotes the number of fibers contained in $B(Y_{k},R_{al})$, and $\omega_{m}$ is the directional vector of fiber $m$. The maximal eigenvalue $\lambda^{+}_{k}$ of $P_{k}$ measures the mean alignment of the fibers in $B(Y_{k},R_{al})$. Its corresponding normalized eigenvector gives the mean direction of the fibers in $B(Y_{k},R_{al})$. A direct computation leads to:

[TABLE]

where

[TABLE]

Note that $\lambda^{+}_{k}=1$ when all the fibers in $B(Y_{k},R_{al})$ have the same direction and $\lambda^{+}_{k}=0$ when the fiber directions are fully random. The intensity of the potential element is then set to:

[TABLE]

where $W_{0}$ and $W_{1}$ are the intensities of the repulsion potential between fiber $k$ and cell $i$, when the local alignment around fiber $k$ is weak or strong respectively.

B.2 Constraints

In order to model adipocyte incompressibility and non overlapping, we assume that the radius of each disk is unaffected whatever mechanical efforts are exerted onto it. The non overlapping constraint between cells $i$ and $j$ is written as an inequality constraint on the following function $\Phi_{ij}$:

[TABLE]

One immediately notes that cells $i$ and $j$ do not overlap if and only if $\Phi_{ij}(X_{i},X_{j})\leq 0$.

To model fiber growth and elongation or conversely rupture, unlinked (resp. linked) intersecting fibers have the possibility to link (resp. unlink) at random times. As long as a pair of linked fibers remains linked, the attachment sites of the two linked fibres are kept at the same point. The maintain of the link between fibers $k$ and $m$ is modeled as equality constraints $\vec{\Psi}_{km}=0$ with:

[TABLE]

where $\ell_{km}$ is the distance of the center of fiber $k$ to its attachment site with fiber $m$ (see Fig. 8) at the moment when the link is created. We use, if $\sin(\theta_{m}-\theta_{k})\neq 0$:

[TABLE]

where $Y^{0}_{k}=(x^{0}_{k},y^{0}_{k})$ are the 2D coordinates of the center of fiber $k$ when the link is created (and similarly for fiber $m$, see Fig. 8).

B.3 Modeling of the main biological phenomena

Pre-adipocyte differentiation: New adipocytes of minimal radius $R_{e}$ are inseminated at random times following a Poisson process of frequency $\nu_{ins}$. The location of the insemination is random with either uniform probability in the domain or with a bias resulting in a higher insemination probability at locations where existing adipocytes are already present. In this last case, the probability of inseminating at a random point $X$, $\mathcal{P}(X,R)$ is a function of the cell density computed in the ball of center $X$ and of radius $R$ and normalized by the maximal possible density in this ball. This models a quorum-sensing process around point $X$ where a pre-adipocyte senses the adipocyte density $\chi(X,R)$, up to a sensing distance $R$. This normalized density reads:

[TABLE]

where $N_{R}$ is the maximal number of cells of radius $R_{\max}$ contained in $B(X,R)$ (for instance $N_{R}=7$ for $R=2R_{\max}$). Note that $\chi\in[0,1]$. Then, the probability of inseminating at $X$, $\mathcal{P}(X,R)$, reads:

[TABLE]

where $\alpha>0$ is the biasing parameter. By denoting $t_{ins}$ the time needed to inseminate 1 cell, we define the characteristic time of the insemination process, $t_{e}$, as the mean time needed to inseminate $N_{A}$ cells:

[TABLE]

Adipocyte growth The volumes of the cells are supposed to grow linearly with time. Given a cell $i$ at time $t$, the radius of cell $i$ at time $t+\Delta t$ reads:

[TABLE]

where $\eta$ is a random number chosen uniformly in $[0,1]$ and $K_{g}$, $\rho_{g}$ are two parameters such that $\frac{K_{g}}{\Delta t}$ is the mean volumic cell growth per unit of time and $\frac{K_{g}\rho_{g}}{\Delta t}$ is related to the standard deviation of the volumic cell growth per unit of time. The characteristic time of cell growth $t_{g}$ is defined as the mean time needed for a cell to reach its maximal radius $R_{\max}$ and reads:

[TABLE]

Fiber growth or rupture: Fiber elongation is modeled by giving fibers the ability to attach to each other. Pairs of unlinked intersecting fibers link together at random times following a Poisson process of frequency $\nu_{\ell}$. To model fiber rupture, two linked fibers unlink following a Poisson process of frequency $\nu_{d}$. Let $(k,m)\in[1,N_{f}]^{2}$ and define $p^{t}_{km}$ such that $p^{t}_{km}=1$ if fibers $k$ and $m$ are linked, 0 otherwise. The time evolution of $p^{t}_{km}$ is given by:

[TABLE]

where $\ell_{km}$ and $\ell_{mk}$ are given by Eq. (11). $\mathbb{P}(p^{t+\Delta t}_{km}=1|p^{t}_{km}=0)$ describes the transition probability for a transition of $p^{t}_{km}$ from [math] to $1$ during the time interval $[t,t+\Delta t]$ while $\mathbb{P}(p^{t+\Delta t}_{km}=0|p^{t}_{km}=1)$ refers to the transition probability for the reverse process. We define $t_{\ell}$ and $t_{d}$ as the characteristic times of the linking and unlinking of fibers: respectively:

[TABLE]

In order to analyse the fiber linking/unlinking process, we define the ratio $\chi_{\ell}$:

[TABLE]

Note that $\chi_{\ell}$ is directly correlated to the fraction of linked fibers (among pairs of intersecting fibers) at equilibrium if the linking-unlinking process was acting alone and will be referred to as the ’linked fiber fraction’ for short. For each unlinking frequency $\nu_{d}$ and each linked fibers fraction $\chi_{\ell}$, the linking frequency is set to $\nu_{f}=\frac{\chi_{\ell}}{1-\chi_{\ell}}\nu_{d}$.

Time scales The reference time is chosen to be ten times the insemination time $t_{ins}$ (inverse of insemination rate) i.e $t_{ref}=10t_{ins}$, and the dimensionless time step is chosen to be $\Delta t=0.1t_{ref}=t_{ens}$. As explained in the main text, new cells are inseminated until reaching a cell volume fraction of 50%. This corresponds to a total number of adipocytes $N_{A}=180$. Therefore the dimensionless time after which insemination stops is $t_{e}=18t_{ref}$. The dimensionless mean time for a cell to reach its maximal radius is $t_{g}=\frac{R_{\max}^{3}\Delta t}{K_{g}}\approx 18t_{ref}$. Finally, the fiber dimensionless mean unlinking time $t_{d}=\frac{1}{\nu_{d}}$ and the ratio $\chi_{\ell}$ which measures the proportion of pairs of linked fibers determine the fiber linking frequency $\nu_{\ell}$ as follows: $\nu_{\ell}=\frac{\chi_{\ell}}{1-\chi_{\ell}}\nu_{d}$ and the dimensionless time for cross-fiber linking $t_{\ell}=\frac{1}{\nu_{\ell}}=\frac{1-\chi_{\ell}}{\chi_{\ell}}t_{d}$. In Fig. 10, we show an example of the different time scales used for the simulations of the main text, for $\chi_{\ell}=0.35$ and $\nu_{d}=10^{-1}$.

B.4 Uzawa algorithm

The constrained minimization problem (4) is solved with an Arrow-Hurwicz-Uzawa type algorithm [42]. To perform the minimization, we introduce a Lagrangian:

[TABLE]

where \lambda=\big{(}\lambda_{ij}\big{)}_{(i,j)\in[1,N_{A}]^{2}}, $\lambda_{ij}>0$ and \vec{\mu}=\big{(}\vec{\mu}_{km}\big{)}_{(k,m)\in\mathcal{N}_{f}}, $\vec{\mu}_{km}\in\mathbb{R}^{2}$ are the sets of Lagrange multipliers of the constraints.

Given a configuration $(C(t_{n}),F(t_{n}))$ at time $t_{n}=n\Delta t$ the configuration $(C(t_{n+1}),F(t_{n+1}))$ at time $t_{n+1}=(n+1)\Delta t$ is defined as the limit as $p\rightarrow\infty$ of the iterative sequence $(C^{p},F^{p})$ where $(C^{p},F^{p},\lambda^{p},\vec{\mu}^{p})$ is defined for all $(i,j)\in[1,N_{a}]^{2}$ and all $(k,m)\in\mathcal{N}_{f}$ by:

[TABLE]

with initial condition $(C^{0},F^{0},\lambda^{0},\vec{\mu}^{0})=(C(t_{n}),F(t_{n}),\lambda^{0},\vec{\mu}^{0})$ and $\lambda_{ij}^{0}=0$, $\vec{\mu}_{km}^{0}=0$, for all $(i,j)\in[1,N_{A}]^{2}$ and all $(k,m)\in\mathcal{N}_{f}$. The parameters $\lambda_{1}$ and $\mu_{2}$ control the actualization of the constraints and their choice is detailed in the next section. The minimization steps $\alpha^{i}_{a}$, $\alpha^{k}_{f}$ and $\alpha^{k}_{\theta}$ control the elementary motion of cell $i$ and fiber $k$ and the elementary rotation of fiber $k$ respectively. Their computation is detailed in the next section. The convergence test of the algorithm reads:

[TABLE]

for a chosen $\epsilon_{r}>0$. Here, $\mathcal{L}^{p}$ is the value of the Lagrangian at iteration $p$. Because of the non convexity of the minimization problem, the uniqueness of the solution to (4) is not ensured and a configuration at each time step corresponds to a local minimum of the minimization problem.

B.5 Choice of the numerical parameters

In this section, the numerical parameters $\alpha^{i}_{a},\alpha^{k}_{f}$ and $\alpha^{k}_{\theta}$ of (B.4) are chosen such that the amplitude of the change of each variable remains controlled. Given three bounds $\delta_{a}$, $\delta_{f}$ and $\delta_{\theta}$, the goal is to ensure $|X_{i}^{p+1}-X_{i}^{p}|\leq\delta_{a}$, $|Y_{k}^{p+1}-Y_{k}^{p}|\leq\delta_{f}$ and $|\theta_{k}^{p+1}-\theta_{k}^{p}|\leq\delta_{\theta}$, for each cell $i$ and fiber $k$. Using (12) and (B.4), the following expressions hold:

[TABLE]

The parameters $\alpha^{i}_{a}$, $\alpha^{k}_{f}$ and $\alpha^{k}_{\theta}$ are consequently set such that:

[TABLE]

The gradients of the potential $\mathcal{W}$ have the following upper bounds for all $i\in[1,N_{a}]$ and $k\in[1,N_{f}]$ (see (1)-(6)):

[TABLE]

where $n^{f}_{i}$, $n^{a}_{k}$ and $\ell^{f}_{k}$ denote the number of fibers interacting with cell $i$, the number of cells interacting with fiber $k$ and the number of fibers linked to fiber $k$ respectively. Here, $d_{0}$ is the value of $d_{0,i}$ evaluated with $R_{i}=R_{\max}$. The following upper bounds for the gradients of the constraint functions (see Eqs. (9)-(11)) are chosen:

[TABLE]

These three gradient bounds are estimated at each iteration of the minimization algorithm and are included into Eqs. (15) to compute the values of the numerical steps. We now turn towards the determination of $\lambda_{1}$ and $\mu_{2}$ of Eqs. (B.4). Dimensionnally, the following estimations can be set from the expression of the Lagrangian Eq. (12) (for all pairs $(i,j)$ and $(k,m)$):

[TABLE]

From the actualization of the constraints given by iterations of Eqs. (B.4), we set:

[TABLE]

Then, parameters $\lambda_{1}$ and $\mu_{2}$ are set to:

[TABLE]

The values of the parameters $\delta_{a}$, $\delta_{f}$ and $\delta_{\theta}$ are taken of the order of $10^{-3}$ and the convergence test tolerance to $\epsilon_{r}=10^{-5}$ (the complete set of the numerical parameters can be found in Table 11).

B.6 Decreasing the computational time

The simulations are performed on a 2D-domain $\Omega=[x_{\min},x_{\max}]\times[y_{\min},y_{\max}]$. In order to reduce the computational time, the domain is divided into sub-squares whose side length $L_{s}$ is a measure of the maximal distance of the agent interactions. The goal is to compute each interaction potential element with the agents located in neighboring sub-squares of the domain only. The procedure is classical and details are omitted. Periodic boundary conditions are set by creating ghost numerical boxes of length $L_{s}$ at each boundary of the domain.

Appendix C Statistical quantifiers

This section is devoted to the computation of statistical quantifiers used to describe cell and fiber structures in both numerical simulations and biological images. A cell cluster is defined as a set of cells almost in contact. Let $\sim_{a}$ be the reflexive and symmetric relation:

[TABLE]

where $\mathcal{N}_{i}$ is the set of cell $i$ neighbors:

[TABLE]

where $\epsilon_{a}$ is the maximal allowed distance up to which two cells not in contact are defined as neighbors and is set to $50\%\max(R_{i},R_{j})$. The equivalence relation $\sim_{A}$ then reads:

[TABLE]

Cells $i$ and $j$ belong to the same cluster if and only if $i\sim_{A}j$. Fig. 12 shows an example of cell cluster separation.

The statistical quantifier $N_{C}$ is defined as the total number of cell clusters which have more than 5 adipocytes per 100 adipocytes.

The statistical quantifier $E$ measures the mean elongation of the cell clusters, and is defined as the number of cells at the boundary of the clusters normalized by the total number of cells in the clusters. As the parameter $E$ is irrelevant for clusters with less than 5 cells, its computation is restricted for clusters $c$ such that $n_{c}>5$ and reads:

[TABLE]

Here, $\mathcal{C}_{c}$ is the set of indices of the cells belonging to cluster $c$, $n_{c}$ is the number of cells in cluster $c$ and $\mathcal{R}$ is the set of indices of all cells with less than 5 neighbors:

[TABLE]

where $\mathcal{N}_{i}$ is defined by Eq. (16).

Finally, in order to determine if the cell clusters have anisotropic shape with a preferred direction, we define the SQ $\Theta^{c}$ as the angle of cluster $c$ shape anisotropy direction. For this purpose, let $X^{c}$ be the center-of-mass of cluster $c$, i.e:

[TABLE]

Then, we define $P^{c}$ as the mean of the projection matrices on the vectors $X_{i}-X^{c}$, for all $i$ in cluster $c$:

[TABLE]

The maximal eigenvalue $\lambda^{+}_{c}$ of $P^{c}$ gives a measure of the shape anisotropy of cell cluster $c$ and its associated eigenvector $u^{c}=(u_{1}^{c},u_{2}^{c})$ gives the shape anisotropy direction. Then, $\Theta^{c}$ is defined as:

[TABLE]

Note that $\Theta^{c}\in[-\frac{\pi}{2},\frac{\pi}{2}]$. The SQ $\Theta$ is then defined as the circular standard deviation of all the angles $\Theta^{c}$ for all clusters:

[TABLE]

where $\bar{R}$ reads:

[TABLE]

Finally, the mean $\bar{\Theta}$ of $\Theta^{c}$ over all the cell clusters reads:

[TABLE]

which ensures that $\bar{\Theta}\in[-\frac{\pi}{2},\frac{\pi}{2}]$. Note that large $\Theta$ corresponds to fully isotropic cell cluster organization, while small $\Theta$ indicates that cell clusters have a preferential direction.

In order to describe the fiber structures, we define a fiber cluster as a set of neighboring quasi-aligned fiber elements and $\Lambda$ as an estimate of the average curvilinear length of such fiber clusters. Finally, $A$ measures the mean alignment of the fibers of a cluster. For this purpose, let us define $\mathcal{M}_{k}$ as the set of neighbors of fiber $k$, quasi-aligned with fiber $k$. Then:

[TABLE]

where $d_{f}(Y_{k},Y_{m})$ reads:

[TABLE]

Here, $Y^{\pm}_{m}=Y_{m}\pm\frac{L_{f}}{2}\omega_{m}$ and $d(X,Y)$ is the distance of point $X$ to point $Y$. Note that $d_{f}(Y_{k},Y_{m})\leq 0$ (resp. $d_{f}(Y_{m},Y_{k})\leq 0$) if the center of fiber $k$ (resp. $m$) is contained in the ellipse of focii $Y^{\pm}_{m}$ (resp. $Y^{\pm}_{k}$) with semi minor axis of length $\tau_{f}L_{f}$ and semi major axis of length $\sqrt{(\frac{L_{f}}{2})^{2}+(\tau_{f}L_{f})^{2}}$. We chose $\tau_{f}=\frac{1}{3}$, which means that a fiber detects a neighboring fiber up to a distance $\frac{L_{f}}{3}$ in its orthogonal direction. This allows us to define fiber clusters as sets of quasi-aligned neighboring fibers. Let us define the reflexive and symmetric relation $\sim_{f}$ by:

[TABLE]

Define the equivalence relation $\sim_{F}$ such that:

[TABLE]

Then, we say that fibers $k$ and $m$ belong to the same cluster if and only if $m\sim_{F}k$. Fig. 13 shows the results of fiber cluster detection from the numerical image displayed on Fig. 12.

The mean elongation of fiber clusters is estimated by $\Lambda$. Given a fiber cluster $c_{f}$ and a division of the simulation domain into squares of side length $L_{f}$, the length of $c_{f}$ is estimated by $L_{f}\Lambda^{c_{f}}_{F}$, where $\Lambda^{c_{f}}_{F}$ is the number of squares which contain the centers of at least one fiber of $c_{f}$. Then, the dimensionless mean fiber cluster elongation $\Lambda$ is defined as the mean of $L_{s}\Lambda^{c_{f}}_{F}$ over all the fiber clusters, normalized by the maximal cell diameter:

[TABLE]

where $N_{T_{f}}$ is the total number of fiber clusters. The longer the fiber clusters, the larger the $\Lambda$.

Finally, we define the SQ $A$ to quantify the mean alignment of the fibers of a cluster. Given a fiber cluster $c_{f}$, the mean alignment of its fibers is defined as the maximal eigenvalue $\lambda^{+}_{c_{f}}$ of the mean projection matrix defined by Eq. (8), where the set $B(V_{k},R_{al})$ is replaced by the set of all fibers of cluster $c_{f}$. Then, $A$ is defined as the mean of the fiber cluster alignment, weighted by the number of fibers in the cluster:

[TABLE]

where $n_{c_{f}}$ is the number of fibers in cluster $c_{f}$.

Appendix D Image processing

This section is devoted to the algorithms and results of the image processing. The goal is to develop segmentation techniques for (a) cell detection and (b) cell cluster detection, in order to compute the SQs on biological images and compare them to those of numerical simulations. It is noteworthy that adipocytes only are visualized in biological images at hand, therefore SQs $E$, $N_{C}$ and $\Theta$ only are accessible from biological images.

(a) Detection/separation of cells First, a fully-automatic method for cell detection based on marker-controlled watershed segmentation has been realized.We use Marker-controlled watershed segmentation, according to the following procedure:

• (i)

The biological image is first filtered by a local median filter which associates to each pixel its mean value in its local neighborhood

• (ii)

The segmentation function is defined as the gradient of the transformed image. The gradient is high at the borders of the objects and low inside the objects.

• (iii)

Compute foreground markers. These are connected blobs of pixels within each of the objects. The morphological techniques ’opening-by-reconstruction’ and ’closing-by-reconstruction’ are used to clean up the image. These operations create flat maxima inside each object that can be located using the intrinsic Matlab function imregionalmax.

• (iv)

Compute background markers. These are pixels that are not part of any object. We perform a simple thresholding of the intensity image: each pixel whose intensity is lower than the mean intensity of the image is set to 0.

• (v)

Modify the segmentation function so that it only has minima at the foreground and background marker locations.

• (vi)

Compute the watershed transform of the modified segmentation function.

Object boundaries are located where $W=0$, where $W$ is the watershed transform of the marked image gradient. This method enables the separation of multiple objects. Each object is characterized by a center (center of mass of the detected region) and a radius $R$ (radius of a circle which has the same area $a$ as the object): $R$ is thus computed as $\sqrt{a/\pi}$.

(b)Detection/separation of lobule-like clusters. For cell cluster detection, a semi-automatic method has been developed. Each sub-images (squares occupying $0.1\%$ of the image area) is filtered by median filtering with the intrinsic function medfilt2 of Matlab. Each output pixel contains the median value in the 3-by-3 neighborhood around the corresponding pixel in the input image. A thresholding of the intensity image fixed at $40\%$ of the mean intensity of the subimage is then applied. The connected objects with 8-connectivity are computed using the Matlab intrinsic function bwlabel. Finally, objects containing less than 2000 pixels (noise objects) are suppressed with the intrinsic Matlab function bwareaopen. If neighboring cell clusters are still visually connected at a point, a line is plotted by hand to separate the two clusters. The process of cell cluster detection is semi-automatic in this sense, but this procedure is sufficient for the purpose of this work, given the low number of biological images to be treated.

Appendix E Results and their analysis

E.1 Simulations of the main text

Here, supplementary results on the simulations of the main text are given. In Fig. 14 (A and C), we show a complete diagram for the values of the Statistical Quantifiers $E$ and $A$ as functions of the fiber unlinking frequency $\nu_{d}$ ($x$-axis) and of the linked fiber fraction $\chi_{\ell}\in[0.1,0.35]$ ($y$-axis), with random insemination and flexural modulus $c_{1}=1$ (see Fig. 4 of section 3 of the main text). In Fig. 14 (B and D), we also show the values of $\Lambda$ and $N_{C}$ corresponding to the simulations of the main text. The values of the SQ are averaged over 10 simulations and we refer the reader to section 3 of the main text for the analysis of the quantifiers $E$ and $A$ and focus here on the SQ $N_{C}$ and $\Lambda$.

Figs. 14 (B and D) first reveal that the fiber mean elongation $\Lambda$ and the number of cell clusters $N_{C}$ have plateau values for $\nu_{d}\leq 10^{-3}$. For a slow fiber linking-unlinking process ($\nu_{d}\leq 10^{-3}$), the fibers in the clusters are poorly aligned (low value of $A$, see section 3 of the main text) and the mean fiber cluster elongation $\Lambda$ is fairly large, meaning that the fibers keep their initial entanglement. For large $\chi_{\ell}=0.35$ and as $\nu_{d}$ increases, the mean fiber cluster elongation $\Lambda$ increases until reaching a maximal value for $\nu_{d}\approx 0.005$. Then, $\Lambda$ loses $50\%$ of its value when $\nu_{d}$ increases in the range $[0.005,0.1]$. Indeed, as explained in Main text, if the linking-unlinking process is too fast, the fiber structures easily align but they are also more sensitive to the compression by the cells and consequently, the fiber elements regroup into shorter fiber clusters, thereby decreasing the fiber length. For a value of the linked fiber fraction $\chi_{\ell}=0.1$, no significant change in the mean cluster number $N_{C}$ arise as the unlinking frequency $\nu_{d}$ increases (Fig. 14 (B)). By contrast, for a value of the linked fiber fraction $\chi_{\ell}=0.35$, the mean cluster number $N_{C}$ increases with increasing unlinking frequency $\nu_{d}$ in the range $[10^{-3},0.05]$ and then stays constant for larger values of $\nu_{d}$ (Fig. 14 (B)). These results show the ability of a connected fiber network to encompass well separated cell clusters for a well chosen fiber linking/unlinking dynamics.

E.2 Influence of the flexural modulus $c_{1}$

Here, we perform a statistical analysis of the influence of the fiber flexural modulus $c_{1}$. Fig. 15 (I) shows simulations with random insemination and two different flexural moduli $c_{1}=0.01$ (first row) and $c_{1}=10$ (second row), for $\chi_{\ell}=0.35$ and three different unlinking frequencies $\nu_{d}=10^{-3}$ (Fig. 15 (I A) and (I D)), $\nu_{d}=10^{-2}$ (Fig. 15 (I B) and (I E)) and $\nu_{d}=0.2$ (Fig. 15 (I C) and (I F)). Fig. 15 (II) shows the values of the mean elongation of cell clusters $E$ (Fig. 15 (II A)), cell cluster number $N_{C}$ (Fig. 15 (II B)), mean alignment of fiber clusters $A$ (Fig. 15 (II C)) and mean elongation of fiber clusters $\Lambda$ (Fig. 15 (II D)), averaged over 10 simulations and plotted as functions of the unlinking frequency $\nu_{d}$ for two different values $c_{1}=0.01$ (black curve) and $c_{1}=10$ (red curve) of the flexural modulus.

The first row of Figs. 15 (I) (Fig. 15 (I A), (I B) and (I C)) shows that for a small flexural modulus $c_{1}=0.01$, the cell structures change from well-separated lobule-like cell clusters (Fig. 15 (I A)) to slightly more elongated clusters (Fig. 15 (I C)) and the fiber structures change from a disorganized fiber network (Fig. 15 (I A) or (I B)) to more aligned fiber clusters (Fig. 15 (I C)), as $\nu_{d}$ increases. For a large flexural modulus $c_{1}=10$ (second row of Fig. 15 I, i.e. Fig. 15 (I D), (I E) and (I F)), the fiber structures are more aligned: the model generates organized fiber clusters able to bend around the cell structures (Fig. 15 (I D) and (I E)). For a fast linking-unlinking (Fig. 15 (I F)), a rigid fiber network composed of long fiber threads which fail to surround the cell clusters is observed.

These observations are confirmed by the values of the Statistical Quantifiers $E$ and $N_{C}$ (for cell clusters) and $A$ and $\Lambda$ (for fiber clusters) shown in Figs. 15 (II) as functions of the unlinking frequency $\nu_{d}$, for $c_{1}=0.01$ (black curves) and $c_{1}=10$ (red curves). Fig. 15 (II A) shows that the flexural modulus $c_{1}$ does not seem to significantly change the mean cell cluster elongation. Fig. 15 (II B) reveals that the number of cell clusters $N_{C}$ is significantly lower for $c_{1}=10$ than for $c_{1}=0.01$ for $\nu_{d}>10^{-3}$. For $c_{1}=10$ and $\nu_{d}\in[10^{-3},0.1]$, this is because ECM rigidity is too large and the fibers fail to separate cell structures, compared to the case $c_{1}=1$ (see section 3 of the main text). For $\nu_{d}>0.1$, we recover the previously described case of a fast fiber linking-unlinking dynamics. As fibers fastly self-organize into long and directed rigid threads, they force the cells to group into chord-like unseparated structures.

Figs. 15 (II C) shows that the mean alignment of the fiber clusters increases with $c_{1}$ (compare the black and red curves), and the difference between the values of $A$ for $c_{1}=0.01$ and for $c_{1}=10$ increases with the fiber unlinking frequency $\nu_{d}$. This is because the fibers are more rigidly maintained with a slow linking-unlinking dynamics than with a fast one (see previous section), and are thus less sensitive to alignment. Fig. 15 (II D) shows that the mean fiber cluster elongation $\Lambda$ is a non monotonous function of $\nu_{d}$ for $c_{1}=0.01$ (black curve) and a monotonically decreasing function for $c_{1}=10$ (red curve). For $c_{1}=0.01$, $\Lambda$ first increases to reach a maximal value at $\nu_{d}\approx 0.05$ and then decreases. This reflects the ability that the fibers have to surround the cell clusters when the flexural modulus is small and the unlinking frequency $\nu_{d}$ is moderate, as already seen in the previous section. This ability is lost with a larger flexural modulus and $\Lambda$ becomes just a decreasing function of $\nu_{d}$.

To sum up, large flexural modulus favors fiber alignment and ECM rigidity compared to small $c_{1}$. Moreover, the choice of the flexural modulus has to be carrefully linked to the fiber linking-unlinking dynamics, which also triggers fiber network alignment and rigidity. For a well calibrated fiber linking-unlinking process and increasing values of $c_{1}$, the structures change from (a) compact middle sized cell clusters in a disorganized fiber network ($c_{1}=0.01$), (b) compact middle sized cell clusters in an organized fiber network ($c_{1}=1$) and (c) elongated cell clusters in fewer quantities inside an organized network ($c_{1}=10$). For $c_{1}<1$, ECM alignment is small and the fiber network cannot easily organize. By contrast, when $c_{1}=1$, ECM alignment is moderate and the fibers that are not too constrained can align. However for $c_{1}>1$, ECM rigidity is too large and this results in elongated cell clusters. Experimentally, it is observed that the lobules are more elongated at the periphery of the tissue than inside. Thus, our results suggest that fibers could be more stretched at the periphery. To support this hypothesis, it would be interesting to develop an experimental quantification method to estimate a local stress tensor similar to what has previously been done for adipocyte stiffness [43]. This analysis suggests that the different morphologies observed in adipose tissues of healthy mice according to the location of fat (central or peripheral)can be due to ECM stiffness. Varying its value is sufficient to break the architecture of the tissue.

E.3 Influence of biased insemination

Figs. 16 (I) shows simulation results obtained with random insemination (Figs. 16 (I A), (I B)) and with biased insemination for $\alpha=10^{-1}$ (Figs. 16 (I C), (I D)), for two values of the unlinking frequency $\nu_{d}$: $\nu_{d}=10^{-4}$ (Figs. 16 (I A), (I C)) and $\nu_{d}=0.1$ (Figs. 16 (I B), (I D)). We refer the reader to Fig. 3 of Main text for simulations with random and biased insemination using intermediate values of $\nu_{d}$ (namely $\nu_{d}=10^{-3}$ and $\nu_{d}=210^{-2}$). Fig. 16 (II) shows the values of the mean elongation of cell clusters $E$ (Fig. 16 (II A)), the number of cell clusters $N_{C}$ (Fig. 16 (II B)), the mean alignment of fiber clusters $A$ (Fig. 16 (II C)) and the mean elongation of fiber clusters $\Lambda$ (Fig. 16 (II D)), averaged over 10 simulations, as functions of the unlinking frequency $\nu_{d}$, for random insemination (in black) and for biased insemination with $\alpha=10^{-1}$ (in red). In these simulations, the flexural modulus between linked fibers is $c_{1}=1$ and the linked fiber fraction $\chi_{\ell}=0.35$. The numerical and model parameters can be found in Table 9.

Figs. 16 (I C) and (I D) show that biased insemination leads to the creation of bigger cell clusters than random insemination (compare with Figs. 16 (I A) and (I B)). This is because biased insemination favors insemination of new adipocytes at locations where cell clusters pre-exist. The regroupment of cells into clusters leaves regions devoid of cells (see Figs. 16 (I C)). For cell clusters, Fig. 16 (II A) shows that biased insemination seems to reduce cell cluster elongation. This can be explained by the fact that cell clusters are larger with biased insemination than with random insemination, which results in a decrease of $E$. In this case, the statistical quantifier $E$ does not allow us to conclude on the form of the cell clusters. Fig. 16 (II B) shows that biased insemination do not significantly change the number of cell clusters. Finally, Figs. 16 (II C) and (II D) show that biased insemination does not have a significant influence on the fiber cluster alignment $A$, but seems to slightly reduce the fiber elongation.

Altogether, this analysis demonstrates that biased insemination with $\alpha=0.1$ does not have a significant impact on the cell and fiber structures for a properly chosen fiber linking/unlinking dynamics. This suggests that, in a sufficiently rigid fiber network, cell clusterization is mainly driven by cell-fiber interactions. Cells and fibers self-organize into middle-sized well-separated cell clusters and aligned fiber structures whatever the type of insemination (random or biased with small $\alpha$) is.

E.4 Anisotropic initial fiber network

As discussed in section 3 of the main text, parts of the adipose tissue reveal an anisotropic cell and fiber organization. In order to obtain oriented cell clusters, we studied the properties of the model starting from an initially anisotropic fiber network. For this purpose, we let the initial fiber directional angles $\theta^{0}_{f}$ be randomly chosen according to a uniform distribution in the interval $[\theta_{1}-\theta_{2},\theta_{1}+\theta_{2}]$, where $\theta_{2}$ is related to the standard deviation of the distribution. Note that the smaller $\theta_{2}$, the more aligned the fibers initially are. By contrast, the simulations shown so far correspond to a fully isotropic initial network, i.e. to the case $\theta_{2}=\pi/2$. The initial number of fiber links was carefully adjusted to be independent of the initial value of $\theta_{2}$ throughout the forthcoming simulations. Indeed, the probability that pairs of fibers intersect is much smaller in an aligned network than in a fully isotropic one. Simulations of Fig. 17 have been obtained with random insemination, flexural modulus $c_{1}=1$ and linked fiber fraction $\chi_{\ell}=0.35$, for different unlinking frequencies $\nu_{d}$ and different values of $\theta_{2}$. Three types of structures have obtained according to the values of $\nu_{d}$ and $\theta_{2}$: (a) lobule-like non oriented cell clusters, (b) lobule-like oriented cell clusters, and (c) elongated and oriented cell clusters. In order to quantify the passage from one morphology to another one, we use the following SQ: the mean cell cluster elongation $E$ and the standard deviation of cell cluster shape anisotropy $\Theta$. We identify the threshold values $E^{*}=0.65$ for $E$ (the same value as in section 3 of the main text) and $\Theta^{*}=0.7$ for $\Theta$. Structures of type (a) correspond to $\Theta>\Theta^{*}$ and $E<E^{*}$. Type (b) is described by $\Theta<\Theta^{*}$ and $E<E^{*}$ and finally type (c) by $\Theta<\Theta^{*}$ and $E>E^{*}$. Fig. 17 shows a phase diagram in the $(E,\Theta)$ plane. Each point in this phase diagram corresponds to statistical quantifiers $(E,\Theta)$ averaged over 10 simulations. The red and blue lines correspond to the separatrix between phases (a) and (b) (of equation $\Theta=\Theta^{*}$) and between phases (b) and (c) (of equation $E=E^{*}$) respectively. Fig. 17 also shows a typical simulation result for each phase, and its position on the phase diagram according to the values of $E$ and $\Theta$.

Fig. 17 shows that, when the initial fiber network is anisotropic, the emergence of a shape anisotropy of the cell clusters depends on the fiber linking-unlinking dynamics. For a slow linking-unlinking dynamics ($\nu_{d}=10^{-3}$) the initial orientation of the fibers must be strongly biased to obtain directionality in the cell and fiber final structures. Otherwise (for $\theta_{2}>\frac{\pi}{5}$), the initial orientation of the network is lost, and cell structures without preferential direction are obtained (see Fig. 17 A). This suggests that for this slow linking-unlinking frequency, cells disturb the initial organization of the fiber network so much that this initial organization is lost. For a fast fiber linking-unlinking dynamics ($\nu_{d}=0.2$ see Fig. 17 C), cell structures are elongated due to the rigidity of the fiber network induced by the fast linking frequency and the action of the alignment torque at the created links. In this case, the fiber network imposes its preferred direction to cell cluster growth and we recover elongated cell clusters as in the case of an initially isotropic fiber network (see section 3 of the main text). Fig. 17 (B) shows that there exist a range of values of $\nu_{d}$ and $\theta_{2}$ for which the model is able to generate lobule-like cell clusters having anisotropic shapes and a preferred shape anisotropy direction. These configurations are obtained for $\nu_{d}\in]10^{-3},10^{-2}[$ and for $\theta_{b}\in[\frac{\pi}{8},\frac{\pi}{4}]$.

Appendix F Shorter fibers, smaller cell-fiber interaction range, fiber-fiber repulsion potential

In this section, we consider the combined effects of shortening the fiber element length, decreasing the cell-fiber interaction range and introducing a fiber-fiber repulsion potential. As described at the end of Section 2.2 of the main text, this leads to a more biologically relevant range of parameters and phenomena. On the other hand, computing with short fiber elements necessitates a large number of such elements and is extremely time consuming. Due to these computational constraints, the simulation of Main Text were run with longer fibers, larger cell-fiber interaction range and no fiber-fiber repulsion potential. The goal of this section is to support the relevance of this approach by showing that the results obtained in the present section with a more biologically relevant set of parameters and phenomena are similar to those described in the Main Text. That a large cell-fiber interaction range may have a similar effect as a smaller range combined with fiber-fiber repulsion makes reasonable sense. Fibers group in clusters. If they undergo a repulsion potential, these clusters are wider and affect the cells in a similar way as if the clusters were thinner (which is the case without fiber-fiber repulsion), but the cell-fiber interaction distance is larger, leading to about the same separation distance between cell clusters.

When using shorter fiber elements, we must be careful that initialization matters if we wish to compare them with longer fiber elements. Also, when the fiber element is shortened, its width encoded in the distance $d_{0}$ (for cell-fiber interaction) must be decreased in the same proportion. Therefore, dividing the fiber length by two indeed means replacing a fiber element by four elements connected altogether to form a diamond-shaped structure like in Fig. 18. So, when using short fibers, we will randomly initialize diamonds like in Fig. 18. However, as soon as the simulation is started the fiber element links in the diamonds may be removed as a consequence of the unlinking process. So, the diamond shape structures may be dissolved or persist according to how fast the unlinking process is. Considering an initial fiber network consisting of fibers longer than cell diameters is biologically plausible as collagen fibers may be produced before any implantation of adipocytes. Indeed, in wound healing, it is well documented that the wound is first populated by fibers before the appearance of the first new cell. Fiber elements in our model are not intended to reproduce a whole biological fiber but they rather produce a discretization of such fibers into rigid elements that can be computationally manipulated. At initialization actual biological fibers may be too long to be represented by a single short fiber element and need to be described by structures made of pre-assembled fiber elements. The chosen diamond shape structure is the simplest of such pre-assembled structure and it allows us to compare the results to those using large fibers as it produces similar potential isolines to the latter.

Now, to take into account fiber-fiber repulsion, we proceed as for fiber-cell repulsion as described in section B.1 and we suppose that the fiber repulsion potential $W_{rep}$ consists of the sum of two-particle potential elements

$\mathbf{W}_{rep}(Y_{k},\theta_{k},Y_{f},\theta_{f})$ modeling the mechanical interaction between fiber $f$ and fiber $k$:

[TABLE]

For two given vectors $Y_{1}$, $Y_{2}\in$ $\mathbb{R}^{2}$ and two angles $\theta_{1},\theta_{2}\in[-\frac{\pi}{2},\frac{\pi}{2}]$, $\mathbf{W}_{rep}=\mathbf{W}_{rep}(Y_{1},\theta_{1},Y_{2},\theta_{2})$ reads:

[TABLE]

where $d(Y_{1},Y_{2},\theta_{2})$ is the distance function defined by Eq. (6) and $W_{f}$ is a model parameter. The fiber-fiber repulsion potential element $W_{rep}$ iso-lines are ellipses with locii located at the two ends of the fiber $(Y_{2},\theta_{2})$. The potential vanishes beyond distance $d_{\xi}$ to the fiber center $Y_{2}$. Parameter $d_{\xi}$ is set such that the length of the semi minor axis is $d_{0}$. In this case, a fiber repels an other fiber up to distance $d_{0}$ in its orthogonal direction. A direct computation gives:

[TABLE]

We consider $N_{f}=800$ fibers of length $L_{f}=1.7$ (note that previously $L_{f}=3$). The fiber-cell interaction range is set to $\tau R_{i}=R_{i}+d_{0}$ with $d_{0}=0.17$ (see Eq. (7) and Fig. 2). Note that previously, $d_{0}=1.32$. When fiber-fiber repulsion is activated, fibers repel their centers to a distance $d_{0}=0.17$ in their orthogonal direction, and we use the value $W_{f}=10$. If not otherwise stated, the other parameters correspond to the simulations of the main text. We consider an initial fiber network composed of $N$ groups of $4$ fibers organized in ’diamond’ configurations (see Fig. 18). The center of gravity of the diamonds is chosen randomly with uniform distribution in the computational domain, and the major axis orientation of the diamonds is randomly chosen in $[-\pi,\pi]$ with uniform probability.

In Fig. 19 (I), we show simulation results for $N_{f}=800$, linked fiber fraction $\chi_{\ell}=0.35$, small fiber unlinking frequency $\nu_{d}=10^{-3}$ for two different fiber lengths: $L_{f}=1.7$ (Figs. 19 (I A) and (I B)), $L_{f}=1.2$ (Figs. 19 (I C) and (I D)) with random initial fiber network (Figs. 19 (I A) and (I C)) and with initial diamond configuration (Figs. 19 (I B) and (I D)). In Figs. 19 (II), we show the values of the statistical quantifiers $E$ (cell cluster elongation (II A)), $A$ (fiber alignment (II B)), $\Lambda$ (fiber cluster elongation (II C)) and $N_{C}$ (number of cell clusters (II D)), as functions of $\nu_{d}$ for $L_{f}=1.7$ and $L_{f}=1.2$ (black and red curves), with random initial fiber network (continuous lines) and with diamond initial configuration (dashed lines).

As shown by Figs. 19 (II), the type of initial fiber configuration does not seem to influence the fiber network (compare continuous and dashed lines of Figs. 19 (II B and C)) but has a strong influence on the cell structures. Indeed, for $L_{f}=1.2$, the mean elongation of cell clusters $E$ undergoes a phase transition as $\nu_{d}$ increases when the initial fiber network is composed of diamond structures. This phase transition is not observed when the fibers are initially randomly distributed. This is due to the fact that for small $\nu_{d}$, due to its internal structure, the diamond like network is less compliant than the random one. As a result, the former exerts more pressure on the cells than the latter, favoring the emergence of small and round clusters of cells. For large $\nu_{d}$, i.e fast remodelling of the fiber network, the internal diamond structure is quickly lost, and we recover the same values of the SQ for diamond and random initial networks.

To sum up, the internal diamond structure of the fiber network increases its rigidity compared to a random initial configuration. Moreover, there exists a critical fiber unlinking frequency for which the fiber network imposes local directional constraints to cell cluster growth, favoring cell cluster elongation. With such an initial diamond structure and $L_{f}=1.2$, we recover the results of the main text with a random initial network and $L_{f}=3$. Note however that the fiber network is less organized with small fibers than with larger ones (compare the values of $A$ and $\Lambda$ between Figs. 19, Fig. 2 of the main text and 14). This is due to the fact that the fibers of the diamonds are not strictly aligned, and this effect adds up with the fiber-fiber repulsion potential which disorganizes the fiber network.

Finally, we show the influence of the linked fiber fraction when fiber-fiber repulsion is activated, for $L_{f}=1.2$ and with diamond-like initial configuration for the fiber network. In Fig. 20, we show contour plots of the SQ $E$ (cell cluster elongation, Plot A), $A$ (fiber alignment, Plot B), $\Lambda$ (fiber cluster elongation, Plot C) and $N_{C}$ (number of cell clusters, Plot D), obtained when varying the linked fiber fraction $\chi_{\ell}$ (vertical axis) and the fiber unlinking frequency $\nu_{d}$ (horizontal axis).

As shown in Fig. 20, we recover the morphologies observed in the main text (A) lobule-like cell clusters surrounded by a disorganized (unaligned) fiber network, (B) lobule-like cell structures in an aligned fiber network, and (C) elongated cell structures in a network composed of long and rigid fiber threads. We also recover the phase transitions in $\nu_{d}$ and $\chi_{\ell}$ for the statistical quantifiers $E$ and $A$ that we observed in Fig. 2 of the main text (with no fiber-fiber repulsion, $L_{f}=3$ and random initial network). The most biologically relevant structures, composed of well-separated lobule-like cell clusters in an organized fiber network, are here again identified around $\nu_{d}\in[10^{-2},10^{-3}]$ for $\chi_{\ell}\geq 0.35$. Note however that the fiber network is less aligned than in the simulations of the main text.

To sum up, adding a fiber-fiber repulsion potential does not lead to a breakdown of the results previously obtained without such a repulsion potential. On the contrary, it reinforces them by making them more robust to parameter changes. These results show that the fiber network of the main text can be seen as a fiber network composed of initially connected smaller fibers which repulse each other.

Appendix G Computational time of the model

The individual based model was coded in FORTRAN90 in sequential mode for each simulation, and sets of simulations were run in parallel for the statistical analysis. The computation time was optimized by the use of a numerical grid to reduce the size of the computational domain around each agent, and one simulation took from few hours to several days as functions of the model parameters. In Fig. 21, we show the computational times -in hours- for simulations of Appendix F, i.e when considering fiber-fiber repulsion, for $\chi_{\ell}=0.1$ and three different values of $\nu_{d}$: $\nu_{d}=10^{-3}$ (blue curves), $\nu_{d}=10^{-2}$ (red curves) and $\nu_{d}=0.1$ (yellow curves). The computational time corresponds to the time a simulation takes to reach the steady state computed at $t=80t_{ref}$, and is measured as the difference between the times at which the data files of the final and initial states are written by the program in the data folder. We plot the computational time as a function of the fiber total area (top figure), and of the number of fiber links (bottom figure). The fiber total area $A_{T}$ corresponds to the area the fibers would occupy if it was in optimal (non-overlapping) configuration, i.e it is computed as the sum of all fiber areas:

[TABLE]

where $d_{0}$ is the length of the semin minor axis of the ellipse around fiber $f$, and $\frac{L_{f}}{2}$ the length of its semi major axis (see Appendix F). Simulations of Figs. 21 are produced with varying the number of fibers $N_{f}=\{200,300,400$ $,500,600,800,1000\}$, the fiber length $L_{f}=\{1.2,2.4,3\}$ and $d_{0}=\{0.12,0.24,0.3\}$.

As shown in Figs. 21, the computational time decreases with the fiber total area $A_{T}$ for all values of $\nu_{d}$, while it increases as a function of the number of fiber links in an almost linear way. This shows that the computational time depends non trivially on the model parameters. As the total area occupied by the fibers increases (by increasing either the number of fibers either each fiber repulsion area), the tissue becomes more and more constrained by the fiber network, and displays fewer degrees of freedom: fibers and cells have less space to move. As a result, the system is closer to its optimal state at each discrete time of the simulation, and this can result in a decrease of the number of iterations in the minimization algorithm. We stress the fact that the computational time displayed in Figs. 21 is the total time of a simulation. This includes the time to compute the interactions between each pair of interacting agents, as well as the time to perform the minimization process between each time steps. As a result, this CPU time is not straightforwardly linked to the number of agents as could be expected. Here, the more constrained the system, the less self-organization is possible, resulting in a small computational time. Finally in the bottom figure of Fig. 21, we show the computational time of the model as a function of the total number of fiber links in the model. As expected, the computational time increases with the number of fiber links, in an almost linear way. This reflects the larger time needed to compute the interactions when increasing the number of fiber links. The simulations of the main text or appendix F correspond to a fiber area approximately equal to $50$, with a mean computational time of 50 hours for $800$ fibers of length $L_{f}=1.7$.

Appendix H Supplementary Informations: Movies

H.1 S1 Video

Simulation: Fixed number of fiber links ($\nu_{d}=\nu_{\ell}=0$), moderately constrained fiber network ($40\%$ of links on the population of crossing fibers). Cells (2D spheres) are represented in red, fibers (straigth lines) in blue and fiber links as green crosses. This video shows a simulation with a fixed number of fiber links, for a moderately constrained fiber network. Cells and fibers self-organize to form lobule-like structures of cells in a stiff fiber network (moderately organized). Biologically relevant cell and fiber structures are obtained, with a fiber network slightly more disorganized than in real tissues.

H.2 S2 Video

Simulation: Fixed number of fiber links ($\nu_{d}=\nu_{\ell}=0$), highly constrained fiber network ($70\%$ of links on the population of crossing fibers). Cells (2D spheres) are represented in red, fibers (straigth lines) in blue and fiber links as green crosses. This video shows a simulation with a fixed number of fiber links, for a highly constrained fiber network. Cells are repelled into zones of low fiber density and newly differentiated cells are then attracted into these zones. This results in the formation of well-separated lobule-like structures of cells. However, the fiber network is disorganized because of the large number of fixed fiber links. Biologically relevant cell structures are obtained, in a fiber network which could represent a pathological situation such as fibrosis.

H.3 S3 Video

Simulation: Dynamical number of fiber links ($\nu_{d},\nu_{\ell}\neq 0$), moderately constrained fiber network ($40\%$ of links on the population of crossing fibers initially), moderate linking/unlinking frequency. Cells (2D spheres) are represented in red, fibers (straigth lines) in blue and fiber links as green crosses. This video shows a simulation with a moderate linking/unlinking frequency. Cells and fibers self-organize to form lobule-like structures of cells in an organized fiber network. The dynamical linking/unlinking favors the organization of the fiber network compared to a fixed number of fiber links (to be compared with H.1). Biologically relevant cell and fiber structures are obtained.

H.4 S4 Video

Simulation: Dynamical number of fiber links ($\nu_{d},\nu_{\ell}\neq 0$), moderately constrained fiber network ($40\%$ of links on the population of crossing fibers initially), fast linking/unlinking. Cells (2D spheres) are represented in red, fibers (straigth lines) in blue and fiber links as green crosses. This video shows a simulation with fast fiber linking/unlinking. The fiber network is more rigid in this regime compared with the case of a slow linking/unlinking process. A preferred fiber direction locally emerges and generates elongated cell clusters. Such elongated structures could model other organs, such as muscles.