Combined finite element model of human proximal femur behaviour considering remodeling and fracture

Disponible en ligne sur

www.sciencedirect.com IRBM 34 (2013) 191–195

Original article

Combined finite element model of human proximal femur behaviour considering remodeling and fracture R. Hambli a,∗ , C.-L. Benhamou b , R. Jennane a , E. Lespessailles b , W. Skalli c , S. Laporte c , J.-D. Laredo d , V. Bousson d , J. Zarka e a Prisme Institute/MMH, 8, rue Léonard-de-Vinci, 45072 Orléans cedex 2, France Inserm U658, IPROS, CHR d’Orléans, 1, rue Porte-Madeleine, 45032 Orléans cedex 1, France c Arts et Métiers ParisTech, laboratoire de biomécanique, Paris, France d CNRS, laboratoire de recherches orthopédiques, université Paris 7, UMR CNRS 7052, B2OA, 10, avenue de Verdun, 75010 Paris, France e CADLM, 43, rue du Saule-Trapu, 91300 Massy, France b

Received 14 January 2013; received in revised form 17 January 2013; accepted 17 January 2013 Available online 16 March 2013

Abstract The purpose of this work was to develop a combined remodeling-to-fracture finite element model allowing for the combined simulation of human proximal femur remodeling under a given boundary conditions followed by the simulation of its fracture behaviour under quasi-static load. The combination of remodeling and fracture simulation into one unified model consists in considering that the femur properties resulting from the remodeling simulation correspond to the initial state for the fracture prediction. The remodeling model is based on a coupled strain and fatigue damage stimulus approach. The fracture model is based on continuum damage mechanics in order to predict the progressive fracturing process, which allows to predict the fracture pattern and the complete force-displacement curve under quasi-static load. To investigate the potential of the proposed unified remodeling-to-fracture model, we performed remodeling simulations on a 3D proximal femur model for a duration of 365 days followed by a side fall fracture simulation reproducing. © 2013 Published by Elsevier Masson SAS.

1. Introduction Healthy human femur adapts its strength and mass to its mechanical use since it endures daily loads. Osteoblasts and osteoclasts cells respectively add and resorb bone during the remodeling process as a response to the mechanical stimulus sensed by osteocytes [1–3]. Results of bone remodeling (mass, microarchitecture and bone material properties of both cancellous and cortical bone) control bone stiffness and strength, which can be assessed by the ultimate force provoking fracture at organ level [4,5]. In clinical practice, current fracture prediction methods are based on the mineral content of bone at a particular instant and do not account for its adaptation, which will occur over time. The prediction of fracture risk could therefore be under or over estimated thereby leading to inaccurate prognoses and

∗

Corresponding author. E-mail address: [email protected] (R. Hambli).

1959-0318/$ – see front matter © 2013 Published by Elsevier Masson SAS. http://dx.doi.org/10.1016/j.irbm.2013.01.011

unnecessary costs. Assessing the mechanical properties of human femur due to a given remodeling cycle and thus predicting its failure is still a challenge that can be addressed using the finite element (FE) method. Over the last years, a large number of FE models have been developed to predict separately bone remodeling and bone fracture [6,7]. Cristofolini et al. [8] reported that the development of enhanced multiscale bone fracture prediction models, requires coupling separated models including remodeling and fracture into a full multiscale model which is still a work in progress. Such FE models would allow clinicians to directly predict bone strength, estimate the risk of fracture and implement relevant preventative treatments for patients. The aim of the current work was to develop an combined R-F FE model to simulate the bone remodeling process under given boundary conditions and duration of time followed by its fracture simulation under a given boundary conditions (side fall, onelegged stance load) in the quasi-static regime (low strain rate). To check the validity of the combined model, we performed remodeling simulation on a 3D proximal femur model for a

192

R. Hambli et al. / IRBM 34 (2013) 191–195

duration of 365 days under daily loading conditions followed by a side fall fracture simulation reproducing an experimental test performed on human femur under a side fall configuration and compared predicted and experimental fracture patterns.



2. Material and methods



Dfat = 1 − 1 −

Proposed combined R-F model is built to perform simulations in two steps (Fig. 1): • step 1: FE simulation of remodeling process is performed under given boundary conditions (load amplitude, angle, frequency and duration of time). The net results of the remodeling is the bone adaptation in the form of heterogeneous distribution of bone material properties; • step 2: FE simulation of fracture under excessive load in order to predict the force-displacement curve, the ultimate force at fracture and the fracture pattern corresponding to the femur mechanical and structural states at the end of remodeling simulation of step 1. The initial conditions of the fracture simulation (material properties, mineralization, fatigue damage and density) correspond to those predicted at the end of the remodeling simulation results. 2.1. Remodeling model The bone adaptation approach used in this study is based on Huiskes et al. [9] theory for its simplicity. The continuum behaviour of bone tissue is modelled using the concept of continuum damage mechanics (CDM) as elastic isotropic behaviour coupled to fatigue damage in the form [10]:   σij = 1 − Dfat aijkl εkl

σ ij is the stress, Dfat is the fatigue damage variable, εkl is the strain and aijkl is the isotropic elasticity stiffness tensor. For high cycle fatigue under purely elastic strain, a nonlinear damage model can be expressed by [10–12]:

(1)

N Nf



1 1−γ

1  1+β

(2)

where: γ and β are material parameters and Nf is the cycle at failure which can be obtained as described by Martin et al. [13]: Nfc = 1.479 × 10−21 Δε−10.3 for compressive loads

(3a)

Nft = 3.630 × 10−32 Δε−14.1 for tensile loads

(3b)

where ε is the amplitude of the applied microstrain. The set of equation describing the fully coupled change in the bone density is given by [10,14]: dρ = αR (S − SR ) dt dρ =0 dt

If S < SR

If SR ≤ S ≤ SF

(4a) (4b)

dρ = αF (S − SF ) dt

If SF < S < SD

(4c)

dρ = αD (S − SD ) dt

If D ≥ Dcfat

(4d)

where ρ is the bone density, t is the time and S is the coupled strain-damage stimulus function. αR , αF and αD denote respectively bone resorption rate, bone formation rate and damage resorption rate. SR , SF and SD denote respectively setpoints for bone resorption, formation and fatigue damage repair.

Fig. 1. Overview of the combined remodeling-to-fracture simulation model in 2 steps implemented into Abaqus code (Abaqus 6.11).

R. Hambli et al. / IRBM 34 (2013) 191–195

The coupled damage-strain remodeling mechanical stimulus is expressed [10]:  1 1 − Dfat σij εij (5) Si = 2 ρ The new density value of the bone tissue is approximated using forward Euler method by: ρt+Δt = ρt + Δρ

(6)

The elastic modulus E at each bone location is calculated according to [10]:   (7) E = C 1 − Dfat ρp where C is experimentally derived constants, (2 ≤ p ≤ 3) [15].

⎧ D=0 ⎪ ⎪ ⎨ D = k εneq ⎪ ⎪ ⎩ D = Dc

(8)

; εeq ≤ εy ; εy < εeq < εf

(9)

; εeq ≥ εf

εeq , k, n, εy εf and Dc are respectively the equivalent strain, the damage factor, the damage exponent, yield strain (damagestrain threshold when damage starts), the strain at fracture and the critical damage at fracture. The equivalent strain εeq is expressed by: εeq =

In the quasi-static regime, the isotropic stress-strain relation of elasticity based damage mechanics is expressed by [7,16]: σij = (1 − D) Cijkl εkl

where D denotes the damage variable, σ ij the stress components, εkl the strains and Cijkl are the elasticity tensor components. It have been reported by other studies that the failure process of bone is controlled by strain [17–19]. Based on these results, an experimentally damage law can be expressed in the general form:



2.2. FE simulation of bone fracture initiation and propagation under quasi-static load

193

2 εij εij 3

(10)

εij are the strain tensor components. Material properties for bone remodeling and fracture and their sources are given in Table 1.

Table 1 Material properties and model parameters for bone used for the remodeling and fracture simulation. Parameters

Notation

Trabecular

Cortical

References

Initial density Initial elastic modulus Poisson ratio Density coefficient Density exponent Fatigue parameter Fatigue exponent Critical fatigue damage for bone repair Threshold strain value for damage initiation Critical damage value at fracture in tension Critical damage value at fracture in compression Strain at fracture in tension Strain at fracture in compression

ρ (g/cm3 ) E0 (MPa) v C (g/cm3 ) p γ β fat Dc ε0 (%) pl Dt pl Dc T εf εC f

0.764 5690 0.3 2014 2.5 0.2 0.4 0.4 0.1 0.95 0.5 1.57 2.5

0.764 9882 0.3 1763 3.2 0.2 0.4 0.4 0.1 0.95 0.5 2.5 4

[7] [10] [13] [10] [12] [7] [12] [12] [7] [7]

Fig. 2. Transition procedure from remodeling-to-fracture simulation: (i) remodeling simulation under a given boundary conditions and cycle; (ii) adaptation of bone properties at the end of the remodeling cycle; (iii) side fall fracture simulation considering the end state of the bone properties after remodeling convergence as initial condition for the fracture simulation; (iv) fracture simulation results (force-displacement curve and fracture pattern).

194

R. Hambli et al. / IRBM 34 (2013) 191–195

Fig. 3. Typical predicted force-displacement curve after remodeling phase. At point (A), elastic deformation initiates. At point (B), cracks initiate and ultimate force is reached. Cracks propagate and failure occurs at point (C). At point (D), the proximal femur is completely broken.

3. Results In the current study, we performed remodeling simulations on a 3D proximal femur model for a duration of 365 days under different daily loading condition followed by a side fall fracture simulation reproducing an experimental fracture test under a side fall configuration. 3.1. Remodeling results An example of the combined remodeling-to-fracture prediction is shown in Fig. 2. Density results were similar to regional densities observed clinically [20] and the density distribution was consistent with many features observed in femoral morphology: trabecular bone of varying density in the head and trochanter, and dense cortical bone in the diaphysis and calcar region of the neck. An example of the predicted force-displacement curve (upon remodeling) is provided in Fig. 3. Predicted proximal femur fracture patterns corresponding to different cracks propagation levels in relation with the force-displacement curve are also reported. Comparison between predicted and experimental complete fracture pattern is given in Fig. 4. One can notice that a very good

agreement is obtained between the predicted and experimental patterns. The proposed combined remodeling-to-fracture FE approach is one of the first models which combines two steps (a) remodeling and (b) progressive fracture simulation of human proximal femurs. The present study illustrates the potential of an integrated model combining remodeling simulation under a given loading cycles followed by a fracture prediction to investigate the effects of bone remodeling on bone strength. We can summarize that the proposed bone remodeling-to-fracture models can be considered a first tool to assess the proximal femur strength after a given remodeling conditions (loading, frequency and duration). This can be very useful for clinicians to prevent fracture risk to osteoporotic patients and to assess whether a given patient will undergoes fracture under some loading conditions (side fall, spontaneous. . .) and implement relevant preventative measures in individual patients. Before potential clinical implementation, more sophisticated mechanobiological models taking into account the bone cells activities must be developed and further validation is needed by performing several experiments with different femur samples and boundary conditions (fall, stance. . .).

Fig. 4. An example of experimental and finite element (FE) predicted complete fracture patterns of a proximal femur.

R. Hambli et al. / IRBM 34 (2013) 191–195

Acknowledgements This work has been supported by French National Research Agency (ANR) through TecSan program (Project MoDos, n◦ ANR-09-TECS-018). References [1] Parfitt AM. Targeted and nontargeted bone remodeling: relationship to basic multicellular unit origination and progression. Bone 2002;30:5–7. [2] Hambli R, Rieger R. Physiologically based mathematical model of transduction of mechanobiological signals by osteocytes. Biomech Model Mechanobiol 2011;11:1–2 [p. 83–93]. [3] Rieger R, Hambli R, Jennane R. Modeling of biological doses and mechanical effects on bone transduction. J Theor Biol 2011;274(1):36–42. [4] Rho JY, Kuhn-Spearing L, Zioupos P. Mechanical properties and the hierarchical structure of bone. Med Eng Phys 1998;20:92–102. [5] Tai K, Dao M, Suresh S, Palazoglu A, Ortiz C. Nanoscale heterogeneity promotes energy dissipation in bone. Nat Mater 2007;6:454–62. [6] Christen D, Webster JD, Müller R. Multiscale modelling and nonlinear finite element analysis as clinical tools for the assessment of fracture risk. Philos Trans R Soc A 2010;368:2653–68. [7] Hambli R. A quasi-brittle continuum damage finite element model of the human proximal femur based on element deletion. Med Biol Eng Comput 2013 [On line]. [8] Cristofolini L, Taddei F, Baleani M, Baruffaldi F, Stea S, Viceconti M. Multiscale investigation of the functional properties of the human femur. Philos Transact A Math Phys Eng Sci 2008;28:366 [1879;33:19–41]. [9] Huiskes R, Ruimerman R, van Lenthe GH, Janssen JD. Effects of mechanical forces on maintenance and adaptation of form in trabecular bone. Nature 2000;405:704–6.

195

[10] Hambli R, Soulat D, Gasser A, Benhamou CL. Strain-damage coupled algorithm for cancellous bone mechano-regulation with spatial function influence. Comput Methods Appl Mech Eng 2009;198:33–6 [1;2673–82]. [11] Hambli R. Multiscale prediction of crack density and crack length accumulation in trabecular bone based on neural networks and finite element simulation. Int J Numerical Methods Biomed Eng 2011;27(4):461– 75. [12] Hambli R. Apparent damage accumulation in cancellous bone using neural networks. J Mech Behav Biomed Mater 2011;4(6):868–78. [13] Martin RB, Burr DR, Sharkey NA. Skeletal Tissue Mechanics. New York: Springer; 1998. [14] McNamara LM, Prendergast JP. Bone remodeling algorithms incorporating both strain and microdamage stimuli. J Biomechanics 2007;40:1381–91. [15] Hernandez CJ, Beaupré GS, Carter DR. A model of mechanobiologic and metabolic infuences on bone adaptation. J Rehabil Res Dev 2000;37(2):235–44. [16] Hambli R, Bettamer A, Allaoui S. Finite element prediction of proximal femur fracture pattern based on orthotropic behaviour law coupled to quasibrittle damage. Med Eng Phys 2012;34(2):202–10. [17] Nagaraja S, Couse TL, Guldberg RE. Trabecular bone microdamage and microstructural stresses under uniaxial compression. J Biomech 2005;38:707–16. [18] Nalla K, Kruzic J, Kinney H, Ritchie O. Mechanistic aspects of fracture and R-curve behavior in human cortical bone. Biomaterials 2005;26:217–31. [19] Niebur GL, Feldstein MJ, Yuen JC, Chen TJ, Keaveny TM. High-resolution finite element models with tissue strength asymmetry accurately predict failure of trabecular bone. J Biomech 2000;33:1575–83. [20] Link TM, Vieth V, Langenberg R, Meier N, Lotter A, Newitt D, et al. Structure analysis of high-resolution magnetic resonance imaging of the proximal femur: in vitro correlation with biomechanical strength and BMD. Calcif Tissue Int 2003;72(2):156–65.