Introduction

Regenerating normal tissue structure after injury is a central goal of regenerative medicine in adult mammals. Overall, striking regenerative properties are present in newborn mammals, but this capability rapidly declines and disappears a few days after birth in most mammals1. This loss is a major issue in spontaneous tissue repair since injury in adults usually leads to scarring rather than regeneration. While scarring rapidly blocks bleeding and reconstitutes a protective barrier2, it has detrimental effects on tissue function. These adverse effects are due to the lack of recovery of tissue architecture and function due to inadequate extracellular matrix (ECM) structuring3.

Classically, the identification of a putative therapeutic target is based on the identification of elements that differentiate regenerative from non-regenerative repair. However, these two repair outcomes are usually studied by comparing different species. Adult spiny mice (Acomys Cahirinus), which show regeneration in adults, have been increasingly used to investigate scar-free regeneration compared to Mus Musculus, the most commonly used mammal in research4. Recently, Sinha et al. used the different regenerative capabilities of reindeer antler velvet compared to back skin, which forms fibrotic scars. Their study underscored the significance of the interactions between fibroblast and immune cells in the process of tissue repair5. Despite the power of these approaches, they can lead to context-specific insights rather than identifying factors specific to repair processes across tissues and species. Identification of regenerative medicine treatments based on these comparative studies may, therefore, fail to activate dormant regenerative capacities in adult mammal tissues across species. To overcome these limitations, we developed an original inducible model of adipose tissue (AT) regeneration in adult mammals that provides a convenient model where regeneration and scarring can be investigated in the same tissue and developmental stage in animals with a shared genetic background6,7,8,9. In this model, a large resection of subcutaneous AT spontaneously drives tissue repair toward scar healing that can be switched toward tissue regeneration following treatment with an antagonist of opioid receptors6.

The complete recovery of tissue architecture during tissue regeneration depends on mechanical interactions between cells and ECM fibrillar components, which provide spatial information for cells. ECM mechanical properties condition its organization which is highly dynamic and results from a combination of three mechanisms: its composition/synthesis, its degradation, and post-translational modifications such as temporally dynamic cross-links between fibrillar components. However, the dynamics and nature of these changes in adult tissue after injury requires further systematic investigations. It is noteworthy that the combined actions of the three mechanisms to generate ECM mechanical forces leading to tissue organization make it difficult to determine the respective importance of each mechanism using in vivo experiments.

To address this complex biological question in a simpler way, and understand the regeneration process as a whole, we previously developed a computational agent-based model (ABM), composed of two agents (segments and growing circles modeling ECM fibers and differentiating cells from AT respectively). This model qualitatively reproduced the tissue repair process and allowed the identification of key factors that control the emergence of tissue architecture such as mechanical cues10,11. It integrated three critical parameters corresponding to the three mechanisms involved in ECM organization in vivo: frequency of fiber synthesis, unlinking, and probability of cross-linking. Changing parameters showed that simple mechanical interactions between adipocytes and fibers may support the emergence of a lobular tissue architecture, similar to the one observed in vivo10,12. These results reinforced the relevance of this computational approach to determine the respective importance of each parameter and to generate new hypotheses to identify in vivo key biomechanical mechanisms controlling the emergence of AT architecture.

Here, we took a step further in the analysis of this computational model to use it as a predictive model of tissue repair. This model mimics both regeneration and scar healing-like architectures. We performed high throughput multiparametric simulations (i.e. in silico experiments) in which ECM parameters were independently changed over a wide range of values and used Machine Learning (ML) models to reveal that ECM cross-linking was the most important parameter explaining tissue repair outcome in silico. The computational model predicted that an early and transient decrease of ECM cross-linking after an injury could be necessary and sufficient to drive tissue regeneration. The computational model-based hypotheses were validated once the time calibration had been carried out. Indeed, by demonstrating that regeneration in adult mammals is unlocked by early and transient inhibition of ECM cross-linking, whereas increased cross-linking drives repair towards scar healing, this work positions the computational model as a predictive tool for tissue regeneration.

Results

A computational 2D model for tissue repair

In our previous conceptual model, a tissue lesion was modeled by the removal of segments (ECM fibers) and circles (adipose cells) (Fig. 1a). The subsequent repair process was obtained by seeding new fibers and the growth of new cells. Cells and fibers interact through mechanical repulsive forces, but when two fibers intersect, they can also form cross-links with a defined probability (Pf), and can thus resist the pressure of growing cells. We mimicked ECM remodeling by setting up fiber synthesis (νf), cross-linking (with a probability, Pf) and unlinking (νad), thus controlling both the density and cross-linking of fibers (Fig. 1a right panel and supplementary Fig. 1a). As shown in our previous theoretical study11, the model allowed to spontaneously reconstruct different types of tissue architectures after injury, by enabling the restoration of a fiber network from which both cell clustering and tissue architecture could emerge. To evaluate tissue reconstruction at the end of simulations, a repair index was calculated that captured the numbers and spatial organization of cells and fibers in the reconstructed tissue compared to the initial non-injured tissue (Supplementary Fig. 1b, c). This repair index \(\Gamma\) was defined in Eq. 1 using morphological quantifiers of the cell and fibrous structures before and after injury:

$$\Gamma =-\sqrt{{\left(\frac{\Delta E}{{E}_{0}}\right)}^{2}+{\left(\frac{\Delta A}{{A}_{0}}\right)}^{2}+{\left(\frac{\Delta {N}_{C}}{{N}_{{C}_{0}}}\right)}^{2}+{\left(\frac{\Delta {N}_{E}}{{N}_{{E}_{0}}}\right)}^{2}}$$
(1)

Where \(\varDelta E,\varDelta A,\varDelta {N}_{C},\,\varDelta {N}_{E}\) corresponded to the differences, before injury, and after reconstruction, in cell cluster elongation, fiber network alignment, number of clusters and number of cells, respectively, and \({E}_{0},{A}_{0},{N}_{{C}_{0}},\,{N}_{{E}_{0}}\) were scaling parameters.

Fig. 1: In silico modeling of AT repair outcomes.
figure 1

a Cellular (left panel) and fibrillar (right panel) phenomenons of the agent-based in silico model of AT repair. Rmin and Rmax correspond to the minimal and maximal radii of cells respectively. b Simulation results after injury: scar repair is defined by a repair index smaller than −2 (lower panel), unstructured tissue is defined by a repair index from −0.7 to −2 (middle panel), and regenerative repair is defined by a repair index above −0.7 during tissue repair steps (upper panel).

Full size image

The model gave rise to three distinct classes of repair outcomes11, regenerative repair (repair index from 0 to - 0.7) was characterized by a reconstructed tissue similar to the non-injured tissue especially in terms of number of adipocytes and lobular organization. In contrast, unstructured repair (repair index from −0.7 to −2) did not reproduce lobular organization. Scar repair (repair index from −2 to −4.5) was characterized by a high number of fibers and a few adipocytes (Fig. 1b). In this previous study, a large sensitivity analysis of the parameters and a temporal calibration of ABM were lacking precluding to use it as a predictive model for in vivo tissue repair.

Fiber cross-linking plays a major role in tissue repair outcome prediction

A sensitivity analysis was performed by generating a large panel of parameter combinations, leading to a dataset of 2218 simulations. In order to consider the stochastic processes of the computational model, each combination of the parameters νf (1 to 9), Pf (0.1 to 0.9), and νad (0.001 to 0.1) was simulated six times. These 2218 simulations provided repair index values ranking from 0 to −4.5 (243 of them being shown in 3D space, Fig. 2a and supplementary Fig. 1d). The sensitivity study, represented on the 3D scatter plot, revealed that with a low fiber unlinking frequency (νad), the combination of specific Pf and νf values were associated with the three repair outcomes. In contrast, a high frequency of fiber unlinking was mainly associated with regenerative repair (Fig. 2a).

Fig. 2: In silico modeling: fiber cross-linking plays a major role in tissue repair outcome prediction.
figure 2

a 3D scatter plot of repair index over a wide range of values of three parameters used in in silico modeling (Pf, νf, νad). b ML pipeline for the predictive model. c RF SHAP-values when varying each parameter for a regenerative outcome. The color in each graph represents the value of the associated parameter – red for high values and blue for low values of the parameter. A positive SHAP value means that the parameter contributes positively to the prediction whereas a negative SHAP value means that the parameter contributes negatively to the prediction. The parameter ranking indicates the importance in the class (repair outcome) prediction: the parameter at the top is the most important for predicting the class. d RF SHAP-values according to each parameter for the scarring class. e RF SHAP dependence plots for parameter values in the regenerative class. f RF SHAP dependence plots of each parameter value for the scarring class.

Full size image

Because fiber organization is the result of the combined actions of the three ECM-related parameters, the respective contribution of each parameter is hard to define. To tackle this complex issue, we applied an explainable ML strategy to reveal the contribution of each parameter by itself considering its interaction with each of the other parameters. Simulations were randomly split into training (80%) and test (20%) sets according to classes (regenerative, unstructured, and scar repair) (Fig. 2b). Both sets were balanced among the different classes (supplementary Fig. 2a). To demonstrate the benefit of a more complex ML model that captures complex and non-linear relationships between parameters, a tree-based ML model Random Forest (RF) was compared to a linear model with regularization (ElasticNet). Models were trained and class predictions were compared (regenerative, unstructured, scar). RF predicted the different repair classes with a better accuracy than ElasticNet (79% versus 72%) (supplementary Fig. 2b). Analyzing confusion matrices, RF results were clearer than ElasticNet showing the benefit of using a tree-based algorithm instead of a linear inspired-model (supplementary Fig. 2c, d). To define the contribution of each variable in repair outcomes predictions, the Shapley Additive exPlanations framework (SHAP) was applied to RF. SHAP graphs represent explanations for regenerative (Fig. 2c) and scar (Fig. 2d) repair class predictions. Parameters (νf, Pf and νad) were ranked according to their importance in the prediction, the parameter at the top being the most important to predict the class. RF results showed that Pf, νad and νf values did not all have the same weight in the prediction of tissue repair outcomes (Pf 43%, νad 32% and νf 24%) (supplementary Fig. 2e). These analyses suggested that Pf and νad, both impacting the number of fiber cross-links, were the most important analyzed parameters in predicting tissue repair outcome.

Low values of Pf and νf and high values of νad had a positive impact on regenerative class prediction (Fig. 2c). More precisely, setting Pf from 0.1 to 0.2, νad from 0.05 to 0.1 and νf from 1 to 3 seemed to systematically drive tissue repair toward regenerative repair (Fig. 2e). In contrast, high values of Pf (from 0.7 to 0.9) and νf (from 7 to 9) as well as low values of νad (from 0.001 to 0.01) had a positive impact on the scar repair class prediction (Fig. 2f).

We next evaluated the temporal changes in fiber cross-linking resulting from interactions between Pf and νad. The number of cross-links over time was recorded. This demonstrated that the three repair outcomes followed the same inverted exponential curve. We found that regenerative repair curve was associated with low number of cross-link values. Our analysis indicated that the outcome of tissue repair could be driven by modulating the number of fiber cross-links in a short temporal window since these values plateaued in the earliest time points of simulations (i.e. in the first Epochs) (supplementary Fig. 2f).

Early and transient modulation of fiber cross-linking is critical to drive tissue repair outcomes in silico

We next temporally calibrated the ABM, by using our in vivo model of AT regeneration6. In the ABM, the non-injured tissue is composed of lobular structures filled with adipocytes. In vivo, after an injury, the wound is filled with fibers within 3 days and adipocytes emerge from 10 days (supplementary Fig. 3). We selected simulations for which the combination of parameters (νf, Pf and νad) allowed wound closure and adipose cell emergence in a timeframe that matched with in vivo observations, which constituted the calibrated model. After this calibration, the number of fiber cross-links plateaued at 6 days post-injury (Fig. 3a), suggesting that day 0 to day 6 post-injury might constitute the critical temporal window during which fiber cross-linking could be modulated to drive tissue repair outcome.

Fig. 3: In silico modeling: Identification and validation of an early and transient treatment window where fiber cross-linking modulation can drive tissue repair outcomes.
figure 3

a Total number of cross-links over time upon regenerative repair, unstructured tissue, and scar repair conditions (n = 7 for each condition). b Simulation results after injury for a Pf=0.7 correspond to scar repair. A transient Pf decrease (Pf = 0.1) during the first 6 days leads to regenerative repair. c Simulation results after injury for a Pf=0.2 correspond to regenerative repair. A transient Pf increase (Pf = 0.8) during the first 6 days leads to scarring repair. d Histogram of the repair index when Pf is constant throughout the simulation (scar and regenerative repair) or for a transient Pf decrease (Pf = 0.1) or increase (Pf=0.8) during the first 6 days. N = 6 simulations for each condition. e Total number of cross-links over time upon regenerative and scar repair or for a transient increase or decrease in Pf. Fixed values of νf = 5 and νad=10-3 were used for all simulations. n = 6 for each condition. Data are expressed as mean +/− SEM and analyzed by Tukey’s test. ****p < 0.0001.

Full size image

We thus tested whether an early and transient change in Pf value was sufficient to reverse the final expected outcome of scar or regenerative simulations. To this end, we used the specific Pf values obtained with SHAP analysis (Fig. 2e, f). Decreasing Pf from 0.7 to 0.1 during the first 6 days after injury switched tissue repair trajectories from scar to regenerative repair (Fig. 3b). In contrast, increasing Pf from 0.2 to 0.8 during the first 6 days after injury switched tissue repair trajectories from regenerative to scar repair (Fig. 3c). These results were confirmed by the quantification of the repair index (Fig. 3d). As expected, changes in tissue repair outcome were associated with a switch in the profiles of the curves representing number of fiber cross-links (Fig. 3e). Taken together, these in silico experiments performed on a calibrated computational tissue model show that early and transient modulation of fiber cross-linking could be sufficient to guide the outcome of tissue repair.

Temporally specific modulation of fiber cross-linking is sufficient to guide tissue repair outcomes in vivo

According to the in silico experiments, we then tested whether an early and transient decrease or increase in fiber cross-linking led to tissue regeneration or scarring respectively using the previously validated in vivo model of AT repair6,8. Indeed, as previously described, treatment with Vehicle or Naloxone Methiodide (Nal-M), an antagonist of opioid receptors, after partial resection of subcutaneous AT leads 1-month post-injury to scar healing or regeneration respectively6,8. During the 6 days following AT resection, Vehicle (scarring) and Nal-M (regenerative) treated mice were treated with beta aminopropionitrile (BAPN) or with Genipin, drugs that inhibit and enhance fiber cross-linking, respectively.

One month after resection, AT repair was analyzed by macroscopic evaluation, light transmittance and quantified fluorescence experiments. BAPN significantly induced AT macroscopic regeneration in scarring mice (Fig. 4a) and a decrease in light transmittance compared to vehicle treated mice, suggesting the presence of new cells in reconstructed area (Fig. 4b). In addition, fluorescence of 300 µm thick tissue sections revealed the presence of mature adipocytes organized as lobules and of a new network of blood vessels in the reconstructed area of BAPN-treated mice (Fig. 4c, Supplementary Fig. 4). This was confirmed by a higher percentage of reconstruction in BAPN compared to vehicle treated mice (Fig. 4d). In contrast, Genipin inhibited macroscopic regeneration and induced AT scarring in regenerative mice (Fig. 4a), associated with an increase in light transmittance, suggesting the absence of new cells in the reconstructed area (Fig. 4b). In addition, reconstructed area was characterized by the presence of a dense collagen fiber network as well as the absence of adipocytes and blood vessels (Fig. 4c, Supplementary Fig. 4), as well as a lower reconstruction percentage (Fig. 4d).

Fig. 4: In vivo validation: Early and transient fiber cross-linking modulation guides tissue repair outcomes.
figure 4

a Morphological image of AT after scarring repair (Vehicle), upon transient inhibition of cross-links (BAPN), under regenerative conditions (Nal-M), and upon transient increase of cross-links (Nal-M + Genipin). Images were taken 1-month post-injury. b Box-plot of light transmittance after Vehicle, BAPN, Nal-M, and Nal-M + Genipin treatments. c Representative fluorescence (300 µm thick tissue sections) images of regenerative (Nal-M and BAPN) and scarring conditions (Vehicle and Nal-M + Genipin) 1-month post-injury. The two lines represent the reconstructed area. Red, blue, and green colors correspond to adipocytes, nuclei, and collagen respectively. Scale bars represent 2,000 µm (panel a) and 100 µm (fluorescence images). d Histogram of the reconstructed area (%) based on the ratio of adipocyte area compared to the whole reconstructed area upon regenerative and scarring repair. e Boxplot of fractal dimension under regenerative, scarring, and BAPN conditions 3 days post-injury. Data are expressed as mean +/- SEM. Data analyzed by Tukey’s test. *p < 0.05; **p < 0.01; ****p < 0.0001.

Full size image

To evaluate ECM quantity which reflect both synthesis and degradation, we quantified total collagen content in the reconstructed area of our tissue samples 3 days post-injury in both regenerative and scarring conditions. We did not observe any difference in collagen content between regenerative or scarring conditions (supplementary Fig. 5). To evaluate ECM organization during treatment windows, we then calculated Fractal Dimension (FD) values on Second Harmonic Generation (SHG) images 3 days post-injury. FD describes the amount of space and self-similarity of structure and is highly sensitive to changes in collagen organization13,14,15. FD values were significantly higher in regenerative than scarring mice. Moreover, ECM cross-linking inhibition through BAPN treatment was sufficient to switch FD values of scarring mice close to regenerative ones (Fig. 4e). These data showed that regeneration and scarring were characterized by specific and distinct early ECM organization patterns that reside in collagen cross-linking consistent with the conclusions of our computational model.

Taken together, these data demonstrated that early and transient modulation of fiber cross-linking guides tissue repair outcomes, both in our in vivo and in silico models.

Discussion

Our study provides biological validation of innovative hypotheses derived from a calibrated computational model used as a digital twin. We validated in vivo that regeneration in adult mammals can be induced by early and transient inhibition of ECM fiber cross-linking after injury. The ability to perform high throughput in silico experiments before implementing in vivo experiments for biological validation greatly reduced the time and number of animals needed to explore the experimental space.

In this study, the initial conceptual computational model was deeply analyzed and calibrated to develop a predictive version. In vivo, ECM organization results from the combined actions of the three mechanisms corresponding to the three ECM-related parameters in our computational model. Our computational model allowed us to explore how ECM mechanical properties led to tissue organization by studying the respective importance of each initial parameter over a wide range of values and temporal window. The systematic investigation and the calibration of the conceptual model led to an independent and temporal exploration of the parameters. The calibration was based on morphological observations of tissue repair process (wound closure and adipocyte emergence timeframes) in the previously described in vivo model of AT regeneration6,8, since reliable methods to measure cross-linking and/or crosslinking activity biologically are currently lacking16. It is interesting to note that, although our computational model is in 2D while the tissue organization is in 3D, we were able to validate the computational hypotheses in vivo. In a next step, a 3D model could be envisioned to better investigate the regeneration of 3D complex structures.

Our in silico results showed that early dynamic ECM linking rates are associated with the final repair outcome. In vivo either too much collagen or too many cross-links or both would inhibit tissue regeneration. We did not observe any difference in collagen content between regenerative or scarring conditions and showed that collagen organization is different between regenerative and scaring conditions with a similar organization between our two regenerative conditions (Nal-M and BAPN). These observations suggested that the difference between regeneration and scarring resides more in collagen cross-linking strengthening the conclusions of our computational model. ECM cross-linking is known to condition ECM topology and physical properties and thus facilitates ECM stabilization17. Dynamic remodeling by changing the connectivity of the ECM network can confer either fluid or solid-like properties to tissues, as previously described in the development of embryos and adult tissues18. In all these studies, a weakly-cross-linked network behaves like a fluid and allows cell rearrangement within the ECM. This transition from a solid to fluid-like behavior also occurs at the tissue level during lung branching morphogenesis, where increased local fluidity enables branching to expand19. In contrast, extensive collagen cross-linking is observed in tissue fibrosis which results in an increasingly stiff and less compliant matrix20. The importance of ECM cross-linking in repair processes is also consistent with previous studies showing that uterus regeneration only occurs within an appropriate range of decellularized matrix cross-links number17.

Although this connectivity appears to be fundamental to the repair outcome, the timing of the interaction is also critical. Our computational tissue repair model revealed that early, transient and tight restriction of ECM cross-link numbers within a given range is critical and sufficient to allow cells to self-organize into an optimal architecture and guide tissue repair towards regeneration.

Our results also indicate that tissue reconstruction after injury seems to result from a self-organization of cells from an initial morphological template, as it has recently been described in skin development18. Indeed, tissue reconstruction in adult mammals seems to be a self-organizing/emergent process where an initial scaffold leads spontaneously to scarring. Altering this initial scaffold sets the stage for the self-organization of cells and recovery of the tissue architecture existing before injury.

Previously, we demonstrated that regeneration in adult tissues was observed after transient and early treatment with an antagonist of endogenous opioids during the 3 days following injury6. Similar results were also obtained on mouse pancreas confirming our conclusions21. We also showed that inhibition of resident macrophages and management of the post-injury inflammatory phase could explain the regenerative effect of endogenous opioids inhibition8. The present study suggests that similar regenerative effects can be obtained without directly managing inflammation and opioid signaling. In all cases, our studies strongly suggest that tissue regeneration is possible but inhibited in adult mammals6 and the present study opens new therapeutic approaches targeting ECM cross-linking while preserving pain management.

In conclusion, our discovery relies on three elements: (i) a biological model relevant to analysis of repair processes in adult mammals, (ii) a calibrated computational model of the biological tissue, (iii) dynamic back and forth between in silico and in vivo models. Taken together, these three elements make our model close to the definition of a digital twin to reveal the complex biomechanical cues controlling tissue architecture emergence, impairment, and recovery after injury. The strategy of generating mechanistic models to formulate sets of mechanisms compatible with literature knowledge and implement them on the computer (i.e. build digital twin candidates) has already been developed by many authors and successfully helped to gain biological insights into real systems22,23,24,25,26,27. Finally, our study highlights the management of cross-linking during the first steps after injury as a relevant target for plastic and reconstructive surgery and regenerative medicine. In addition, our work reveals the relevance of computational strategies to delineate fine biological processes with a minimal number of in vivo experiments.

Methods

Animals

All experiments in mice were performed on 7-week-old male mice. C57BL/6 OlaHSD mice were obtained from Envigo Laboratories and housed in a conventional facility at the Faculty of Pharmacy in Toulouse. Animals were group-housed (3 to 5 per cage) in a controlled environment (12-hour light/dark cycles at 21 °C) with unrestricted access to water and a standard chow diet. Animals were maintained in accordance with the ARRIVE guidelines of the European Community Council. Mice were killed by cervical dislocation. All experiments were carried out in compliance with European Community Guidelines (2010/63/UE) and approved by the MESRI (Ministry of Higher Education, Research and Innovation) and the French ethics committee (protocol reference: 40880-2023010612201954 v6).

Subcutaneous AT resection

Control mice used for the baseline control were not subjected to surgery. For mice that underwent unilateral resection of subcutaneous AT, animals were anaesthetized by inhalation of isoflurane 2.5%. A single abdominal incision was then made to access and excise 35% of the right AT between the lymph node and groin, the skin was closed with 3 suture points.

In vivo treatments

Mice were treated with 50 µl of naloxone methiodide (Nal-M) (daily subcutaneous injections, 17 mg/kg, N129, Sigma Aldrich) or NaCl (daily subcutaneous injections) from day 0 to 3 after AT resection. For cross-link regulation, mice were treated with 50 µl of lysyl oxidase inhibitor, β-aminopropionitrile (daily subcutaneous injections, 150 mg/kg, A3134, Sigma Aldrich) from day 0 to 6 after AT resection. Genipin treatment was administrated daily as a co-injection with Nal-M during the first 3 days after AT resection and Genipin only (daily subcutaneous injection, 5 mg/kg, G4796, Sigma Aldrich) from day 4 to 6 after AT resection.

Regeneration evaluation by light transmittance analyses

One-month post-resection, mouse AT were fixed (PFA 4%, 24 h) and placed on a coverglass. The samples were observed through a binocular microscope and pictures were taken using a 48MP camera (Nikon, D5000). All images were opened in ImageJ software and the contrast was normalized using the histogram equalization function. Regions of interest (ROI) were drawn to include the reconstructed area and the mean gray value in this ROI was calculated. Higher mean gray values mean that the tissues were not reconstructed since the light can go through the reconstructed area and lower values were associated with reconstructed tissue where the reconstructed area is darker thanks to the presence of new adipocytes.

Fluorescence and SHG imaging

Mice AT were fixed (PFA 4%, 24 h), embedded in agarose gel (Sigma A0169), and cut into 300 µm tissue sections with a Vibroslice for Campden Instruments. Tissue sections were incubated for 1.5 h in PBS/0.2% triton at room temperature and then 2 h in Bodipy 493 (Invitrogen 03922) in the dark. Tissue sections were incubated for 2 h in Draq5 (Thermofisher, 62251) and 2 h in lectin (Eurobio Scientific, RL1102). The Bodipy concentration was used according and the manufacturer instructions, Draq5 and Lectin were used at 5 mM and 5 µg/ml respectively. Imaging was performed using a Biphotonic Laser Scanning Microscope (LSM880 Carl Zeiss, Jena, Germany) with an objective lens LCI ‘Pan Apochromat’ 10X/0.45 or ‘LD C-Apochromat’ 40x/1.1. Bodipy, lectin, and Draq5 were excited using 488, 561, and 633 nm lasers, and SHG was excited using the 800 nm biphoton laser.

Fractal dimension (FD) analyses

For FD analysis, SHG tiles were individually opened in ImageJ, and background subtractions were performed before binarizing all images from the Z-stack. FD was then estimated using the box-counting method in ImageJ. The software considered box-counting in two dimensions, allowing quantification of pixel distribution in this space. The FD was based on a series of grids with different sizes (boxes) over an image and the recorded data (counting) for each successive box size. The results were expressed as the FD of the object that is DF = Log N/Log r; N, where N is the number of boxes that cover the pattern, and r is the magnification or the inverse of the box size. FD were thereby calculated using the ImageJ software set between 0 and 2, with 0 corresponding to an image without signal (0 pixels) and 2 to an image full of signal.

Picrosirius Red (PSR) Staining

AT were collected from non-injured mice 3 and 10 days after injury and fixed (PFA 4%, 24 h). AT were dehydrated using the following successive alcohol treatment protocol: 10 min in 70% EtOH, 15 min in 95% EtOH, 35 min in 100% EtOH, and 35 min in Bioclear solution. Then AT samples were embedded in paraffin for 24 h. 5 µm sections were cut, deparaffined by successive baths in Bioclear and 100% EtOH, and then rinsed in distilled water. Tissue sections were incubated in Fast Green/Citrate buffer solution (0.04%) for 15 min, rinsed in distilled water, incubated in Fast Green/Picrosirius Red (0.1%) for 15 min, washed in distilled water, and dehydrated in 2 min EtOH and 10 min Bioclear solution and mounted with Eukitt Mounting Medium. Tissue sections were imaged using Lyon Platform Imaging (CIQLE) with a Zeiss AxioScan 7 microscope.

Total collagen assay

AT were collected from Vehicle, BAPN, Nal-M, and Nal-M+Genipin treated mice 3 days after injury. AT were placed in a tube with a stainless bead and frozen. They were pulverized 2 times 2 min at 25 Hz with the Tissue Lyser (Retsch MM300). Then pulverized tissue was digested overnight at 4 °C with pepsine solution (Sigma P7012) 0.1 mg/ml in acetic acid 0,5 M. After centrifugation (10 min, 12,000 rpm, 4 °C), the supernatant containing solubilized collagen was collected. Collagen quantification was performed with the Sircol collagen assay (Biocolor S1000) according to the manufacturer instructions and the absorbance measurement was done at 540 nm with Varioskan LUX (Thermo Scientific).

Machine learning analyses

The dataset was composed of 2,218 simulations from 405 combinations of parameters. The dataset was split into train (80%) and test (20%) sets. One tree-based model (RandomForest Classifier) and one linear-based model with regularization (ElacticNet) were used to predict the final result of simulations (regenerative, unstructured, and scar repair). To better explain the prediction results, a SHAP model was applied to calculate the influences and interactions of each parameter on the output of the predictive model. RF analysis was performed using Scikit-learn, seaborn and shap libraries, and the following RandomForestClassifier hyperparameters (max_depth = 6, min_samples_split = 10, n_estimators = 15, criterion=’gini’).

Statistical analyses

The number of animals used in each study is indicated directly on the figure or in the figure legends. Measurements were taken from distinct samples. Mice were randomly allocated to the different groups and investigators were blinded to analyses. All results are given as means ± SEM or ± SD for the barplot or curve plot and median ± min/max for the boxplot. Statistical differences were measured using an ANOVA test when there were more than two groups and Tukey-HSD post-hoc or Kruskal-Wallis tests (two-sided) were used to determine statistical differences between each group. All statistical analyses were carried out using RStudio software. p < 0.05 was considered as significance level. The following symbols for statistical significance were used throughout the manuscript: *p < 0.05; **p < 0.01; ***p < 0.001, ****p < 0.0001.

Mathematical modeling

The 2D mathematical model features cells described as 2D spheres represented by their centers and radii and fibers described as segments of fixed length. The model was implemented on a 2D square ___domain with periodic boundary conditions. Cells (adipocytes) and fibers were represented as 2D-spheres of time-dependent radius and segments respectively. The agents (cells and fibers) interacted via mechanical repulsive interactions. Cells were modeled as incompressible, i.e., non-overlapping, while fiber-fiber and fiber-cell interactions were modeled via soft repulsion potentials, allowing some interpenetration of the agents. Intersecting fibers were able to link or unlink following random (Poisson) processes of frequencies νlink and νunlink, respectively. In contrast, linked fiber pairs were constrained to maintain the position of their cross-link during motion. Linked fibers were subjected to an alignment force at their junctions. The model was divided into two steps (morphogenesis and reconstruction). The morphogenesis step was started from a random distribution of fibers linked with Pf probability. New cells were randomly inseminated into the ___domain with a minimal radius and allowed to grow with linear volumic growth until a maximal radius was reached. The cell insemination process was stopped when the maximal number of cells was reached. As shown in [9], this model was able to spontaneously generate tissue architectures at equilibrium consisting of lobular cell structures in organized fiber networks which served as a basis for the reconstruction step.

The reconstruction steps were studied starting from morphogenesis simulations, described above, as a starting point. Tissue injury was modeled by removing all the components (cells and fibers) located in a region of the simulation square ___domain. This created a gradient at the border of the wounds that activates the production of an ‘injury signaling chemical’. This chemical was produced at the front of the tissue (borders of the wound), allowed to diffuse in the tissue, and had a finite lifetime. This signal inhibited the insemination of new cells and locally activated the insemination of new fibers modeled as a Poisson process of frequency νf. Upon insemination, new fibers were automatically linked with a proportion (Pf) of their intersecting neighbors, and unlinked with frequency νad. The wound was then filled with new fibers and the injury-signaling chemicals disappeared, enabling insemination of new cells modeled as a Poisson process. The frequency of new cell insemination depended on: (a) The density of injury-signaling molecules: new cells were more likely to appear where the density of injury-signaling chemical is low, (b) the number of fiber links in the ECM: new cells were more likely to appear in regions with low ECM fiber links, and (c) the density of existing agents: new cells were more likely to appear where existing fibers were already present. As in the morphogenesis step, new cells grew according to a linear volumic growth up to a preset maximal radius.