© X. Zhang et al. Published by EDP Sciences 2025

1 Introduction

The healing of fracture is a complex biological process, which is usually classified into direct healing (primary healing) and indirect healing (secondary healing) in medicine [1]. Direct healing refers to the direct connection of the fracture ends through close contact and stable fixation (such as internal fixation with plates and screws) without the formation of obvious callus. Indirect healing refers to the completion of healing through the gradual calcification of fibrocartilaginous callus, which is the main way of natural healing. Secondary healing is typically divided into stages of inflammation, periosteal thickening, soft callus formation, and hard callus formation. Both mechanical stimulation and growth factors can affect fracture healing. Appropriate mechanical stimulation will promote fracture healing, while inappropriate mechanical stimulation can lead to delayed healing or non-healing [2]. Besides mechanical stimulation that affects fracture healing, biological factors also play an important role in the fracture healing process.

Various computer numerical models have been developed to predict tissue differentiation during secondary healing. These models include qualitative and quantitative analyses and can be roughly classified into mechanical regulation healing models, biological regulation healing models, and mechanical-biological regulation healing coupling models.

1.1 Mechanical-regulation healing models

Studies have shown that the local mechanical environment during fracture healing plays an important role in regulating tissue differentiation. Pauwels established the theoretical framework for mechanical stimulation and tissue differentiation for the first time through systematic experiments [3]. On this basis, the Claes team used the elastic finite element model to quantitatively analyze the correlation between strain and hydrostatic pressure and callus tissue differentiation, and proposed a tissue differentiation theory based on local stress and strain state [4,5]. Their study confirmed that the mechanical environmental parameters at the calcification interface are the key factors that determine tissue differentiation, and this conclusion was double-factor validation by finite element simulation and animal experiments.

The Prendergast research group innovatively constructed a biphasic finite element model, and revealed the correspondence between mechanical stimulation intensity and tissue differentiation type by introducing two different mechanical parameters: strain stimulation and solid-liquid relative velocity stimulation [6]. Their study found that high-intensity mechanical stimulation induces fibrous tissue formation, while moderate stimulation is conducive to the production of cartilage and bone tissue. At the same time, the Carter team focused on the regulation of hydrostatic pressure and tensile strain. By predicting the distribution of mechanical parameters at the fracture site, the corresponding relationship between different mechanical environments and the formation of fibrous tissue, cartilage and bone tissue was clarified [7].

In terms of research methods, scholars such as Lacroix have developed a new prediction algorithm, which successfully reproduced the spatial temporal evolution law of intramembrane osteogenesis and intracartilage osteogenesis during fracture healing. Existing research mainly uses strain invariant as mechanical signal indicators, but some scholars also use parameters such as fluid velocity and fluid shear stress for analysis. This study provide a multi-angle theoretical basis for understanding that mechanics stimulates and regulates tissue regeneration [8].

Epari and other scholars compared and analyzed the effects of shear/torsion load and axial compression on fracture healing through finite element simulation, and found that axial movement is more conducive to promoting healing, while shear movement may delay the repair process. This finding is consistent with clinical observation, that is, proper axial fretting can stimulate callus formation, while excessive shear force may interfere with healing. In addition, the study also shows that dynamic hydrostatic pressure can promote the synthesis of extracellular matrix, while static compression can inhibit cartilage formation, further confirming that the time-varying characteristics of mechanical stimulation play a key role in regulating tissue regeneration [9].

Ren et al. [10] developed a mechanical control model based on fuzzy logic, which can simulate the healing process under different fracture shapes and multi-axial loads, and successfully predict the occurrence of delayed healing or even nonunion. This model provides a theoretical basis for clinical assessment of fracture healing risk. In addition, Naveiro et al. [11] and Schwarzen et al. [12] introduced the grid growth algorithm, combined with the spatial neighborhood function to optimize the mechanical model, so that it can more accurately simulate the healing process of different fracture types (such as transverse, oblique or comminuted fractures).

Pietsch et al. [13] put forward a new numerical model based on fluid mechanics interface capture technology, which can describe the dynamic changes of tissue growth more accurately. Similarly, some studies also use multi-scale modeling method to combine the macro-mechanical environment with the biological response at the cell level to reveal the mechanical-biological coupling mechanism more comprehensively.

1.2 Biological-regulation healing models

Bailón-Plaza and Van Der Meulen presented a mathematical model that studied the effects of growth factors (TGF-β1 and BMP-2/4) on fracture healing. The model simulated mesenchymal cell migration, differentiation into chondrocytes and osteoblasts, and extracellular matrix synthesis and degradation. It predicted that osteogenic growth factor production by osteoblasts and the duration of initial growth factor release were critical for complete ossification. The model successfully replicated healing stages and highlighted the regulatory roles of growth factors in bone regeneration, providing a framework for understanding normal and pathological healing processes [14].

Carlier proposed a MOSAIC model, introducing a vascular endothelial growth factor, which enhanced the ability to study the impact of molecular mechanisms on angiogenesis and its relationship with bone formation, and allowed for research to be conducted in a more mechanical manner at different temporal and spatial scales [15].

Geris et al. developed a multi-scale computational model to investigate the role of oxygen in bone fracture healing. The model integrated oxygen's influence on cellular processes such as proliferation, differentiation, and apoptosis, based on experimental data. It combined tissue-level PDEs with cellular-level agent-based descriptions of angiogenesis, including the DII4-Notch signaling pathway. The model simulated normal healing in a small defect (0.5 mm), predicting tissue patterns consistent with experimental observations. It highlighted the importance of oxygen gradients, showing how hypoxia in central callus regions favored chondrogenesis, while higher oxygen near bone ends promoted osteogenesis. A sensitivity analysis confirmed the model's robustness to parameter variations. For a critical-sized defect (3 mm), the model predicted non-union due to severe hypoxia causing cell death and disrupted vascularization, emphasizing the need for timely oxygen delivery. Comparisons with the MOSAIC model (which lacked oxygen dependency) demonstrated that oxygen dynamics were critical for accurate predictions, particularly in impaired healing. The study concluded that oxygen tension is a key determinant of fracture healing, offering insights for improving therapies, such as enhancing vascularization in large defects. Limitations included simplified oxygen delivery assumptions and the need for further validation of diffusion coefficients in healing tissues [16].

Several studies have developed mathematical models to analyze the inflammatory phase of fracture healing. Kojouh [17] introduced a nonlinear ODE-based framework to examine early inflammation, integrating immune cells, tissue cells, and signaling molecules. Expanding on prior research, Trejo et al. [18] incorporated distinct macrophage phenotypes into their model to explore their role in inflammatory regulation. Similarly, Zhang et al. [19] formulated a PDE-based computational model to simulate cellular and cytokine interactions within the fracture callus, specifically assessing how TNF influences healing in both healthy and diabetic scenarios. These studies collectively enhance understanding of inflammation's regulatory mechanisms in bone repair.

1.3 Mechanobiological coupling regulation healing model

The mechanobiological coupling model for bone healing integrates finite element analysis (FEA) with partial differential equations (PDEs) to simulate the healing process. FEA is employed to quantify mechanical stimuli (e.g., stress, strain, and fluid dynamics) within the callus, while PDE systems model cellular behaviors such as migration, proliferation, and differentiation. These mechanical parameters are either directly incorporated into the biological equations or combined with molecular factors as inputs for fuzzy logic algorithms to predict tissue development. The modeling framework typically initiates with the fracture site's granulation tissue phase. FEA-derived mechanical data is then coupled with biological PDE systems to analyze how mechanical cues influence cellular activities. Pioneering work by Lacroix and Prendergast established the foundation for such biomechanical healing models, which later evolved to incorporate growth factor effects and vascular remodeling. For instance, Simon et al. introduced blood perfusion as a spatiotemporal variable in vascular network modeling, while Amnet and Hofer applied fuzzy logic to linear elastic FEA simulations for dynamic healing prediction. These approaches collectively enhance the understanding of cortical and trabecular bone regeneration mechanisms [20].

Isaksson et al. [21] developed a mechano-regulatory bone-healing model incorporating cell-phenotype-specific activities to improve predictions of tissue regeneration. It coupled cellular mechanisms—proliferation, differentiation, migration, and matrix production—with mechanical stimuli using finite element analysis. The model simulated four cell types and three extracellular matrices, updating tissue properties iteratively. Compared to phenomenological models, it more accurately predicted normal fracture healing, excessive mechanical loading effects, periosteal stripping, and impaired cartilage remodeling. The results demonstrated the importance of cell-specific rates and mechanistic approaches in understanding bone repair, offering potential for clinical applications and further research in mechanobiology.

Wang et al. proposed a 3D fracture-healing model combining finite element analysis and fuzzy logic to simulate tissue regeneration. It considered mechanical stimuli and blood supply, using 21 fuzzy rules to predict tissue differentiation. A transverse fracture model was developed, and interfragmentary movement (IFM) was analyzed under stable and unstable conditions. The results aligned with experimental data, showing decreased IFM over time as callus stiffness increased. Compared to prior 2D models, this approach improved accuracy, with average errors of 0.029 mm and 0.075 mm for stable and unstable cases, respectively. The model demonstrated the influence of mechanical and biological factors on healing but omitted growth factors. The research was implemented via Visual Studio 2012 and validated against experimental trends [22].

Frame and Schmitt developed a novel multi-tissue computational framework incorporating advanced mechanobiological principles, including tissue-specific differentiation, strain-mediated growth modulation, and mechanical regulation of tissue maturation. Their unified approach simultaneously simulated bone remodeling and fracture repair processes under shared mechanical regulatory mechanisms, demonstrating applicability in implant therapy. Separately, Vavva et al established a biologically driven predictive model for bone regeneration, emphasizing ultrasound-enhanced vascularization. Later collaborative work by Vavva and Carlier integrated ultrasonic stimulation into mechanobiological bone healing models, revealing how mesenchymal stem cell sources influence both healing trajectories and temporal dynamics while confirming ultrasound's therapeutic potential in osseous repair [23].

A key challenge in coupled mechanobiological modeling lies in the computational complexity of solving multiple partial differential equations and executing finite element analyses across iterative cycles. Research indicates that elevated mechanical stimuli typically promote fibrous tissue development, whereas moderated loading conditions favor osteogenesis. During each computational cycle, the model dynamically updates tissue distributions, subsequently recalibrating material properties and geometric configurations to inform mechanical parameter calculations for subsequent iterations [24].

Rousseau proposed a peri-implant healing model and presented a set of taxis-diffusion-reaction equations that comprehensively consider the combined influence of mechanical and biochemical factors [25]. Milan, through porous bone conduction biomaterials, imitated the osteogenic characteristics of the periosteum, determined the optimal loading for promoting bone formation and differentiation, and proposed that the rehabilitation plan for patients should be adapted to reproduce the best mechanical stimulation [26].

This paper presents a coupled mechanical-biological regulation model and a numerical modeling method for achieving fracture healing, revealing the intrinsic roles of various factors influencing fracture healing and providing assistance for the effective treatment of patients.

2 Materials and methods

A biomechanical model was developed to simulate fracture healing, incorporating the activities of mesenchymal stem cells, chondrocytes, and osteoblasts. The activities of these three cell types depend on mechanical stimuli, generating cartilage and bone tissues, which cause spatiotemporal changes in tissue distribution. Under specific mechanical, biochemical, and cellular conditions, individual cells may exhibit various behaviors, including migration, proliferation, differentiation, apoptosis, matrix production, or matrix degradation those based on their microenvironment and interactions with neighboring cells. The coupled model uses two separate modules: a fluid-solid phase coupling module and a bio-module described by partial differential equations for cell activities.

2.1 Mechanical module: finite element model to calculate mechanical stimulation

An equivalent simulation model of a tibia with a Taylor frame, as shown in Figure 1, was established for the mechanical module. This model is based on previous experiments [22,27]. The fixator is designed to allow appropriate axial movement of bone fragments while maintaining rigidity in other directions. According to the axial nonlinear force-deflection behavior previously described and illustrated in Figure 2, the external fixator is modeled as a spring system installed parallel to the callus [28]. The parameters of this spring system include a basic stiffness of K_bas = 4600 N/mm and a pretension F_pre ranging from 100 N to 800 N, as shown in Figure 3.

As shown in Figure 3, the proximal end of the fracture is subjected to force in the axial direction, while the distal end is fixed. The proximal and distal ends of the fracture are connected by an external fixator through surgery, thereby achieving mechanical balance.

The finite element model features a transverse bone fracture with a 4 mm gap, surrounded by an external callus. Fragment displacement is initialized at 1 mm, reflecting realistic fracture biomechanics and tissue patterns consistent with Claes and Heigele's observations. Leveraging axial symmetry, the computational model reduces the geometry to an eighth-section representation incorporating cortical bone, marrow, and callus tissues. Simulation boundary conditions define impermeable surfaces at: (1) the callus periphery, (2) cortical bone extremities, and (3) the medullary cavity, mimicking fascial encapsulation [28]. The mechanical regulation framework [27] governs cell fate determination, where mesenchymal stem cell differentiation (osteogenic/chondrogenic/fibrogenic) is controlled by a mechanical stimulus index S (S = γ/a+v/b). Here, γ represents octahedral shear strain in the solid matrix, v denotes interstitial fluid velocity, with empirical constants a = 0.0375 and b = 3 μm s⁻¹ defining the differentiation thresholds. In the early stages of fracture healing, values of S in the range of 0.267 ≤ S < 1 stimulate the differentiation into osteoblasts and bone; 0.01 ≤ S < 0.267 stimulates osteoblast differentiation and bone maturation; 1 ≤ S < 3 stimulates chondrocytes and cartilage formation, and S ≥ 3 stimulate fibroblasts and fibrous tissue formation, as shown in Figure 4.

Simulating physiologically relevant cyclic loading patterns [29], as shown in Figure 5, a cyclic load of 300 N is applied to the cortical bone at the top of the model. In the research of REN T [30], it was reported that the delay of fracture healing was caused by torsional instability. Therefore, it is required that the axial stiffness of the fracture fixation be sufficient to reduce the influence of torsion.

Biophysical stimuli are calculated under maximum load conditions. The mechanical analysis is based on standard poroelastic analysis conducted in ABAQUS. All finite elements are assigned linear elastic material properties. The cortical components comprise dense bone, while the callus comprises a varying tissue mixture. Each callus element has its material properties, updated at every time step of the simulation. The detailed process of these updates and how they are implemented will be described in the simulation implementation section.

Fig. 1 Three-dimensional model of tibial fracture fixed by Taylor frame.

Fig. 2 Nonlinear force-displacement relationship characterizing the external fixator's mechanical behavior.

Fig. 3 Two-dimensional axisymmetric finite element model of standardized fracture callus fixed with Taylor nonlinear frame.

Fig. 4 Mechano-regulation concept regulating cell differentiation to form fibrous tissue, cartilage, or bone. The two types of mechanical stimulation used are tissue strain and interstitial fluid flow. The arrow curve shows the differentiation pathway under biophysical stimulation.

Fig. 5 Loading mode of physiological load.

2.2 Biological part: PDE model to determine(describe) the cell process

Secondary fracture healing is typically divided into stages of inflammation, periosteal thickening, soft callus formation, and hard callus formation. To describe the tissue changes during the healing process more conveniently, the fracture area in the simulation model is divided into the following regions, as shown in Figure 6.

The fracture region is divided into five zones: the external callus zone (I), cortical bone callus zone (II), periosteal callus zone (III), cortical zone (IV), and bone marrow zone (V). In this study, we focused on cell differentiation and tissue formation in zones I-III, where mesenchymal stem cells may originate from bone marrow, periosteum or surrounding tissues, depending on the choice of clinical treatment surgery. These cells differentiate into osteoblasts, chondrocytes or fibroblasts under the biomechanical stimulation of the local environment, thereby generating cartilage and bone tissue in the callus and completing the healing process.

Due to its axisymmetry, the spatiotemporal changes in tissue distribution during healing are displayed in one-eighth of the area. The r-z coordinates are rotated for ease of observation, as shown in Figure 7. In this model, both the geometric and load conditions satisfy symmetry. In the simulation, when the mesh density is sufficient to capture the key features, the stress/strain distribution at the symmetrical interface and the error compared to the full model are usually less than 5%. This can be verified through a rough comparison of the 3D mesh.

The mathematical model proposed in this study builds upon the model by Bailón-Plaza and Van Der Meulen [14,31]. This model describes fracture healing as spatiotemporal changes in the densities of seven variables: mesenchymal stem cells (c_m), chondrocytes (c_c), osteoblasts (c_b), composite fibrous/cartilage extracellular matrix (m_c), bone cell-extracellular matrix (m_b), and general osteogenic (g_b) and chondrogenic (g_c) growth factors. Figure 8 provides a schematic diagram of tissue differentiation in this biological model.

The mathematical model describing the interactions of these variables is a slip-diffusion-reaction type partial differential equation, equations (1)–(7). During the process of fracture healing, the Slip-diffusion-reaction type partial differential equation can be used to describe the dynamic coupling process of cell migration, growth factor diffusion and bone tissue regeneration. The model has the ability to optimize fracture fixation schemes (such as mechanical loading strategies) and predict the risk of delayed healing.

Cell density, matrix density, and growth factor concentration are the main variables. These variables are denoted as c_i, m_i, and g_i, respectively, where the subscripts m, c, and b represent mesenchymal cells, chondrocytes, and osteoblasts, respectively. These subscripts are consistent in all variables, which are functions of two spatial variables (r, z) and one temporal variable, t. Cell density (number of cells per ml) is targeted at mesenchymal stem cells c_m = c_m(r,z,t), chondrocytes c_c = c_c(r,z,t), and osteoblasts c_b = c_b(r,z,t). The total matrix density (g ml⁻¹), m = m(r,z,t), is composed of cartilage matrix density m_c = m_c(r,z,t) and bone matrix density m_b = m_b(r,z,t). Finally, the concentrations of chondrogenic growth factor g_c = g_c(r,z,t) and osteogenic growth factor g_b = g_b(r,z,t) are modeled (ng ml⁻¹).

$$\frac{\partial c_m}{\partial t}=\nabla \cdot [D_m\nabla c_m-C_m{c}_m\nabla m]+c_m[P_m-B_m{c}_m]-F_b{c}_m-F_c{c}_m$$(1)

$$\frac{\partial c_c}{\partial t}=c_c[P_c-B_c{c}_c]+F_b{c}_m-F_{bc}{c}_c$$(2)

$$\frac{\partial c_b}{\partial t}=c_b[P_b-B_b{c}_b]+F_b{c}_m+F_{bc}{c}_c-d_b{c}_b$$(3)

$$\frac{\partial m_c}{\partial t}=(P_{cs}-Q_{cd1}{m}_c)\times (c_m+c_c)-Q_{cd2}{m}_c{c}_b$$(4)

$$\frac{\partial m_b}{\partial t}=(P_{bs}-Q_{bd}{m}_b)c_b$$(5)

$$\frac{\partial g_c}{\partial t}=\nabla [D_{gc}\nabla g_c]+E_{gc}{c}_c-d_{gc}{g}_c$$(6)

$$\frac{\partial g_b}{\partial t}=\nabla [D_{gb}\nabla g_b]+E_{gb}{c}_b-d_{gb}{g}_b$$(7)

where C_m and D_m represent the corresponding diffusion coefficients for chemotactic and haptotactic cells, respectively, it is assumed that the diffusion coefficients depend on the volume fractions of bone (m_b) and cartilage (m_c), as shown below: D_m = D_m0(1-m_c-m_b), where D_m0 is the initial diffusion coefficient. The total cell density is given by c_tot = c_m+c_c+c_b, where P_m, P_c, and P_b are the proliferation rates for each cell type. The coefficients B_m, B_c, and B_b are determined by the limit cell densities K_lm, K_lc, and K_lb, such that B_m = P_m/K_lm, B_c = P_c/K_lc, and B_b = P_b/K_lb. Functions F_b, F_c, and F_bc associate cell differentiation with growth factor concentration. P_m, P_c, P_b, F_b, F_c, and F_bc depend on the mechanical stimulus S. Using estimates calculated by Olsen et al, the dimensionless values of Pbs and Qbd in this case range from 0.2 to 2 [32].

The migration of mesenchymal stem cells consists of random and directed movement. Random movement is modeled as a haptokinetic process which refers to the random migration behavior of cells under certain immobilized stimulation. The form of the haptokinetic coefficient is based on Olsen et al. [32]. Random movement is affected by the total matrix density (defined as m = m_c+m_b), which causes cell movement to cease to a certain extent, thereby avoiding the excessive production of extracellular matrix.

$$D_m=\frac{D_h}{(K_h^2+\text{m})}m$$(8)

$$C_m=\frac{C_k}{{(K_k+m)}^2}.$$(9)

The proliferation of mesenchymal stem cells and other cell types is modeled using a logistic growth function, with the proliferation rate dependent on the surrounding matrix density [33].

$$P_i=\frac{P_{io}}{(K_i^2+m^2)}m\begin{array}{ccc}\hfill \mathrm{}\hfill & \hfill for\hfill & \hfill i=m,b,c\hfill \end{array}.$$(10)

Mesenchymal stem cell differentiation into osteoblasts or chondrocytes is mediated by osteogenic and chondrogenic growth factors [23].

$$F_b=\frac{Y_1}{H_1+{\text{g}}_b}{g}_b$$(11)

$$F_c=\frac{Y_2}{H_2+{\text{g}}_c}{g}_c.$$(12)

Finally, the replacement of chondrocytes by endochondral ossification depends on the density of connective/chondrocyte ECM and the concentration of osteogenic growth factors [34].

$$F_{bc}=\frac{m_c^6}{(B_{ec}^6+m_c^6)}\times \frac{Y_3}{(H_3+g_b)}{g}_b.$$(13)

Fig. 6 Diagram of the fracture cross-section depicting various zones within the injury site.

Fig. 7 Area showing temporal and spatial changes in cell differentiation at the fracture site (Coordinate rotation).Points a, b, and c respectively represent the cortical bone gap, the medullary cavity, and the external callus.

Fig. 8 (A) Diagram of cell proliferation, migration, and differentiation (B) Diagram of different stages in the fracture healing process. Proliferation (circular arrow), differentiation (vertical arrow), production of growth factors (gb: osteogenic growth factor, gc: cartilage growth factor), and extracellular matrix (mb: bone extracellular matrix, mc: cartilage extracellular matrix, m: the total matrix).

2.3 Geometry, initial and boundary conditions

The cell model defines initial conditions with mesenchymal stem cells (MSCs) originating from the periosteum, bone marrow interface, callus periphery, and scattered within the callus. Initially, the callus consists of granulation tissue, with most cell densities set to zero except for MSC influx from adjacent soft tissues, periosteal cambium, and marrow. Material properties initially reflect granulation tissue. MSC contributions are prioritized, with the periosteal layer playing a dominant role over surrounding soft tissues (Fig. 6). The initial density of the extracellular matrix is low, with two initial sources of osteogenic and chondrogenic growth factors corresponding to the damaged periosteum and the cortex beneath the hematoma, as shown in Figure 9. The initial concentrations of these growth factors are set as boundary conditions corresponding to certain healing days. The dimensionless boundary conditions used are c_m-init = 0.5 and 1, representing the MSC fronts from the surrounding soft tissue and the periosteum, respectively; m_init = 0.1, g_c-init = 20, and g_b-init = 20, corresponding to inflammatory conditions (Appendix). In the mechanical model, axial loads are applied at the top of the proximal fragment of cortical bone and supported axially at the distal fragment. Load is applied in the axial direction with a force of 300 N, representing the magnitude of load during normal walking post-operation. Nodes at the axis of rotational symmetry are fixed radially, and symmetry boundary conditions are applied to the model's sides. The fluid exchange between the callus and bone marrow is allowed through shared nodes, with flux across the external boundary being constrained.

Fig. 9 Geometric model specification. (A) Initial inflammatory hematoma with sparse extracellular matrix, and (B) chondrogenic and osteogenic growth factor origins which localized at the injury site and subperiosteal cortex, respectively.

2.4 Numerical implementation: iterative healing simulation

In this simulation framework, all cell densities are initialized to zero except those defined by boundary conditions. The process begins by running the cell model, which is called the mechanical model, utilizing the initial material properties of the callus. Once the mechanical simulation is completed, a subroutine transfers the mechanical stimuli from each element and integration point to the cell model. The cell model then simulates one day of the healing process, evaluating various active cellular processes. Subsequently, another subroutine transfers the updated material properties back to the mechanical model in preparation for the next iteration. This cycle is repeated until the fracture healing process is complete, as detailed in the flowchart shown in Figure 10. Throughout the simulation, the material properties of cortical bone and bone marrow remain unchanged, and all material properties are derived from the literature [27], as documented in Table 1.

The flowchart illustrates the differentiation of tissue over time. Biophysical stimuli values are calculated in the simulation by applying the mechanoregulation rules. Once the new tissue phenotype is determined, the “material properties” are adjusted using a mixing rule. The unit elastic modulus of the fracture site is the sum of the products of the elastic modulus of each pure tissue and the cube of the corresponding tissue concentration, as shown in equation (14); the Poisson's ratio is the sum of the products of the Poisson's ratio of each pure tissue and the corresponding tissue concentration, as shown in equation (15).

$$E_{ele}=E_{bone}{c}_{bone}^3+E_{cart}{c}_{cart}^3+E_{conn}{c}_{conn}^3$$(14)

$$v_{ele}=v_{bone}{c}_{bone}+v_{cart}{c}_{cart}+v_{conn}{c}_{conn}.$$(15)

The change in the healing state during the fracture healing process is equal to the current healing state plus the change in the healing state brought about by cell migration and differentiation. By taking the partial derivative of the overall healing state function with respect to time during the fracture healing process, the state change at any moment during the healing process can be obtained.

$$T_{i+1}=\frac{\partial T}{\partial t}\mathrm{\Delta }t+T_i.$$(16)

Therefore, the fracture healing process can be regarded as a cyclic iterative process with respect to time. The tissue state variable at time i + 1 is equal to the tissue state variable at time i plus the tissue change from time i to time i + 1. Cell differentiation is mainly influenced by the current mechanical environment and biological factors.

The simulation program for the cell model is written in Python, and the model's predicted healing response is presented by reporting the spatial distributions of cell density, extracellular matrix density, and growth factor concentration within the callus at different time points. The model's seven variables, c(t,x) (Eqs. 91–7)), must remain non-negative to prevent numerical instability. Besides ensuring the non-negativity of the model's numerical solution, the algorithm must also comply with the conservation of mass. The finite volume technique is employed due to its inherent mass-conserving characteristics. The spatial domain is covered by an equally spaced computational grid (with grid size x = y = 0.02 cm), and the nondimensional partial differential equations are solved on this grid using an alternating direction implicit (ADI) and explicit finite difference method.

This method features an in-built automatic step size control function to ensure that the error caused by each time step (0.24 hours) remains below the user-specified tolerance of 10⁻⁴m for displacement and 10⁻¹ pa for pressure, while minimizing computational costs as much as possible.

Fig. 10 Flow chart of dynamic fracture healing model predicting tissue differentiation over time.

3 Results

3.1 Tissue distribution

The results demonstrate that bone regeneration during fracture healing can be achieved using the proposed coupled model through iterative data exchange via subroutines. The sequence of tissue differentiation occurs in the same order as observed experimentally. The initial healing stage features granulation tissue as the primary component of the callus. Computational modeling forecasts bone formation via two pathways: intramembranous ossification and endochondral ossification through cartilage. The spatial arrangement of differentiated cells and tissues suggests this outcome stems from combined biophysical cues. The model simulates the fracture healing process described in experiments (as shown in Fig. 11). Mesenchymal stem cells (MSCs) diffuse into the fracture area from surrounding tissues such as muscle, periosteum, and bone marrow, forming the initial soft callus. By day 10, MSCs have spread across most of the fracture injury area (Fig. 11C), and the callus formation becomes prominent, completing by day 15. MSCs differentiate into various cells, such as chondrocytes and osteoblasts near the fracture site and at the cortical regions (Figs. 11A and 11B). This process reflects how cellular behavior and tissue formation during bone healing are influenced by the spatial and temporal distribution of mechanical and biological stimuli, leading to the organized development of different types of tissue essential for fracture repair.

Starting from day 5, chondrocytes gradually appear at the fracture site. This timing aligns with the arrival of mesenchymal stem cells in areas enriched with diffusible cartilage growth factors. By day 30, many proliferating chondrocytes are present near the fracture, and this condition persists until day 40. During this period, osteoblasts appear and converge along the cortical periosteal regions. There is a notable consistency between the locations of chondrocyte formation and cartilage growth factor diffusion points (as shown in Figs. 11A and 12A). This entire chondrocyte generation takes approximately 40 days, after which the endochondral ossification process commences. This reflects the transition from cartilage to bone, where the initially formed cartilage tissue is gradually replaced by bone tissue, signifying a critical phase in the fracture healing process. This sequence and timing are crucial for ensuring proper alignment and stability of the newly formed bone structure.

Starting from day 5, osteoblasts gradually begin to appear in the cortical periosteal regions. This coincides with the arrival of mesenchymal stem cells in areas enriched with osteogenic growth factors. By day 30, many proliferating osteoblasts are present near the fracture and continue until day 40. These cells advance along the fracture site, aligning with osteogenic growth factors' diffusion points (as shown in Figs. 11B and 12B). The generation of chondrocytes takes approximately 40 days, followed by the process of intramembranous ossification. After 40 days, the entire callus is filled with extracellular matrix from bone cells. This phase marks a crucial transition toward the maturation and stabilization of the fracture site, as the initial soft callus becomes more rigid and mineralized, contributing to the structural integrity necessary for complete bone healing. The consistency between osteoblast formation and the diffusion of growth factors highlights the crucial role of well-orchestrated biophysical stimuli in successful fracture repair.

Both endochondral and intramembranous ossification processes rely on their respective tissue formation constants, determining the rates of new bone formation. In the model, endochondral ossification is particularly sensitive to osteogenic growth factor production changes. Low levels of these growth factors can result in slow and incomplete endochondral ossification in the external callus, potentially leading to non-union of the fracture. Endochondral ossification predominantly occurs throughout most of the fracture callus area (as shown in Fig. 13A), the concentration of cartilage matrix in the fracture line area is approximately 75% in about 40 days. In contrast, intramembranous ossification takes place in the cortical periosteal regions (illustrated in Fig. 13B) and gradually approaches the fracture gap. In about 40 days, the concentration of bone matrix near the periosteal cortical bone reaches approximately 90%. Through the combined processes of endochondral and intramembranous ossification, bone tissue ultimately occupies the entire fracture callus (depicted in Fig. 13C). Following this phase, the bone enters the remodeling stage, culminating in the completion of bone healing. This progression underscores the importance of adequate growth factor presence and the orchestrated sequence of ossification processes, which ensure successful fracture repair by transforming the initial soft callus into a fully integrated, structurally sound bone.

The predicted data of our model was compared with the work of Pratik et al [35]. As can be seen from Table 2, the data trends predicted by both are basically consistent.

Fig. 11 Numerical solution of cell density after fracture. (A) chondrocytes (days 1, 5, 10, 20, 30, 40) (B) osteoblasts (days 1, 5, 10, 20, 30, 40) (C) mesenchymal stem cells (days 1, 5, 10, 15).

Fig. 12 Numerical solution of growth factor concentration after fracture. (A) chondrogenesis in callus (days 1, 5, 10, 20, 30, 40) and (B) osteogenesis in callus (days 1, 5, 10, 20, 30, 40) growth factor concentrations.

Fig. 13 Numerical solution of extracellular matrix density after fracture. (A) Cartilage extracellular matrix in callus (day 1, 5, 10, 20, 30, 40), (B) extracellular matrix in callus (day 1, 5, 10, 20, 30, 40) density and (C) total extracellular matrix in callus (day 1, 5, 10, 20, 30, 40).

3.2 Mechanical stimulation

The fluid flow patterns during bone fracture healing are quite complex. In the hematoma phase of healing, the callus consists of granulation tissue. It is subjected to large strains at the fracture gap, resulting in the highest fluid flow, which decreases as cartilage forms, As shown in Figure 14. When the bone begins to differentiate, as shown in Figures 14A and 14B, fluid flow increases because the woven bone is more permeable than cartilage. In the external callus, fluid flow is generally low. However, as mesenchymal stem cells differentiate extensively into chondrocytes and fill the entire callus tissue, this increase in solid-phase hardening reduces fluid flow capacity. Initially, fluid flow in the medullary cavity is predicted to be very low but relatively stable. Due to the difficulty in defining the permeability of various differentiated tissues, the calculated fluid flow here might deviate from reality. However, it does not affect the conclusion that fluid distribution within the callus tissue is uneven.

Perren (1979) suggested that as new tissue differentiates in the fracture gap during the healing process, the strain in the regenerative tissue gradually decreases. This study's numerical simulation confirmed the internal fragment strain theory (As shown in Fig. 15), showing that strain amplitude is greatest at the cortical bone gap (point b) at the start of healing. At the same time, it is relatively stable in the medullary cavity (point a) and the external callus (point c). As bone tissue forms, the strain gradually decreases. The simulation predicts that the reduction in strain across the fracture depends on the combined effects of external fixation and tissue differentiation. In the early stages of fracture healing, the external fixator provides primary support, and clinical practice aims to keep the strain below 20%. As bone tissue grows, adjusting the parameters of the external fixator removes stress shielding, allowing physiological loads to be applied to the bone tissue gradually.

During the period of fracture healing, the stimulation direction should be axial. Twisting is not conducive to fracture healing. Regarding the influence on mechanical strength, absolute stability (strain <2%) is required during the inflammatory stage to avoid interference with the process of hematoma organization. In the cartilage callus stage, the optimal strain is 2%−10% to stimulate endochondral ossification. Shear stress >10% is prone to cause fibrous tissue hyperplasia and delay mineralization. Gradually reduce the strain to <5% to promote lamellar bone replacement of woven bone.

Fig. 14 Prediction of fluid flow. Three characteristic positions, a, b, and c, were selected respectively in the medullary cavity, fracture gap, and external callus area.

Fig. 15 Prediction of octahedral shear strain. Three characteristic positions a, b, and c were selected respectively in the medullary cavity, fracture gap, and external callus region.

4 Discussion

In this study, a mechanobiological tissue differentiation model incorporating mixed-cell activity was developed to simulate various processes in fracture regeneration. A key feature of the simulation algorithm is its inclusion of processes such as cell diffusion, proliferation, differentiation, and matrix production in the callus, which are predicted to have a significant impact on healing patterns and rates (as shown in Fig. 11). This is combined with a pre-existing mechanoregulation algorithm, where the magnitudes of deviatoric strain and fluid flow velocity govern cell and tissue differentiation. The algorithm was chosen because it has been shown to provide more accurate predictions than other methods, aligning closely with experimental observations.

The model utilizes a Taylor frame for external fixation, which is minimally invasive and avoids the impact on mesenchymal stem cell sources seen with treatments like intramedullary nails and fixation plates. Despite this, the model remains widely applicable, capable of predicting changes in healing patterns due to compromised periosteal stripping and cartilage remodeling, thus closely mimicking a realistic fracture healing environment. Altering the extent of periosteal stripping has been a reliable experimental method for inducing delayed or non-union healing. The model can differentiate between various stem cell sources by setting different initial cell concentration values . For example, in treatments using intramedullary nails or fixation plates, where tissue differentiation in the fracture healing process relies on cells from the surrounding callus, external callus healing may become the primary healing mode, as illustrated in Figure 13.

This research investigates the regulatory dynamics between cellular activities and mechanical cues during bone repair using a combined mathematical and physiological approach. The computational framework effectively captures the temporal evolution of healing, demonstrating that mesenchymal stem cells require adequate time to disperse within the hematoma before achieving localized expression. Subsequently, these cells differentiate into cartilage and intramembranous bone tissue within the callus under growth factor mediation, while endochondral ossification progresses concurrently with cartilage maturation. As depicted in Figure 13A, the ossification front advances at a rate influenced by cartilage turnover and bone deposition, with the cartilaginous matrix completing mineralization first. The model further reveals that endochondral ossification persists over an extended duration.

When subjected to physiological loading, the framework accurately replicates both endochondral and intramembranous ossification processes, with mechanical response predictions consistent with existing studies. Variations in interfragmentary motion, clinically and experimentally observed, are modulated by fixator stiffness and preload. Overloading may lead to prolonged cartilage retention, causing delayed union or non-union, a phenomenon successfully simulated here. Comparable tissue alterations have been noted in prior tissue-level models under elevated mechanical stimuli.

Findings indicate that only integrative models incorporating biomechanical and biological variables can reliably predict healing outcomes across stable and unstable fixation scenarios. This iterative modeling approach facilitates the examination of diverse mechanical environments and cellular mechanisms, enhancing understanding of impaired healing. Such insights may refine clinical fixation strategies and reduce reliance on animal testing.

The internal fragment strain theory during the fracture healing process holds that the mechanical strain (tissue deformation) at the fracture ends directly determines the type and quality of the healing tissue. Strain within 2% promotes direct bone healing (primary healing), where osteoblasts can directly form lamellar bone without the need for cartilage intermediaries. Strain between 2% and 10% induces cartilage formation (indirect healing), generating callus through endochondral ossification. Strain above 10% leads to fibrous tissue proliferation, potentially causing delayed healing or non-healing. The scope of this study is indirect healing, while direct healing and abnormal healing will be the focus of our next stage of work.

5 Conclusion

In summary, the model proposed in this work successfully simulates tissue differentiation during the fracture healing process and provides a deeper understanding of the regulatory role of growth factors and their interactions with the mechanical environment. The results indicate that the tissue formation patterns predicted by the model are similar to those observed experimentally. An additional advantage of the model is its ability to evaluate relative timing, individual cell activity, concentration, and the impact of each parameter on the tissue regeneration process, allowing comparison with experimental results. Furthermore, the model can predict changes in healing patterns due to compromised periosteal stripping and cartilage remodeling.

During the normal indirect healing process of a fracture, an external fixator (such as the Taylor frame) provides initial stability during the hematoma stage, allowing mesenchymal stem cells to diffuse to the fracture site. Under the support of the external fixator and axial force, micro-movement at the fracture site stimulates the differentiation of mesenchymal stem cells into chondrocytes and osteoblasts. As the external fixator is dynamically adjusted, the mechanical stimulation is optimized, enabling the biological response to complete the remodeling of the bone structure.

The model includes components such as cartilage extracellular matrix, bone extracellular matrix, chondrocytes, osteoblasts, and mesenchymal stem cells. Although some literature includes fibroblast extracellular matrix, fibroblasts, oxygen levels, vascular endothelial cells, and intracellular factor activities in the model, these additions increase complexity. In normal healing processes, fibroblasts and the fibrous matrix do not play significant roles. In contrast, endothelial cells and oxygen levels can be considered variables, which will be accounted for in future research. Other limitations to be noted include: (1) Limitations of structural finite element analysis: This is constrained by assumptions regarding geometry, loading, and material properties, with material attributes sourced from relevant literature. (2) Simplified model assumptions: An axisymmetric finite element model simulates completely symmetric transverse fractures with stable external callus without considering the natural growth of the callus. The fracture gap size is maintained constant through external fixation. (3) Load type limitation: Only axial loads are considered, while the effects of torsional loads are ignored. These conclusions provide a valuable foundation for further research while highlighting the need to explore more complexities and real-world impacts in healing processes.

Funding

This work was supported by the Zhejiang Science and Technology Plan Project—National Key Research and Development Program of China (Grant No. 2019C03075).

Conflicts of interest

All authors certify that they have no affiliations with or involvement in any organization or entity with any financial or non-financial interest in the subject matter or materials discussed in this manuscript.

Data availability statement

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Author contribution statement

Qinghua Yang, Yi Xun, and Xinxing Zhang contributed to the conception of the study. Qinghua Yang, Xinxing Zhang, and Jun Qian performed the experiment and analyzed the data. Qinghua Yang obtained grant funding. Xinxing Zhang and Yi Xun carried out software programming. Xinxing Zhang and Yunsheng Mao wrote the original draft. Qinghua Yang and Yi Xun performed the review and editing. All authors read and approved the final manuscript.

References

All Tables

All Figures

	Fig. 1 Three-dimensional model of tibial fracture fixed by Taylor frame.
In the text

	Fig. 2 Nonlinear force-displacement relationship characterizing the external fixator's mechanical behavior.
In the text

	Fig. 3 Two-dimensional axisymmetric finite element model of standardized fracture callus fixed with Taylor nonlinear frame.
In the text

	Fig. 4 Mechano-regulation concept regulating cell differentiation to form fibrous tissue, cartilage, or bone. The two types of mechanical stimulation used are tissue strain and interstitial fluid flow. The arrow curve shows the differentiation pathway under biophysical stimulation.
In the text

	Fig. 5 Loading mode of physiological load.
In the text

	Fig. 6 Diagram of the fracture cross-section depicting various zones within the injury site.
In the text

	Fig. 7 Area showing temporal and spatial changes in cell differentiation at the fracture site (Coordinate rotation).Points a, b, and c respectively represent the cortical bone gap, the medullary cavity, and the external callus.
In the text

Fig. 8 (A) Diagram of cell proliferation, migration, and differentiation (B) Diagram of different stages in the fracture healing process. Proliferation (circular arrow), differentiation (vertical arrow), production of growth factors (gb: osteogenic growth factor, gc: cartilage growth factor), and extracellular matrix (mb: bone extracellular matrix, mc: cartilage extracellular matrix, m: the total matrix).

In the text

	Fig. 9 Geometric model specification. (A) Initial inflammatory hematoma with sparse extracellular matrix, and (B) chondrogenic and osteogenic growth factor origins which localized at the injury site and subperiosteal cortex, respectively.
In the text

	Fig. 10 Flow chart of dynamic fracture healing model predicting tissue differentiation over time.
In the text

	Fig. 11 Numerical solution of cell density after fracture. (A) chondrocytes (days 1, 5, 10, 20, 30, 40) (B) osteoblasts (days 1, 5, 10, 20, 30, 40) (C) mesenchymal stem cells (days 1, 5, 10, 15).
In the text

	Fig. 12 Numerical solution of growth factor concentration after fracture. (A) chondrogenesis in callus (days 1, 5, 10, 20, 30, 40) and (B) osteogenesis in callus (days 1, 5, 10, 20, 30, 40) growth factor concentrations.
In the text

	Fig. 13 Numerical solution of extracellular matrix density after fracture. (A) Cartilage extracellular matrix in callus (day 1, 5, 10, 20, 30, 40), (B) extracellular matrix in callus (day 1, 5, 10, 20, 30, 40) density and (C) total extracellular matrix in callus (day 1, 5, 10, 20, 30, 40).
In the text

	Fig. 14 Prediction of fluid flow. Three characteristic positions, a, b, and c, were selected respectively in the medullary cavity, fracture gap, and external callus area.
In the text

	Fig. 15 Prediction of octahedral shear strain. Three characteristic positions a, b, and c were selected respectively in the medullary cavity, fracture gap, and external callus region.
In the text