Skip to main content

Advertisement

Computational modeling of human bone fracture healing affected by different conditions of initial healing stage

Abstract

Background

Bone healing process includes four phases: inflammatory response, soft callus formation, hard callus development, and remodeling. Mechanobiological models have been used to investigate the role of various mechanical and biological factors on bone healing. However, the effects of initial healing phase, which includes the inflammatory stage, the granulation tissue formation, and the initial callus formation during the first few days post-fracture, are generally neglected in such studies.

Methods

In this study, we developed a finite-element-based model to simulate different levels of diffusion coefficient for mesenchymal stem cell (MSC) migration, Young’s modulus of granulation tissue, callus thickness and interfragmentary gap size to understand the modulatory effects of these initial phase parameters on bone healing.

Results

The results quantified how faster MSC migration, stiffer granulation tissue, thicker callus, and smaller interfragmentary gap enhanced healing to some extent. However, after a certain threshold, a state of saturation was reached for MSC migration rate, granulation tissue stiffness, and callus thickness. Therefore, a parametric study was performed to verify that the callus formed at the initial phase, in agreement with experimental observations, has an ideal range of geometry and material properties to have the most efficient healing time.

Conclusions

Findings from this paper quantified the effects of the initial healing phase on healing outcome to better understand the biological and mechanobiological mechanisms and their utilization in the design and optimization of treatment strategies. It is also demonstrated through a simulation that for fractures, where bone segments are in close proximity, callus development is not required. This finding is consistent with the concepts of primary and secondary bone healing.

Peer Review reports

Background

Bone healing is a complex four-phase process, which starts with an inflammatory response and hematoma formation, resulting in granulation tissue development at 3–7 days post-fracture. Following this initial phase, a cartilaginous soft callus is formed from the granulation tissue in 2–4 weeks. After 2–4 months, this formation develops into a hard bony callus that surrounds the fracture site. The ossified callus is restructured for several months to years until the final bone structure is achieved, which generally resembles the original (pre-fracture) morphology of the bone [1, 2]. While the bone healing process has been experimentally studied for several decades [3,4,5,6,7], mechanobiological models have been used more recently to study the effects of both mechanical loading and biological factors on cellular activities and tissue formation following fracture [1, 8]. Such models can be used to study different factors that impact the healing process; to predict outcomes under different mechanical or biological conditions; and in response to new treatment strategies [9,10,11].

In mechanobiological modeling, mechanical factors such as strain or stress in fracture sites are typically estimated using finite element (FE) analysis. Mechanical stimuli, biological factors, and chemical stimuli influence biological processes and cellular activities, such as mesenchymal stem cell (MSC) migration, tissue differentiation, angiogenesis, and growth factor secretion, which in turn influence and regulate the bone healing process [1, 12,13,14,15,16,17,18,19,20]. Most mechanobiological models of bone healing consider a predefined callus with an ideal fixed geometry and predefined material properties [12,13,14, 21], where they neglect the initial phases of healing (i.e. the inflammatory stage, hematoma evolution to form granulation tissue and initial callus development during the first few days post-fracture) [1]. However, few studies have accounted for callus geometry development in their simulations by assuming that it is similar to volume expansion due to the application of thermal loading [22,23,24] or swelling pressure [9, 25]. These numerical mechanisms (i.e. thermal expansion and swelling pressure) are regulated through mechanobiological rules and should be considered as an improvement in accounting for callus geometry development; however, they may not simulate the actual mechanism of callus geometry development, especially during the initial phase of healing [1]. Another limitation of the current studies is characterization of the material properties of the hematoma and granulation tissue during the initial phase [1, 21, 26, 27].

On the other hand, a growing body of experimental studies has highlighted the critical role of initial phases of healing on the bone healing process and outcome [2]. For instance, inhibiting the initial post-fracture inflammatory response through anti-inflammatory treatment has been reported to impair granulation tissue formation and callus development, consequently delaying or preventing healing [28, 29]. Moreover, interfragmentary gap size and initial stability of the fracture site (i.e., fixation level of interfragmentary motion) are critical factors, which specify the form of healing (i.e., primary or secondary healing) and the recovery time. In primary bone healing, where the distance between bone fracture surfaces is very small and is completely constrained by fixation, no callus is formed. Secondary bone healing involves callus formation, where callus size partially depends on the interfragmentary motion levels conducive to healing [23, 30,31,32,33,34]. Moreover, the callus geometry is shown to be an optimal shape to endure the mechanical loading during the healing process [35,36,37]..

Therefore, we hypothesize that the initial phase has a contributory mechanobiological effect on the overall bone healing process, resulting in the formation of an initial callus with an ideal range of geometry and material properties to achieve the most efficient healing time. To that end, we utilized a pre-developed finite element-based model by Lacroix & Prendergast (2002) [30] to simulate the bone healing process in models with different diffusion coefficients of MSC migration, granulation tissue Young’s moduli, callus geometries, and interfragmentary gap sizes. These parameters modulate the outcome of bone healing during its initial phase, which involves inflammatory stage, hematoma evolution to form granulation tissue and initial callus development during the first few days post-fracture. The diffusion coefficient can specify local levels of MSC density, especially during the initial post-fracture days [23]. The elastic modulus of granulation tissue determines the mechanical response level of the fracture site during the initial phase [1]. The mechanical response of fracture sites and MSC density depend on callus thickness during the healing process, including the initial phase [30]. Interfragmentary gap size and mechanical stability of fracture site can alter the callus thickness especially at the inflammatory stage and soft callus phase [32, 38]. In this parametric study, we aim to investigate how these factors and the callus developed at the initial healing phase influence healing time and healing pattern.

Methods

The mechanobiological regulation outlined by Prendergast et al. (1997) [17] was utilized to determine tissue differentiation type under applied mechanical loading (Fig. 1-A). As a general expression, high levels of mechanical stimuli result in fibrous tissue formation, intermediate levels promote cartilaginous tissue formation, and lower levels lead to bone formation. This mechanobiological regulation was smoothed and modified based on Sapotnick and Nackenhorst’s work [39], in order to prevent abrupt changes in tissue differentiation categories (Fig. 1-B) [39].

Fig. 1
figure1

a Mechanobiological regulation by Prendergast et al. (1997) [17]. b Smoothed mechanobiological regulation based on Sapotnick and Nackenhorst (2015) [39]. c Left: Callus geometry dimensions, including thickness (d) and interfragmentary half gap size (h). Right: FE mesh and boundary conditions of stress analysis where the blue elements are marrow, green elements are bone, and red elements are callus

A human bone shaft was modeled as a hollow cylinder with a transverse cut perpendicular to the cylindrical axis. An axisymmetric biphasic finite-element analysis of the bone was developed using linear poroelastic material properties for the involved tissues according to the model presented by Lacroix & Prendergast (2002) [30]. The FE model was made of the 4-node quadrilateral, bilinear displacement, and bilinear pore pressure elements (Fig. 1-C, right). For the base model with a 4 mm callus thickness (i.e. d = 4 mm in Fig. 1-C left) and a 3 mm interfragmentary gap size (i.e. h = 1.5 mm in Fig. 1-C left), there were 311 elements in the marrow, 366 elements in the bone fragment and 2034 elements in the callus (Fig. 1-C). Boundary conditions were applied at the bottom and left borders of the model as shown in Fig. 1-C, left. Bone, bone marrow, cartilage, and fibrous tissue were modeled as linear poroelastic biphasic materials [40,41,42], with material properties shown in Table 1 [21, 30]. The bone healing process was simulated for up to 120 iterations (days), with results obtained for each day using an iterative process. The iterative simulation of healing process was stopped either when 120 iterations were completed or sooner when a complete bony callus was achieved (i.e. a complete bony callus is achieved when every element of callus gains Young’s modulus higher than 2 GPa). In each iteration, an axial load was applied to the top end of the bone and was increased linearly from 0 to 500 N in 1 s, similar to the model presented by Lacroix & Prendergast (2002) [30], to calculate fluid flow and octahedral shear strain for each element (ABAQUS version 6.13–2, Simulia, Providence, RI, USA). Through a separate finite element-based diffusion analysis, MSC migration was simulated for every iteration of the simulation to determine the spatial and temporal MSCs distribution using \( \frac{\partial c}{\partial t}=D{\nabla}^2c \), where c is the MSC density, D is the diffusion coefficient of MSC migration and t is time. For the base model, a value of 0.5 mm2/day was considered as the MSC diffusion coefficient. Bone marrow and periosteal surface of the bone and soft tissues surrounding the callus were considered as the MSC migration sources. As MSCs migration initiates during the initial phase of healing, we adjusted the initial MSC distribution accordingly. Therefore, we first performed a preliminary MSC diffusion analysis to calculate the local MSC density in day 7 post-fracture (i.e. by the end of the initial healing phase). Then we started the first iteration of the healing simulation by implementing the preliminary MSC density. A mesh convergence study was performed for the finite element analysis of the base model and the models with different values of MSC diffusion coefficients to eliminate any mesh dependence in the final results.

Table 1 Material properties [21, 30]

Cells within each callus element differentiated into tissues or matrices such as bone, cartilage, or fibrous tissue as a result of the local state of mechanical parameters and MSC density. Following the rule of mixtures, the average material properties of the newly formed tissue and those from the nine previous days were calculated at each step of calculation to update each element’s material properties [21]. The updated material properties were used in FE analyses of the next iteration.

We repeated the numerical simulation for models with a wide range of diffusion coefficients of MSC migration, granulation tissue elastic moduli (denoted by Eg), callus thicknesses (denoted by d) and interfragmentary half gap sizes (denoted by h). To specify an appropriate range of variation for each parameter, we considered a base model [30] with normal values of 0.5 mm2/day, 1 MPa, 4 mm and 1.5 mm for MSC diffusion coefficient, granulation tissue Young’s modulus, callus thickness, and interfragmentary half gap size, respectively. For the upper bound of MSC diffusion coefficient range, it was increased until a state of saturation was observed and for the lower bound, it was reduced until nonunion or delayed healing was observed. For other parameters, a similar approach was conducted to determine the upper and the lower bounds. However, we stopped at 2 MPa for the upper bound of granulation tissue Young’s modulus, since values higher than 2 MPa are even stiffer than fibrous tissue or bone marrow, which is not probable for a relatively fresh blood clot [43]. As a result, the following domains of variables have been specified (please see requisite scripts in Additional files 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31 and 32):

  • [0.001, 0.01, 0.1, 0.5, 1, 10, 100] mm2/day for MSC diffusion coefficient

  • [0.01, 0.05, 0.1, 0.2, 0.5, 1, 2] MPa for Young’s modulus of granulation tissue

  • [1,2,3,4,5,6,7,8] mm for callus thickness

  • [0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4] mm for interfragmentary half gap size.

We considered the time associated with complete development of the following structures as possible healing indices: 1) cartilaginous callus (CC), 2) bony bridging (BB), and 3) bony callus (BC) [31, 44]. It was assumed that a cartilaginous callus is developed when a cartilaginous connection is formed between two bone fragments (i.e., a sequence of elements exists with Young’s modulus higher than 10 MPa to connect the bone fragment with bottom border of the callus) [44]. Bony bridging is achieved when a bony connection forms between the two bone fragments (i.e. a sequence of elements exists with Young’s modulus higher than 2 GPa to connect the bone fragment with bottom border of the callus). Finally, a bony callus is achieved when every element of the whole callus has Young’s modulus greater than 2 GPa [21].

Results

The simulation results for models with different levels of diffusion coefficients varying from 0.001 mm2/day to 100 mm2/day are outlined in Fig. 2. At the start of the simulation, MSCs migrate from the three sources mentioned above into the fracture site. For the cases with diffusion coefficient of 0.5 mm2/day, level of MSCs density within the whole callus was greater than 50% of the maximum allowed cell density at day 5. When diffusion coefficient increased to 10 mm2/day, level of MSCs density became greater than 50% of the maximum allowed cell density at the end of day 1. However, when diffusion coefficient decreased to 0.1 mm2/day, level of MSCs was higher than 50% of the maximum allowed cell density after 40 days, and when it decreased to 0.01 mm2/day, level of MSCs never reached the greater than 50% of maximum allowed cell density threshold in 120 days. In models with a small diffusion coefficient value (i.e. 0.001 to 0.01 mm2/day), a long delay in healing was predicted, resulting in the formation of an incomplete bony callus after 120 days. Models with a diffusion coefficient in the range of 0.1 to 1 mm2/day predicted a normal healing process with the formation of a complete bony callus within 120 days of simulation. Further increase in the diffusion coefficient affected neither the healing process nor the timeline. Moreover, interfragmentary strain reduced by 0, 5 and 10% at day 1, day 3 and day 7, respectively, and maximum fluid flow reduced by 0, 0 and 5% at the same days, respectively, when MSC diffusion coefficient increased from 0.5 mm2/day to 100 mm2/day. On the other hand, interfragmentary strain increased 0, 7 and 50% at day 1, day 3 and day 7, respectively, and maximum fluid flow increased 0, 0 and 20% at the same days, respectively, when MSC diffusion coefficient reduced from 0.5 mm2/day to 0.01 mm2/day.

Fig. 2
figure2

a Healing pattern at different days during the healing process. The days are selected to show the onset of cartilaginous callus (CC), bony bridging (BB) and bony callus (BC) formation in models with different diffusion coefficients, D. b Effect of the diffusion coefficient on the healing duration (i.e. No. of days) associated with the onset of cartilaginous callus, bony bridging and bony callus formation. In this set of simulations, Eg = 1 MPa, d = 4 mm, and h = 1.5 mm

The role of granulation tissue Young’s modulus on the healing process is demonstrated in Fig. 3. No considerable changes were observed in the healing outcome for elastic modulus values ranging from 0.01 to 0.2 MPa since cartilaginous callus occurred at day 23 to 25, bony bridging occurred at day 46 to 48, and bony callus occurred at day 66 to 70. However, by increasing the elastic modulus from 0.2 MPa to 2 MPa, cartilaginous callus was formed 10 days earlier, while bony bridging occurred 16 days earlier, followed by the development of bony callus 24 days earlier. Also, interfragmentary strain reduced by 33, 37 and 45% at day 1, day 3 and day 7, respectively, and maximum fluid flow reduced by 0, 14 and 36% at the same days, respectively, when Young’s modulus of granulation tissue increased from 1 MPa to 2 MPa. On the other hand, interfragmentary strain increased 306, 257 and 144% at day 1, day 3 and day 7, respectively, and maximum fluid flow increased 237, 212 and 190% at the same days, respectively, when granulation tissue Young’s modulus reduced from 1 MPa to 0.1 MPa.

Fig. 3
figure3

a Healing pattern at different days during the healing process. The days are selected to show the onset of cartilaginous callus (CC), bony bridging (BB) and bony callus (BC) formation in models with different elastic moduli of granulation tissue, Eg. b Effect of granulation tissue’s elastic modulus on the healing duration (i.e., No. of days) associated with the onset of cartilaginous callus, bony bridging, and bony callus formation. In this set of simulations, D = 0.5 mm2/day, d = 4 mm, and h = 1.5 mm

The modeling results for different sizes of callus thickness are exhibited in Fig. 4. An extremely small callus thickness (1 mm) was predicted to develop into a fibrous callus and nonunion. A small 2 mm callus thickness progressed to a cartilaginous callus in 2 months, a bony bridge in 3 months and bony callus in 4 months. A callus thickness range from 3 to 6 mm led to a cartilaginous callus in 2–3 weeks, bony bridge in 4–6 weeks and complete bony callus in 6–10 weeks. Callus thicknesses greater than 6 mm enhanced the speed of bone healing, as the bony callus was completed within 1 month for thicknesses ranging from 7 to 8 mm. Moreover, interfragmentary strain reduced by 3, 31 and 75% at day 1, day 3 and day 7, respectively, and maximum fluid flow reduced by 3, 32 and 71% at the same days, respectively, when callus thickness increased from 4 mm to 8 mm. On the other hand, interfragmentary strain increased 13, 20 and 52% at day 1, day 3 and day 7, respectively, and maximum fluid flow increased 78, 0 and 36% at the same days, respectively, when callus thickness decreased from 4 mm to 1 mm.

Fig. 4
figure4

a Healing pattern at different days during the healing process. The days are selected to show the onset of cartilaginous callus (CC), bony bridging (BB), bony callus (BC), and fibrous callus (FC) formation in models with different callus thicknesses, d. b Effect of callus thicknesses on the healing duration (i.e., No. of days) associated with the onset of cartilaginous callus, bony bridging, and bony callus formation. In this set of simulations, D = 0.5 mm2/day, Eg = 1 MPa, and h = 1.5 mm

The effect of interfragmentary half gap size on bone healing, where h is varied between 0.5 mm to 4 mm is shown in Fig. 5. For a 0.5 mm interfragmentary half gap size, a cartilaginous callus was predicted at day 4, bony bridging occurred at day 13, and complete bony callus occurred in 33 days. For a 4 mm interfragmentary half gap size, cartilaginous callus was achieved in 1 month, bony bridging occurred in 2 months, and complete bony callus occurred in 3 months. An increase in interfragmentary half gap size from 0.5 mm to 4 mm consistently delays the bone healing process, increasing the healing time. Also, interfragmentary strain reduced by 40, 62 and 81% at day 1, day 3 and day 7, respectively, and maximum fluid flow reduced by 0, 25 and 59% at the same days, respectively, when interfragmentary half gap size reduced from 1.5 mm to 0.5 mm. On the other hand, interfragmentary strain increased 35, 42 and 84% at day 1, day 3 and day 7, respectively, and maximum fluid flow increased 144, 129 and 217% at the same days, respectively, when interfragmentary half gap size increased from 1.5 mm to 4 mm.

Fig. 5
figure5

a Healing pattern at different days during the healing process. The days are selected to show the onset of cartilaginous callus (CC), bony bridging (BB) and bony callus (BC) formation in models with different interfragmentary half gap sizes, h. b Effect of interfragmentary half gap sizes on the healing duration (i.e., No. of days) associated with the onset of cartilaginous callus, bony bridging, and bony callus formation. In this set of simulations, D = 0.5 mm2/day, Eg = 1 MPa, and d = 4 mm

The day corresponding to the onset of bony bridging for three different callus thicknesses (d = 3, 5, and 7 mm) is shown in Fig. 6, where the MSC diffusion coefficient is varied between 0.01 and 10 mm2/day. The results are presented for three different values of granulation tissue Young’s modulus (Eg = 0.1, 1 and 2 MPa). It should be noted that for the callus thickness of 1 mm, boney bridging does not occur in 120 days in the simulations, regardless of the level of MSC diffusion coefficient and granulation tissue Young’s modulus considered in this set of simulations. Thus, no results are shown for the callus thickness of 1 mm. In general, the onset of bridging occurs quicker for the models with a thicker callus. Faster MSC migration and a stiffer granulation tissue also expedite the healing, resulting in a quicker formation of boney bridging.

Fig. 6
figure6

Onset of bony bridging in models with different callus thicknesses, MSC diffusion coefficient, and Young’s modulus of granulation tissue

Discussion

We used a well-established model of the bone healing process presented by Lacroix & Prendergast (2002) [30] to design a parametric study in order to computationally quantify effects of the initial phase of healing on the healing outcome. We reviewed the performance of our simulation approach and FE model to see whether the results are compatible with previous computational studies and experimental observations. In the base model, our numerical simulations predict that cartilaginous callus is achieved in 2–3 weeks from the start of the healing process, bony bridging occurs in 1 month, and complete bony callus is developed in less than 2 months. This development timeline matches fairly well with clinical observations, as well as the results presented in previous numerical investigations [5, 21, 30]. In addition to the timeline, pattern of tissue formation in our simulation is comparable with other studies [21, 30, 45]. Bone formation initially begins in the external region of the initial callus, far from the fracture site where mechanical stimuli are at their lowest local values [45, 46]. Gradually this initial bone formation provides mechanical support for the fracture site and thus reduces mechanical stimuli and initiates bone formation at other regions of callus such as near the bone marrow and fracture gap [30, 46].

As the simulation results outlined, models with a larger value of diffusion coefficient for MSC migration, stiffer granulation tissue, and a thicker callus thickness predict lower level of mechanical stimuli and faster healing process. An increase in the diffusion coefficient for MSC migration means that the MSCs can be distributed more rapidly across the callus area for differentiation. In our simulations, when the diffusion coefficient was less than 0.1 mm2/day, our simulations indicated insufficient supply of MSCs to support differentiation and tissue formation, which subsequently delayed healing or resulted in nonunion. This is consistent with the nonunion results predicted by Geris et al. [11], when the MSC sources of migration were removed. On the other hand, no considerable change in the healing process was observed by increasing the diffusion coefficient to values greater than 1 mm2/day. At this stage, MSCs are present in high volume in the callus, and thus the healing time is rather limited by MSC differentiation or tissue formation rates. In other words, MSCs are readily available throughout the callus, but no improvement in healing occurs, as MSCs cannot differentiate and form tissue at a faster rate [13, 47].

In addition, callus development serves to support mechanical loading and provide the desired stability for bone healing [8, 35, 36]. Also, a significant increase in fluid flow was observed for small and soft calluses which may shift the mechanical stimuli to the domain where mostly fibrous tissue can be formed. Hence, calluses with small thicknesses or those made of soft granulation tissue are not able to support the applied mechanical loading and provide a suitable environment for the proper tissue type formation. Based on this study, callus thicknesses smaller than 3 mm or granulation tissue softer than 0.5 MPa resulted in delayed healing or nonunion. On the other hand, a callus thicker than 6 mm does not result in improvements in healing. Larger callus size results in prolonged resorption and remodeling [48]. Granulation tissues with an elastic modulus higher than 2 MPa are even stiffer than fibrous tissue or bone marrow, which is not probable for a relatively fresh blood clot [43]. Therefore, after a certain level, there is no need for a larger or stiffer callus to support mechanical loading and stabilize the fracture site. According to the findings of this study, there is an ideal range that has also been observed in experimental studies [5, 49,50,51,52] (Fig. 7) for each initial phase parameter (i.e. 0.1–1 mm2/day for migration rate, 1–2 MPa for Young’s modulus of granulation tissue, 3–6 mm for callus thickness). As shown in Fig. 7 A, MSCs mostly spread out over the callus during the first week in our simulations with the ideal range of migration rate. On the other hand, experimental observations by Iwaki et al. [49] and Wang et al. [53] showed that MSCs mostly spread out over the rat callus during the day 2 to day 7 post-fracture. As shown in Fig. 7 B, the granulation tissue indentation modulus, measured by Leong et al. [26] in a rat (i.e. 0.99 MPa), completely matches the ideal range of granulation tissue Young’s modulus predicted in this paper (i.e. 1–2 MPa). As shown in Fig. 7 C, the predicted ideal range of callus geometry and gap size was also in agreement with the experimental observations made by De Bastiani et al. [52] in humans, Augat et al. [54], and Epari et al. [46] in sheep, and Boer et al. [55] in goats. Thus, simulation results interestingly outline that the formed callus at the initial phase of healing (i.e. normal healing that is observed in experimental studies and clinical environments) contains optimal geometry and material properties to have the most efficient healing time.

Fig. 7
figure7

Comparison of the simulation results (i.e., the optimal ranges for the initial healing phase parameters) with the experimental studies for a: MSC diffusion, b: Granulation tissue elastic stiffness, and c: Callus size. (with permission to reuse from the publishers)

As indicated by our results, increasing the interfragmentary gap size delays bone healing, and shrinking the gap expedites it [30, 38]. This was seen in simulations with a 0.5 mm interfragmentary half gap size, where bony bridging and complete bony callus formation occurred in 2 weeks and 1 month, respectively. The remarkable impact of smaller interfragmentary gap size motivated us to investigate its effects on the smallest callus sample with the thickness of 1 mm (i.e., the sample where no sign of healing was seen in 120 days when combined with a 1.5 mm interfragmentary half gap size) (Fig. 4). Interestingly, a normal pattern of healing was observed when a very small 0.25 mm interfragmentary half gap size was paired with a very small 1 mm-thick callus (Fig. 8). The results matched experimental and clinical observations [32, 33, 38] and emphasized that a larger callus is necessary, when the interfragmentary gap is enlarged, to have a normal pattern of healing. Figure 8 also indicated that if bone fragments were tightly positioned with respect to one another in the secondary form of bone healing, almost no callus development would be needed which was in agreement with the concept of primary bone healing [4, 56]. These findings highlight the potential capability of bone healing models in understanding the basis and plausible mechanisms behind clinical observations [10].

Fig. 8
figure8

Healing pattern at different days during the healing process. The days are selected to show the onset of cartilaginous callus (CC), bony bridging (BB), bony callus (BC), and fibrous callus (FC) formation in models with different interfragmentary half gap sizes, h. In this set of simulations, D = 0.5 mm2/day, Eg = 1 MPa, and d = 1 mm

The quality of cartilaginous callus, position of bony bridging, and pattern of healing can also be affected by changes in the initial phase of healing. An increase in diffusion coefficient shifts the bony bridging position from the exterior of the callus towards the middle, and increase the average stiffness of the cartilaginous callus. Increasing the interfragmentary gap size also changes the position of bony bridging from the exterior of the callus to the inside. However, in some cases, it is not entirely clear how the initial phase affects the healing pattern. For instance, no differences were observed in the bony bridging position or cartilaginous callus stiffness, following the change in callus thickness or elastic modulus of granulation tissue.

As one of the limitations of this study, we only focused on material properties and geometrical factors of the initial callus as the outcome of the initial phase of healing. Other factors such as angiogenesis, effects of growth factors, oxygen tension, or type of loading were not directly investigated since complementary experimental studies are needed to provide reliable data to include them in the simulation. Also, material properties of the granulation tissue in the initial phase of healing, including elastic modulus or diffusion coefficient for MSC migration, have not been studied and analyzed well under different conditions of healing [1, 21]. Therefore, a precise range of material properties is not available for the initial callus formed at the initial phase to compare with our simulation results. However, some estimates have been conducted in previous simulation studies for the material properties of granulation tissue, which are in agreement with our reported optimal range [13, 47]. Similar to the relevant computational studies [13, 21], we modeled the involved tissues by homogenous material properties which is a concern during the initial healing phase as more heterogeneity is expected there. Thus, further experimental investigations on the material properties of the newly formed tissues in the initial healing phase are required to establish the required material properties of the initial callus properly. Moreover, we assumed that the callus size was fixed after the initial phase of healing. This assumption is consistent with clinical observations, where the callus geometry develops during the initial phase of healing and is resorbed during the remodeling phase [5, 32, 38]. Similar to the model presented by Lacroix and Prendergast (2002) [30], we assumed that MSC migration is mainly governed by diffusion [57,58,59] and we considered the effects of fluid flow as a mechanical stimulus for MSC activities [21].

Conclusions

In conclusion, we have outlined the importance of the initial phase of healing, resulting in the formation of the initial callus with a range of geometry and material properties for optimal healing time. Findings from this work quantified the effects of the four important initial phase parameters on healing outcome. Consequently, there are well-established models to simulate soft callus formation, hard callus development, and remodeling phases of healing; however, one part is missing to complete the puzzle, and that is the initial phase of healing. This study emphasizes that the initial phase of healing should not be ignored in modeling of the healing process. Results from this study also raise questions about the clinical applications and the mechanisms of the initial healing phase such as how can we regulate these parameters at the initial healing phase to achieve the most efficient healing time? And how do micro-motions at fracture site, biological factors, and immune system response influence callus size and the level of granulation tissue formation at the initial phase of healing? As a future direction, a comprehensive model is required to simulate bone healing from the initial phase of healing to the end, considering both biology and mechanics. There are well-established models to simulate soft callus formation, hard callus development and remodeling phases of healing. However, they lack modeling of the initial phase of healing. This study illustrates the potential of addressing the initial phase of healing in a comprehensive simulation. Hence, further experimental investigations on the biological and mechanical factors in early stage of healing are required to develop more robust and predictive models that can simulate healing from the beginning to the end, and to better understand how clinicians can control and modulate the initial phase with its parameters.

Availability of data and materials

For this research, we developed finite element models in ABAQUS for stress analysis of fractured bone and mass diffusion of mesenchymal stem cells migration. They are all in ABAQUS input file format where handled by python code. All of the input files and python codes are attached in the supplementary material.

Abbreviations

BB:

Bony bridging

BC:

Bony callus

c:

Mesenchymal stem cells density

CC:

Cartilaginous callus

d:

Callus thickness

D:

Diffusion coefficient of mesenchymal stem cells

E:

Elastic modulus

Eg:

Elastic modulus of granulation tissue

FC:

Fibrous callus

FE:

Finite Element

h:

interfragmentary half gap size

MSC:

Mesenchymal Stem Cells

t:

time

References

  1. 1.

    Ghiasi MS, Chen J, Vaziri A, Rodriguez EK, Nazarian A. Bone fracture healing in mechanobiological modeling: a review of principles and methods. Bone Reports. 2017.

  2. 2.

    Claes L, Recknagel S, Ignatius A. Fracture healing under healthy and inflammatory conditions. Nat Rev Rheumatol. 2012;8(3):133–43.

  3. 3.

    Einhorn TA, Gerstenfeld LC. Fracture healing: mechanisms and interventions. Nat Rev Rheumatol. 2015;11(1):45–54.

  4. 4.

    Perren SM. Physical and biological aspects of fracture healing with special reference to internal fixation. Clin Orthop Relat Res. 1978;138:22.

  5. 5.

    McKibbin B. The biology of fracture healing in long bones. J Bone Joint Surg Br. 1978;60-B(2):150–62.

  6. 6.

    Friedenberg Z, FRENCH G. The effects of known compression forces on fracture healing. Surg Gynecol Obstet. 1952;94(6):743–8.

  7. 7.

    Morgan EF, Salisbury Palomares KT, Gleason RE, Bellin DL, Chien KB, Unnikrishnan GU, Leong PL. Correlations between local strains and tissue phenotypes in an experimental model of skeletal healing. J Biomech. 2010;43(12):2418–24.

  8. 8.

    Pivonka P, Dunstan CR. Role of mathematical modeling in bone fracture healing. BoneKEy Rep. 2012;1.

  9. 9.

    Ament C, Hofer E. A fuzzy logic model of fracture healing. J Biomech. 2000;33(8):961–8.

  10. 10.

    Carlier A, Geris L, Lammens J, Van Oosterwyck H. Bringing computational models of bone regeneration to the clinic. Wiley Interdiscip Rev Syst Biol Med. 2015;7(4):183–94.

  11. 11.

    Geris L, Reed AAC, Vander Sloten J, Simpson AHRW, Van Oosterwyck H. Occurrence and treatment of bone atrophic non-unions investigated by an integrative approach. PLoS Comput Biol. 2010;6(9):e1000915.

  12. 12.

    Geris L, Sloten JV, Oosterwyck HV. Connecting biology and mechanics in fracture healing: an integrated mathematical modeling framework for the study of nonunions. Biomech Model Mechanobiol. 2010;9(6):713–24.

  13. 13.

    Lacroix D, Prendergast PJ, Li G, Marsh D. Biomechanical model to simulate tissue differentiation and bone regeneration: application to fracture healing. Med Biol Eng Comput. 2002;40(1):14–21.

  14. 14.

    Bailón-Plaza A, van der Meulen MCH. Beneficial effects of moderate, early loading and adverse effects of delayed or excessive loading on bone healing. J Biomech. 2003;36(8):1069–77.

  15. 15.

    Carter DR, Beaupré GS, Giori NJ, Helms JA. Mechanobiology of skeletal regeneration. Clin Orthop Relat Res. 1998;355:S41–55.

  16. 16.

    Claes LE, Heigele CA. Magnitudes of local stress and strain along bony surfaces predict the course and type of fracture healing. J Biomech. 1999;32(3):255–66.

  17. 17.

    Prendergast PJ, Huiskes R, Søballe K. Biophysical stimuli on cells during tissue differentiation at implant interfaces. J Biomech. 1997;30(6):539–48.

  18. 18.

    Boerckel JD, Uhrig BA, Willett NJ, Huebsch N, Guldberg RE. Mechanical regulation of vascular growth and tissue regeneration in vivo. Proc Natl Acad Sci. 2011;108(37):E674–80.

  19. 19.

    Epari DR, Taylor WR, Heller MO, Duda GN. Mechanical conditions in the initial phase of bone healing. Clin Biomech. 2006;21(6):646–55.

  20. 20.

    Zhang L, Richardson M, Mendis P. Role of chemical and mechanical stimuli in mediating bone fracture healing. Clin Exp Pharmacol Physiol. 2012;39(8):706–10.

  21. 21.

    Isaksson H, Wilson W, van Donkelaar CC, Huiskes R, Ito K. Comparison of biophysical stimuli for mechano-regulation of tissue differentiation during fracture healing. J Biomech. 2006;39(8):1507–16.

  22. 22.

    Anderson DD, Thomas TP, Campos Marin A, Elkins JM, Lack WD, Lacroix D: Computational techniques for the assessment of fracture repair. Injury 2014, 45, Supplement 2:S23-S31.

  23. 23.

    Gómez-Benito MJ, García-Aznar JM, Kuiper JH, Doblaré M. Influence of fracture gap size on the pattern of long bone healing: a computational study. J Theor Biol. 2005;235(1):105–19.

  24. 24.

    Gómez-Benito MJ, García-Aznar JM, Kuiper JH, Doblaré M. A 3D computational simulation of fracture callus formation: influence of the stiffness of the external fixator. J Biomech Eng. 2005;128(3):290–9.

  25. 25.

    Isaksson H, Comas O, van Donkelaar CC, Mediavilla J, Wilson W, Huiskes R, Ito K. Bone regeneration during distraction osteogenesis: Mechano-regulation by shear strain and fluid velocity. J Biomech. 2007;40(9):2002–11.

  26. 26.

    Leong P, Morgan E. Measurement of fracture callus material properties via nanoindentation. Acta Biomater. 2008;4(5):1569–75.

  27. 27.

    Manjubala I, Liu Y, Epari DR, Roschger P, Schell H, Fratzl P, Duda G. Spatial and temporal variations of mechanical properties and mineral content of the external callus during bone healing. Bone. 2009;45(2):185–92.

  28. 28.

    Barry S. Non-steroidal anti-inflammatory drugs inhibit bone healing: a review. Vet Comp Orthop Traumatol. 2010;23(6):385.

  29. 29.

    Allen HL, Wase A, Bear W. Indomethacin and aspirin: effect of nonsteroidal anti-inflammatory agents on the rate of fracture repair in the rat. Acta Orthop Scand. 1980;51(1–6):595–600.

  30. 30.

    Lacroix D, Prendergast PJ. A mechano-regulation model for tissue differentiation during fracture healing: analysis of gap size and loading. J Biomech. 2002;35(9):1163–71.

  31. 31.

    Carlier A, van Gastel N, Geris L, Carmeliet G, Van Oosterwyck H. Size does matter: an integrative in vivo-in Silico approach for the treatment of critical size bone defects. PLoS Comput Biol. 2014;10(11):e1003888.

  32. 32.

    Augat P, Margevicius K, Simon J, Wolf S, Suger G, Claes L. Local tissue properties in bone healing: influence of size and stability of the osteotomy gap. J Orthop Res. 1998;16(4):475–81.

  33. 33.

    Claes LE, Heigele CA, Neidlinger-Wilke C, Kaspar D, Seidl W, Margevicius KJ, Augat P. Effects of mechanical factors on the fracture healing process. Clin Orthop Relat Res. 1998;355:S132–47.

  34. 34.

    Lu Y, Lekszycki T. Modelling of bone fracture healing: influence of gap size and angiogenesis into bioresorbable bone substitute. Mathematics and Mechanics of Solids. 2016;1081286516653272.

  35. 35.

    Ribeiro FO, Folgado J, Garcia-Aznar JM, Gómez-Benito MJ, Fernandes PR. Is the callus shape an optimal response to a mechanobiological stimulus? Med Eng Phys. 2014;36(11):1508–14.

  36. 36.

    Comiskey D, MacDonald B, McCartney W, Synnott K, O'Byrne J. Predicting the external formation of a bone fracture callus: an optimisation approach. Comput Methods Biomech Biomed Engin. 2012;15(7):779–85.

  37. 37.

    Garland T, Huey RB. Testing symmorphosis: does structure match functional requirements? Evolution. 1987;41(6):1404–9.

  38. 38.

    Claes L, Augat P, Suger G, Wilke H-J. Influence of size and stability of the osteotomy gap on the success of fracture healing. J Orthop Res. 1997;15(4):577–84.

  39. 39.

    Sapotnick A, Nackenhorst U: A Mechanically Stimulated Fracture Healing Model Using a Finite Element Framework. In: Biomedical Technology. Volume 74, edn. Edited by Lenarz T, Wriggers P. Berlin: Springer-Verlag Berlin; 2015: 41–53.

  40. 40.

    Miramini S, Zhang L, Richardson M, Pirpiris M, Mendis P, Oloyede K, Edwards G. Computational simulation of the early stage of bone healing under different configurations of locking compression plates. Comput Methods Biomech Biomed Engin. 2015;18(8):900–13.

  41. 41.

    Miramini S, Zhang L, Richardson M, Mendis P, Oloyede A, Ebeling P. The relationship between interfragmentary movement and cell differentiation in early fracture healing under locking plate fixation. Australas Phys Eng Sci Med. 2016;39(1):123–33.

  42. 42.

    Miramini S, Zhang L, Richardson M, Mendis P, Ebeling PR. Influence of fracture geometry on bone healing under locking plate fixations: a comparison between oblique and transverse tibial fractures. Med Eng Phys. 2016;38(10):1100–8.

  43. 43.

    Hori R, Lewis J. Mechanical properties of the fibrous tissue found at the bone-cement interface following total joint replacement. J Biomed Mater Res. 1982;16(6):911–27.

  44. 44.

    Repp F, Vetter A, Duda G, Weinkamer R. The connection between cellular mechanoregulation and tissue patterns during bone healing. Med Biol Eng Comput. 2015;53(9):829–42.

  45. 45.

    Wilson C, Schuetz M, Epari D. Effects of strain artefacts arising from a pre-defined callus domain in models of bone healing mechanobiology. Biomech Model Mechanobiol. 2015;14(5):1129–41.

  46. 46.

    Epari DR, Schell H, Bail HJ, Duda GN. Instability prolongs the chondral phase during bone healing in sheep. Bone. 2006;38(6):864–70.

  47. 47.

    BailÓN-Plaza A, Van Der Meulen MCH. A mathematical framework to study the effects of growth factor influences on fracture healing. J Theor Biol. 2001;212(2):191–209.

  48. 48.

    Peter C, Cook WO, Nunamaker DM, Provost MT, Seedor JG, Rodan GA. Effect of alendronate on fracture healing and bone remodeling in dogs. J Orthop Res. 1996;14(1):74–9.

  49. 49.

    Iwaki A, Jingushi S, Oda Y, Izumi T, Shida JI, Tsuneyoshi M, Sugioka Y. Localization and quantification of proliferating cells during rat fracture repair: detection of proliferating cell nuclear antigen by immunohistochemistry. J Bone Miner Res. 1997;12(1):96–102.

  50. 50.

    Goodship A, Cunningham J: Pathophysiology of functional adaptation of bone in remodelling and repair in-vivo. 2001.

  51. 51.

    Kim W, Kohles SS. Optical acquisition and polar decomposition of the full-field deformation gradient tensor within a fracture callus. J Biomech. 2009;42(13):2026–32.

  52. 52.

    De Bastiani G, Aldegheri R, Renzi Brivio L. The treatment of fractures with a dynamic axial fixator. J Bone Joint Surg Br. 1984;66(4):538–45.

  53. 53.

    Wang F-S, Yang K, Kuo Y-R, Wang C-J, Sheen-Chen S-M, Huang H-C, Chen Y-J. Temporal and spatial expression of bone morphogenetic proteins in extracorporeal shock wave-promoted healing of segmental defect. Bone. 2003;32(4):387–96.

  54. 54.

    Augat P, Merk J, Genant H, Claes L. Quantitative assessment of experimental fracture repair by peripheral computed tomography. Calcif Tissue Int. 1997;60(2):194–9.

  55. 55.

    Den Boer F, Bramer J, Patka P, Bakker F, Barentsen R, Feilzer A, De Lange E, Haarman H. Quantification of fracture healing with three-dimensional computed tomography. Arch Orthop Trauma Surg. 1998;117(6–7):345–50.

  56. 56.

    Cowin SC: Bone mechanics handbook: CRC press; 2001.

  57. 57.

    Web-based Injury Statistics Query and Reporting System (WISQARS) [http://www.cdc.gov/injury/wisqars/].

  58. 58.

    Ghimire S, Miramini S, Richardson M, Mendis P, Zhang L. Effects of dynamic loading on fracture healing under different locking compression plate configurations: a finite element study. J Mech Behav Biomed Mater. 2019.

  59. 59.

    Ganadhiepan G, Miramini S, Patel M, Mendis P, Zhang L. Bone fracture healing under Ilizarov fixator: influence of fixator configuration, fracture geometry and loading. Int J Numer Method Biomed Eng. 2019:e3199.

Download references

Acknowledgments

The authors would like to acknowledge Professor Hashemi from Northeastern University for fruitful discussions on this work.

Funding

This work was supported in part by a grant from the United States National Science Foundation’s Civil, Mechanical, and Manufacturing Innovation (Grant No. 1634560), and by the Department of Mechanical and Industrial Engineering at Northeastern University. The funding source had no role in study design, data analysis, or data interpretation.

Author information

MG conducted the computational modeling, helped with the conceptualization, and wrote the manuscript. AN, AV, and ER contributed to conceptualization, design of the study, and review of the manuscript. JC helped with the analysis of the results and writing the manuscript. All authors have read and approved the final version of this manuscript.

Correspondence to Ashkan Vaziri or Ara Nazarian.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare no other competing interests. Ara Nazarian is a Associate Editor for BMC Musculoskeletal Disorders.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Additional file 1. 2D-9-d = 1.py: a python script to handle models when callus thickness was 1 mm. (PY 21kb)

Additional file 2. 2D-9-d = 2.py: a python script to handle models when callus thickness was 2 mm. (PY 21kb)

Additional file 3. 2D-9-d = 3.py: a python script to handle models when callus thickness was 3 mm. (PY 21kb)

Additional file 4. 2D-9-d = 4.py: a python script to handle models when callus thickness was 4 mm. (PY 21kb)

Additional file 5. 2D-9-d = 5.py: a python script to handle models when callus thickness was 5 mm. (PY 21kb)

Additional file 6. 2D-9-d = 6.py: a python script to handle models when callus thickness was 6 mm. (PY 21kb)

Additional file 7. 2D-9-d = 7.py: a python script to handle models when callus thickness was 7 mm. (PY 21kb)

Additional file 8. 2D-9-d = 8.py: a python script to handle models when callus thickness was 8 mm. (PY 21kb)

Additional file 9. Job-2D-stress-d = 1–0.inp: an ABAQUS input file for stress-strain analysis when callus thickness was 1 mm. (INP 306kb)

Additional file 10. Job-2D-stress-d = 2–0.inp: an ABAQUS input file for stress-strain analysis when callus thickness was 2 mm. (INP 416kb)

Additional file 11. Job-2D-stress-d = 3–0.inp: an ABAQUS input file for stress-strain analysis when callus thickness was 3 mm. (INP 573kb)

Additional file 12. Job-2D-stress-d = 4–0.inp: an ABAQUS input file for stress-strain analysis when callus thickness was 4 mm. (INP 747kb)

Additional file 13. Job-2D-stress-d = 5–0.inp: an ABAQUS input file for stress-strain analysis when callus thickness was 5 mm. (INP 964kb)

Additional file 14. Job-2D-stress-d = 6–0.inp: an ABAQUS input file for stress-strain analysis when callus thickness was 6 mm. (INP 1200kb)

Additional file 15. Job-2D-stress-d = 7–0.inp: an ABAQUS input file for stress-strain analysis when callus thickness was 7 mm. (INP 1450kb)

Additional file 16. Job-2D-stress-d = 8–0.inp: an ABAQUS input file for stress-strain analysis when callus thickness was 8 mm. (INP 1640kb)

Additional file 17. Job-2D-diffusion-d = 1–0.inp: an ABAQUS input file for mass diffusion analysis when callus thickness was 1 mm. (INP 280kb)

Additional file 18. Job-2D-diffusion-d = 2–0.inp: an ABAQUS input file for mass diffusion analysis when callus thickness was 2 mm. (INP 383kb)

Additional file 19. Job-2D-diffusion-d = 3–0.inp: an ABAQUS input file for mass diffusion analysis when callus thickness was 3 mm. (INP 529kb)

Additional file 20. Job-2D-diffusion-d = 4–0.inp: an ABAQUS input file for mass diffusion analysis when callus thickness was 4 mm. (INP 695kb)

Additional file 21. Job-2D-diffusion-d = 5–0.inp: an ABAQUS input file for mass diffusion analysis when callus thickness was 5 mm. (INP 898kb)

Additional file 22. Job-2D-diffusion-d = 6–0.inp: an ABAQUS input file for mass diffusion analysis when callus thickness was 6 mm. (INP 1120kb)

Additional file 23. Job-2D-diffusion-d = 7–0.inp: an ABAQUS input file for mass diffusion analysis when callus thickness was 7 mm. (INP 1360kb)

Additional file 24. Job-2D-diffusion-d = 8–0.inp: an ABAQUS input file for mass diffusion analysis when callus thickness was 8 mm. (INP 1540kb)

Additional file 25. Job-2D-diffusion-d = 1-E.inp: an ABAQUS input file for final results presentation when callus thickness was 1 mm. (INP 279kb)

Additional file 26. Job-2D-diffusion-d = 2-E.inp: an ABAQUS input file for final results presentation when callus thickness was 2 mm. (INP 383kb)

Additional file 27. Job-2D-diffusion-d = 3-E.inp: an ABAQUS input file for final results presentation when callus thickness was 3 mm. (INP 528kb)

Additional file 28. Job-2D-diffusion-d = 4-E.inp: an ABAQUS input file for final results presentation when callus thickness was 4 mm. (INP 694kb)

Additional file 29. Job-2D-diffusion-d = 5-E.inp: an ABAQUS input file for final results presentation when callus thickness was 5 mm. (INP 898kb)

Additional file 30. Job-2D-diffusion-d = 6-E.inp: an ABAQUS input file for final results presentation when callus thickness was 6 mm. (INP 1120kb)

Additional file 31. Job-2D-diffusion-d = 7-E.inp: an ABAQUS input file for final results presentation when callus thickness was 7 mm. (INP 1360kb)

Additional file 32. Job-2D-diffusion-d = 8-E.inp: an ABAQUS input file for final results presentation when callus thickness was 8 mm. (INP 1540kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Ghiasi, M.S., Chen, J.E., Rodriguez, E.K. et al. Computational modeling of human bone fracture healing affected by different conditions of initial healing stage. BMC Musculoskelet Disord 20, 562 (2019). https://doi.org/10.1186/s12891-019-2854-z

Download citation

Keywords

  • Bone fracture healing
  • Inflammatory stage
  • Initial callus size
  • Granulation tissue material properties
  • Migration rate
  • Mechanobiological modeling
  • Finite element analysis