Mathematical modeling of fracture healing in mice: comparison between experimental data and numerical simulation results

Abstract

The combined use of experimental and mathematical models can lead to a better understanding of fracture healing. In this study, a mathematical model, which was originally established by Bailón-Plaza and van der Meulen (J Theor Biol 212:191–209, 2001), was applied to an experimental model of a semi-stabilized murine tibial fracture. The mathematical model was implemented in a custom finite volumes code, specialized in dealing with the model’s requirements of mass conservation and non-negativity of the variables. A qualitative agreement between the experimentally measured and numerically simulated evolution in the cartilage and bone content was observed. Additionally, an extensive parametric study was conducted to assess the influence of the model parameters on the simulation outcome. Finally, a case of pathological fracture healing and its treatment by administration of growth factors was modeled to demonstrate the potential therapeutic value of this mathematical model.

Appendix A

This appendix contains the non-dimensionalized equations and parameters of the mathematical model, the temporal and spatial parameters of the model, the scaling values for the seven dependent variables and the values of boundary and initial conditions. The seven dependent variables are the concentration of mesenchymal stem cells (c _m ), chondrocytes (c _c ) and osteoblasts (c _b ), the density of the fibrous/cartilaginous (m _c ) and bone (m _b ) ECM (total matrix density m=m _c +m _b ) and the concentration of a chondrogenic (g _c ) and an osteogenic (g _b ) growth factor. The meaning of the model parameters is explained in Bailón-Plaza and van der Meulen [4]. Most of the parameter values were derived from experimental results reported in literature but a few of them were estimated. Some of these estimated parameter values were altered (boldface) in this study, in order to obtain a better correspondence between simulation results and experimental data.

• Non-dimensionalized equations and parameters

  $$\frac{\partial c_m}{\partial t} = \nabla \left[D_{cm} \nabla c_m - Cc_m \nabla m \right] + A_{m} c_{m} \left[1 - \alpha_{m} c_{m} \right] - F_{1} c_{m} - F_{2} c_{m} $$

  (2)
  $$\frac{\partial c_c}{\partial t} = A_c c_c \left[1 - \alpha_c c_c\right] + F_2 c_m - F_3 c_c$$

  (3)
  $$\frac{\partial c_b}{\partial t} = A_b c_b \left[1 - \alpha_b c_b\right] + F_1 c_m + F_3 c_c - d_b c_b$$

  (4)
  $$\frac{\partial m_c}{\partial t} = P_{cs} \left(1 - \kappa_c m_c\right)\times\left(c_m + c_c \right) - Q_{cd2} m_c c_b$$

  (5)
  $$\frac{\partial m_b}{\partial t} = P_{bs} \left(1 - \kappa_b m_b\right)c_b $$

  (6)
  $$\frac{\partial g_c}{\partial t} = \nabla \left[D_{gc} \nabla g_c \right] + E_{gc} c_c - d_{gc} g_c $$

  (7)
  $$\frac{\partial g_b}{\partial t} = \nabla \left[D_{gb} \nabla g_b \right] + E_{gb} c_b - d_{gb} g_b $$

  (8)
  $$\begin{aligned} D_{cm} &= \frac{D_h}{\left(K_h^2 + m^2\right)}m \quad F_1 = \frac{Y_1}{\left(H_1 + g_b\right)} g_b \\ C &= \frac{C_k}{\left(K_k + m\right)^2} \quad F_2 = \frac{Y_2}{\left(H_2 + g_c\right)} g_c \\ A_m &= \frac{A_{mo}}{\left(K_m^2 + m^2\right)}m \quad F_3 = \frac{m_c^6}{\left(B_{ec}^6 + m_c^6\right)} \times \frac{Y_3}{\left(H_3 + g_b\right)} g_b \\ A_c &= \frac{A_{co}}{\left(K_c^2 + m^2\right)}m \quad E_{gc} = \frac{G_{gc} g_c}{\left(H_{gc} + g_{c}\right)} \times \frac{m}{\left(K_{gc}^3 + m^3\right)} \\ A_b &= \frac{A_{bo}}{\left(K_b^2 + m^2\right)}m \quad E_{gb} = \frac{G_{gb}g_b}{\left(H_{gb} + g_b\right)}\\ \end{aligned}$$

  D _h =0.014 [23], K _h =0.25 [23], C _k =0.0034 [23], K _k =0.5 [23], A _mo =1.01, K _m =0.1, α _m =1, A _co =0.101, K _c =0.1, α _c =1, A _bo =0.202 [17], K _b =0.1, α _b =1, Y ₁=10, H ₁=0.1, Y ₂=50, H ₂=0.1, Y ₃ =500, B _ec =1.5, H ₃=0.1, d _b =0.1, P _cs =0.2 [29], κ _c =1 [29], Q _cd2 =1.5, P _bs =2, κ _b =1, D _gc =0.005 [26, 38], D _gb =0.005 [26, 38], G _gc =50, H _gc =1, K _gc =0.1, G _gb =350, H _gb =1, d _gc =100 [9, 13, 15], d _gb =100 [9, 13, 15].

• Scaling parameters: T=1 day, L=3.5 mm, c ₀=10⁶ cells ml⁻¹, m ₀=0.1 g ml⁻¹ [29], g ₀=100 ng ml⁻¹ [26].

• Initial and boundary conditions: $m_{c\_ini} = 0.1,$ ${\bf c}_{\bf m\_bc} = {\bf 0.02}$ (during the first 3 days), $g_{c\_bc} = g_{b\_bc} = 2$ (during 20 days).