In silico dynamic characterization of the femur: Physiological versus mechanical boundary conditions

1 In silico dynamic characterization of the femur: physiological versus 1 mechanical boundary conditions 2 E. Reina-Romo, J. Rodríguez-Vallés, J.A. Sanz-Herrera 3 School of Engineering, University of Seville, 41092 Seville (Spain) 4 5 Abstract 6 It is established that bone tissue adapts and responds to mechanical loading. Several 7 studies have suggested an existence of positive influence of vibration on the bone 8 mass maintenance. Thus, some bone regeneration therapies are based on vibration of 9 bone tissue under circumstances of disease to stimulate its formation. Frequency of 10 loading should be properly selected and therefore a correct characterization of the 11 dynamic properties of this tissue may be critical for the success of such orthopedic 12 techniques. On the other hand, many studies implement vibration techniques with in 13 silico models. Numerical results are exclusively dependent on properties of bone 14 tissue, i.e. geometry, density distribution and stiffness, as well as boundary conditions. 15 In the present study, the influence of boundary conditions and material properties on 16 the dynamic characteristics of bone tissue was explored in a human femur. Bone shape 17 and density were directly reconstructed from computer tomographies, whereas 18 natural frequencies and modes of vibration were obtained for different boundary 19 conditions including physiological and mechanical ones. Results of this study show the 20 moderate effect of material properties compared to the much substantial effect of 21 boundary conditions. A factor of 2 in the natural frequency was obtained depending on 22 imposed boundary conditions, highlighting the importance in the selection of 23 appropriate conditions in the analysis of the bone organ. 24 Keywords: Bone Mechanics, Finite Element Method, Natural Frequency, Modal 25 Analysis. 26 27 1. Introduction 28 It is established that bone tissue adapts and responds to mechanical loading. The 29 characteristics of the mechanical loads, such as peak magnitude of ground reaction 30 force, peak rate of force production, and repetition rate of these loads, are known to 31 be important in regulating bone regeneration [1]. According to experimental studies, 32 maximum bone formation directly depends on loading frequencies and targeted 33 locations of bones [2-6]. Zhao et al. [7] have demonstrated in an in vivo study 34 performed on mice that the application of loads with frequencies near the natural 35 frequency of bone favors bone formation. Rubin et al. [8] have shown that very-low36 This is an Accepted Manuscript of an article published by Elsevier in Medical Engineering and Physics on August 2018, available at: https://doi.org/10.1016/j.medengphy.2018.06.001 Copyright 2018 Elsevier. En idUS: Licencia Creative Commons CC BY-NC-ND 2 magnitude, frequency vibration of 30Hz, increases significantly (by 34.2 %) the density 37 of the trabecular bone in the proximal human femur compared to controls. Therefore, 38 a possible application to increase bone mass under circumstances of disease would be 39 to apply vibration loading stimulation: whole-body vibration, bone bending, axial 40 loading, and joint loading [10-13]. Whole-body vibration seems to prevent/reverse 41 sarcopenia and possibly osteoporosis [9]. Vibration analysis techniques have been 42 used successfully to determine mechanical properties of human bone [14], to monitor 43 fractures [15-16], to quantify the stability of dental implants [17-18] and to detect 44 forms of femoral prosthesis loosening [19]. In the orthopaedic field, this technique 45 may have potential due to apparent benefits it offers. 46 The vibration characteristics of bone have been analyzed both experimentally [10,2047 21] and numerically [7,14,19,22-23]. Although Weiss et al. [24] have shown the 48 importance of the boundary conditions in vibration studies in a hip endoprothesis 49 system, the most commonly used boundary conditions are the free ones [14,22], 50 which are known to be far from the in vivo environment. For example, Campoli et al. 51 [22] already analyzed by means of finite element analysis how the natural frequencies 52 changed by varying the density and shape of the human femur in free boundary 53 conditions. The behavior of the femoral head prosthesis has also been analyzed in 54 silico with modal analyses assuming free boundary conditions and constant mechanical 55 properties [23]. Pérez & Seral-García [19] used other boundary conditions to simulate 56 numerically the change in the resonance frequency during the osseointegration 57 process of a cementless human hip system. In that study mechanical properties were 58 assigned based on the level of Hounsfield Unit (HU). 59 A number of studies have already shown that geometry and material properties at 60 boundaries and spatiotemporal distributions within bone have a profound effect on 61 bone mechanobiology variables such as fluid velocities and pore pressures [25-26]. 62 However, as far as the authors know, there is no vibration study that has numerically 63 analyzed the influence of the properties of bone tissue, i.e. geometry, density 64 distribution and stiffness, as well as boundary conditions on the dynamic 65 characteristics of bone. Therefore, the aim of this study is to analyze in silico the 66 dynamic characteristics of a human femur under different boundary conditions and 67 assumptions found in the literature. Different numerical set-ups of the femur, which 68 differed in the material properties assignment and boundary conditions, have been 69 analyzed. On the one hand, either constant material properties or a technique based 70 on HU was used. On the other hand, four different boundary conditions are 71 implemented and analyzed: free, diaphysis, condyle and physiological boundary 72 conditions with the purpose of determining the dynamic characteristics of the femur 73 as closely as possible as they can be found in vivo. Using the information provided in 74 this study, a proper determination of bone resonance frequencies could be conducted, 75 which in turn may be of special relevance to calibrate existing vibration therapies 76 3 performed clinically and to help to improve the treatment of some diseases such as 77 the osteoporosis. 78 2. Materials and Methods 79 80 2.1. Geometry 81 A set of 664 computer axial tomographies (CATs) of a cadaveric adult male human 82 femur were selected. CAT resolution includes a pixel size of 0.78x0.78 mm and a slice 83 thickness of 0.8 mm. Geometry was built by means of the software MIMICS 10.0® as 84 follows: outer (cortical) and inner (marrow) regions were selected manually at 85 diaphyseal region attending to the grey scale level of pictures, as well as spongy bone 86 at proximal areas. Other regions of non interest in the study, such as hip, knee or tibia 87 were not included in the geometry. Then, 3D geometry of the right femur was 88 reconstructed by using the above mentioned software. 89 2.2 Boundary conditions 90 Four different boundary conditions are implemented and analyzed in the paper, 91 namely, free, diaphysis (fig. 1a), condyle (fig. 1b) and physiological boundary 92 conditions (fig. 1c) [27]. They are explained next. 93 Free boundary conditions 94 In this case, femur was free of loading and displacement restrictions. This is the most 95 suitable boundary conditions used both in experimental modal analysis [14,20] as well 96 as numerical modal analysis of the bone tissue [14,22-23], due to both reliability and 97 simplicity at the laboratory level and finite element analysis. 98 Diaphysis boundary conditions 99 Displacements are prescribed at nodes in the mid-diaphysis where the cut was 100 performed (see fig. 1a) [27,28]. Both nodes at exterior (cortical) and interior (marrow) 101 regions were considered. 102 Condyle boundary conditions 103 All the displacements are prescribed at the three nodes shown in fig. 1b in the distal 104 condyles. 105 4 (a) (b) (c) (d) (e) (f) 106 Figure 1. Finite element model of different implemented boundary conditions: (a) diaphysis, 107 (b) condyle and (c) physiological. Tetrahedral finite element mesh at (d) the exterior (cortical), 108 (e) the interior (marrow and spongy bone) regions; (f) finite element mesh for the diaphysis 109 boundary conditions analysis. 110 Physiological boundary conditions 111 The so-called physiological conditions are the ones proposed by Speirs et al. [27]. In 112 that work, authors conclude that this set of boundary conditions is the one that best 113 represents the femur deformations at real conditions. They are implemented as 114 follows (fig. 1c). First, the three translational degrees of freedom at a node placed in 115 the knee centre are prescribed, specifically at the joint of the tibia with the femur. 116 5 Then, two displacements at a node in the femur head are also prescribed. At this node, 117 only the displacement in the direction along an axis towards the knee center is 118 allowed. Finally, the sixth degree of freedom was constrained at a node placed on the 119 distal lateral epicondyle in order to avoid rigid body motions (see fig. 1c). 120 2.3 Finite element mesh 121 Femur geometry was exported into the finite element software ANSYS 14.5®. 122 Tetrahedral 8 node 3D quadratic elements (SOLID 185) were used including 60227 123 elements and 12020 nodes in the femur mesh (see figs. 1d and e). 124 For diaphysis boundary conditions, a different finite element model was used (see fig. 125 1f). In this case a cut was performed at the mid-diaphysis [27] (see fig. 1f). 126 2.4 Mechanical properties of bone tissue 127 Linear, elastic and isotropic mechanical properties were considered in this study 128 [14,22]. Material properties were characterized by Young’s modulus and Poisson’s 129 ratio. Nonetheless, femur domain was considered to be heterogeneous, such that each 130 finite element in the model was featured by different mechanical parameters. 131 Mechanical properties have been assigned with two different approaches: depending 132 on the level of HU, regardless its location and, manually, considering three different 133 material sets. 134 Firstly, the technique based on HU was implemented. Bone density was linearly 135 correlated with HU from 0,5 g cm3 ⁄ of the spongy bone to 1,952 g cm3 ⁄ of the cortical 136 bone [22]. Mathematically, this relation attends to the following curve: 137 ρ (g cm3 ⁄) = 0,000968 ×HU + 0,5 (1) 138 where 𝜌 is the bone tissue density and HU is the Hounsfield Unit, which varied from 0 139 to 1500. Bone density is related to the mechanical properties following Beaupre et al. 140 [29] as: 141 If ρ ≤ 1,2 g cm3 ⁄ ; E = 2014ρ2,5 (MPa), ν = 0,2 (2) 142 If ρ > 1,2 g cm3 ⁄ ; E = 1763ρ3,2 (MPa), ν = 0,32 143 where E and 𝜈 are the Young’s modulus and Poisson’s ratio, respectively. A set of 10 144 values of density, uniformly distributed from lower (0,5 g cm3 ⁄) to upper (1,952 145 g cm3 ⁄) bounds, respectively, were considered. Consequently, a number of 10 146 different pairs (E, 𝜈) were estimated for the mechanical properties along the bone 147 tissue through Eq. (2) according to its estimated density. This was included as different 148 materials in the finite element mesh by using utilities of MIMICS 10.0® 149 and ANSYS ICEM 14.5® softwares. 150 6 Secondly, the finite element model was tested assuming constant mechanical 151 properties. Three different material sets were distinguished: cortical bone, bone 152 marrow and spongy bone (fig. 1e). In this case, mechanical properties are given in 153 table 1. 154 Material Density (𝐤𝐠 𝐦𝟑 ⁄) Young’s modulus (GPa) Poisson’s ratio (-) Cortical bone [22] 1800 13 0,3 Bone marrow [29] 1060 0,001 0.5 Spongy bone [22] 500 0,6 0,12 155 Table 1. Mechanical properties of finite element model of the femur when considered as 156 constant. 157 2.5 Analysis type 158 A linear modal testing analysis was used on the FE model to assess the vibration 159 characteristics of the human femur. The modal analysis package of ANSYS software 160 was used. To implement the numerical procedure, first the mode-extraction method to 161 be used for the modal analysis is chosen. In this study, the Block Lanczos method is 162 selected since it is a fast, efficient and robust algorithm to perform modal analysis. 163 Next the frequency range has to be defined. To include the first five modes of 164 vibrations, a range from 0 to 1000 Hz is taken. 165 A harmonic response analysis was also performed to obtain the vibration amplitude of 166 the femur at a specific frequency range under a vertical load applied at the femoral 167 head for the three different cases analyzed (condyle, diaphysis and physiological 168 boundary conditions). This linear analysis has also been performed using the 169 commercial software ANSYS. The harmonic response analysis method chosen is the 170 Full method within a frequency range from 0 to 2000 Hz for all the boundary 171 conditions analyzed. 172 The numerical analyses were run in a laptop PC Intel 1.8 GHz (1 core) with 8GB RAM. 173 CPU time of the analyses was estimated in 8 minutes. 174 3 Results 175 Table 2 summarizes the dynamic characteristics of the femur analyzed. The first five 176 natural frequencies as well as their corresponding modes of vibration are detailed for 177 the different boundary conditions (BC) analyzed (free BC; diaphysis BC; condyle BC; 178 physiological BC) and the two different material properties set-ups performed 179 (constant mechanical properties and Hounsfield Units). The vibration plane of the 180 7 modes of vibration is specified: frontal bending, sagittal bending, combined bending 181 (both frontal and sagittal bending) and torsion (T). 182 Natural frequencies range from 245 to 814 Hz for the free boundary conditions, 108 to 183 887 Hz for the diaphysis boundary conditions, 222 to 803 Hz for the condyle boundary 184 conditions and 107 to 782 Hz for the physiological boundary conditions. In order to 185 compare qualitatively the four analyzed cases, the normalized natural frequencies for 186 the different set-ups analyzed have been represented (fig. 2). Normalization has been 187 performed with respect to the lowest frequency. It can be observed the moderate 188 effect of material properties compared to the much substantial effect of boundary 189 conditions. In particular, a high difference can be found from case to case of boundary 190 conditions, being a factor higher than 2 in some frequencies. The biggest differences 191 are found between free and condyle boundary conditions versus physiological ones. 192 On the other hand, differences due to the fact of considering constant mechanical 193 properties along the bone tissue are not so acute although important: a difference in 194 the range 5-25% can be found for different boundary conditions across natural 195 frequencies. Higher variations are found for the case of the condyle boundary 196 conditions, whereas considering constant mechanical properties for the remaining 197 cases is an assumable hypothesis (fig. 2). 198 To illustrate the modes of vibration for one of the set-ups analyzed, fig. 3 shows the 199 deformed shape of the bone for the first five natural frequencies in the case of the 200 physiological boundary conditions and constant mechanical properties assignment. It 201 may be observed that sagittal bending modes are given at the first and fifth natural 202 frequencies (fig. 3 and Table 2), frontal bending modes at the second and fourth 203 natural frequencies and torsional mode at the third natural frequency. In the case of 204 free boundary conditions, frontal bending modes are shown at the fourth natural 205 frequency, sagittal bending modes at the fifth frequency, torsional mode at the third 206 natural frequency and combined bending modes for the first and second natural 207 frequencies (Table 2). For the diaphysis BC, frontal bending is shown for modes 1 and 208 4, sagittal bending for modes 2 and 5 and torsional bending for mode 3. On the other 209 hand, in the case of the condyle boundary conditions, Table 2 shows front sagittal 210 modes at natural frequencies 1 and 4, frontal bending modes at 2 and 5 and torsional 211 mode at the third natural frequency. 212 The results of the harmonic analysis are shown in fig. 4. It shows the dimensionless 213 dynamic amplification factor, i.e. the vertical displacement vector, versus its 214 corresponding value in a static analysis with the same loading conditions, in a spectral 215 plot, for the different analyzed boundary conditions in the case in which material 216 properties are assigned based on the HU method. As can be observed, the natural 217 frequencies differed in the three cases as well as the amplitude of vibration pointing 218 out the importance (quantitatively and qualitatively) of the boundary conditions 219 chosen on the dynamical response obtained (fig. 4). 220 8 221 222 Free BC Diaphysis BC Condyle BC Physiological BC Const HU Const HU Const HU Const HU 1 ω1 (Hz) 271 245 111 108 276 222 146 107 φ1 CB CB FB FB SB SB SB SB 2 ω2 (Hz) 310 274 122 119 398 283 269 165 φ2 CB CB SB SB FB FB FB FB 3 ω3 (Hz) 660 567 592 531 468 378 330 240 φ3 T T T T T T T T 4 ω4 (Hz) 790 703 826 738 775 643 702 483 φ4 FB FB FB FB SB SB FB FB 5 ω5 (Hz) 814 760 887 861 803 716 782 619 φ5 SB SB SB SB FB FB SB SB Table 2. Summary of the dynamic characteristics of the femur: natural frequencies ω (Hz) and 223 associated vibration modes φ. Results are shown for the different analyzed boundary conditions 224 (BC): free BC; diaphysis BC; condyle BC; physiological BC and for the two different material 225 properties set-ups performed: constant mechanical properties (Const) and Hounsfield units (HU). 226 For the modes of vibration, the corresponding vibration plane is specified: frontal bending (FB), 227 sagittal bending (SB), combined bending (CB) and torsion (T). 228 229 9 230 Figure 2.First 5 normalized natural frequencies (ω) of the femur under different analyzed 231 boundary conditions (BC): free, diaphysis, condyle and physiological for the two mechanical 232 properties assignment set-ups: constant mechanical properties (C) and Hounsfield Unit (HU). 233 Normalization has been performed with respect to the lowest frequency for each set of 234 analyzed boundary conditions and each natural frequency. 235 236 0 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 C HU C HU C HU C HU C HU Free BC Diaphysis BC Condyle BC Physiological BC 16 [39] Kiiski J, Heinonen A, Jarvinen TL, Kannus P, Sievanen H. Transmission of vertical 423 whole body vibration to the human body. J Bone Miner Res 2008; 23:1318–1325. 424 [40] Matsumoto Y, Griffin M. Dynamic response of the standing human body exposed 425 to vertical vibration: influence of posture and vibration magnitude. J Sound Vibr 1998; 426 212:85–107. 427 [41] Crewther B, Cronin J, Keogh J. Gravitational forces and whole body vibration: 428 implications for prescription of vibratory stimulation. Phys Ther Sport 2004; 5:37–43. 429 [42] Pasqualini M, Lavet C, Elbadaoui M, Vanden-Bossche A, Laroche N, Gnyubkin V, 430 Vico L. Skeletal site-specific effects of whole body vibration in mature rats: from 431 deleterious to beneficial frequency dependent effects. Bone 2013; 55:69–77. 432 [43] Dionello CF, Sá-Caputo D, Pereira HV, Sousa-Gonçalves CR, Maiworm AI, Morel DS, 433 Moreira-Marconi E, Paineiras-Domingos LL, Bemben D, Bernardo-Filho M. Effects of 434 whole body vibration exercises on bone mineral density of women with 435 postmenopausal osteoporosis without medications: Novel findings and literature 436 review. J Musculoskelet Neuronal Interac 2016; 16:193-203. 437 [44] Nyman JS, Granke M, Singleton RC, Pharr GM. Tissue-Level Mechanical Properties 438 of bone contributing to fracture risk. Curr Osteoporos Rep 2016; 14:138-150. 439 [45] Matcuk GR, Mahanty SR, Skalski MR, Patel DB, White EA, Gottsegen CJ. Stress 440 fractures: Pathophysiology, clinical presentation, imaging features, and treatment 441 options. Emerg Radiol 2016; 23:365-375. 442 443