journal of the mechanical behavior of biomedical materials 29 (2014) 328 –338
Available online at www.sciencedirect.com
www.elsevier.com/locate/jmbbm
Research Paper
An interface finite element model can be used to predict healing outcome of bone fractures J.A. Aliertaa, M.A. Pe´rezb, J.M. Garcı´a-Aznarb,n a
Escuela Politécnica Superior del Ejército, Ministry of Defense, Madrid, Spain Mechanical Engineering, Aragón Institute of Engineering Research (I3A), University of Zaragoza, Ed. Betancourt-Campus Río Ebro, C/María de Luna, 50018 Zaragoza, Spain b
ar t ic l e in f o
abs tra ct
Article history:
After fractures, bone can experience different potential outcomes: successful bone con-
Received 16 May 2013
solidation, non-union and bone failure. Although, there are a lot of factors that influence
Received in revised form
fracture healing, experimental studies have shown that the interfragmentary movement
20 September 2013
(IFM) is one of the main regulators for the course of bone healing. In this sense,
Accepted 23 September 2013
computational models may help to improve the development of mechanical-based
Available online 8 October 2013
treatments for bone fracture healing. Hence, based on this fact, we propose a combined
Keywords:
repair-failure mechanistic computational model to describe bone fracture healing. Despite
Bone fracture healing
being a simple model, it is able to correctly estimate the time course evolution of the IFM
Finite element prediction
compared to in vivo measurements under different mechanical conditions. Therefore, this
Design of fracture fixators
mathematical approach is especially suitable for modeling the healing response of bone to
Interfragmentary movement
fractures treated with different mechanical fixators, simulating realistic clinical conditions. This model will be a useful tool to identify factors and define targets for patient specific therapeutics interventions. & 2013 Elsevier Ltd. All rights reserved.
1.
Introduction
Fractures are a common orthopedic problem and the fracture healing is a natural process that regenerates bone to its original state and function. Several factors influence these bone healing events, such as genetic, cellular and biochemical factors, blood supply, neural and hormonal regulation, age, the type of fracture interfragmentary motion and fracture geometry (Einhorn, 2005; Goodship et al.,1993; Hadjiargyrou et al., 1998; Jagodzinski and Krettek, 2007; Marsell and Einhorn, 2011). However, the most common orthopedic treatments consist on the mechanical stabilization of the bone fracture gap, regulating the interfragmentary movement. n
Corresponding author. Tel.: þ34 976 76 27 96; fax: þ34 976 76 26 70. E-mail address: [email protected] (J.M. García-Aznar).
1751-6161/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jmbbm.2013.09.023
The role of the mechanical stabilization on the design of fracture fixators has evolved very much in the last years, updating from the concept of “open reduction and internal fixation (ORIF)” to “biological fixation” (Perren, 2002). In both cases, differences are presented from a biological and mechanical point of view. In fact, the concept of “biological fixation” is based on the application of the fixator as a minimally invasive percutaneous osteosynthesis (MIPO). So, the contact of the implant with bone is reduced at maximum, avoiding the damage to the blood supply. In addition to biological differences, there is also a different concept in the role of the mechanical stabilization. In the “biological fixation” primary ossification is avoided, promoting the occurrence of a secondary
journal of the mechanical behavior of biomedical materials 29 (2014) 328 –338
ossification mechanism that induces the formation of bone callus (Manjubala et al., 2009; Vetter et al., 2010). Therefore, the mechanical design concept is updated from an absolute stability to a controlled mechanical instability that favors the formation of bone callus. Hence, the mechanical stabilization due to the stiffness of the fixator has acquired a relevant role for achieving a successful bone healing (Bishop et al., 2003; Chao et al., 1989; Draper et al., 1997; Goodship et al., 1993; Richardson et al., 1994; Wehner et al., 2011). Many studies have examined the role of different local mechanical conditions on bone fracture healing: the influence of the magnitude of interfragmentary movement (IFM) and the initial size of the fracture gap (Claes et al., 1997; Gómez-Benito et al., 2011; Goodship and Kenwright, 1985), the stiffness of the fixation (Schell et al., 2005) and the type of movement of the gap (Augat et al., 2003; Bishop et al., 2006). Computational models of the fracture healing process may prove useful in the determination of optimal mechanical treatments. In fact, several mechano-biological models based on finite element simulations studied the influence of local mechanical conditions on biological events that regulate the temporal and spatial evolution of the different ossification mechanisms that occur during healing (Andreykiv et al., 2008; Checa and Prendergast, 2009; Isaksson et al., 2008, 2009a; Lacroix and Prendergast, 2002b; Loboa et al., 2001; ReinaRomo et al., 2011; Shefelbine et al., 2005; Simon et al., 2011; Wehner et al., 2010). Most of these models have focused on how different mechanical stimuli (such as: fluid flow, octahedral shear strain, deviatoric strain, strain energy density, etc.) are able to predict spatial pattern distribution of tissues regulated by intramembranous and endochondral ossification (Checa and Prendergast, 2009; Lacroix and Prendergast, 2002a; Wehner et al., 2010). Other models have combined tissue differentiation rules with callus growth (Comiskey et al., 2013; García-Aznar et al., 2007; Gómez-Benito et al., 2005) or callus resorption (Ament and Hofer, 2000; Byrne et al., 2011; Isaksson et al., 2009b; Lacroix and Prendergast, 2002b) depending on the temporal stage of bone healing that they studied. Although these mathematical models are very useful to understand the fundamental cellular mechanisms that locally regulate tissue differentiation and callus growth/ resorption, they are not really helpful for the simulation of realistic bone fractures, where a whole-organ analysis is required. In fact, these models have only analyzed simple tranversal or obliques (Comiskey et al., 2013; Loboa et al., 2001) bone fractures with two fragments. However, fractures are much more complexes with very different shapes of the fracture line, complicated anatomical locations and a different number of fragments. Therefore, the main aim of the present study is to propose a phenomenological computational model able to simulate the temporal recovery of mechanical properties of the fracture zone during the healing process, which is regulated by the mechanical stability. If this approach proves feasible, it offers the possibility of using computer simulations in the clinical treatment of complex fractures with multiple fragments, complicated geometries and different anatomical locations and in other orthopedic applications where bone regeneration occurs.
2.
329
Material and methods
The bone fracture gap was modeled through the incorporation of interface elements that connected the two fracture ends simulating the discontinuity in the displacement field between fragments. A mathematical model is here proposed to simulate the temporal evolution of the separation between both fragments, regulated by the mechanical conditions existent in the fracture gap.
2.1.
Fracture gap/interface mechanical behavior
To model the fracture gap behavior, 6-nodes and 8-nodes cohesive elements (Fig. 1) are used to connect the fracture ends (García-Aznar et al., 2009). The thickness of these elements, which is the dimension of the gap fracture, is thin enough to consider it negligible with respect to the overall dimensions of the bone fracture. Accordingly, the behavior of these elements is directly established in terms of one traction versus one separation law. The nominal traction stress vector, t, consists of three components: tn, ts, tt, which represents the normal (tn) (along
Fig. 1 – Constitutive model: (a) cohesive elements, (b) shear traction and (c) normal traction.
330
journal of the mechanical behavior of biomedical materials 29 (2014) 328 –338
the local 3-direction) and the two shear tractions (ts, tt) (along the local 1- and 2-directions), respectively (Fig. 1). The corresponding separations are denoted by δn, δs and δt. Thus, the nominal strains can be defined as εn ¼
δn δs δt ; εs ¼ ; εt ¼ h h h
ð1Þ
where h is the original thickness of the cohesive element (fracture gap thickness). The purpose of the present study is to simulate the temporal evolution of the bone healing process. Therefore, to quantify the degree of healing, union degree, we use α as the main state variable. This variable is normalized, corresponding α¼ 1 to a totally successful bone healing and α¼0 to a completely non-union or malunion. The union degree α, may be also related with the so-called damage variable, d (Moreo et al., 2007), through the following expression: d ¼ 1 α
ð2Þ
The net effect of the bone healing process on the macroscopic interface mechanical properties is what it is represented by the growth of α. In addition, the model allows the possibility of reducing the union degree α, which means mechanical failure of the gap, because it is not able to bear the local demands. For each element of the fracture gap, the union degree, α, can be defined as K0i α ¼ Ki
i ¼ n; s; t
ð3Þ
where K0i defines the linear stiffness in a completely healed fracture gap along the direction i (the maximum fracture gap stiffness available), and Ki is the current interface stiffness in the i direction (Fig. 1b and c). The three components of the nominal traction stress vector, t, can be calculated, using the above defined variables, through the following equation: ti ¼ K0i αεi
2.2.
i ¼ n; s; t
ð4Þ
Mechanical damage: α decrease
If mechanical demands in the fracture gap are important, a reduction of the degree of union may occur as consequence of its mechanical failure. For the shear traction, we assume the degree of union is reduced by the following law (Fig. 1b): ( ðε ε Þ i ci ðε0i εci Þ ; if εαi oεi oεci i ¼ s; t ð5Þ αi ¼ 0; if εi Zεci where, ε0i maximum shear strain at the linear region for α¼ 1 (perfectly healed interface), εαi maximum shear strain at the linear region for each value of α and εci the maximum shear strains allowed. If εci are exceeded, the fracture gap would have no stiffness at all. For the compression traction (εn o0) we assume a linear behavior with respect to the associated strain until εn reaches the value 1, which means that the fracture ends are in contact. If lower values of εn are reached the simulation will be aborted because no penetration is allowed.
For the tensile traction (εn Z0) we assume an initially linear behavior with respect to the associated normal strain, followed by a linear decay (Fig. 1c): ( ðε ε Þ n cn ; if εαn oεn oεcn ð6Þ αn ¼ ðε0n εcn Þ 0; if εn Zεcn where ε0n, εαn and εcn have an equivalent meaning that those defined for shear traction. In Fig. 1b and c the dashed line represents the behavior law (X¼ε (strain); Y ¼t (stress)) of a completely healed bone gap (α¼1). The continuous line represents the behavior law (t–ε) of a particular bone gap through the healing process (αo1). In these figures, the parts of the graphs with negative slope represent the Eq. (5) (Fig. 1b) and 6 (Fig. 1c). In this part of the graph, α decreases (α_ o0) because the limit strain (εα) is reached (for example, due to a very high load). If this fact occurs during the healing process, the healing will be impaired or at least delayed. Therefore, we finally select as union degree, α, the minimum value associated to each direction: α ¼ min fαs ; αt ; αn g
2.3.
ð7Þ
Biological union: growth of union degree α
The union degree, α, can increase due to bone healing, and therefore it is calculated as the addition of two contributions: α ¼ αc þ αb
ð8Þ
where αc reflects the recovery of the mechanical properties of the fracture gap mainly due to the formation of cartilage, and αb reflects the recovery of the mechanical properties of the fracture gap due to bone formation. In order to estimate the temporal evolution of α ðα_ ¼ α_ c þ α_ b Þ, we propose a regulatory law based on three different healing zones (adapted from Claes and Heigele, 1999) which are established as of ffi the qffiffiffiffiffiffiffiffiffiffiffiffiffiffi a function compression, εn, of each interface and the shear strains εshear ¼ ε2s þ ε2t element (Fig. 2a) with the following characteristics: Zone 1 non-union: The compression and shear strains are so high that there is a high interfragmentary movement (IFM) which impairs a successful fracture healing. Only fibrous tissue is formed and there is no change in αc nor in αb ðα_ c ¼ α_ b ¼ 0Þ Zone 2 cartilage formation: The compression and shear strains are moderate, in the same way that the IFM. We consider that in those elements, whose strains belong to this zone, cartilage will be formed and therefore αc will grow. However, these elements should be under this strain threshold during a certain period of time (Mc) before the growing of αc, this time represents the maturation time needed for the cartilage to appear close to the fracture gap. It is established that αc will have an exponential evolution with the time (Fig. 2b), therefore once the maturation time for cartilage (Mc) is reached, the value of αc will grow exponentially, providing that the strains keep in zone 2, during a fixed period of time, tc (cartilage formation time), until it reaches the maximum value (αcmax). If during the healing process, αc is
journal of the mechanical behavior of biomedical materials 29 (2014) 328 –338
growing and the strains reach the zone 3, αc will continue growing with the same rate until the maximum value (αcmax). Zone 3 bone formation: In this zone strains and IFM are lower than in zone 2, as a result, bone is formed, so αb will grow. We also define, Mb, as the maturation time for the formation of bone. It is established that αb will have a linear evolution with the time (Fig. 2c), therefore once the maturation time for cartilage (Mb) is reached, the value of αb will
331
grow linearly, providing that the strains keep in the zone 3, until it reaches the maximum value (αbmax) at the total healing time, th. If αb begins to grow because the conditions for bone formation are achieved and αc is still growing, a simultaneous growing of both contributions will happen. According to the definition of α, its maximum value is 1 (corresponding to a totally successful bone healing), therefore the limits for αc (αcmax) and αb (αbmax) must satisfy the following equation: αmax ¼ αb max þ αc max ¼ 1
ð9Þ
sh The model parameters values represented in Fig. 2 (Lsh c , Lb , comp comp Lc , Lb , αcmax, αbmax, Mc, Mb, tc, th) will be later defined by tuning the model to experimental results. The model was numerically implemented in an Abaqus user subroutine (UMAT) and all the finite element analyses (FEA) were carried out in ABAQUS v.6.11.
2.4.
Examples of application for the calibration process
In order to determine the value of the parameters that define the model, we made a quantitative comparison between the numerical results (obtained from our model) and the experimental ones found in the literature. With this idea in mind, we simulated two animal experiments (compression and shear) from the literature that are detailed below and summarized in Table 1.
2.4.1.
Fig. 2 – (a) Healing zones, (b) temporal evolution of αc and (c) temporal evolution of αb.
Compression experiments
To estimate the compression parameters (Lcomp , Lcomp , αcmax, c b αbmax, Mc, Mb, tc, th), we simulated the compression experiments performed by Claes et al. (1997). They determined the influence of the size of the osteotomy gap and interfragmentary strain on callus formation in the healing of the sheep metatarsus (Claes et al., 1997). Therefore, a cylinder of external diameter of 17 mm was simulated (Fig. 3a). The fracture gap was subjected to an external compression load of 330 N (Claes et al., 1995). The gap size was 2 mm and the initial interfragmentary strain (IFS¼IFM/h) was approximately 31%, that was achieved with the use of a fixator, simulated by four linear springs with a total stiffness of 348 N/mm. In order to avoid rigidbody motions the nodes at the base of the model were fully constrained. Bone was considered as a linear elastic material for all analyses where the mechanical properties are shown in Table 2.
Table 1 – Examples of model application. Goal
Experiment
Bone
Load condition
Reference
Calibration process
Compression Shear
Sheep metatarsus Sheep tibia
Compression Shear
Claes et al. (1997) Bishop et al. (2006)
Evaluation of the model predictive capacity
Compression Combined load case Comminuted fracture
Sheep metatarsus Human tibia Sheep metatarsus
Compression Combined Combined
Claes et al. (1997) Wehner et al. (2010) –
332
journal of the mechanical behavior of biomedical materials 29 (2014) 328 –338
Fig. 3 – FE model for the simulation of (a) compression and (b) shear experiments in order to calibrate the corresponding parameters of the model.
Table 2 – Material properties of the tissue considered for all analyses.
Table 3 – Fixation stiffness for each experiment simulated in the compression predictive simulation.
Tissue type
Elastic modulus in MPa
Poisson's ratio
Experiment
Stiffness (N/mm)
Cortical bone Trabecular bone
10,000a 200a
0.2a 0.2a
Initial gap ¼1 mm; IFS ¼ 7% Initial gap ¼1 mm; IFS ¼ 31% Initial gap ¼2 mm; IFS ¼ 7%
3300 764 1650
a
Hernandez et al. (2001).
2.4.2.
Shear experiments
In a similar way, it is necessary that the parameters of the model allow also to simulate correctly the healing process under shear conditions. However, no experiments in the literature have been found in which, the time course of the interfragmentary movement (IFM) during the healing is evaluated under a pure shear load case. Nevertheless Bishop et al. (2006) performed experiments studying the healing of a transverse tibial osteotomy in sheep where interfragmentary torsional shear (torsion) was applied under carefully controlled conditions and some data of this experiment are used for the estimation of the rest of the parameters of our sh model (Lsh c , Lb ). We simulated these experiments taking the geometry and the loads from Bishop et al. (2006), (Fig. 3b). A cylinder of external diameter of 20 mm was simulated. A torsional rotation (torsion) was applied over a 2.4 mm osteotomy with a maximum interfragmentary principal strain magnitude of 25%. This corresponded to a 7.21 rotation calculated in the outer cortical radius of 10 mm. A load limit of 1670 N mm was imposed to allow interfragmentary strain to decrease with increasing callus stiffness. In order to avoid rigid-body motions the nodes at the base of the model were constrained.
To calibrate the model two data sets in the related article (Bishop et al., 2006) were used: – The stiffness increased steadily during the experiment with the load limit achieved within the fourth week. – After eight weeks of healing, the bending stiffness was approximately the 80% of the intact bone. The union degree α, is the variable that represents the improvement of the stiffness and then is directly related to the bending stiffness.
2.5. Examples of application for the evaluation of the model predictive capacity In order to evaluate the predictive capacity of the model, three additional experiments (compression, combined load case and comminuted fracture) that are summarized in Table 1 were simulated using the parameters previously obtained in the calibration process.
2.5.1.
Compression experiments
Three compression experiments, similar to the described in the model calibration section (Fig. 3a), were performed by
journal of the mechanical behavior of biomedical materials 29 (2014) 328 –338
Claes et al. (1997). The differences among them were the gap dimension and the initial IFS. In fact, gaps of 1 mm and 2 mm and initial IFS of 7% and 31% were fixed in these experiments. We simulated these three experiments with the same FE model used in the model calibration section, changing the gap dimension and the stiffness of the fixator (Table 3). The same loading and boundary conditions were used.
2.5.2.
Combined load case experiment
We found no more data in the literature of pure shear experiments of fracture healing than those used in the calibration (Bishop et al., 2006). That fact prevented us from validating our model exclusively with this type of loading condition. However, Wehner et al. (2010) performed a clinical experiment in a human patient with a fracture in the diaphyseal tibia. A combined compression-shear load case was applied, and at several time points during healing, the IFM was determined during a single leg stance. Thanks to these available data, simulating this experiment we were able to evaluate the potential of our model under a complex combined load case. A human tibia was scanned and the images stored in Dicom format. MIMICS V. 15.0 was used to reconstruct the geometry (Fig. 4a) through the segmentation process. Afterwards, using this software, an osteotomy was made in the diaphyseal part. In addition, a unilateral external fixator was modeled according to the specifications shown in Wehner et al. (2010). This fixator, once imported in MIMICS, was implanted in the tibia, and the assembly was meshed with tetrahedral elements. We used ABAQUS to create cohesive elements in the fracture gap. The final FE model had a total number of nodes of 19,499 and a total number of elements of 91,558 (118 linear wedge cohesive elements and 91,440 linear tetrahedral elements). The element size used was inside the asymptotic region of convergence and represents a good trade off between numerical accuracy and computational cost (Fig. 4c).
333
To simulate the experiment, we applied a vertical load of 700N (corresponding to a single leg stance of a patient of approximately 70 kg (Wehner et al.,2010)) at the proximal condyles of the tibia (Fig. 4b). Due to this load, a combined load state (compression-shear) was generated in the fracture gap zone. The nodes of the base were constrained so that, the rigid-body motions are avoided.
2.5.3.
Comminuted fracture
A complex fracture was simplified and simulated in a sheep metatarsus. The comminuted fracture modeled consisted of an oblique fracture with three fragments (Fig. 5). The loads, fixation system and bone geometry were the same that those used in the Section 2.4.1. Due to the complex geometry of the fracture zone the vertical load applied generated combined mechanical demands in the different fracture gaps. Four different simulations were made changing the stiffness of the external fixator in order to get the best fixation
Fig. 5 – FE model for the simulation of the comminuted fracture.
1.6 Simulated
1.4
Experimental (Claes,1997)
IFM (mm)
1.2 1 0.8 0.6 0.4 0.2 0 0
10
20
30
40
50
60
Healing time (days)
Fig. 4 – Predictive simulation using combined load case: (a) geometry, (b) applied loads and (c) FE of the fracture gap zone.
Fig. 6 – Compression calibration results: evolution of the interfragmentary motion in the computational simulation compared with the experimental data in a 2-mm gap of a sheep (Claes et al., 1997).
334
journal of the mechanical behavior of biomedical materials 29 (2014) 328 –338
Table 4 – Model parameters. Parameter
Value
Mc—Maturation time for cartilage (days) (see Fig. 2b) Mb—Maturation time for bone (days) (see Fig. 2c) αcmax—Maximum value of αc (see Fig. 2b) αbmax—Maximum value of αb (see Fig. 2c) tc—Cartilage formation time (days) (see Fig. 2b) th—Total healing time (days) (see Fig. 2c) Lcomp —Compression strain limit for cartilage formation (see Fig. 2a) c
4a 10b 0.3b 0.7b 36b 60c 0.6b
comp
Lb
0.25b
—Compression strain limit for bone formation (see Fig. 2a)
Lsh c —Shear Lsh b —Shear
0.25b
strain limit for cartilage formation (see Fig. 2a)
0.15b
strain limit for bone formation (see Fig. 2a) K0i (i ¼n,s,t)—Initial linear stiffness (MPa) ε0i (i¼ n,s,t)—Maximum strain at the linear region εci (i¼n,s,t)—Maximum allowed strain
50d 0.3b 1b
a
Park et al. (1998). Estimated. c Klein et al. (2004). d Bishop et al. (2003). b
system for this fracture. The stiffness of the fixator used for each simulation was 400, 448, 500 and 560 N/mm.
3.
Results
Here, in first step, we discuss the results of the calibration process, and next, we also explain additional experiments to evaluate the predictive capacity of the model.
3.1.
Calibration process results
3.1.1.
Compression experiments results
Evolution of the IFM obtained experimentally (Claes et al., 1997) and computationally has been represented in Fig. 6 for a 2-mm gap of a sheep. The experimental variability has been also represented and the computational results are within previous range of variation. Initially, there is no fracture healing, because the IFM is relatively high (1 mm), as the simulation evolves, the IFM is reduced. After 40 days of healing, the IFM has been reduced by a ratio of 80% (see Fig. 6).
3.1.2.
Shear experiments results
The temporal evolution of the union degree, α, of the fracture gap has been estimated under a torsional load case. We have been able to correctly reproduce the experimental results obtained by Bishop et al. (2006). Experimentally the torsion load limit of 1670N mm is achieved within the fourth week, using our model, we estimate1 that this load limit is achieved the 32nd day.
3.1.3.
Parameters of the model
Using the results obtained in the calibration process, the parameters that best fit the simulated and the experimental 1 We determine the torsion load limit as the resultant torsion moment, which is computed by means of the integration of the shear traction stress vector over the bone section.
results are shown in Table 4. Some of the model parameters come from previous studies of the literature and the rest are estimated with the results of Sections 3.1.1and 3.1.2.
3.2.
Evaluation of the model predictive capacity—results
3.2.1.
Compression experiments results
The experimental results obtained by Claes et al. (1997) for different fracture gap sizes and IFS have been represented in Fig. 7, where we can also see the experimental range of variation. In addition, we include in this figure the numerical results corresponding to the computational simulations. As we can see, a good agreement is found between numerical and experimental results for all the studied conditions.
3.2.2.
Combined load case experiment results
Experimental results from Wehner et al. (2010) obtained for a combined load case have been represented and compared with the computational predicted ones in Fig. 8. Although there were no more experimental data in order to represent its range of variation, our computational results are very close to the experimental ones. The IFM compression is very high at the initial healing stages, and as the healing process evolves the IFM is reduced both experimentally and computationally due to the combined load case. However, in the case of shear loads, the shear IFM slightly decreases from the first days until the last days of healing (Fig. 8).
3.2.3.
Comminuted fracture results
The temporal evolution of the vertical displacement of the fixation points has been represented for four different fixators (Fig. 9). It can be observed that the stiffness value that leads to better results is K ¼ 500 N/mm, because the final interfragmentary displacements are lower, which means better healing conditions. We have also represented the temporal evolution of the bone reaction (Fig. 10) i.e. the load that the bone is able to
335
journal of the mechanical behavior of biomedical materials 29 (2014) 328 –338
0.6 0.25 0.5
Simulated
0.4
Experimental (Claes,1997) Gap=2 mm; Initial IFS=7%
Simulated Experimental (Claes,1997) Gap=1mm; Initial IFS=7%
0.15
IFM (mm)
IFM (mm)
0.2
0.1
0.3 0.2
0.05
0.1
0
0 0
10
20
30
40
50
60
0
10
20
Healing time (days)
30
40
50
60
Healing time (days)
0.3 Simulated 0.25
Experimental (Claes,1997) Gap=2 mm; Initial IFS=7%
IFM (mm)
0.2
0.15
0.1
0.05
0 0
10
20
30
40
50
60
Healing time (days)
Fig. 7 – Evolution of the interfragmentary motion in the computational simulation compared with the experimental data (Claes et al., 1997): (a) gap¼ 1 mm; initial IFS¼ 7%, (b) gap¼ 1mm; initial IFS ¼31%, and (c) gap¼ 2 mm; initial IFS¼7%. 1.6
Experimental IFM compression (Wehner et al., 2010)
1.4
Experimental IFM shear (Wehner et al., 2010) Simulated IFM compression
IFM (mm)
1.2 1
Simulated IFM shear
0.8
If we compare the different stiffness able to heal a simple transverse fracture (K¼ 348 N/mm see Section 2.4.1) and a comminuted fracture (K¼ 500 N/mm) under the same loading conditions, we can observe that the second one requires a higher stiffness to achieve successful healing.
0.6 0.4
4.
Discussion
0.2 0 0
10
20
30
40
50
60
70
Healing time (days)
Fig. 8 – Predicted time course of the interfragmentary movement (IFM) due to compression and shear stresses compared to in vivo measurements of the patients (Wehner et al., 2010).
bear during the healing process. The best result is for a stiffness of K ¼500 N/mm, however it can be seen that K ¼560 N/mm (more rigid fixation) gives better results in the first healing stages probably due to the primary bone formation. The more flexible fixation systems give poorer results probably because the strains are too high and the healing process is not successful.
In this study we propose a coupled computational model able to simulate the normal and the impaired (mechanical failure) bone healing regulated by the mechanical stability or the interfragmentary movement in the fracture gap. Indeed, the model predicted the temporal recovery/failure of mechanical properties of the fracture gap during the healing. The aim of this model is to identify the role of different and combined degrees of stabilization for fracture fixation. The model here presented is based on the main mechanisms that regulate the biology of fracture healing. In fact, bone heals by either direct intramembranous ossification or indirect fracture healing, which consists of both intramembranous and endochondral ossification (Marsell and Einhorn, 2011). Therefore, this model combines both approaches, modeling both direct and indirect ossification or equivalently intramembranous and endochondral ossification.
336
journal of the mechanical behavior of biomedical materials 29 (2014) 328 –338
Fig. 9 – Comminuted fracture results: temporal evolution of the vertical displacement of the fixation points (see position of these points in Fig. 5) for different fixation stiffness.
Fig. 10 – Temporal evolution of bone reaction (N) for different values of stiffness of the fixator. Primary healing (intramembranous ossification) requires extremely rigid fixation systems and negligible spaces between bone fragments. The quality of regenerated bone in primary healing cannot be always guaranteed and therefore a new fracture may occur. In those parts of the fracture gap in which, only primary healing occurs, we assume that there will be only the temporal evolution of one part (αb) of the union degree (α) (contribution of α mainly due to the bone formation). On contrast, secondary healing (endochondral ossification) takes place under more flexible fixation that allows
a higher, but controlled, instability of the fracture gap region. In this second step, bone formation is involved restoring the mechanical stability and recovering the mechanical properties in the bone fracture gap. In any case, secondary healing is the most common in clinical treatments, since the high quality of the regenerated tissue and also the probability of a successful healing is higher than in primary healing (Marsell and Einhorn, 2011). Therefore, when secondary healing occurs we assume a double coupled pathway: first we consider that there will be evolution of αc (contribution of α due mainly to the cartilage formation) secondly this cartilage will be replaced by bone (increasing αb). We have to keep in mind that our model is based on some simplifications, which present relevant implications for the conclusions and the future applications. First, due to the cohesive elements used in its implementation, it can only be used with fractures where the dimension of the gap is thin enough to consider it negligible with respect to the overall dimensions of the bone fracture. Therefore, we have to keep in mind that this model is not valid for large bone defects. Secondly, different mechano-biological models (Andreykiv et al., 2008; Checa and Prendergast, 2009; Isaksson et al., 2008, 2009a; Lacroix and Prendergast, 2002b; Loboa et al., 2001; Reina-Romo et al., 2011; Shefelbine et al., 2005; Simon et al., 2011; Wehner et al., 2010) have been proposed where a detailed study of the whole callus is performed, analyzing
journal of the mechanical behavior of biomedical materials 29 (2014) 328 –338
how different mechanical stimuli are spatially distributed in the callus in order to understand the fundamental cellular mechanisms that locally regulate tissue differentiation and callus growth. However, in this work, we do not study how the mechanical stimulus is distributed in the whole callus dominium, but we take into account the temporal progression of the callus stiffness in the fracture gap regulated by the mechanical environment in this gap. Thirdly, we have considered the following time-dependent assumptions for the temporal evolution of the union degree α: (1) an exponential evolution of the stiffness during the formation of the cartilaginous structure (αc) and (2) a linear evolution of the gap rigidity due to bone formation (αb). Both assumptions are based on previous experiments that we observed in-vivo for sheep animals in tibia fractures (Gómez-Benito et al., 2011). Fourthly, we have estimated the parameters of the model (calibration process, Section 2.4.) using experiments of bone healing in sheep and, afterwards, we have used these parameters in the predictive simulations (Section 2.5.) to simulate the experiments of bone healing both in human and sheep, obtaining good agreement in both cases. However, in the clinical application of this model to simulate fracture healing in human fractures, a finest estimation of the parameters should be done taking into account different aspects like sex, age, weight, etc. Finally, this model is only regulated by mechanical factors; however, as it is well known that bone healing is also regulated by other factors (chemical, genetic, biological etc.). This model could be considered as a useful first approach that could be improved in future works including these additional factors. Therefore, this work presents a simple model able to describe the temporal evolution of the bone fracture mechanical properties under different mechanical stability conditions. In addition, this model is very versatile for the simulation of realistic bone fracture geometries (oblique, comminuted, spiral and compound) allowing a whole-organ analysis. Therefore, this model can be used to study and improve healing of bone fragments stabilized by different fixation systems, such as, locking plates, nails, external fixators and intramedullary screws.
Acknowledgments The research leading to these results has received funding from the (European Commission) Seventh Framework Programme (FP7/2007-2013) under grant agreement nº286179.
r e f e r e n c e s
Ament, C., Hofer, E.P., 2000. A fuzzy logic model of fracture healing. Journal of Biomechanics 33, 961–968. Andreykiv, A., van Keulen, F.V., Prendergast, P.J., 2008. Simulation of fracture healing incorporating mechanoregulation of tissue differentiation and dispersal/proliferation of cells. Biomechanics and Modeling in Mechanobiology 7, 443–461. Augat, P., Burger, J., Schorlemmer, S., Henke, T., Peraus, M., Claes, L., 2003. Shear movement at the fracture site delays healing in a
337
diaphyseal fracture model. Journal of Orthopaedic Research 21, 1011–1017. Bishop, N.E., Schneider, E., Ito, K., 2003. An experimental two degrees-of-freedom actuated external fixator for in vivo investigation of fracture healing. Medical Engineering and Physics 25, 335–340. Bishop, N.E., Rhinjn, M.v., Tami, I., Corveleijn, R., Schneider, E., Ito, K., 2006. Shear does not necessarily inhibit bone healing. Clinical Orthopaedics and Related Research 443, 307–314. Byrne, D.P., Lacroix, D., Prendergast, P.J., 2011. Simulation of fracture healing in the tibia: mechanoregulation of cell activity using a lattice modeling approach. Journal of Orthopaedic Research, October, 1496–1503. Claes, L.E., Wilke, H.-J., Augat, P., Rübenacker, S., Margevicius, K.J., 1995. Effect of dynamization on gap healing of diaphyseal fractures under external fixation. Clinical Biomechanics 10 (5), 227–234. Claes, L.E., Augat, P., Suger, G., Wilke., H.-J., 1997. Influence of size and stability of the osteotomy gap on the success of fracture healing. Journal of Orthopaedic Research 15, 577–584. Claes, L.E., Heigele, C.A., 1999. Magnitudes of local stress and strain along bony surfaces predict the course and type of fracture healing. Journal of Biomechanics 32, 255–266. Comiskey, D., MacDonald, B.J., McCartney, W.T., Synnott, K., O’Byrne, J., 2013. Predicting the external formation of callus tissues in oblique bone fractures: idealised and clinical case studies. Biomechanics and Modeling in Mechanobiology http: //dxdoi.org/10.1007/s10237-012-0468-6. Chao, E., Aro, H., DGLewallen, Kelly, P., 1989. The effect of rigidity on fracture healing in external fixation. Clinical Orthopaedics and Related Research 241, 24–35. Checa, S., Prendergast, P., 2009. A mechanobiological model for tissue differentiation that includes angiogenesis: a latticebased modeling approach. Annals of Biomedical Engineering 37 (1), 129–145. Draper, E.R.C., Strachan, R.K., Hughes, S.P.F., Nicol, A.C., Paul, J.P., 1997. The design and performance of an experimental external fixator with variable axial stiffness and a compressive force transducer. Medical Engineering and Physics 19 (8), 660–695. Einhorn, T.A., 2005. The science of fracture healing. Journal of Orthopaedic Trauma 19 (10), S4–S6. García-Aznar, J.M., Kuiper, J.H., Gómez-Benito, M.J., Doblaré, M., Richardson, J.B., 2007. Computational simulation of fracture healing: Influence of interfragmentary movement on the callus growth. Journal of Biomechanics 40, 1467–1476. García-Aznar, J.M., Pérez, M.A., Moreo, P., 2009. Modelling of interfaces in biomechanics and mechanobiology. Computer Modeling in Engineering & Sciences: CMES 48 (3), 271–301. Gómez-Benito, M.J., García-Aznar, J.M., Kuiper, J.H., Doblaré, M., 2005. Influence of fracture gap size on the pattern of long bone healing: a computational study. Journal of Theoretical Biology 235, 105–119. Gómez-Benito, M.J., González-Torres, L.A., Reina-Romo, E., Grasa, J., Seral, B., García-Aznar, J.M., 2011. Influence of high-frequency cyclical stimulation on the bone fracture-healing process: mathematical and experimental models. Philosophical Transactions Series A, Mathematical, Physical, and Engineering Sciences 369, 4278–4294. Goodship, A.E., Kenwright, J., 1985. The influence of induced micromovement upon the healing of experimental tibial fractures. Journal of Bone and Joint Surgery (British Volume) 67-b (4), 650–655. Goodship, A.E., Watkins, P.E., Rigby, H.S., Kenwright, J., 1993. The role of fixator frame stiffness in the control of fracture healing. An experimental study. Journal of Biomechanics 26 (9), 1027–1035. Hadjiargyrou, M., McLeod, K., Ryaby, J., Rubin, C., 1998. Enhancement of fracture healing by low intensity ultrasound. Clinical Orthopaedics and Related Research 355 (S), 216–229.
338
journal of the mechanical behavior of biomedical materials 29 (2014) 328 –338
Hernandez, C.J., Beaupre, G.S., Keller, T.S., Carter, D.R., 2001. The influence of bone volume fraction and ash fraction on bone strength and modulus. Bone 29 (1), 74–78. Isaksson, H., Donkelaar, C.C.v., Huiskes, R., Ito, K., 2008. A mechano-regulatory bone-healing model incorporating cellphenotype specific activity. Journal of Theoretical Biology 252, 230–246. Isaksson, H., Donkelaar, C.C.v., Ito, K., 2009a. Sensitivity of tissue differentiation and bone healing predictions to tissue properties. Journal of Biomechanics 42, 555–564. Isaksson, H., Gröngröft, I., Wilson, W., Donkelaar, C.C.v., Rietbergen, B.v., Tami, A., Huiskes, R., Ito, K., 2009b. Remodeling of fracture callus in mice is consistent with mechanical loading and bone remodeling theory. Journal of Orthopaedic Research 27 (5), 664–672. Jagodzinski, M., Krettek, C., 2007. Effect of mechanical stability on fracture healing – an update. Injury 38 (S1), S3–S10. Klein, P., Opitz, M., Schell, H., Taylor, W.R., Heller, M.O., Kassi, J.-P., Kandziora, F., Duda, G.N., 2004. Comparison of unreamed nailing and external fixation of tibial diastases – mechanical conditions during healing and biological outcome. Journal of Orthopaedic Research 22, 1072–1078. Lacroix, D., Prendergast, P.J., 2002a. Three-dimensional simulation of fracture repair in the human tibia. Computer Methods in Biomechanics and Biomedical Engineering 5 (5), 369–376. Lacroix, D., Prendergast, P.J., 2002b. A mechano-regulation model for tissue differentiation during fracture healing: analysis of gap size and loading. Journal of Biomechanics 35, 1163–1171. Loboa, E.G., Beaupré, G.S., Carter, D.R., 2001. Mechanology of initial pseudarthrosis formation with oblique fractures. Journal of Orthopaedic Research 19, 1067–1072. Manjubala, I., Liu, Y., Epari, D.R., Roschger, P., Schell, H., Fratzl, P., Duda, G.N., 2009. Spatial and temporal variations of mechanical properties and mineral content of the external callus during bone healing. Bone 45, 185–192. Marsell, R., Einhorn, T., 2011. The biology of fracture healing. Injury 42 (6), 551–555. Moreo, P., Pérez, M.A., García-Aznar, J.M., Doblaré, M., 2007. Modelling the mechanical behaviour of living bony interfaces.
Computer Methods in Applied Mechanics and Engineering 196, 3300–3314. Park, S.-H., O’Connor, K., McKellop, H., Sarmiento, A., 1998. The influence of active shear or compressive motion on fracture healing. Journal of Bone and Joint Surgery (American Volume) 80-A (6), 868–878. Perren, S.M., 2002. Evolution of the internal fixation of long bone fractures. Journal of Bone and Joint Surgery (British Volume) 84-B (8), 1093–1110. Reina-Romo, E., Gómez-Benito, M.J., Domínguez, J., Niemeyer, F., Wehner, T., Simon, U., Claes, L.E., 2011. Effect of the fixator stiffness on the young regenerate bone after bone transport: computational approach. Journal of Biomechanics 44, 917–923. Richardson, J.B., Cunningham, J.L., Goodship, A.E., O’Connor, B.T., Kenwright, J., 1994. Measuring stiffness can define healing of tibial fractures. Journal of Bone and Joint Surgery (British Volume) 76-B (3), 389–394. Schell, H., Epari, D.R., Kassi, J.P., Bragulla, H., Bail, H.J., Duda, G.N., 2005. The course of bone healing is influenced by the initial shear fixation stability. Journal of Orthopaedic Research 23, 1022–1028. Shefelbine, S.J., Augat, P., Claes, L., Simon, U., 2005. Trabecular bone fracture healing simulation with finite element analysis and fuzzy logic. Journal of Biomechanics 38, 2440–2450. Simon, U., Augat, P., Claes, L., Utz, M., 2011. A numerical model of the fracture healing process that describes tissue development and revascularisation. Computer Methods in Biomechanics and Biomedical Engineering 14 (1), 79–93. Vetter, A., Epari, D.R., Robin Seidel, H.S., Fratzl, P., Duda, G.N., Weinkamer, R., 2010. Temporal tissue patterns in bone healing of sheep. Journal of Orthopaedic Research 28 (11), 1440–1447. Wehner, T., Claes, L., Niemeyer, F., Nolte, D., Simon, U., 2010. Influence of the fixation stability on the healing time—a numerical study of a patient-specific fracture healing process. Clinical Biomechanics 25, 606–612. Wehner, T., Penzkofer, R., Augat, P., Claes, L., Simon, U., 2011. Improvement of the shear fixation stability of intramedullary nailing. Clinical Biomechanics 26, 147–151.
Available online at www.sciencedirect.com
www.elsevier.com/locate/jmbbm
Research Paper
An interface finite element model can be used to predict healing outcome of bone fractures J.A. Aliertaa, M.A. Pe´rezb, J.M. Garcı´a-Aznarb,n a
Escuela Politécnica Superior del Ejército, Ministry of Defense, Madrid, Spain Mechanical Engineering, Aragón Institute of Engineering Research (I3A), University of Zaragoza, Ed. Betancourt-Campus Río Ebro, C/María de Luna, 50018 Zaragoza, Spain b
ar t ic l e in f o
abs tra ct
Article history:
After fractures, bone can experience different potential outcomes: successful bone con-
Received 16 May 2013
solidation, non-union and bone failure. Although, there are a lot of factors that influence
Received in revised form
fracture healing, experimental studies have shown that the interfragmentary movement
20 September 2013
(IFM) is one of the main regulators for the course of bone healing. In this sense,
Accepted 23 September 2013
computational models may help to improve the development of mechanical-based
Available online 8 October 2013
treatments for bone fracture healing. Hence, based on this fact, we propose a combined
Keywords:
repair-failure mechanistic computational model to describe bone fracture healing. Despite
Bone fracture healing
being a simple model, it is able to correctly estimate the time course evolution of the IFM
Finite element prediction
compared to in vivo measurements under different mechanical conditions. Therefore, this
Design of fracture fixators
mathematical approach is especially suitable for modeling the healing response of bone to
Interfragmentary movement
fractures treated with different mechanical fixators, simulating realistic clinical conditions. This model will be a useful tool to identify factors and define targets for patient specific therapeutics interventions. & 2013 Elsevier Ltd. All rights reserved.
1.
Introduction
Fractures are a common orthopedic problem and the fracture healing is a natural process that regenerates bone to its original state and function. Several factors influence these bone healing events, such as genetic, cellular and biochemical factors, blood supply, neural and hormonal regulation, age, the type of fracture interfragmentary motion and fracture geometry (Einhorn, 2005; Goodship et al.,1993; Hadjiargyrou et al., 1998; Jagodzinski and Krettek, 2007; Marsell and Einhorn, 2011). However, the most common orthopedic treatments consist on the mechanical stabilization of the bone fracture gap, regulating the interfragmentary movement. n
Corresponding author. Tel.: þ34 976 76 27 96; fax: þ34 976 76 26 70. E-mail address: [email protected] (J.M. García-Aznar).
1751-6161/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jmbbm.2013.09.023
The role of the mechanical stabilization on the design of fracture fixators has evolved very much in the last years, updating from the concept of “open reduction and internal fixation (ORIF)” to “biological fixation” (Perren, 2002). In both cases, differences are presented from a biological and mechanical point of view. In fact, the concept of “biological fixation” is based on the application of the fixator as a minimally invasive percutaneous osteosynthesis (MIPO). So, the contact of the implant with bone is reduced at maximum, avoiding the damage to the blood supply. In addition to biological differences, there is also a different concept in the role of the mechanical stabilization. In the “biological fixation” primary ossification is avoided, promoting the occurrence of a secondary
journal of the mechanical behavior of biomedical materials 29 (2014) 328 –338
ossification mechanism that induces the formation of bone callus (Manjubala et al., 2009; Vetter et al., 2010). Therefore, the mechanical design concept is updated from an absolute stability to a controlled mechanical instability that favors the formation of bone callus. Hence, the mechanical stabilization due to the stiffness of the fixator has acquired a relevant role for achieving a successful bone healing (Bishop et al., 2003; Chao et al., 1989; Draper et al., 1997; Goodship et al., 1993; Richardson et al., 1994; Wehner et al., 2011). Many studies have examined the role of different local mechanical conditions on bone fracture healing: the influence of the magnitude of interfragmentary movement (IFM) and the initial size of the fracture gap (Claes et al., 1997; Gómez-Benito et al., 2011; Goodship and Kenwright, 1985), the stiffness of the fixation (Schell et al., 2005) and the type of movement of the gap (Augat et al., 2003; Bishop et al., 2006). Computational models of the fracture healing process may prove useful in the determination of optimal mechanical treatments. In fact, several mechano-biological models based on finite element simulations studied the influence of local mechanical conditions on biological events that regulate the temporal and spatial evolution of the different ossification mechanisms that occur during healing (Andreykiv et al., 2008; Checa and Prendergast, 2009; Isaksson et al., 2008, 2009a; Lacroix and Prendergast, 2002b; Loboa et al., 2001; ReinaRomo et al., 2011; Shefelbine et al., 2005; Simon et al., 2011; Wehner et al., 2010). Most of these models have focused on how different mechanical stimuli (such as: fluid flow, octahedral shear strain, deviatoric strain, strain energy density, etc.) are able to predict spatial pattern distribution of tissues regulated by intramembranous and endochondral ossification (Checa and Prendergast, 2009; Lacroix and Prendergast, 2002a; Wehner et al., 2010). Other models have combined tissue differentiation rules with callus growth (Comiskey et al., 2013; García-Aznar et al., 2007; Gómez-Benito et al., 2005) or callus resorption (Ament and Hofer, 2000; Byrne et al., 2011; Isaksson et al., 2009b; Lacroix and Prendergast, 2002b) depending on the temporal stage of bone healing that they studied. Although these mathematical models are very useful to understand the fundamental cellular mechanisms that locally regulate tissue differentiation and callus growth/ resorption, they are not really helpful for the simulation of realistic bone fractures, where a whole-organ analysis is required. In fact, these models have only analyzed simple tranversal or obliques (Comiskey et al., 2013; Loboa et al., 2001) bone fractures with two fragments. However, fractures are much more complexes with very different shapes of the fracture line, complicated anatomical locations and a different number of fragments. Therefore, the main aim of the present study is to propose a phenomenological computational model able to simulate the temporal recovery of mechanical properties of the fracture zone during the healing process, which is regulated by the mechanical stability. If this approach proves feasible, it offers the possibility of using computer simulations in the clinical treatment of complex fractures with multiple fragments, complicated geometries and different anatomical locations and in other orthopedic applications where bone regeneration occurs.
2.
329
Material and methods
The bone fracture gap was modeled through the incorporation of interface elements that connected the two fracture ends simulating the discontinuity in the displacement field between fragments. A mathematical model is here proposed to simulate the temporal evolution of the separation between both fragments, regulated by the mechanical conditions existent in the fracture gap.
2.1.
Fracture gap/interface mechanical behavior
To model the fracture gap behavior, 6-nodes and 8-nodes cohesive elements (Fig. 1) are used to connect the fracture ends (García-Aznar et al., 2009). The thickness of these elements, which is the dimension of the gap fracture, is thin enough to consider it negligible with respect to the overall dimensions of the bone fracture. Accordingly, the behavior of these elements is directly established in terms of one traction versus one separation law. The nominal traction stress vector, t, consists of three components: tn, ts, tt, which represents the normal (tn) (along
Fig. 1 – Constitutive model: (a) cohesive elements, (b) shear traction and (c) normal traction.
330
journal of the mechanical behavior of biomedical materials 29 (2014) 328 –338
the local 3-direction) and the two shear tractions (ts, tt) (along the local 1- and 2-directions), respectively (Fig. 1). The corresponding separations are denoted by δn, δs and δt. Thus, the nominal strains can be defined as εn ¼
δn δs δt ; εs ¼ ; εt ¼ h h h
ð1Þ
where h is the original thickness of the cohesive element (fracture gap thickness). The purpose of the present study is to simulate the temporal evolution of the bone healing process. Therefore, to quantify the degree of healing, union degree, we use α as the main state variable. This variable is normalized, corresponding α¼ 1 to a totally successful bone healing and α¼0 to a completely non-union or malunion. The union degree α, may be also related with the so-called damage variable, d (Moreo et al., 2007), through the following expression: d ¼ 1 α
ð2Þ
The net effect of the bone healing process on the macroscopic interface mechanical properties is what it is represented by the growth of α. In addition, the model allows the possibility of reducing the union degree α, which means mechanical failure of the gap, because it is not able to bear the local demands. For each element of the fracture gap, the union degree, α, can be defined as K0i α ¼ Ki
i ¼ n; s; t
ð3Þ
where K0i defines the linear stiffness in a completely healed fracture gap along the direction i (the maximum fracture gap stiffness available), and Ki is the current interface stiffness in the i direction (Fig. 1b and c). The three components of the nominal traction stress vector, t, can be calculated, using the above defined variables, through the following equation: ti ¼ K0i αεi
2.2.
i ¼ n; s; t
ð4Þ
Mechanical damage: α decrease
If mechanical demands in the fracture gap are important, a reduction of the degree of union may occur as consequence of its mechanical failure. For the shear traction, we assume the degree of union is reduced by the following law (Fig. 1b): ( ðε ε Þ i ci ðε0i εci Þ ; if εαi oεi oεci i ¼ s; t ð5Þ αi ¼ 0; if εi Zεci where, ε0i maximum shear strain at the linear region for α¼ 1 (perfectly healed interface), εαi maximum shear strain at the linear region for each value of α and εci the maximum shear strains allowed. If εci are exceeded, the fracture gap would have no stiffness at all. For the compression traction (εn o0) we assume a linear behavior with respect to the associated strain until εn reaches the value 1, which means that the fracture ends are in contact. If lower values of εn are reached the simulation will be aborted because no penetration is allowed.
For the tensile traction (εn Z0) we assume an initially linear behavior with respect to the associated normal strain, followed by a linear decay (Fig. 1c): ( ðε ε Þ n cn ; if εαn oεn oεcn ð6Þ αn ¼ ðε0n εcn Þ 0; if εn Zεcn where ε0n, εαn and εcn have an equivalent meaning that those defined for shear traction. In Fig. 1b and c the dashed line represents the behavior law (X¼ε (strain); Y ¼t (stress)) of a completely healed bone gap (α¼1). The continuous line represents the behavior law (t–ε) of a particular bone gap through the healing process (αo1). In these figures, the parts of the graphs with negative slope represent the Eq. (5) (Fig. 1b) and 6 (Fig. 1c). In this part of the graph, α decreases (α_ o0) because the limit strain (εα) is reached (for example, due to a very high load). If this fact occurs during the healing process, the healing will be impaired or at least delayed. Therefore, we finally select as union degree, α, the minimum value associated to each direction: α ¼ min fαs ; αt ; αn g
2.3.
ð7Þ
Biological union: growth of union degree α
The union degree, α, can increase due to bone healing, and therefore it is calculated as the addition of two contributions: α ¼ αc þ αb
ð8Þ
where αc reflects the recovery of the mechanical properties of the fracture gap mainly due to the formation of cartilage, and αb reflects the recovery of the mechanical properties of the fracture gap due to bone formation. In order to estimate the temporal evolution of α ðα_ ¼ α_ c þ α_ b Þ, we propose a regulatory law based on three different healing zones (adapted from Claes and Heigele, 1999) which are established as of ffi the qffiffiffiffiffiffiffiffiffiffiffiffiffiffi a function compression, εn, of each interface and the shear strains εshear ¼ ε2s þ ε2t element (Fig. 2a) with the following characteristics: Zone 1 non-union: The compression and shear strains are so high that there is a high interfragmentary movement (IFM) which impairs a successful fracture healing. Only fibrous tissue is formed and there is no change in αc nor in αb ðα_ c ¼ α_ b ¼ 0Þ Zone 2 cartilage formation: The compression and shear strains are moderate, in the same way that the IFM. We consider that in those elements, whose strains belong to this zone, cartilage will be formed and therefore αc will grow. However, these elements should be under this strain threshold during a certain period of time (Mc) before the growing of αc, this time represents the maturation time needed for the cartilage to appear close to the fracture gap. It is established that αc will have an exponential evolution with the time (Fig. 2b), therefore once the maturation time for cartilage (Mc) is reached, the value of αc will grow exponentially, providing that the strains keep in zone 2, during a fixed period of time, tc (cartilage formation time), until it reaches the maximum value (αcmax). If during the healing process, αc is
journal of the mechanical behavior of biomedical materials 29 (2014) 328 –338
growing and the strains reach the zone 3, αc will continue growing with the same rate until the maximum value (αcmax). Zone 3 bone formation: In this zone strains and IFM are lower than in zone 2, as a result, bone is formed, so αb will grow. We also define, Mb, as the maturation time for the formation of bone. It is established that αb will have a linear evolution with the time (Fig. 2c), therefore once the maturation time for cartilage (Mb) is reached, the value of αb will
331
grow linearly, providing that the strains keep in the zone 3, until it reaches the maximum value (αbmax) at the total healing time, th. If αb begins to grow because the conditions for bone formation are achieved and αc is still growing, a simultaneous growing of both contributions will happen. According to the definition of α, its maximum value is 1 (corresponding to a totally successful bone healing), therefore the limits for αc (αcmax) and αb (αbmax) must satisfy the following equation: αmax ¼ αb max þ αc max ¼ 1
ð9Þ
sh The model parameters values represented in Fig. 2 (Lsh c , Lb , comp comp Lc , Lb , αcmax, αbmax, Mc, Mb, tc, th) will be later defined by tuning the model to experimental results. The model was numerically implemented in an Abaqus user subroutine (UMAT) and all the finite element analyses (FEA) were carried out in ABAQUS v.6.11.
2.4.
Examples of application for the calibration process
In order to determine the value of the parameters that define the model, we made a quantitative comparison between the numerical results (obtained from our model) and the experimental ones found in the literature. With this idea in mind, we simulated two animal experiments (compression and shear) from the literature that are detailed below and summarized in Table 1.
2.4.1.
Fig. 2 – (a) Healing zones, (b) temporal evolution of αc and (c) temporal evolution of αb.
Compression experiments
To estimate the compression parameters (Lcomp , Lcomp , αcmax, c b αbmax, Mc, Mb, tc, th), we simulated the compression experiments performed by Claes et al. (1997). They determined the influence of the size of the osteotomy gap and interfragmentary strain on callus formation in the healing of the sheep metatarsus (Claes et al., 1997). Therefore, a cylinder of external diameter of 17 mm was simulated (Fig. 3a). The fracture gap was subjected to an external compression load of 330 N (Claes et al., 1995). The gap size was 2 mm and the initial interfragmentary strain (IFS¼IFM/h) was approximately 31%, that was achieved with the use of a fixator, simulated by four linear springs with a total stiffness of 348 N/mm. In order to avoid rigidbody motions the nodes at the base of the model were fully constrained. Bone was considered as a linear elastic material for all analyses where the mechanical properties are shown in Table 2.
Table 1 – Examples of model application. Goal
Experiment
Bone
Load condition
Reference
Calibration process
Compression Shear
Sheep metatarsus Sheep tibia
Compression Shear
Claes et al. (1997) Bishop et al. (2006)
Evaluation of the model predictive capacity
Compression Combined load case Comminuted fracture
Sheep metatarsus Human tibia Sheep metatarsus
Compression Combined Combined
Claes et al. (1997) Wehner et al. (2010) –
332
journal of the mechanical behavior of biomedical materials 29 (2014) 328 –338
Fig. 3 – FE model for the simulation of (a) compression and (b) shear experiments in order to calibrate the corresponding parameters of the model.
Table 2 – Material properties of the tissue considered for all analyses.
Table 3 – Fixation stiffness for each experiment simulated in the compression predictive simulation.
Tissue type
Elastic modulus in MPa
Poisson's ratio
Experiment
Stiffness (N/mm)
Cortical bone Trabecular bone
10,000a 200a
0.2a 0.2a
Initial gap ¼1 mm; IFS ¼ 7% Initial gap ¼1 mm; IFS ¼ 31% Initial gap ¼2 mm; IFS ¼ 7%
3300 764 1650
a
Hernandez et al. (2001).
2.4.2.
Shear experiments
In a similar way, it is necessary that the parameters of the model allow also to simulate correctly the healing process under shear conditions. However, no experiments in the literature have been found in which, the time course of the interfragmentary movement (IFM) during the healing is evaluated under a pure shear load case. Nevertheless Bishop et al. (2006) performed experiments studying the healing of a transverse tibial osteotomy in sheep where interfragmentary torsional shear (torsion) was applied under carefully controlled conditions and some data of this experiment are used for the estimation of the rest of the parameters of our sh model (Lsh c , Lb ). We simulated these experiments taking the geometry and the loads from Bishop et al. (2006), (Fig. 3b). A cylinder of external diameter of 20 mm was simulated. A torsional rotation (torsion) was applied over a 2.4 mm osteotomy with a maximum interfragmentary principal strain magnitude of 25%. This corresponded to a 7.21 rotation calculated in the outer cortical radius of 10 mm. A load limit of 1670 N mm was imposed to allow interfragmentary strain to decrease with increasing callus stiffness. In order to avoid rigid-body motions the nodes at the base of the model were constrained.
To calibrate the model two data sets in the related article (Bishop et al., 2006) were used: – The stiffness increased steadily during the experiment with the load limit achieved within the fourth week. – After eight weeks of healing, the bending stiffness was approximately the 80% of the intact bone. The union degree α, is the variable that represents the improvement of the stiffness and then is directly related to the bending stiffness.
2.5. Examples of application for the evaluation of the model predictive capacity In order to evaluate the predictive capacity of the model, three additional experiments (compression, combined load case and comminuted fracture) that are summarized in Table 1 were simulated using the parameters previously obtained in the calibration process.
2.5.1.
Compression experiments
Three compression experiments, similar to the described in the model calibration section (Fig. 3a), were performed by
journal of the mechanical behavior of biomedical materials 29 (2014) 328 –338
Claes et al. (1997). The differences among them were the gap dimension and the initial IFS. In fact, gaps of 1 mm and 2 mm and initial IFS of 7% and 31% were fixed in these experiments. We simulated these three experiments with the same FE model used in the model calibration section, changing the gap dimension and the stiffness of the fixator (Table 3). The same loading and boundary conditions were used.
2.5.2.
Combined load case experiment
We found no more data in the literature of pure shear experiments of fracture healing than those used in the calibration (Bishop et al., 2006). That fact prevented us from validating our model exclusively with this type of loading condition. However, Wehner et al. (2010) performed a clinical experiment in a human patient with a fracture in the diaphyseal tibia. A combined compression-shear load case was applied, and at several time points during healing, the IFM was determined during a single leg stance. Thanks to these available data, simulating this experiment we were able to evaluate the potential of our model under a complex combined load case. A human tibia was scanned and the images stored in Dicom format. MIMICS V. 15.0 was used to reconstruct the geometry (Fig. 4a) through the segmentation process. Afterwards, using this software, an osteotomy was made in the diaphyseal part. In addition, a unilateral external fixator was modeled according to the specifications shown in Wehner et al. (2010). This fixator, once imported in MIMICS, was implanted in the tibia, and the assembly was meshed with tetrahedral elements. We used ABAQUS to create cohesive elements in the fracture gap. The final FE model had a total number of nodes of 19,499 and a total number of elements of 91,558 (118 linear wedge cohesive elements and 91,440 linear tetrahedral elements). The element size used was inside the asymptotic region of convergence and represents a good trade off between numerical accuracy and computational cost (Fig. 4c).
333
To simulate the experiment, we applied a vertical load of 700N (corresponding to a single leg stance of a patient of approximately 70 kg (Wehner et al.,2010)) at the proximal condyles of the tibia (Fig. 4b). Due to this load, a combined load state (compression-shear) was generated in the fracture gap zone. The nodes of the base were constrained so that, the rigid-body motions are avoided.
2.5.3.
Comminuted fracture
A complex fracture was simplified and simulated in a sheep metatarsus. The comminuted fracture modeled consisted of an oblique fracture with three fragments (Fig. 5). The loads, fixation system and bone geometry were the same that those used in the Section 2.4.1. Due to the complex geometry of the fracture zone the vertical load applied generated combined mechanical demands in the different fracture gaps. Four different simulations were made changing the stiffness of the external fixator in order to get the best fixation
Fig. 5 – FE model for the simulation of the comminuted fracture.
1.6 Simulated
1.4
Experimental (Claes,1997)
IFM (mm)
1.2 1 0.8 0.6 0.4 0.2 0 0
10
20
30
40
50
60
Healing time (days)
Fig. 4 – Predictive simulation using combined load case: (a) geometry, (b) applied loads and (c) FE of the fracture gap zone.
Fig. 6 – Compression calibration results: evolution of the interfragmentary motion in the computational simulation compared with the experimental data in a 2-mm gap of a sheep (Claes et al., 1997).
334
journal of the mechanical behavior of biomedical materials 29 (2014) 328 –338
Table 4 – Model parameters. Parameter
Value
Mc—Maturation time for cartilage (days) (see Fig. 2b) Mb—Maturation time for bone (days) (see Fig. 2c) αcmax—Maximum value of αc (see Fig. 2b) αbmax—Maximum value of αb (see Fig. 2c) tc—Cartilage formation time (days) (see Fig. 2b) th—Total healing time (days) (see Fig. 2c) Lcomp —Compression strain limit for cartilage formation (see Fig. 2a) c
4a 10b 0.3b 0.7b 36b 60c 0.6b
comp
Lb
0.25b
—Compression strain limit for bone formation (see Fig. 2a)
Lsh c —Shear Lsh b —Shear
0.25b
strain limit for cartilage formation (see Fig. 2a)
0.15b
strain limit for bone formation (see Fig. 2a) K0i (i ¼n,s,t)—Initial linear stiffness (MPa) ε0i (i¼ n,s,t)—Maximum strain at the linear region εci (i¼n,s,t)—Maximum allowed strain
50d 0.3b 1b
a
Park et al. (1998). Estimated. c Klein et al. (2004). d Bishop et al. (2003). b
system for this fracture. The stiffness of the fixator used for each simulation was 400, 448, 500 and 560 N/mm.
3.
Results
Here, in first step, we discuss the results of the calibration process, and next, we also explain additional experiments to evaluate the predictive capacity of the model.
3.1.
Calibration process results
3.1.1.
Compression experiments results
Evolution of the IFM obtained experimentally (Claes et al., 1997) and computationally has been represented in Fig. 6 for a 2-mm gap of a sheep. The experimental variability has been also represented and the computational results are within previous range of variation. Initially, there is no fracture healing, because the IFM is relatively high (1 mm), as the simulation evolves, the IFM is reduced. After 40 days of healing, the IFM has been reduced by a ratio of 80% (see Fig. 6).
3.1.2.
Shear experiments results
The temporal evolution of the union degree, α, of the fracture gap has been estimated under a torsional load case. We have been able to correctly reproduce the experimental results obtained by Bishop et al. (2006). Experimentally the torsion load limit of 1670N mm is achieved within the fourth week, using our model, we estimate1 that this load limit is achieved the 32nd day.
3.1.3.
Parameters of the model
Using the results obtained in the calibration process, the parameters that best fit the simulated and the experimental 1 We determine the torsion load limit as the resultant torsion moment, which is computed by means of the integration of the shear traction stress vector over the bone section.
results are shown in Table 4. Some of the model parameters come from previous studies of the literature and the rest are estimated with the results of Sections 3.1.1and 3.1.2.
3.2.
Evaluation of the model predictive capacity—results
3.2.1.
Compression experiments results
The experimental results obtained by Claes et al. (1997) for different fracture gap sizes and IFS have been represented in Fig. 7, where we can also see the experimental range of variation. In addition, we include in this figure the numerical results corresponding to the computational simulations. As we can see, a good agreement is found between numerical and experimental results for all the studied conditions.
3.2.2.
Combined load case experiment results
Experimental results from Wehner et al. (2010) obtained for a combined load case have been represented and compared with the computational predicted ones in Fig. 8. Although there were no more experimental data in order to represent its range of variation, our computational results are very close to the experimental ones. The IFM compression is very high at the initial healing stages, and as the healing process evolves the IFM is reduced both experimentally and computationally due to the combined load case. However, in the case of shear loads, the shear IFM slightly decreases from the first days until the last days of healing (Fig. 8).
3.2.3.
Comminuted fracture results
The temporal evolution of the vertical displacement of the fixation points has been represented for four different fixators (Fig. 9). It can be observed that the stiffness value that leads to better results is K ¼ 500 N/mm, because the final interfragmentary displacements are lower, which means better healing conditions. We have also represented the temporal evolution of the bone reaction (Fig. 10) i.e. the load that the bone is able to
335
journal of the mechanical behavior of biomedical materials 29 (2014) 328 –338
0.6 0.25 0.5
Simulated
0.4
Experimental (Claes,1997) Gap=2 mm; Initial IFS=7%
Simulated Experimental (Claes,1997) Gap=1mm; Initial IFS=7%
0.15
IFM (mm)
IFM (mm)
0.2
0.1
0.3 0.2
0.05
0.1
0
0 0
10
20
30
40
50
60
0
10
20
Healing time (days)
30
40
50
60
Healing time (days)
0.3 Simulated 0.25
Experimental (Claes,1997) Gap=2 mm; Initial IFS=7%
IFM (mm)
0.2
0.15
0.1
0.05
0 0
10
20
30
40
50
60
Healing time (days)
Fig. 7 – Evolution of the interfragmentary motion in the computational simulation compared with the experimental data (Claes et al., 1997): (a) gap¼ 1 mm; initial IFS¼ 7%, (b) gap¼ 1mm; initial IFS ¼31%, and (c) gap¼ 2 mm; initial IFS¼7%. 1.6
Experimental IFM compression (Wehner et al., 2010)
1.4
Experimental IFM shear (Wehner et al., 2010) Simulated IFM compression
IFM (mm)
1.2 1
Simulated IFM shear
0.8
If we compare the different stiffness able to heal a simple transverse fracture (K¼ 348 N/mm see Section 2.4.1) and a comminuted fracture (K¼ 500 N/mm) under the same loading conditions, we can observe that the second one requires a higher stiffness to achieve successful healing.
0.6 0.4
4.
Discussion
0.2 0 0
10
20
30
40
50
60
70
Healing time (days)
Fig. 8 – Predicted time course of the interfragmentary movement (IFM) due to compression and shear stresses compared to in vivo measurements of the patients (Wehner et al., 2010).
bear during the healing process. The best result is for a stiffness of K ¼500 N/mm, however it can be seen that K ¼560 N/mm (more rigid fixation) gives better results in the first healing stages probably due to the primary bone formation. The more flexible fixation systems give poorer results probably because the strains are too high and the healing process is not successful.
In this study we propose a coupled computational model able to simulate the normal and the impaired (mechanical failure) bone healing regulated by the mechanical stability or the interfragmentary movement in the fracture gap. Indeed, the model predicted the temporal recovery/failure of mechanical properties of the fracture gap during the healing. The aim of this model is to identify the role of different and combined degrees of stabilization for fracture fixation. The model here presented is based on the main mechanisms that regulate the biology of fracture healing. In fact, bone heals by either direct intramembranous ossification or indirect fracture healing, which consists of both intramembranous and endochondral ossification (Marsell and Einhorn, 2011). Therefore, this model combines both approaches, modeling both direct and indirect ossification or equivalently intramembranous and endochondral ossification.
336
journal of the mechanical behavior of biomedical materials 29 (2014) 328 –338
Fig. 9 – Comminuted fracture results: temporal evolution of the vertical displacement of the fixation points (see position of these points in Fig. 5) for different fixation stiffness.
Fig. 10 – Temporal evolution of bone reaction (N) for different values of stiffness of the fixator. Primary healing (intramembranous ossification) requires extremely rigid fixation systems and negligible spaces between bone fragments. The quality of regenerated bone in primary healing cannot be always guaranteed and therefore a new fracture may occur. In those parts of the fracture gap in which, only primary healing occurs, we assume that there will be only the temporal evolution of one part (αb) of the union degree (α) (contribution of α mainly due to the bone formation). On contrast, secondary healing (endochondral ossification) takes place under more flexible fixation that allows
a higher, but controlled, instability of the fracture gap region. In this second step, bone formation is involved restoring the mechanical stability and recovering the mechanical properties in the bone fracture gap. In any case, secondary healing is the most common in clinical treatments, since the high quality of the regenerated tissue and also the probability of a successful healing is higher than in primary healing (Marsell and Einhorn, 2011). Therefore, when secondary healing occurs we assume a double coupled pathway: first we consider that there will be evolution of αc (contribution of α due mainly to the cartilage formation) secondly this cartilage will be replaced by bone (increasing αb). We have to keep in mind that our model is based on some simplifications, which present relevant implications for the conclusions and the future applications. First, due to the cohesive elements used in its implementation, it can only be used with fractures where the dimension of the gap is thin enough to consider it negligible with respect to the overall dimensions of the bone fracture. Therefore, we have to keep in mind that this model is not valid for large bone defects. Secondly, different mechano-biological models (Andreykiv et al., 2008; Checa and Prendergast, 2009; Isaksson et al., 2008, 2009a; Lacroix and Prendergast, 2002b; Loboa et al., 2001; Reina-Romo et al., 2011; Shefelbine et al., 2005; Simon et al., 2011; Wehner et al., 2010) have been proposed where a detailed study of the whole callus is performed, analyzing
journal of the mechanical behavior of biomedical materials 29 (2014) 328 –338
how different mechanical stimuli are spatially distributed in the callus in order to understand the fundamental cellular mechanisms that locally regulate tissue differentiation and callus growth. However, in this work, we do not study how the mechanical stimulus is distributed in the whole callus dominium, but we take into account the temporal progression of the callus stiffness in the fracture gap regulated by the mechanical environment in this gap. Thirdly, we have considered the following time-dependent assumptions for the temporal evolution of the union degree α: (1) an exponential evolution of the stiffness during the formation of the cartilaginous structure (αc) and (2) a linear evolution of the gap rigidity due to bone formation (αb). Both assumptions are based on previous experiments that we observed in-vivo for sheep animals in tibia fractures (Gómez-Benito et al., 2011). Fourthly, we have estimated the parameters of the model (calibration process, Section 2.4.) using experiments of bone healing in sheep and, afterwards, we have used these parameters in the predictive simulations (Section 2.5.) to simulate the experiments of bone healing both in human and sheep, obtaining good agreement in both cases. However, in the clinical application of this model to simulate fracture healing in human fractures, a finest estimation of the parameters should be done taking into account different aspects like sex, age, weight, etc. Finally, this model is only regulated by mechanical factors; however, as it is well known that bone healing is also regulated by other factors (chemical, genetic, biological etc.). This model could be considered as a useful first approach that could be improved in future works including these additional factors. Therefore, this work presents a simple model able to describe the temporal evolution of the bone fracture mechanical properties under different mechanical stability conditions. In addition, this model is very versatile for the simulation of realistic bone fracture geometries (oblique, comminuted, spiral and compound) allowing a whole-organ analysis. Therefore, this model can be used to study and improve healing of bone fragments stabilized by different fixation systems, such as, locking plates, nails, external fixators and intramedullary screws.
Acknowledgments The research leading to these results has received funding from the (European Commission) Seventh Framework Programme (FP7/2007-2013) under grant agreement nº286179.
r e f e r e n c e s
Ament, C., Hofer, E.P., 2000. A fuzzy logic model of fracture healing. Journal of Biomechanics 33, 961–968. Andreykiv, A., van Keulen, F.V., Prendergast, P.J., 2008. Simulation of fracture healing incorporating mechanoregulation of tissue differentiation and dispersal/proliferation of cells. Biomechanics and Modeling in Mechanobiology 7, 443–461. Augat, P., Burger, J., Schorlemmer, S., Henke, T., Peraus, M., Claes, L., 2003. Shear movement at the fracture site delays healing in a
337
diaphyseal fracture model. Journal of Orthopaedic Research 21, 1011–1017. Bishop, N.E., Schneider, E., Ito, K., 2003. An experimental two degrees-of-freedom actuated external fixator for in vivo investigation of fracture healing. Medical Engineering and Physics 25, 335–340. Bishop, N.E., Rhinjn, M.v., Tami, I., Corveleijn, R., Schneider, E., Ito, K., 2006. Shear does not necessarily inhibit bone healing. Clinical Orthopaedics and Related Research 443, 307–314. Byrne, D.P., Lacroix, D., Prendergast, P.J., 2011. Simulation of fracture healing in the tibia: mechanoregulation of cell activity using a lattice modeling approach. Journal of Orthopaedic Research, October, 1496–1503. Claes, L.E., Wilke, H.-J., Augat, P., Rübenacker, S., Margevicius, K.J., 1995. Effect of dynamization on gap healing of diaphyseal fractures under external fixation. Clinical Biomechanics 10 (5), 227–234. Claes, L.E., Augat, P., Suger, G., Wilke., H.-J., 1997. Influence of size and stability of the osteotomy gap on the success of fracture healing. Journal of Orthopaedic Research 15, 577–584. Claes, L.E., Heigele, C.A., 1999. Magnitudes of local stress and strain along bony surfaces predict the course and type of fracture healing. Journal of Biomechanics 32, 255–266. Comiskey, D., MacDonald, B.J., McCartney, W.T., Synnott, K., O’Byrne, J., 2013. Predicting the external formation of callus tissues in oblique bone fractures: idealised and clinical case studies. Biomechanics and Modeling in Mechanobiology http: //dxdoi.org/10.1007/s10237-012-0468-6. Chao, E., Aro, H., DGLewallen, Kelly, P., 1989. The effect of rigidity on fracture healing in external fixation. Clinical Orthopaedics and Related Research 241, 24–35. Checa, S., Prendergast, P., 2009. A mechanobiological model for tissue differentiation that includes angiogenesis: a latticebased modeling approach. Annals of Biomedical Engineering 37 (1), 129–145. Draper, E.R.C., Strachan, R.K., Hughes, S.P.F., Nicol, A.C., Paul, J.P., 1997. The design and performance of an experimental external fixator with variable axial stiffness and a compressive force transducer. Medical Engineering and Physics 19 (8), 660–695. Einhorn, T.A., 2005. The science of fracture healing. Journal of Orthopaedic Trauma 19 (10), S4–S6. García-Aznar, J.M., Kuiper, J.H., Gómez-Benito, M.J., Doblaré, M., Richardson, J.B., 2007. Computational simulation of fracture healing: Influence of interfragmentary movement on the callus growth. Journal of Biomechanics 40, 1467–1476. García-Aznar, J.M., Pérez, M.A., Moreo, P., 2009. Modelling of interfaces in biomechanics and mechanobiology. Computer Modeling in Engineering & Sciences: CMES 48 (3), 271–301. Gómez-Benito, M.J., García-Aznar, J.M., Kuiper, J.H., Doblaré, M., 2005. Influence of fracture gap size on the pattern of long bone healing: a computational study. Journal of Theoretical Biology 235, 105–119. Gómez-Benito, M.J., González-Torres, L.A., Reina-Romo, E., Grasa, J., Seral, B., García-Aznar, J.M., 2011. Influence of high-frequency cyclical stimulation on the bone fracture-healing process: mathematical and experimental models. Philosophical Transactions Series A, Mathematical, Physical, and Engineering Sciences 369, 4278–4294. Goodship, A.E., Kenwright, J., 1985. The influence of induced micromovement upon the healing of experimental tibial fractures. Journal of Bone and Joint Surgery (British Volume) 67-b (4), 650–655. Goodship, A.E., Watkins, P.E., Rigby, H.S., Kenwright, J., 1993. The role of fixator frame stiffness in the control of fracture healing. An experimental study. Journal of Biomechanics 26 (9), 1027–1035. Hadjiargyrou, M., McLeod, K., Ryaby, J., Rubin, C., 1998. Enhancement of fracture healing by low intensity ultrasound. Clinical Orthopaedics and Related Research 355 (S), 216–229.
338
journal of the mechanical behavior of biomedical materials 29 (2014) 328 –338
Hernandez, C.J., Beaupre, G.S., Keller, T.S., Carter, D.R., 2001. The influence of bone volume fraction and ash fraction on bone strength and modulus. Bone 29 (1), 74–78. Isaksson, H., Donkelaar, C.C.v., Huiskes, R., Ito, K., 2008. A mechano-regulatory bone-healing model incorporating cellphenotype specific activity. Journal of Theoretical Biology 252, 230–246. Isaksson, H., Donkelaar, C.C.v., Ito, K., 2009a. Sensitivity of tissue differentiation and bone healing predictions to tissue properties. Journal of Biomechanics 42, 555–564. Isaksson, H., Gröngröft, I., Wilson, W., Donkelaar, C.C.v., Rietbergen, B.v., Tami, A., Huiskes, R., Ito, K., 2009b. Remodeling of fracture callus in mice is consistent with mechanical loading and bone remodeling theory. Journal of Orthopaedic Research 27 (5), 664–672. Jagodzinski, M., Krettek, C., 2007. Effect of mechanical stability on fracture healing – an update. Injury 38 (S1), S3–S10. Klein, P., Opitz, M., Schell, H., Taylor, W.R., Heller, M.O., Kassi, J.-P., Kandziora, F., Duda, G.N., 2004. Comparison of unreamed nailing and external fixation of tibial diastases – mechanical conditions during healing and biological outcome. Journal of Orthopaedic Research 22, 1072–1078. Lacroix, D., Prendergast, P.J., 2002a. Three-dimensional simulation of fracture repair in the human tibia. Computer Methods in Biomechanics and Biomedical Engineering 5 (5), 369–376. Lacroix, D., Prendergast, P.J., 2002b. A mechano-regulation model for tissue differentiation during fracture healing: analysis of gap size and loading. Journal of Biomechanics 35, 1163–1171. Loboa, E.G., Beaupré, G.S., Carter, D.R., 2001. Mechanology of initial pseudarthrosis formation with oblique fractures. Journal of Orthopaedic Research 19, 1067–1072. Manjubala, I., Liu, Y., Epari, D.R., Roschger, P., Schell, H., Fratzl, P., Duda, G.N., 2009. Spatial and temporal variations of mechanical properties and mineral content of the external callus during bone healing. Bone 45, 185–192. Marsell, R., Einhorn, T., 2011. The biology of fracture healing. Injury 42 (6), 551–555. Moreo, P., Pérez, M.A., García-Aznar, J.M., Doblaré, M., 2007. Modelling the mechanical behaviour of living bony interfaces.
Computer Methods in Applied Mechanics and Engineering 196, 3300–3314. Park, S.-H., O’Connor, K., McKellop, H., Sarmiento, A., 1998. The influence of active shear or compressive motion on fracture healing. Journal of Bone and Joint Surgery (American Volume) 80-A (6), 868–878. Perren, S.M., 2002. Evolution of the internal fixation of long bone fractures. Journal of Bone and Joint Surgery (British Volume) 84-B (8), 1093–1110. Reina-Romo, E., Gómez-Benito, M.J., Domínguez, J., Niemeyer, F., Wehner, T., Simon, U., Claes, L.E., 2011. Effect of the fixator stiffness on the young regenerate bone after bone transport: computational approach. Journal of Biomechanics 44, 917–923. Richardson, J.B., Cunningham, J.L., Goodship, A.E., O’Connor, B.T., Kenwright, J., 1994. Measuring stiffness can define healing of tibial fractures. Journal of Bone and Joint Surgery (British Volume) 76-B (3), 389–394. Schell, H., Epari, D.R., Kassi, J.P., Bragulla, H., Bail, H.J., Duda, G.N., 2005. The course of bone healing is influenced by the initial shear fixation stability. Journal of Orthopaedic Research 23, 1022–1028. Shefelbine, S.J., Augat, P., Claes, L., Simon, U., 2005. Trabecular bone fracture healing simulation with finite element analysis and fuzzy logic. Journal of Biomechanics 38, 2440–2450. Simon, U., Augat, P., Claes, L., Utz, M., 2011. A numerical model of the fracture healing process that describes tissue development and revascularisation. Computer Methods in Biomechanics and Biomedical Engineering 14 (1), 79–93. Vetter, A., Epari, D.R., Robin Seidel, H.S., Fratzl, P., Duda, G.N., Weinkamer, R., 2010. Temporal tissue patterns in bone healing of sheep. Journal of Orthopaedic Research 28 (11), 1440–1447. Wehner, T., Claes, L., Niemeyer, F., Nolte, D., Simon, U., 2010. Influence of the fixation stability on the healing time—a numerical study of a patient-specific fracture healing process. Clinical Biomechanics 25, 606–612. Wehner, T., Penzkofer, R., Augat, P., Claes, L., Simon, U., 2011. Improvement of the shear fixation stability of intramedullary nailing. Clinical Biomechanics 26, 147–151.
