Mathematical Modeling of Bone Growth and Regeneration After Fractures Using Differential Equations: Qualitative Analysis and Numerical Simulation

Abstract

Bone fracture healing is a coordinated, multistage process governed by interactions among cell populations, biochemical signaling, and the local mechanical environment. This paper develops a four-dimensional nonlinear autonomous system of ordinary differential equations describing the early and intermediate phases of secondary fracture repair. The model tracks the densities of mesenchymal stem cells ( ), chondrocytes ( ), osteoblasts ( ), and a generic osteoinductive growth factor ( ), representing morphogens such as transforming growth factor-  (TGF- ) or bone morphogenetic proteins (BMPs). We model cell differentiation with saturating Holling type-II (Michaelis–Menten) functional responses rather than simple linear rates, reflecting the finite number of membrane receptors available for growth-factor binding. We show that all solutions with non-negative initial data remain non-negative and are uniformly ultimately bounded, so that the dynamics are confined to a compact invariant region . Three equilibria are identified: a trivial non-healing state , a fibrotic (cartilage/bone-free) state , and an interior healed state . Local stability is analyzed through the Jacobian matrix and the Routh–Hurwitz criteria, and a sufficient condition for global asymptotic stability of  is obtained using a Volterra-type Lyapunov function together with LaSalle’s invariance principle. A fourth-order Runge–Kutta simulation illustrates the expected qualitative sequence of events: stem-cell recruitment, transient soft-callus formation, and its gradual replacement by bone.