To investigate the inflammatory stage of bone fracture healing, a multiscale in silico model has been developed. The model, named COMMBINI (COmputational Model of Macrophages dynamics in the Bone INjury innate Immunoresponse), aims to simulate the biological and mechanical environment during the progression of the healing of a long bone fracture. To date, only the cellular and molecular modules of COMMBINI have been developed with the support of PhySiCell [Ghaffarizadeh et al. (39)], an open-source software that simulates the cells as single entities within an agent-based model. These virtual cells perform phenotype-specific activities (e.g. migration, proliferation) and regulate the molecular level (e.g. consumption and production of cytokines).

During the inflammatory stage, the cellular level plays a major role, as the innate immune cells actively contribute to initiating the healing response. To simulate this cellular level, an agent-based model has been developed. With this approach, each cell was simulated independently and not as a passive component of a cell population, thereby providing stochasticity to the investigation and guaranteeing the spatio-temporal variability that characterizes biological systems [Wehrens et al. (40); Allen et al. (41)]. The simulation was performed within a geometrical domain that represents the shape of a murine tibia fracture over a virtual period of 3 days, encompassing the early inflammatory stage. For the current study, a fracture opening in the center of the bone was simulated. The size of the fracture gap depends on the specific case study under investigation (cfr. § 2.3, 2.7). The model geometry was created by assuming a hollow cylinder as a simplified shape for the bone and a spheroid shape for the callus domain (healing region), following the same assumptions as previous studies [Wang and Yang (42); Borgiani et al. (43); Perier-Metz et al. (44)]. The healing region is the spatial domain where cell activities and molecular dynamics are simulated. Boundary conditions are imposed on the surfaces of the healing region. Bone marrow is simulated as a reservoir of non-polarized macrophages: they are recruited from the bone marrow compartment to invade the healing region. Furthermore, once the inflammation is over, the macrophages leave the region and emigrate back to the marrow compartment. The same conditions are imposed on the curved surface of the healing region, which simulates the periosteal boundaries. A zero-flux condition is imposed on the surface of the bone cortex as it is assumed that cells cannot migrate and cytokines cannot diffuse through it. The 2D model is generated by an intersecting plane along the middle axis ( Figure 1A ).

The iterative nature of the model allowed the investigation of the cellular environment evolution with a time resolution of Δt _cell = 1 min. In each iteration, virtual cells within the Region of Interest (ROI) perform specific actions based on phenotype-specific ratios, and the cellular environment is updated accordingly. Four different cell phenotypes are described in this computer model: non-polarized macrophages (M0), pro-inflammatory macrophages (M1), anti-inflammatory macrophages (M2) and polymorphonuclear neutrophils (PMN) ( Figure 1B ). The PMNs are the only cell type simulated within the healing region at the initial time-point. They are uniformly distributed within the region with an initial concentration [PMN]₀ . Macrophages start to appear from the first iterations of the simulation onward. At the molecular level, an initial concentration of fracture debris (Db₀ ) is homogeneously distributed within the healing region. This initial condition is crucial as the presence of debris chemotactically promotes the invasion of the healing region by the immune cells. No inflammatory cytokines are simulated within the region at the initial time-point but they start being secreted from the first iteration onwards. The M0 recruitment from the marrow cavity and tissues surrounding the healing region is stimulated by the presence of debris. The PMNs and macrophages phagocytose the debris, leading to a decrease in its concentration and recruitment capacity as healing progresses. To simulate this behavior in the computational model, the M0 recruitment ratio follows a dynamic pattern that decreases along with the physiological reduction of debris concentration within the healing region [Trejo et al. (38)]:

Equation (1) halts macrophage recruitment when no more debris needs to be removed or if the maximal concentration of macrophages is reached within the healing region. The dynamic is regulated by two parameters, whose variation can lead to faster or slower recruitment of macrophages: k_R(M0) is the maximum non-polarized macrophage recruitment ratio and [MΦ] _max is the maximal macrophage concentration allowed within the healing region. To create a realistic evolution of the cellular environment during the inflammatory stage of bone healing, the macrophages and PMNs perform additional activities, i.e. they migrate, proliferate, polarize and are subject to apoptosis ( Figure 1C ).

In COMMBINI, cellular migration is stochastically simulated as a sequence of “jumps” in random directions to create a movement pathway in the 2D space [Allen et al. (41)]. The span of each jump is defined by the migration speed (k_v ) associated with each cell phenotype. Cellular proliferation is simulated by generating a daughter cell with identical characteristics to its mother cell in one of the neighboring positions. The proliferation ratio (k_p ) associated with each cell phenotype determines the frequency of cell division within each iteration. Apoptosis is simulated as the removal of cells by programmed cell death. The apoptosis ratio (k_a ) of a cell increases with the accumulation of phagocyted debris [Bratton and Henson (45)] and the number of other cells in its vicinity, mimicking the consumption of essential nutrients for survival. Furthermore, macrophages have the ability to change their phenotype in response to the surrounding inflammatory environment as perceived at the molecular environment (more details in § 2.2). The M0 macrophages can, under specific molecular conditions, polarize into either an M1 or M2 phenotype ( Figure 1B )[Yunna et al. (46)]. In our model, this process is simulated as the change of the phenotype flag associated with the macrophage. Following the phenotype switch, the virtual macrophage adjusts its behavior by modifying its parameter values and algorithm dynamics according to the characteristics assigned to the new phenotype. Moreover, although infrequent, interpolarization can occur between M1 and M2 phenotypes, depending on whether pro-inflammatory macrophages reside in an anti-inflammatory environment, and vice versa ( Figure 1B )[Yunna et al. (46)]. Inter-polarization into pro-inflammatory macrophages is rare compared to the interpolarization into the anti-inflammatory phenotype due to the natural progression of the bone healing process. To avoid unnecessary complexity, this model does not include further subdivisions within the M2 subtypes. However, the model can be readily expanded to include such subdivisions if the scope is extended beyond the inflammatory stage to incorporate the repair phase.

In the discrete agent-based model, the apoptosis, proliferation and polarization conditions are reported as probability values for the respective event to occur within the iteration period Δt _cell . Therefore, during each iteration, a random floating point value between 0 and 1 (precision 10⁻⁶) is assigned to each cell for each event. If the value exceeds the probability value, the event does not get triggered (N paths in Figure 1C ). Conversely, if the value is lower than the probability value, the cell is removed, generates a daughter cell or changes its phenotype (Y paths in Figure 1C ). For cell proliferation, the position of the daughter cell is randomly selected from the four adjacent positions that are not occupied by other cells. Migration is performed at every iteration by allowing the cell to jump multiple times to adjacent positions based on their migration speed and the spatial and temporal resolution of the model. In this model, assuming a spatial resolution of 1 µm (cellular model spatial resolution) and an iteration period of Δt _cell = 1 min (temporal resolution), a PMN (k_v = 5.00 µm min⁻¹) will perform five jumps during each iteration. The direction of each jump is randomly chosen among the four surrounding positions that are not occupied by other cells, when chemotaxis is not involved. However, a large part of the phagocytic cells included in this work is driven by the fracture debris gradient. Chemotaxis is incorporated into the model by directing cell movement according to the gradient of the chemotactic agent concentration ( Figure 1C ).

While macrophages are recruited, PMNs promote the onset of the inflammatory response. In the first version of COMMBINI, PMNs are the only non-macrophage population considered at the cellular level. At the start of the simulation, PMNs are uniformly distributed within the healing region with an initial concentration [PMN]₀. Through the course of the inflammation, PMNs are recruited from the surrounding tissues by following a dynamic analogous to (1). PMNs are short-lived cells that tend to disappear from the healing region after triggering the initial inflammatory signal and its amplification [Summers et al. (47)]. Therefore, the proliferation of PMNs is not included in the model ( Figure 1B ). To simulate the natural behavior of neutrophils, PMNs simulated in COMMBINI release pro-inflammatory cytokines and clear debris from their surroundings to generate a pro-inflammatory environment [Kovtun et al. (48, 49)].

The cellular level has a mutual regulatory relationship with the molecular level. Consequently, we simulated the molecular model within the same agent-based model that simulates the cellular environment. The dynamics of cytokine concentration at the molecular level are simulated using partial differential equations (PDE) with function descriptions obtained from the literature ( Supplementary Table 2 ). The equations were solved using the BioFVM solver [Ghaffarizadeh et al. (50)] on a 2000 µm x 2000 µm square 2D grid within the healing region, with a resolution of 10 µm. The concentration of each inflammatory cytokine is evaluated in each grid element. This setup enables multiscale interactions, as each element in the molecular model shares its position with one or more cells in the cellular environment, according to the common coordinate system. The activities of the cells within the same element are regulated by the cytokine concentration within it ( Figure 1C ). Conversely, the presence of cells within each element regulates the intrinsic variation of cytokine concentration, reproducing phenotype-specific dynamics ( Figure 1C ). Macrophage polarization is regulated by the molecular level as the macrophages simulated at the cellular level polarize according to the cytokine concentration predicted in the same spatial location of the healing region ( Figure 1C ). Tumor Necrosis Factor alpha (TNFα) and Interleukin 10 (IL10) have been chosen for this model to respectively represent pro- and anti-inflammatory cytokines at the molecular level. Therefore, we described the macrophage polarization rules as probability functions, which are regulated by the concentration of those cytokines [Trejo et al. (38)]:

In equations (2 - 5), the parameters k represent the macrophage polarization ratios and the parameters a represent the cytokine half-saturation for macrophage polarization.

The molecular environment is, in turn, regulated by the immune cells ( Figure 1B ). These release pro- and anti-inflammatory cytokines, according to the dynamics included in the model. In addition to TNFα and IL10, the model includes Transforming Growth Factor beta (TGFβ) and Interferon gamma (IFNγ), as they regulate cell activity in the healing region: e.g. TGFβ lowers secretion of pro-inflammatory cytokines by M1 macrophages [Nagaraja et al. (51)], and IFNγ downregulates macrophage proliferation ( Figure 1B ). All the cell-specific cytokine secretion dynamics simulated in this model are reported in Table 1 .

The cytokines diffuse through the molecular level by following Fick’s law of diffusion, with a specific diffusion coefficient (D) associated with each cytokine. Neumann boundary conditions (no-flux) have been assigned to the borders of the healing region and bone cortex. Additionally, decay rates (d) have been set for each cytokine to simulate their enzymatic degradation, leading to a decrease in concentration. A temporal resolution of Δt _mol = 1 s was assigned to iteratively simulate the dynamics within the molecular level. To coordinate the temporal dynamics between the two levels, which are characterized by different temporal resolutions, the cellular environment updates every 60 iterations of the molecular level.

Additionally, the molecular level simulates the dynamical spatio-temporal variation of the concentrations of debris within the healing region. In this study, the term debris is used to define the agglomerate of dead cell bodies and necrotic tissue pieces resulting from the bone fracture. The presence of debris elicits the release of Damage Associated Molecule Pattern (DAMP) inflammatory stimuli. The distribution of debris concentration is included at the molecular level as a biological variable capable of influencing the inflammatory stage development [Chow et al. (52)]. In COMMBINI, the macrophages follow the debris concentration gradient at the molecular level to orient their migration towards the zones of the healing region characterized by a higher concentration of debris. Phagocytosis has been implemented in the model as the capacity of macrophages and PMNs to remove debris in their spatial surroundings, hence clearing the healing region. An engulfment ratio k_e was defined to quantify the debris phagocyted by those cells within the iteration period.

The model parameters at both cellular and molecular levels were obtained from previously published in vitro works that investigated the biological characteristics of macrophages and cytokines ( Supplementary Tables 1, 2 ). Afterward, a parameter calibration was performed to minimize the differences between the simulation outcomes and the experimental results from dedicated in vivo studies through the use of immunofluorescent imaging of macrophage populations. The in vivo experiments have received approval from the Ethical Committee for Animal Experimentation of the KU Leuven (approval number 020/2022). Tibial osteotomies (1 mm) were created in male C57BL/6 mice, fixated with an external Ilizarov fixator as previously described [van Gastel et al. (53)]. Three samples were obtained from the animals at 3 days post-fracture and prepared for immunohistology. The samples were fixated in formalin overnight at 4°C and decalcified with an edetic acid (EDTA) solution. The decalcified fracture samples were embedded in paraffin and 5 µm thick sections were mounted on glass slides. One slide from the center of each sample has been selected for immunofluorescence staining, obtaining n = 3 ex vivo images to use for calibration. The slides were deparaffinized with Histo-Clear (National Diagnostics, cat. no. HS-202) and dehydrated, followed by enzymatic antigen retrieval using 1 mg mL⁻¹ Pepsin in 0.02M HCl. The samples were blocked with 5% bovine serum albumin (BSA) in phosphate-buffered saline (PBS) with 0.1% Tween20 (Merck, cat. no. P1379) and 0.01% Tergitol (AppliChem, cat. no. A9780) for 45 minutes at room temperature. The samples were stained with immunofluorescent markers for macrophages and their specific subtypes. DAPI identifies all the nuclei and the Cluster of Differentiation 68 (CD68) is a general marker for macrophages [Schlundt et al. (26)]. Co-expression of CD68 and CD80 is specific for pro-inflammatory macrophages, while co-expression of CD68 and CD206 identifies anti-inflammatory macrophages [Schlundt et al. (26)]. The samples have been incubated overnight at 4°C with a 1:500 dilution of anti-CD68 antibody (ThermoFisher, cat. no. 14-0681-82) and a 1:100 dilution of anti-CD80 (ThermoFisher, cat. no. PA5-85913) or anti-CD206 antibody (ThermoFisher, cat. no. PA5-101657) in blocking buffer. On the second day, the samples were incubated for 4 hours at room temperature with Goat anti-rat IgG, Alexa Fluor Green 488 antibody with a 1:500 dilution (ThermoFisher, cat. no. A-11006) and Donkey anti-Rabbit IgG, Alexa Fluor Red 594 with a 1:1000 dilution (ThermoFisher, cat. no. A32754) in blocking buffer. Since bone is autofluorescent, the Vector TrueVIEW autofluorescence quenching kit (Vector, cat. no. SP-8400-15) was used. Finally, a counterstaining was performed with 5 µg mL⁻¹ DAPI for 10 minutes. The samples were dried and mounted in VECTASHIELD Vibrance Antifade (Vector, cat. no. H-1700). Samples were imaged using the Olympus IX83 inverted microscope within 48 hours. The sections were conserved at -20°C for additional image acquisitions.

A custom Python script was developed to analyze the immunofluorescent images and extract quantitative information at the cellular level. The outcome of the pipeline generates a fully segmented image with spatial information about macrophage distribution. Whole-cell segmentation was performed by Mesmer (DeepCell), a deep-learning tool trained on an extensive database of tissue image data and validated by experts [Greenwald et al. (54)]. Dimension filtering is applied to the images and the elements with a surface area below 80 µm² or larger than 200 µm² are not classified as cells [Cannon and Swanson (55)]. An ROI is chosen on the immunofluorescent image by selecting the fracture region, avoiding the bone cortex and staining artifacts. All cells within the ROI are labeled according to phenotype and quantified. For each macrophage phenotype, concentrations are calculated by dividing the number of cells by the ROI area. This data is compared with the macrophage concentrations simulated by the cellular level in the agent-based model. To perform a more direct qualitative comparison, in silico immunofluorescence was generated as output of the computational model by assigning the same color-coded pattern to the virtual cells as the immunofluorescent images. For example, the bright green fluorescence assigned to the CD68 channel was used to paint the cytoplasm of all the virtual macrophages, as they are supposed to express that marker ( Figure 2 ). Co-marking is represented by the chromatic combination of the two markers: i.e. M1 cells that are co-marked by CD68 (green) and CD80 (red) are represented in silico with a yellow color ( Figure 2 ). Additionally, this novel computational technique delivers a dataframe that contains information about all the cells identified within the ROI of the immunofluorescence image. Each cell is categorized in detail according to its size, the 2D position of its centroid and marker positivity.

The parameter calibration of the computer model was performed by following an optimization pathway to reduce the difference between quantified experimental and simulation outcomes. The calibration process can be time-exhaustive when many model parameters are included. Therefore, a sensitivity analysis was performed to determine the model parameters that most strongly influence the quantitative outcome of COMMBINI. The model was run multiple times with different combinations of parameter values. Reduction of the number of simulation runs was possible by cutting non-necessary repetitions with the support of Taguchi’s orthogonal arrays [Kacker et al. (56)]. This strategy is convenient when many parameters have to be analyzed: the model is regulated by 36 parameters and a 2-level sensitivity analysis would have required 2³⁶ simulation repetitions to analyze all the parameter combinations (full factorial). With Taguchi’s orthogonal array, we reduced this number to 72, drastically dropping the estimated runtime of the analysis. An analysis of variance (ANOVA) was performed on the model outputs to evaluate the percentage of the total sum of squares (%TSS) for each parameter [Isaksson et al. (57)]. The absolute value of this percentage represents how sensitive the output is to variation of the parameter value: higher %TSS expresses a more significant influence on the output. The sign associated with the %TSS indicates the influence on the output variation: if positive, an increase in the parameter value results in an increase in the output value and vice versa. For each output, the four most influential parameters were selected according to the highest %TSS absolute value.

Once the most significant parameters were identified by the sensitivity analysis, we calibrated them with the support of a Genetic Algorithm (GA) [McCall (58)]. A fitness function was generated employing data from experimental images, with the aim of reducing the quantitative differences between the in silico model and ex vivo immunofluorescence images. Numerical differences between experimental data (e.g. concentration of macrophages) and the corresponding quantitative output from the agent-based model of the cellular level were employed as the fitness function. The GA follows an evolutionary approach based on subsequent generations, aiming to minimize the fitness function. If a combination of parameter values did not reduce the function, it was removed by the algorithm in the following generation, allowing it to keep only the most promising ones. The selection of the most promising values and their cross-combination with the other components of the population minimized, generation after generation, the fitness function until a predetermined threshold was met. A more detailed explanation of the GA methodology employed to calibrate this model is reported in Supplementary Materials .

Validation of the results was performed on a different dataset of experimental immunofluorescent images (n = 2), previously reported by Schlundt et al. (26). Differently from the dataset that was used for calibration, the model of the validation set is characterized by a smaller fracture gap size (0.7 mm), in a different bone (femur) from female mice. The mouse strain (C57BL/6) was analogous to our in-house experiment and the same immunofluorescent staining markers were used to investigate the macrophage distribution in ex vivo images. The model domain was adapted to match the dimensions of the validation experiment’s bone and fracture gap. The biological parameter values obtained from the GA calibration process were validated by quantitatively comparing the macrophage populations concentrations simulated on this new domain and the ones measured from ex vivo immunofluorescent images. The success of the validation process supports the claim that confirms the assertion that the additional calibration step using data from in vivo experiments is important and leads to a more accurate representation of the inflammatory phase of fracture healing in murine long bones than when using parameter values derived from in vitro experiments reported in the literature.

Due to the involvement of the discrete agent-based framework, the multiscale model has a stochastic nature. The variability is shown by the mean and standard deviation of multiple repetitions (n = 5) of the simulation under the same investigative conditions and different initial random seeds. One-tailed student’s T-test was performed to investigate the differences between the calibrated and non-calibrated models.

When in vitro experiments reported in the literature are used to parametrize the model, the simulation results show a concentration of macrophages within the healing region of 346.4 ± 9.3 mm⁻² after 1 day, followed by an average increase of 12.7% between day 1 and day 3. Specifically, at day 1 the M0 concentration is 207.5 ± 8.1 mm⁻², the M1 concentration is 99.2 ± 7.4 mm⁻² and the M2 concentration is 39.7 ± 7.1 mm⁻². As the inflammation progresses, the concentrations vary between day 1 and day 3: M0 decreases by 84.4 ± 4.1%, M1 and M2 increase 2.2-fold (± 0.3) and 3.2-fold (± 0.8) respectively ( Figure 3A ). At the molecular level, the cellular engulfment leads to a reduction in fracture debris over time, resulting in the complete clearance of debris from the healing region within 3 days ( Figure 3B ). Pro- and anti-inflammatory cytokines secreted by immune cells exhibit analogous dynamics throughout the onset of bone healing, though pro-inflammatory cytokine secretion is more intense during the early stage of healing ( Figure 3C ), followed by a delayed anti-inflammatory wave ( Figures 3C, D ).

When considering the total macrophage concentration output, the ANOVA test revealed that the in silico model exhibited the highest sensitivity to the macrophage recruitment ratio (k_{R ( M 0)}) during the early stage of inflammation (%TSS = 46.9% at day 1, reduced to %TSS = 3.3% at day 3). In the later stage, it was observed that the initial concentration of PMNs ([PMN]₀) had the largest impact, although with a negative trend (%TSS = -37.5% at day 3). Additionally, the non-polarized macrophage proliferation ratio (k_{p ( M 0)}) influenced the results at day 1 (%TSS = 12.9%), while the pro-inflammatory macrophage proliferation ratio (k_{p ( M 1)}) had a greater effect on the output at day 3 (%TSS = 15.9%). Furthermore, the debris engulfment ratio associated with PMNs (k_{e (PMN)}) exhibited an influence on the predicted macrophage concentration at day 3, with a negative trend (%TSS = -13.4%). The complete list of %TSS associated with each parameter at day 1 and 3 is reported in Supplementary Table 3 .

By minimizing the fitness function, defined as the difference in the macrophage concentration within the healing region between values obtained from computer simulations and experiments on day 3 postfracture, the GA identified the optimal combination of values for the most influential parameters at that time-point ([PMN]₀, k_{R ( M 0)}, k_{p ( M 1)}, k_{e (PMN)}). The algorithm converged after nine generations for the parameters ( Figure 4 ), and it resulted in a clear tendency for higher macrophage proliferation rates (k_{p ( M 1)}) to better capture the experimental data (1.07 10⁻³ min⁻¹, +28.5% compared to literature value). Calibrated values for macrophage recruitment and neutrophil engulfment ratios showed smaller yet still considerable divergence from literature-based values (k_{R ( M 0)} = 2.33 10⁻² h⁻¹, +10.9%; k_{e (PMN)} = 2.71 10⁻³ min⁻¹, -18.6%) and the initial PMN population tended to maintain the concentration value found in the literature ([PMN]₀ = 984.38 µm⁻³, -1.6%). Throughout the iteration of the GA, the average difference between in silico output and ex vivo immunofluorescent image quantification decreased from 240.9 mm⁻² to 107.1 mm⁻², resulting in a 56.5% reduction of the fitness function ( Figure 4 ).

When the model was run with the optimized parameters, the M0 concentration peaked around day 1 (213.1 ± 17.4 mm⁻²) and decreased with the progression of the inflammation (35.7 ± 6.5 mm⁻² on day 3). Pro- and anti-inflammatory macrophage concentrations increased from day 1 (M1: 122.1 ± 14.6 mm⁻², M2: 43.2 ± 8.3 mm⁻²) to day 3 (M1: 281.9 ± 26.6 mm⁻², M2: 140.0 ± 25.4 mm⁻²) ( Figure 5A ). The M1 concentration showed a significant influence of the calibration on day 1 (p = 0.016) and day 3 (p = 0.008), in contrast to both M0 and M2 concentrations where no significant influence was observed (p > 0.05).

When comparing the results of the model using literature-based values ( Figure 3A ) with those obtained with the calibrated model, we observed that the qualitative dynamics of macrophage concentration during the inflammation processes remained unaltered for all the subtypes. However, there was an increase in the number of cells within the healing region. The pro-inflammatory macrophage concentration in particular increased (+31.0%) due to the GA-driven increment of the proliferation ratio. This observation aligns with the fitness objective of the calibration to reduce the difference in macrophage concentration on day 3 between ex vivo immunofluorescence (518.8 ± 8.3 mm⁻², identified as CD68+ cells) and computer model results (non-calibrated: 389.3 ± 36.5 mm⁻²; GA optimized: 457.6 ± 51.5 mm⁻², p = 0.033) ( Figure 5B ).

In the case of the 0.7 mm fracture dataset, we observed a decrease in macrophage concentration within the healing region. The ex vivo immunofluorescence data at day 3 showed a concentration of 414.0 ± 22.6 mm⁻² (CD68+). This reduction in macrophage concentration is also predicted by the GA-calibrated in silico model for the 0.7 mm gap. Specifically, the in silico model predicted a concentration of 415.2 ± 28.8 mm⁻² at day 3 post-fracture ( Figure 5B ).