In order to develop a predictive model which takes into account the mechanical parameters involved in the fracture process, a correlation between these magnitudes and those ones measured in clinical terms is firstly required. To this effect, Carter and Hayes [31
] established a direct relationship between Young's modulus and bone density for low strain rates (0.01):
On the other hand, the relationship between the BMD value and the apparent density is adjusted, according to experimental results, as:
being λ a parameter that depends on the sample data. For the present study a value of 9/25 was assigned for λ, in order to fit the actual data presented in [32
BMD is the current standard for diagnosis of osteoporosis. Over 100 worldwide published papers assessing hip BMD evolution, both in natural conditions and in patients under drug therapy were selected for analysis. Among them, three treatments have been selected for the comparative study: alendronate 10 mg per day [33
], oral ibandronate 2.5 mg per day [34
] and PTH 1–84 100 mg per day [35
]. These drugs have proven to be effective, and their BMD evolution curves show the required regularity in the analyzed time period.
Bisphosphonates settle in the bone tissue and its effect persists during some time after its administration, and hence it is completely accepted in the clinical practice the use of intermittent or discontinuous treatments. In most of the studies a maximum five years follow-up is done, although a treatment continuous for ten years is also accepted. Moreover, there are several studies concerning the safety of long term treatments (ten years) [36
]. With respect to the use of combined therapies (teriparatide and bisphosphonates), it is a usual clinical practice in the osteoporosis treatment [43
Regarding the natural history of BMD, the average curve published by Mazess [48
] was chosen as a reference. Briefly, the natural curve of BMD linked with age was compared to the curves of BMD in patients treated with bisphosphonates or derivatives of parathyroid hormone.
Since the standard adjustment techniques do not provide enough accuracy and reliability to be applied to the predictive model, higher level, more complex, adjustments have been set out to obtain continuous curves of regression. These make it possible to obtain continuous curves for the BMD evolution. The evolution trends of the different selected cases in the study can be extrapolated from those continuous curves, ensuring a consistent behavior. The following adjustments were proposed:
These adjustments (Equations 3, 4 and 5) have been applied to the four examined curves (natural evolution and the three therapies). Several choices have been made during this adjustment: the lowest mean square error, the closest to unity correlation coefficient R2
, and a 10-year standardized follow-up period with the possibility of extrapolation to longer periods up to 15 years (Figure ). That extrapolation is intended to be consistent in mathematical terms, independently of the actual duration of the considered study. This condition will permit in the future an easy incorporation of more long term studies when available. Standardization of BMD measurements included in each published paper was required in order to make comparisons among them [49
], because de data were obtained from different densitometric equipment depending on the study (Hologic, Lunar and Norland).
s Damage Mechanics model [50
] was selected for the simulation of the degenerative process. This model defines the mechanical damage, D
, as a function of the equivalent strain, εc
the equivalent strain, where εI, εII, εIII are the principal components of the strain tensor. This mechanical damage leads to a decrease of the bone stiffness, according with the relationship:
corresponds to the Young’s modulus of healthy bone and E
is the actual value of Young’s modulus of the bone with cumulative mechanical damage. The κ
constants depend on the critical damage value, the critical strain and the strain threshold, ε0
, below which no damage is produced (εc
corresponds to the equivalent strain value that produces the critical damage, Dcri
. Once the damage model is defined, a relationship between the damage level and the probability of fracture must be set. Obviously several factors are involved in the equation and may difficult the calculations: damage location, amount of damage, range and type of load cycles, and so on. However, since osteoporosis is a generalized disease affecting bone mass extensively, and it occurs in older people with little variation in their life habits, a simplified model could be used. In this work, a model of fracture probability based on the law of Paris [51
], which explains the stable crack growth under monotonic loading, was developed. The number of loading cycles needed to increase damage from the Di
value to the Df
value was expressed as follows:
are material parameters associated with the law of Paris, γ
is a parameter related to the stress intensity factor defined in Linear Elastic Fracture Mechanics, which depends on geometric aspects and stress distribution, Δσ
is the variation of stress in each loading cycle and ω
is a parameter which relates the mechanical damage to the crack size a
By normalizing the probability of fracture, assigning value 1 to critical damage (D
and value 0 to no damage (D
), we finally obtain the probability of fracture as a function of both the number of cycles and the damage:
In Equation (11), the magnitude Nmax
represents the maximum number of load cycles necessary for a mechanical damage equivalent to the critical value, Dcri
. A value 5
has been given to the β
coefficient for cortical bone, according to Taylor [52
], and, in accordance with Kargarnovin [50
], a critical damage, Dcri
, of 0
and a critical strain, εcri
, of 0
have been considered, with a strain threshold, ε0
, of 0
(no damage is produced by strains below 0
Progression curves are basic references and provide only general information. In order to apply the model to specific patients we must consider, in addition to the trend, the reference density value of the patient. According to the law of interpolated natural evolution, the density matching with the age of the patient is given by:
being ρN(t0) the average density measured in the patient and ρ0 the matching value from the reference progression curve, that is, the actual patient’s density doesn’t necessarily be equal to the value corresponding to the average curve for the considered population. There is an offset that should be added to the reference progression curve of density, providing a translation of the average curve, allowing an adaptation for each individual patient, so:
When a drug therapy is applied to the same patient, a similar correction is required because a new offset arises. So, the progression curve for this patient under a treatment would be:
In the Equations (12) to (14), the subscripts or superscripts N and T represent natural evolution or evolution with treatment, respectively. All these adjustments provide the estimated BMD value for any type of patient, at any age, and under any prescribed therapy. From this value, mechanical properties of bone can be calculated. It must be noticed that the considered curves represent the mean evolutionary curves for the population, and an individual patient could not follow the curve exactly but in an approximate way.
Densitometric data of the healthy femur have been taken as the starting point for this study, based on a previous work of our group on the biomechanical behavior of a femoral stem [32
]. In that survey, densitometric data were correlated both with the apparent volumetric density, and with the Young moduli of each of the Gruen zones. Table depicts the BMD, the apparent volumetric density, and the Young moduli correlation corresponding to a healthy femur.
Table 1 BMD, apparent density and Young’s modulus for the data corresponding to the study of Herrera
From all the previous calculations, an evolutionary algorithm has been implemented (Figure ), which has been used combined with a finite element model of the femur [32
] using the Abaqus software [53
]. The procedure was based on a finite element model of the upper half of the femur, made up of tetrahedral elements with quadratic approximation (C3D10 in Abaqus nomenclature), as shown in Figure and b. To generate the model a 3D laser scanner Roland PIZCA was used for the outer geometry and thirty transverse direction tomographic cross-sections and eight longitudinal direction cross-sections were taken using CAT (General Electric Brightspeed Elite) to determine the geometry of the cancellous bone. Moreover a 3D reconstruction was made to obtain the shape of the medullar cavity [32
Evolutionary algorithm for the prediction of fracture probability.
3D FE model of proximal femur: a) Adaptation of Gruen zones to the healthy femur; b) Boundary conditions on the proximal femur.
In the final mesh, consisting of 408094 elements and 75223 nodes, cortical bone, cancellous bone and bone marrow have been differentiated (229931 elements for cortical bone, 166220 elements for cancellous bone and 11943 elements for bone marrow). To guarantee the accuracy in the FE results, a sensitivity analysis was performed with a mesh refinement in order to achieve a convergence towards a minimum of the potential energy, with a tolerance of 1% between consecutive meshes.
In addition, the femur has been divided into 7 zones of Gruen [54
], since the available information on BMD was related to these zones (Figure ). The mechanical properties obtained from the evolutionary algorithm have been assigned to each of the Gruen zones. For this purpose, the apparent density is calculated for each Gruen zone from its BMD measure according Equation (2), and then, the Young’ modulus is obtained from Equation (1). The same process was followed for cancellous bone. Any case, cortical bone stresses are slightly influenced by the mechanical properties of cancellous bone and bone marrow since their stiffness is much lower than the corresponding to cortical bone. The function of cancellous bone, in mechanical terms, consists in providing stability to the thin walled cortical bone. A linear elastic behavior has been defined for all the materials. The final collection of results was focused on the proximal femur, mainly on the area located within the edge of the femoral head and the subtrochanteric region.
As a boundary condition, the middle third of the femoral diaphysis was clamped (Figure ), since this area is far enough from the proximal femur to avoid significant perturbations in the stress distribution. Thus, computational cost was reduced if compared with the entire femur model.
In the finite element model construction, as important as the reaction strength on the femoral head due to the body weight, is the inclusion of the muscle forces to be considered in the simulation. In our model only the abductor muscle forces were included, in line with several authors [55
]. As a rule of thumb, abductor muscles cause a reaction on the femoral head of 2.75 times the body weight. However, load increases as much as 4 times the body weight when the heel impacts the ground, and during the double support stage of the gait [57
]. The latter situation was considered, as the worst case, in order to set the boundary conditions. According to data from the study [32
] a 79.3 kg body-weight was set as a reference. So, two load conditions were imposed (Figure ):
• Reaction strength on the femoral head due to the body weight (3110 N).
• Load due to the abductor muscles, applied to the proximal area of the greater trochanter (1360 N).
The above values are slightly higher than those included in recent studies [58
]. This model has been used in predicting the evolution of femoral fracture probability, by comparing the natural history and the expected evolution under different therapies.