Integrated computational and in vivo models reveal Key Insights into macrophage behavior during bone healing

Myeloid-derived monocyte and macrophages are key cells in the bone that contribute to remodeling and injury repair. However, their temporal polarization status and control of bone-resorbing osteoclasts and bone-forming osteoblasts responses is largely unknown. In this study, we focused on two aspects of monocyte/macrophage dynamics and polarization states over time: 1) the injury-triggered pro- and anti-inflammatory monocytes/macrophages temporal profiles, 2) the contributions of pro- versus anti-inflammatory monocytes/macrophages in coordinating healing response. Bone healing is a complex multicellular dynamic process. While traditional in vitro and in vivo experimentation may capture the behavior of select populations with high resolution, they cannot simultaneously track the behavior of multiple populations. To address this, we have used an integrated coupled ordinary differential equations (ODEs)-based framework describing multiple cellular species to in vivo bone injury data in order to identify and test various hypotheses regarding bone cell populations dynamics. Our approach allowed us to infer several biological insights including, but not limited to,: 1) anti-inflammatory macrophages are key for early osteoclast inhibition and pro-inflammatory macrophage suppression, 2) pro-inflammatory macrophages are involved in osteoclast bone resorptive activity, whereas osteoblasts promote osteoclast differentiation, 3) Pro-inflammatory monocytes/macrophages rise during two expansion waves, which can be explained by the anti-inflammatory macrophages-mediated inhibition phase between the two waves. In addition, we further tested the robustness of the mathematical model by comparing simulation results to an independent experimental dataset. Taken together, this novel comprehensive mathematical framework allowed us to identify biological mechanisms that best recapitulate bone injury data and that explain the coupled cellular population dynamics involved in the process. Furthermore, our hypothesis testing methodology could be used in other contexts to decipher mechanisms in complex multicellular processes.


Introduction
The tightly-coupled relationship between bone-forming osteoblasts and bone-resorbing osteoclasts in bone remodeling and healing is well established [1]. Bone remodeling is initiated by osteoclastic turnover of aged and compromised bone tissue. Molecular cues derived from bone resorption subsequently drive mesenchymal precursor expansion and differentiation into osteoblasts for formation of new bone [1]. Bone healing on the other hand begins with osteoblastic bone callus deposition that is subsequently remodeled by osteoclasts [1]. Beyond this classic paradigm of the bone modeling unit (BMU), studies are increasingly identifying other cellular populations and factors that also contribute to the maintenance of bone. Macrophages of the myeloid lineage play critical roles in inflammation, wound healing and cancer progression [2]. Recent studies have also shed light on their contribution to bone biology. While osteoclasts have traditionally been known as the tissue resident macrophage of the bone, more recent studies identified a novel population of bone-resident macrophages, osteomacs, which facilitate osteoblast bone formation [3,4]. Additionally, in the context of bone healing, macrophages have been documented to rapidly infiltrate sites of bone injury to clear cellular debris in a process called efferocytosis and elicit subsequent inflammatory response and mineralized callus formation [1]. Monocytes and macrophages are major components of the bone immune infiltrate following injury [1,[5][6][7]. Recent studies using genetic or pharmacological depletion of macrophages demonstrated significantly delayed time to bone repair [4,5,[8][9][10]. The diversity of macrophage function owes to its versatility in polarizing and responding to environmental cues [1,[6][7][8]. These critical functions ensure the right temporal sequence of events necessary for healthy and timely bone repair after injury. For instance, IL-4 and TNFα have been shown to promote different macrophage polarization states and impact bone healing [11][12]. As an example, acute pro-inflammatory factors such as TNFα can improve bone repair while prolonged administration has the opposite effect [11,12]. There are, however, a number of gaps in our understanding of monocyte and macrophage population and polarization behavior, including but not limited to: 1) the contributions of pro-versus anti-inflammatory macrophages in coordinating bone injury response, 2) whether macrophages are directly involved in control of osteoclasts and osteoblasts populations and activities during bone injury, and 3) the main mechanisms that govern pro-and anti-inflammatory macrophages population dynamics.
While in vitro and in vivo experimentation techniques can capture the behavior of individual populations with high resolution, they do not allow for understanding the simultaneous interplay between multiple cell types whose numbers change over time. This obstacle can be overcome with the integration of experimental data to computational approaches in order to model the interactions occurring during bone injury repair. This type of approach has already been applied to other disease contexts like cancer [13][14][15][16][17][18][19][20][21][22]. Amongst the possible types of modeling approaches, agent-based models, such as discrete-continuum Hybrid Cellular Automata can examine mechanisms at the cellular scale leading to emergence of nontrivial macroscopic patterns [23]. One advantage of such an approach, is the possibility to inform the model with experimentally measured parameters. However, these parameters, such as macrophage polarization rate for example, can sometimes be exceptionally difficult to measure in vivo or in vitro. On the other hand, systems of Ordinary Differential Equations (ODEs) model individual populations over time under a well-mixed assumption and are often used to estimate in vivo parameters. While they do not describe cellular mechanisms as finely as agent-based models, their relative computational simplicity make them a convenient tool to identify key parameters through data fitting [24][25][26]. Multiple mathematical approaches have been used to study bone remodeling and repair [27][28][29][30][31]. The vast majority use systems of ODEs to model bone cell populations in homeostatic bone remodeling and bone disease such as osteoporosis and multiple myeloma [32][33][34][35][36][37]. Bone remodeling is a physiological program that is tightly regulated spatially and temporally. Other groups have considered the role of space in the process and represented cell population either as continuous spatial field, describing the dynamics by a set of partial differential equations (PDE) [38,39], or as individual agents by an agent-based model approach [22]. These models have largely focused on the interaction between bone-building osteoblasts and bone-resorbing osteoclasts, mostly ignoring the role of immune and inflammatory cells. Although these models have addressed biologically and clinically relevant questions, very few studies, one of which is from our group [40] have quantitatively compared predictions of bone injury dynamics to longitudinal biological data. Some studies have included the role of inflammatory cells like macrophages, but they remain theoretical and have not been experimentally validated [41,42]. The role of inflammation, and that of macrophages in particular, is recognized as being key for coordinating the bone injury response in vivo but, to date, how monocyte/macrophage populations coordinate it and interact directly with osteoblasts and osteoclasts (and vice versa) during bone remodeling has not been thoroughly examined [1,6]. Here we use experimental, in combination with published data, to integrate osteoblasts, osteoclasts, bone, naïve, pro-and anti-inflammatory monocytes and macrophages into a coupled ODE model of the bone ecosystem. This approach allowed for the interrogation of key hypotheses that explain the bone healing program, such as the polarization and clearance dynamics of monocytes/macrophages, interactions between anti-inflammatory macrophages and pro-inflammatory monocytes/macrophages, and how pro-and anti-inflammatory monocytes/macrophages modulate osteoclast and osteoblast behaviors. We posit this integrated approach can be used to uncover mechanisms driving bone injury repair dynamics and to identify key strategies aimed at shortening bone healing times.

Quantitative data of cell populations during bone injury repair dynamics
During non-critical bone injury healing, the following sequence of steps occur: early inflammation and hematoma formation, direct intramembranous bone deposition into a mineralized callus by osteoblasts [43][44][45], and callus remodeling by osteoclasts [1,3,9,46]. Bone-forming osteoblasts and bone-resorbing osteoclasts are critical mediators of these steps, and their numbers shift accordingly during each phase of repair. In order to temporally quantitate bone cells during injury, we extracted multi-cellular longitudinal data from an experimental model of bone injury repair whereby non-critical epiphyseal fracture was generated in mice by direct intratibial injection [40,46-49] (Fig 1A). In the injured tibias of the mice, bone volume and cell populations numbers were quantified at baseline (day 0) and at day 1, 2, 3, 7 and 14 (n = 5/ time point) following injury. High-resolution μCT analysis of the site of bone injury demonstrated changes in bone volume (BV/TV) subsequent to bone injury ( Fig 1B). Bone volume remained diminished over a 48-hour period prior to a five-day long expansion, beyond  baseline levels. By day 14, the bone volume returned toward homeostasis. Consistent with other published observations, the overall bone volume dynamics were accompanied by corresponding sequential waves of osteoblast and osteoclast numbers [49-51] (Fig 1B). Interestingly, the overlaid data reveal alternating waves of osteoblasts and osteoclasts. In the same longitudinal study, the contralateral tibia from each mouse was additionally subject to flow cytometry to derive dynamics of total and polarized myeloid populations [2,4,5,7, (Fig  1B). The myeloid dataset shows that pro-inflammatory monocytes and macrophages spiked within the first 48 hours while anti-inflammatory macrophages were observed between 24 and 72 hours (Fig 1B). Importantly, a fainter but prolonged secondary wave of pro-inflammatory monocytes-macrophages was noted, an observation which is in line with past studies in other inflammatory contexts [51, [106][107][108].

Mathematical modeling reveals key insights into myeloid behaviors during bone injury repair
In order to shed light on monocyte-macrophage dynamics, we interrogated hypotheses regarding population dynamics, differentiation, lifespan and plasticity. To this end we built a coupled ordinary differential equations-based framework describing seven cell populations as well as the bone volume temporal dynamics. The cell populations we considered in the model were bone-building osteoblasts, bone-resorbing osteoclasts, naïve monocytes, pro-inflammatory monocytes, naïve macrophages, pro and anti-inflammatory macrophages. To properly integrate cell population temporal data into this framework, we first curated common literature observations and hypotheses regarding osteoblasts, osteoclasts, monocyte and macrophage behavior during tissue injury healing (Table 1). In the model, the initial dynamics (osteoblast expansion, osteoclast decrease, macrophage polarization, monocyte infiltration) are triggered by injury factors [106]. We assume that the amount of factors released from an injury are proportional to the bone damage induced and are the primary driver of myeloid response. Myeloid cells are known to infiltrate the bone and polarize into pro-inflammatory status to clear cellular debris when exposed to injury-associated factors (Table 1) [1]. Therefore, an injury variable was included in the model that drives the initial pro-inflammatory response by monocytes/macrophages. In the model, this injury variable is being depleted by a decay rate term that is proportional to the number of pro-inflammatory cells (Equations on S1 and S3 Figs). Table 1. Reference sources of predominant mechanisms. Established biological behaviors and functions of bone cell populations. Framework for a comprehensive and coupled 9-population ODE model is constructed based off of summarizing known published interactions between each population. Inclusion of select hypotheses for each ambiguous aspect of myeloid biology is based on the prevalence of their corresponding publications (at least seven supporting references for each).

PLOS COMPUTATIONAL BIOLOGY
Integrated computational framework deciphers macrophage biology during bone healing c. Injury factors drive resident pro-inflammatory macrophages polarization; pro-inflammatory macrophages repolarize into anti-inflammatory macrophages when phagocytosing cellular debris; Anti-inflammatory macrophages naturally repolarize into proinflammatory macrophages in absence of stimulus (plasticity); Pro-inflammatory macrophages drive pro-inflammatory monocyte polarization.
Whereas evidence for all these mechanisms have been shown in vitro or in vivo, our goal here is to identify the ones that can recapitulate in vivo dynamics, in order to define main mechanisms that drive bone injury dynamics. We focus here on parsimonious hypothesis combinations, where for a given cell dynamics aspect (a, b or c), one mechanism only is considered, e.g c1 as opposed to c1&2. We described the hypotheses using ordinary differential equations and integrated them into a mathematical model framework to assess the ability of each hypothesis combination to recapitulate the experimental data. This resulted in 18 ODE models, each describing a unique combination of hypotheses. In each permutation, we assessed how well the model fitted to the experimental data ( Table 2). The models were ranked based on their goodness of fit, which was measured by the Akaike Information Criterion (AIC, Table 2), and number of residuals lower than 1 ( Table 3). The fits to experimental data were obtained in two different ways, with the choice of two different functionals to minimize, J 2 and J 1 . Here we present the results obtained with J 1 but the conclusions remained the same with  Table 2. Akaike information criterion (AIC) for comprehensive ODE of all 18 combinations of hypotheses. Akaike information criterion (AIC) for comprehensive ODE of all 18 combinations of hypotheses. Left columns denote the hypotheses from each of three mechanisms tested. The AIC scores resulting from J 2 and J 1 minimization for each model are shown on the right and vary dramatically across models. Comparing AICs reveal one best combination (boxed in red) and the worst fitting model is highlighted in blue lines. the remaining combinations were substantially lower ( Table 2). Looking to another metric of goodness of fit, the number of residuals lower than 1, a3-b2-c2 is clearly the best combination, with 25 residuals lower than one, and the rest of the models is far apart. With 15, 13 and 14 residuals lower than one, respectively, a3-b1-c2, a3-b1-c1 and a3-b2-c1 rank pretty low regarding to the residuals metric (Table 3). Interestingly, some hypothesis combinations do better than a3-b1-c2, a3-b2-c1 and a3-b1-c1 in term of number of residuals lower than 1, but worst in term of AIC. In conclusion, combination a3-b2-c2 does substantially better than all other combinations for both goodness of fit criteria (and for both J 2 and J 1 optimizations). Of note, the best-fitting hypothesis a3-b2-c2 (Fig 3) assumes osteoblasts are the main osteoclastogenesis driver (Fig 2 Hypothesis a3), and that anti-inflammatory macrophages play an important role in suppressing osteoclasts and pro-inflammatory macrophages [91] (Fig 2 Hypothesis c2). It also suggests that initial osteoblast expansion is driven by factors associated with the onset of bone injury [106]. The coupled mathematical model also allows for the estimation of polarization rates over time for pro-and anti-inflammatory macrophages in these two scenarios (Table 4). By comparison, some hypotheses combinations such as a2-b2-c3 yielded significantly poorer fits (AIC of 113 for a2-b2-c3, the worst fitting one,

Model simulations are consistent with independent published experimental data
Analysis of the literature reveals several factors that are important regulators of bone injury repair, such as tumor necrosis factor alpha (TNFα), interleukin-4 (IL-4), interferon-γ (IFNγ) and oncostatin M (OSM) [69,98,[109][110][111][112][113][114][115][116]. For example, studies in mice genetically deficient for OSM exhibited reduced bone formation and osteoblasts numbers at the non-critical bone injury site [49,117]. OSM is produced by anti-inflammatory macrophages and promotes osteoblast expansion and activity [49, [117][118][119]. To assess the robustness of our bone injury repair mathematical model we simulated the effect of OSM depletion on osteoblast number and determined if the model would recapitulate the qualitative temporal dynamics of osteoblast, osteoclast population and bone volume as shown in an independent experimental dataset [49]. The best fitting hypothesis combination model integrates hypotheses a3, b2 and c2 (red boxes in a-c). a mechanism a3 assumes that osteoblasts and anti-inflammatory macrophages promote and inhibit osteoclast formation, respectively. b mechanism b2 assumes that injury factors promote osteoblast expansion. c mechanism c2 assumes that injury factors promote pro-inflammatory monocytes/macrophages polarization. Pro-inflammatory monocytes/macrophages promote anti-inflammatory macrophages polarization, which in return drive depolarization of monocytes/macrophages back to the naive state. d schematic representation of the model using a3-b2-c2 hypothesis combination. Arrows represent positive (green) or negative (red) types of cellular interactions. e Temporal plots and corresponding goodness of fit metrics (AIC and R2s) across all populations, obtained through J 1 minimization. https://doi.org/10.1371/journal.pcbi.1009839.g003 We found that reducing the effect of anti-inflammatory macrophages on osteoblast expansion by 50%, mineralization activity by 50% and osteoclast inhibition by 80% yielded similar osteoblast, osteoclast and bone dynamics to those obtained from the OSM-deficient mice ( Fig 5). Of note, osteoblast and bone levels are below baseline osteoclast number and remain largely unchanged between treatment and control in both the experimental data and model predictions. While not examined in vivo by the independent study, our mathematical model additionally generated corresponding predictions of the effect of OSM depletion on monocyte/ macrophage dynamics ( Fig 5). Interestingly, OSM depletion increased anti-inflammatory macrophage population and transiently decreased pro-inflammatory populations. Collectively, our model predictions are in qualitative accordance with this independent experimental dataset. This suggests that the model can be used for understanding the roles of myeloid cells in the bone ecosystem during bone injury healing and for developing therapies to accelerate and improve the process.

Discussion
Much remains to be discovered about how cells in the bone ecosystem collectively orchestrate the bone injury repair. Pharmacological and genetic experimental approaches can provide Table 4. Parameter values of the best fitting model a3b2c2. Table showing biological description of each mathematical variable with data-derived initial conditions and units. First column is the variable notation used in the equations for each parameter. Second column is the biological meaning of the variable. Third column is the initial condition for each variable, typically an initial cell population level. Fourth column is the variable unit. information as to the importance of key populations of cells, such as macrophages, but these approaches seldomly address the direct and indirect effects of other cell types involved in bone injury repair. Mathematical modeling has the advantage of being able to consider complex biological processes resulting from the interactions between several cellular populations, but their relevance is limited by the availability of biological parameters and validation data. Here, by combining experimental and mathematical models, we have investigated the interactions between cell populations in bone that synchronously orchestrate the bone injury repair program. To do so, we built a mathematical model that captures the dynamics of seven cell populations and the bone mass. Importantly, the system of equations is coupled so that each cell Other hypotheses combinations fail to recapitulate in vivo data. The example of the worst fitting hypothesis a2-b2-c1 fails to recapitulate in vivo data. a mechanism a2 assumes that pro-inflammatory monocytes/macrophages and osteoblasts promote and inhibit osteoclast formation, respectively. b mechanism b2 assumes that injury factors promote osteoblast expansion. c mechanism c1 assumes that injury factors promote pro-inflammatory monocytes/ macrophages and anti-inflammatory macrophages polarization. Anti-inflammatory macrophages drive depolarization of monocytes/macrophages back to the naive state. d schematic representation of the model using a2-b2-c1 hypothesis combination. Arrows represent positive (green) or negative (red) types of cellular interactions. e temporal plots and corresponding goodness of fit metrics (AIC and R2s) across all populations, obtained through J 1 minimization. type can regulate the activity of other cell types. This interplay between the different populations poses a challenge regarding the reconciliation between model dynamics and experimental data but gives credence to the novel insights it has allowed us to uncover. These include 1) Anti-inflammatory macrophages drive early osteoclast inhibition and pro-inflammatory phenotype suppression 2) pro-inflammatory macrophages are involved in osteoclast activation (bone resorptive activity), whereas osteoblastic cells promote osteoclast differentiation 3) Proinflammatory monocytes/macrophages rise during two expansion waves, which can be explained by the anti-inflammatory macrophages-mediated inhibition phase between the two waves.

Parameter
Experimentally, as described in our previous study [40], we observe a rapid expansion of pro-inflammatory monocytes and macrophages in the first 24 hours with anti-inflammatory macrophages emerging shortly thereafter and persisting for up 48 hours. Of the hypotheses tested by the model, a3-b2-c2 provided the best fit of model simulations to experimental data. Under this set of assumptions, anti-inflammatory macrophages cause retraction of the pro-inflammatory population and facilitates osteoblast expansion and mineralization/stabilization of the injury site. With the natural depletion of the anti-inflammatory population, the remaining injury-associated factors causes a second expansion of pro-inflammatory macrophages and monocytes that in turn enhance osteoclast formation and activity. This increased activity is essential for the resorption of the mineralized callus at the site of injury and the return to bone homeostasis in the given time frame. To our knowledge very few reports have proposed a role and mechanism for this second inflammatory wave during bone healing. A role of MSCs and osteoblastic cells has been proposed for a second increase in inflammatory cytokines like TNF [56]. However, this two waves pattern has been observed in larger spectrum of inflammatory contexts, not only in bone [51, [106][107][108], suggesting that this temporal profile is not bone specific. Importantly, to our knowledge the present study is the first to propose a mathematical framework of bone healing where bone and immune cell populations are fully coupled and to inform such a model with experimental longitudinal data of all these populations.
Our hypothesis combination approach allowed us to explore the polarization properties of monocytes/macrophages that can be difficult to determine in vivo. For example, the best fitting ODE model, a3-b2-c2 allowed us to estimate the rates of pro-and anti-inflammatory macrophage polarization and indicates that pro-inflammatory macrophages do not re-polarize into an anti-inflammatory phenotype given the time frame, which goes against studies suggesting macrophage plasticity and reprogramming, at least in the context of bone injury repair. While not disputing the possibility that macrophages can repolarize, our results suggest that, based on the timing of the acquired experimental timing points, repolarization does not appear to be the main mechanism that recapitulates macrophage polarization dynamics. Additional insights provided by the ODE model include estimations on macrophage lifespan during the healing process and the contributory roles of pro-inflammatory macrophages and monocytes to the process. This information can be critical for therapies that target specific myeloid populations during bone injury repair in a bid to accelerate bone healing.
Another important result of our study is that the hypothesis combination that fit the best, a3-b2-c2, implies that osteoblast and pro-inflammatory monocytes-macrophages have distinct roles in osteoclast biology. According to this model, osteoblast drives osteoclast differentiation, whereas pro-inflammatory monocytes-macrophages drive osteoclast resorptive activity. In most studies, this distinction is not made and both cell types are assumed to contribute to both osteoclast differentiation and resorptive activity. The data we present here suggest a distinction in osteoclast supportive functions. This hypothesis combination has been also used to simulate OSM depletion bone injury process. Model simulations were in accordance with experimental data from an independent published study.
An important aspect of mathematical modeling is that it allows us to distill the key cell species and molecules driving the bone dynamics. Of note, our model does not consider the potential roles of other cell types in the bone ecosystem that could contribute, such as T cells.
Our results suggest that that integrating myeloid populations into the model provides enough resolution to explain the process of non-critical bone injury repair. Our mathematical framework is flexible enough, however, that the effects of other immune cells such as T cells could be included. Additionally, we are aware the hypotheses we have identified throughout this study, while the most common, do not cover exhaustively all myeloid behaviors described in literature. However, the hypothesis testing pipeline we have devised enables us to efficiently adapt our model to reflect any additional hypotheses. Another important aspect that our model currently neglects is that of bone quality which requires a different set of data acquisition techniques and further refinement of the mathematical approach. Moreover, bone healing is a spatially regulated process and having this aspect included in the model would be an exciting refinement in order to explore further mechanistic aspects of bone structure and regeneration.
Through our hypothesis combination approach, we have integrated established biology into a mathematical framework describing cell population dynamics during non-critical bone injury repair. One potential application of our framework is to investigate how time to healing subsequent to bone injury can be reduced. Existing studies have shown that bone healing times can be impacted in modulating pro-and anti-inflammatory macrophages [10, 109,110,117]. While the model is parameterized with mouse data, there is much overlap between mice and humans with respect to the phases of the bone injury repair program. Thus, using our existing workflow, we can conceivably re-parameterize our model with human patient-derived data to further its potential as a relevant prospective tool for the clinic.
In conclusion, we have developed a coupled ordinary differential equation (ODE) system of the bone ecosystem that models the interplay between 8 key cellular populations during bone injury repair. The model yields several novel findings regarding macrophage dynamics and macrophage impact on osteoblasts and osteoclasts dynamics. Further, the model can also provide novel insights into phenomena that are hard to measure in vivo such as rate of pro-or anti-inflammatory polarization over time. A better understanding of bone healing will have clinical translatability allowing, for instance, to accelerate the process and improve patient outcomes. The model accounts for coupling between these population and will be useful in developing therapeutic strategies/interventions that shorten healing times. Further, the model has broad applicability and can be used as a platform to examine other bone diseases such as osteoarthritis and skeletal malignancies such as bone-metastatic cancer.

Ethics statement
Data from in vivo bone injury experiment were derived from our previous work [40]. In this study, all animal studies were performed in accordance with Guidelines for the Care and Use of Laboratory Animals published by the National Institutes of Health, and approved by the Animal Care and Use Committee at the University of South Florida, under IACUC Protocol R5857 (CCL). Male C57BL/6 mice (5-6 weeks old) were purchased from Jackson Laboratory. Mice (n = 30) were subject to tibial bone injury by penetration of a 28-gauge (0.3062mm diameter) syringe through the knee epiphysis to mid-shaft.

Intratibial bone injury model
Data from in vivo bone injury experiment were derived from our previous work [40]. In this study, all animal studies were performed in accordance with Guidelines for the Care and Use of Laboratory Animals published by the National Institutes of Health, and approved by the Animal Care and Use Committee at the University of South Florida, under IACUC Protocol R5857 (CCL). Male C57BL/6 mice (5-6 weeks old) were purchased from Jackson Laboratory. Mice (n = 30) were subject to tibial bone injury by penetration of a 28-gauge (0.3062mm diameter) syringe through the knee epiphysis to mid-shaft. Mice tibias at baseline and at days 1, 2, 3, 7 and 14 (n = 5/time point) were collected for analysis. Temporal population data was used to parameterize subsequent mathematical models.

Micro-computed topography
Bone volume data was derived from formalin-fixed tibias by micro-computed topography (μCT) scanning using Scanco μ35 scanner. Endosteal trabecular bone volume was analyzed 100μm away from the tip of growth plate to clear the dense bone nature of the growth plate. 1000μm along the midshaft of each bone was then scanned and analyzed using built-in functions (n = 30 bones; 5/time point).

TRAcP staining
After μCT analysis, tibia bones were decalcified using 14% EDTA for 3 weeks for further staining quantitation and analyses. Decalcified bones were sectioned at 4μm thickness. Sections were enzymatically stained for tartrate-resistant acid phosphatase (TRAcP) for osteoclast numbers based on manufacturer's protocol [120]. Stained slides were imaged using the Evos Auto microscope to capture 20X photos which included injury site and its immediate periphery. All TRAcP positive (red) multinucleated osteoclasts within 5μm radius from injury were counted, and mathematically converted to osteoclasts / bone marrow volume (#OCL/μm 3 ) for each bone at each time point.

Immunofluorescence staining and quantitation
FFPE tibia bones were further sequentially sectioned and baked at 56˚C in preparation for immunofluorescence staining of osteoblast (RUNX2 at 1:500; Abcam Cat. No. ab81357) and nuclear staining (DAPI). Deparaffined and rehydrated slides were subject to heat-induced antigen retrieval method. Sections were then blocked and incubated in primary antibodies diluted in 10% normal goat serum in TBS overnight at 4˚C. Subsequently, slides were stained with secondary Alexa Fluor 568-conjugated antibody at 1:1000 at room temperature for 1 hour under light-proof conditions. Stained slides were stained with DAPI for nuclear contrast and mounted for imaging at 20X using Zeiss upright fluorescent microscope to include the injury site as well as the immediate peripheral tissue. All runx2 positive cells (red staining colocalizing with DAPI) within 5μm radius from injury were counted and mathematically converted to osteoblasts / bone marrow volume (#OBL/μm 3 ) for each bone at each time point.

Flow cytometry and analysis
Harvested contralateral injured tibias (n = 30; 5/time point) had tips removed and were subjected to centrifugation at 16,000g for 5 seconds for isolation of whole bone marrow for flow cytometry staining and analysis. Red blood cells were lysed using RBC Lysis Buffer from Sigma Aldrich

Mathematical and computational methods
Comprehensive model structure. In order to efficiently describe the process of model building regarding the different hypothesis combinations (Fig 2), we begin with a generic set of coupled equations, using a formulation that is valid for all hypotheses combinations (S3 Fig). For all populations, homeostasis was described by a replenishment term and a clearance term. For polarized monocytes/ macrophages, no replenishment was considered at homeostasis, as the baseline measured by flow cytometry was close to zero. Here is the detailed description, equation by equation :  Equation 1: Naive monocytes.
Naive monocytes are assumed to be replenished at a constant rate H Mo , and to die at a rate δ Mo Mo. δ Mo , the lifespan parameter, was retrieved from literature (Table 4), and H Mo was estimated so that monocyte level at homeostasis would match experimentally measured monocyte baseline. The term I(M 1 , Mo 1 , D) corresponds to the number of monocytes infiltrating the bone marrow per unit of time due to injury factors and pro-inflammatory cells. It is equal to I 1 (M 1 + Mo 1 ) + I 2 D. As indicated by the mathematical formulation, inflammation-associated monocyte recruitment is driven by injury signals on one hand, and pro-inflammatory monocytes/macrophages on the other hand. This reflects the fact that cellular debris and pro-inflammatory cells produce factors that help recruiting additional monocytes [1].
Pro-inflammatory macrophages are assumed to be absent of the bone marrow under homeostatic conditions, as experimental al baseline level did not 0.3% across the replicates.
Anti-inflammatory macrophages are assumed to be absent of the bone marrow under homeostatic conditions, as experimental al baseline level did not exceed 0.07% across the replicates.
Pro-inflammatory monocytes are assumed to be absent of the bone marrow under homeostatic conditions, as experimental al baseline level did not exceed 0.16% across the replicates. The term p 3 (D, M 1 , Mo 1 )Mo represents the number of pro-inflammatory monocytes generated (from naive pool) per unit of time, as a function of injury factors and pro-inflammatory cells. The term depol 3 (M 2 )Mo 1 represents the number of pro-inflammatory monocytes that revert back to a naive state per unit of time, as a function of anti-inflammatory macrophages. The term depol 22 M 2 represents the number of anti-inflammatory macrophages that reprogram into a pro-inflammatory phenotype. The term p 22 (D)M 1 represents the number of pro-inflammatory macrophages that reprogram into an anti-inflammatory phenotype, as a function of injury factors.

Pro-inflammatory macrophages/monocytes and anti-inflammatory macrophages: Hypothesis c1.
Both CD11b+Ly6C+ monocytes and CD11b+Ly6C-macrophages can polarize into a proinflammatory phenotype [1], in response to injury signals or to factors produced by already present pro-inflammatory cells [1]. In assumptions c1 and c2, monocytes and macrophages polarize into pro-inflammatory monocytes and macrophages respectively through two terms. The first is proportional to the amount of injury signals present and the second is proportional to the amount of pro-inflammatory cells present. In the assumption c3, injury signals polarize local resident macrophages, and those in turn promote polarization of pro-inflammatory monocytes. In c3, pro-inflammatory macrophages repolarize into anti-inflammatory macrophages by the uptake of cellular debris/injury signals [121]. In this scenario, by plasticity, antiinflammatory macrophages naturally return to a pro-inflammatory phenotype in absence of signals [122]. Polarization into anti-inflammatory macrophages was assumed to be proportional to the amount of injury signals for assumption c1, proportional to the amount of proinflammatory cells for assumption c2, and a transition term from pro-inflammatory macrophages for c3.
Both CD11b+Ly6C+ monocytes and CD11b+Ly6C-macrophages can polarize into a proinflammatory phenotype [1], in response to injury signals or to factors produced by already present pro-inflammatory cells [1]. In assumptions c1 and c2, monocytes and macrophages polarize into pro-inflammatory monocytes and macrophages respectively through two terms. The first is proportional to the amount of injury signals present and the second is proportional to the amount of pro-inflammatory cells present. In the assumption c3, injury signals polarize local resident macrophages, and those in turn promote polarization of pro-inflammatory monocytes. In c3, pro-inflammatory macrophages repolarize into anti-inflammatory macrophages by the uptake of cellular debris/injury signals [121]. In this scenario, by plasticity, anti-inflammatory macrophages naturally return to a pro-inflammatory phenotype in absence of signals [122]. Polarization into anti-inflammatory macrophages was assumed to be proportional to the amount of injury signals for assumption c1, proportional to the amount of proinflammatory cells for assumption c2, and a transition term from pro-inflammatory macrophages for c3.
Both CD11b+Ly6C+ monocytes and CD11b+Ly6C-macrophages can polarize into a proinflammatory phenotype [1], in response to injury signals or to factors produced by already present pro-inflammatory cells [1]. In assumptions c1 and c2, monocytes and macrophages polarize into pro-inflammatory monocytes and macrophages respectively through two terms. The first is proportional to the amount of injury signals present and the second is proportional to the amount of pro-inflammatory cells present. In the assumption c3, injury signals polarize local resident macrophages, and those in turn promote polarization of pro-inflammatory monocytes. In c3, pro-inflammatory macrophages repolarize into anti-inflammatory macrophages by the uptake of cellular debris/injury signals [121]. In this scenario, by plasticity, antiinflammatory macrophages naturally return to a pro-inflammatory phenotype in absence of signals [122]. Polarization into anti-inflammatory macrophages was assumed to be proportional to the amount of injury signals for assumption c1, proportional to the amount of proinflammatory cells for assumption c2, and a transition term from pro-inflammatory macrophages for c3. Equation 6: Osteoblasts.
Osteoblast are assumed to be replenished at a rate H OB OC, proportional to osteoclasts. This reflects the ability of osteoclasts to produce osteogenic signals like transforming growth factor β (TGFβ) and bone morphogenetic proteins (BMPs) [123]. Similar assumptions are considered in published works of homeostatic bone remodeling [32][33][34]. Osteoblasts are assumed to die at a rate δ OB OBB. δ OB , the lifespan parameter, was retrieved from literature (Table 4), and H OB was estimated so that osteoblast level at homeostasis would match experimentally measured osteoblast baseline. The term γ OB (D, M 2 ) represents the number of osteoblasts generated per unit of time, as a function of injury factors and anti-inflammatory factors. Osteoblast dynamics: Hypothesis b1.
The osteoblast clearance term was assumed to be proportional to the bone volume, in order to account for osteoblast differentiation into osteocytes, when resorbing bone matrix [123]. A similar assumption is made in the model developed by Ryser et al. that describes bone remodeling as a spatial evolutionary game [124]. During injury, an extra term for osteoblast expansion is present, driven by anti-inflammatory macrophages (hypothesis b1) or injury factors (hypothesis b2), both supported by literature [1,3,6,7,125]. Osteoblast dynamics: Hypothesis b2.
The osteoblast clearance term was assumed to be proportional to the bone volume, in order to account for osteoblast differentiation into osteocytes, when resorbing bone matrix [123]. A similar assumption is made in the model developed by Ryser et al. that describes bone remodeling as a spatial evolutionary game [124]. During injury, an extra term for osteoblast expansion is present, driven by anti-inflammatory macrophages (hypothesis b1) or injury factors (hypothesis b2), both supported by literature [1,3,6,7,125].
Osteoclasts are assumed to be replenished at a rate d OC (OB, M 1 , Mo 1 , M 2 )M, which reflects differentiation of macrophages into osteoclasts, as a function of osteoblasts, pro-inflammatory macrophages, pro-inflammatory monocytes, anti-inflammatory macrophages [126][127][128]. This reflects osteoclastic factors produced by osteoblasts (RANKL) and pro-inflammatory monocytes/macrophages (IL-1, TNF), as well as anti-osteoclastic factors produced by osteoblasts (OPG) and anti-inflammatory macrophages (transforming growth factor β (TGFβ), IL-10). The term δ OC (M 2 , OB)OC represents the number of osteoclasts dying per unit of time, as a function of anti-inflammatory macrophages and osteoblasts. This reflects factors produced by anti-inflammatory macrophages (IL-10, TGFβ) and osteoblasts (OPG) that reduce osteoclast lifespan.
This osteoclast formation term was assumed to be proportional to osteoblasts for assumption a3, reflecting the ability of osteoblastic cells to produce RANKL, which is an essential mediator of osteoclast formation [123]. For the other assumptions, homeostatic osteoclast replenishment was assumed to be constant. This term had an additional contribution from pro-inflammatory monocytes/macrophages for assumptions a1 and a2, representing the ability of proinflammatory monocytes/macrophages to produce factors like IL-1 and TNF that favor osteoclast formation [123,129]. Osteoclast formation was divided by an inhibitory term, a linear function of anti-inflammatory macrophages for assumptions a1 and a3, and a linear function of osteoblasts for assumption a2. The first assumption reflects factors produced by anti-inflammatory macrophages, like IL-10, that disrupt osteoclast formation [86], whereas the second reflects the ability of osteoblasts to produce osteoprotegerin (OPG), a RANKL decoy receptor [123]. Moreover, this inhibition affects not only the ability of monocytes-macrophages to fuse and form osteoclasts, but also their life span. Indeed RANKL is necessary for osteoclast survival since OPG produced by osteoblasts reduces their life span [130]. Similarly, anti-inflammatory macrophages produce TGFβ, which is known to drive osteoclast apoptosis [131]. Osteoclast dynamics: Hypothesis a2.
This osteoclast formation term was assumed to be proportional to osteoblasts for assumption a3, reflecting the ability of osteoblastic cells to produce RANKL, which is an essential mediator of osteoclast formation [123]. For the other assumptions, homeostatic osteoclast replenishment was assumed to be constant. This term had an additional contribution from pro-inflammatory monocytes/macrophages for assumptions a1 and a2, representing the ability of proinflammatory monocytes/macrophages to produce factors like IL-1 and TNF that favor osteoclast formation [123,129]. Osteoclast formation was divided by an inhibitory term, a linear function of anti-inflammatory macrophages for assumptions a1 and a3, and a linear function of osteoblasts for assumption a2. The first assumption reflects factors produced by anti-inflammatory macrophages, like IL-10, that disrupt osteoclast formation [86], whereas the second reflects the ability of osteoblasts to produce osteoprotegerin (OPG), a RANKL decoy receptor [123]. Moreover, this inhibition affects not only the ability of monocytes-macrophages to fuse and form osteoclasts, but also their life span. Indeed RANKL is necessary for osteoclast survival since OPG produced by osteoblasts reduces their life span [130]. Similarly, anti-inflammatory macrophages produce TGFβ, which is known to drive osteoclast apoptosis [131]. Osteoclast dynamics: Hypothesis a3.
This osteoclast formation term was assumed to be proportional to osteoblasts for assumption a3, reflecting the ability of osteoblastic cells to produce RANKL, which is an essential mediator of osteoclast formation [123]. For the other assumptions, homeostatic osteoclast replenishment was assumed to be constant. This term had an additional contribution from pro-inflammatory monocytes/macrophages for assumptions a1 and a2, representing the ability of proinflammatory monocytes/macrophages to produce factors like IL-1 and TNF that favor osteoclast formation [123,129]. Osteoclast formation was divided by an inhibitory term, a linear function of anti-inflammatory macrophages for assumptions a1 and a3, and a linear function of osteoblasts for assumption a2. The first assumption reflects factors produced by anti-inflammatory macrophages, like IL-10, that disrupt osteoclast formation [86], whereas the second reflects the ability of osteoblasts to produce osteoprotegerin (OPG), a RANKL decoy receptor [123]. Moreover, this inhibition affects not only the ability of monocytes-macrophages to fuse and form osteoclasts, but also their life span. Indeed RANKL is necessary for osteoclast survival since OPG produced by osteoblasts reduces their life span [130]. Similarly, anti-inflammatory macrophages produce TGFβ, which is known to drive osteoclast apoptosis [131]. Equation 8: Bone volume.
Bone dynamics is described by two terms: a bone resorption term, which is the volume of bone resorbed per unit of time and is assumed to be proportional to the number of osteoclasts, and a bone formation term, which is the volume of bone formed per unit of time and is assumed to be proportional to the number of osteoblasts. Such assumptions have been broadly used across a large variety of modeling studies [32-34]. δ B (1 + α)OCB is the bone resorption term and is the sum of the two terms δ B OCB (homeostatic resorption) and α(M 1 + Mo 1 ) δ B OCB (pro-inflammatory monocyte/macrophage mediated resorption). As indicated by this mathematical formulation, bone resorption was assumed to also be proportional to the bone mass. This reflects the fact that more bone volume increases the likelihood for bone resorption. Furthermore, this formulation ensures bone mass stays strictly positive in the model. P B (1 + β)OB is the bone formation term and is the sum of the two terms P B OB (homeostatic bone formation) and βM 2 P B OB (anti-inflammatory macrophage mediated bone formation). Under homeostasis, bone is formed at rate P B OB and resorbed at rate δ B OCB. Resorption rate δ B , as well as as parameters α and β, were calibrated on bone volume dynamics, and bone formation rate P B was then imposed by the relation P B OB 0 = δ B OC 0 B 0 , which ensures that bone volume remains at homeostasis when osteoblasts and osteoclasts are at homeostatic levels. As indicated by mathematical formulations, bone formation and resorption terms were assumed to linearly increase with respect to anti-inflammatory and pro-inflammatory monocytes/ macrophages, respectively. This accounts for the fact that anti-inflammatory macrophages typically produce osteogenic factors like TGFβ or OSM, that are known to promote osteoblast expansion and bone mineralization [3], and for the fact that pro-inflammatory monocytes/macrophages typically produce osteolytic factors like TNF and IL-1, that are known to promote osteoclast resorptive activity [132]. Equation 9: Injury factors.
Injury factors dynamics consists in an exponential type of decay, with a decay rate δ D (M 1 + Mo 1 ) proportional to pro-inflammatory monocytes/macrophages number, which represents how pro-inflammatory monocytes/macrophages uptake cellular debris, which in return reduces pro-inflammatory signals. Population homeostasis. In order to estimate the homeostatic cell replenishment parameters, we set them equal to the clearance term (lifespan) which was either based on literature values or calibrated directly from experimental data (S1 Fig and Tables 1, 4 and 5).
ODE solver. The ODE45 function of Matlab was used to solve the differential equation system. The experimental baseline values (time 0) were used as initial conditions. Parameter estimation method. To estimate parameters for goodness of fit, we defined the following objective function: Where i represents the time point index and j the variable index, α represents the parameter set used to evaluate the model function f, D ij represents the experimental data of variable j at time point i, σ i represents the experimental standard deviation (over all the animals of a given time point), and N represents the number of time points. The functional J 2 corresponds to the weighted least squares criterion. The functional J 1 is the Tchebychev approximation, which considers the maximal residual instead of the sum of the residuals. In both cases, the choice of the max over the observed variables of the sum of the squares of the residuals was motivated to make sure that every variable was fitted with equal relative importance. Indeed, in the case of the minimization of the sum of the squares over all the variables, it is sometimes possible to find an optimum in minimizing certain variables at the detriment of others. This way, we ensure that all variables are equally well fitted. The reason for considering the objective function J 1 in addition to the classical criterion J 2 is to avoid neglecting any time point in the fit. The least squares functional allows sometimes to find an optimum optimizing certain time points at the detriment of others. In this current study, this is a potentially big issue, as the time sampling is not homogenous across the time points, meaning that biological dynamics of importance might be ignored, while still producing a good fit under the least squares metric.
In order to minimize this function representing the error estimate between data and model, we used the Matlab function fminsearch with a penalization term to stay in a parameter range set with reasonable boundaries.
In order to rank the models in term of goodness of fit, we used AIC, that is defined as follows: where p is the number of parameters, and the functional J is either J 2 or J 1 .
OSM knockout data. In order to retrieve data from the OSM knockout independent dataset, we used webplotdigitizer to collect datapoint from the plot presented in Guihard et al [49].  in dMo/dt). Parameters with no reference publication were estimated to obtain best possible fits to temporal dynamics data (parameters in red) and are listed in Table 3. In all equations, black terms correspond to homeostatic dynamics, whereas red terms correspond to injury dynamics.  inflammatory macrophages promote osteoblast expansion. c Mechanism c2 assumes that injury factors promote pro-inflammatory monocytes/macrophages polarization. Pro-inflammatory monocytes/macrophages promote anti-inflammatory macrophages polarization, which in return drive depolarization of pro-inflammatory monocytes/macrophages back to the naive state. d Schematic representation of the model using a3b1c2 hypothesis combination.

Supporting information
Arrows