A 3D in Silico Multi-Tissue Evolution Model Highlights the Relevance of Local Strain Accumulation in Bone Fracture Remodeling

Since 5–10% of all bone fractures result in non-healing situations, a thorough understanding of the various bone fracture healing phases is necessary to propose adequate therapeutic strategies. In silico models have greatly contributed to the understanding of the influence of mechanics on tissue formation and resorption during the soft and hard callus phases. However, the late-stage remodeling phase has not been investigated from a mechanobiological viewpoint so far. Here, we propose an in silico multi-tissue evolution model based on mechanical strain accumulation to investigate the mechanobiological regulation of bone remodeling during the late phase of healing. Computer model predictions are compared to histological data of two different pre-clinical studies of bone healing. The model predicted the bone marrow cavity re-opening and the resorption of the external callus. Our results suggest that the local strain accumulation can explain the fracture remodeling process and that this mechanobiological response is conserved among different mammal species. Our study paves the way for further understanding of non-healing situations that could help adapting therapeutic strategies to foster bone healing.


1
Model parameter sensitivity analysis 1.1 Methods 56 different parameters were used in the multi-tissue evolution model. To identify their relative contribution to the simulation outcome, a parameter sensitivity analysis was performed with three levels for each parameter: -1 (divided by 2), 0 (baseline value) and +1 (multiplied by 2). This analysis was performed on a simplified geometry: a 3D beam with a square section of side 30mm and length 90mm subjected to a loading of 1 MPa applied at the free end of the beam, the other one being clamped (Supplementary Figure 1). Initial tissue fractions were defined homogeneously in the beam as follows: 40% mature bone and 10% of each other immature and mature tissue fractions. Tetrahedral elements of average size 3.5 mm were used. 3 different ROIs were defined: lower, middle and upper 10-mm horizontal slice in the beam (Supplementary Figure 1). Based on a Taguchi array, 243 simulation were run, from which the individual contribution of each parameter could be derived by calculating the sum of squared deviations (SSD) on the output value of interest: the average Young's modulus in the upper, middle and lower third ROI after 2500 and 5000h simulation time. Taking the Young's modulus into account allowed to integrate the various tissue types into one unique output value.

Results: Most influential parameters after 5000h
The 7 most influential parameters for the Young's modulus average values in the different ROIs after 5000h were: the bone yield strain, b B,III and c B,III (coefficients relating the bone mechano-response to the third principal strain), b B,I and c B,I (coefficients relating the bone mechano-response to the first principal strain), m B M (a bone maturation function coefficient), m B R (a bone resorption function coefficient). The bone yield strain was derived from an experimental study (Bayraktar et al., 2004) and was therefore assumed to be known.
To further investigate the effect of the other 6 parameters, 2 simulations were run for 10000h: a "lowlevel" simulation (1) where all 6 parameters had their value leading to reduced Young's moduli and a "high-level" simulation (3) where all 6 parameters had their value leading to increased Young's moduli (Supplementary Table 1). Those 2 simulations were compared to the baseline case (2) (1) to (3) (first row) and (2,4,5) (second row).

Results: Most influential parameters after 2500h
The 6 most influential parameters for the Young's modulus average values in the different ROIs after 2500h were: m B M , l B M , k B M (the 3 coefficients of the bone maturation function), k B R (a bone resorption function coefficient), the bone formation rate α B and the bone maturation rate γ B .
To further investigate the effect of those 6 parameters, 2 simulations were run for 10000h: a "lowlevel" simulation (4) where all 6 parameters had their value leading to reduced Young's moduli and a "high-level" simulation (5) where all 6 parameters had their value leading to increased Young's moduli (Supplementary Table 1). Those 2 simulations were compared to the baseline case (2) (all parameters taken as defined in the baseline computer model). In those cases, longer-term predictions did not differ much between the baseline case and the high-or low-level simulations (Supplementary Figure 2d-f); only the dynamics of the prediction was different (Supplementary Figure 3).

Conclusions
The parameter sensitivity analysis revealed that the simulations are mostly sensitive to a few bone mechano-response and maturation or resorption coefficients at the longer-term ( (Frost, 1996;Prendergast et al., 1997;Turner, 1998;Claes and Heigele, 1999;Martin and Seeman, 2008). The values for m B M and m B R were determined based on preliminary examples to achieve consistent results (Frame et al., 2017). The same values were therefore employed for the study described here.