Numerical estimation of bone density and elastic constants distribution in a human mandible

Numerical estimation of bone density and elastic constants distribution in a human mandible J. M. Reina a,∗J. M. Garc´ ıa-Aznar bJ. Dom´ ınguez aM. Doblar´ eb 24th February 2006 aDepartment of Mechanical Engineering, University of Seville, Escuela Superior de Ingenieros, Camino de los Descubrimientos s/n E-41092 Sevilla, Spain bGroup of Structural Mechanics and Materials Modelling, Aragon Institute of Engineering Research (I3A), University of Zaragoza, Mar´ ıa de Luna, 3, E-50018 Zaragoza, Spain ∗Corresponding author. Tel. +34-954487311 Fax. +34-954487295 Word count: 3200 1 1 INTRODUCTION 1 Introduction Many analyses may be found in the literature that use the FEM (Siegele and Solt´ esz, 1989; del Valle et al., 1997; Meijer et al. 1992; Meijer et al., 1994) to determine the stress level in dental implants and the surrounding bone. In those models, it is necessary to establish the mechanical properties of the materials involved (i.e., bone and the material of the implant). The most widely used implants are made with metallic materials, with well-known elastic properties. This is not the case of bone material. Its complex behaviour has been a subject of intense research for long (Beaupr´ e et al., 1990; Cowin and Hegedus, 1976; Doblar´ e and Garc´ ıa, 2002; Huiskes et al., 1987; Hazelwood et al., 2001; Jacobs, 1994). The difficulties arise from its heterogeneity and anisotropy, apart from the important fact that, as a living tissue, its microstructure and mechanical properties evolve with time. The problem of heterogeneity is traditionally solved using macroscopic models with averaged mechanical properties (Siegele and Solt´ esz, 1989; del Valle et al., 1997; Meijer et al. 1992; Meijer et al., 1994). Some of them also distinguish between different areas where mechanical properties are different, including the anisotropic behaviour (e.g. Korioth et al., 1992). The evolution of the microstructure and mechanical properties with time is related to bone remodelling. This phenomenon was studied during the second half of the 19th Century, by Wolff (1986), but it was not formulated mathematically until 1976, by Cowin and Hegedus (1976). Many bone remodelling models have been formulated since, taking as starting point the ideas established by these authors. These models have been traditionally used to predict density distributions in various bones, but mainly in the femur. Many models are able to predict the bone density, but only a few can predict the anisotropy distribution with reasonable accuracy. One of these latter was developed by Doblar´ e and Garc´ ıa (2002) and applied to the proximal femur. Starting from an arbitrary initial situation (uniform density ρ =0.5g/cm3and isotropic behaviour), and applying the normal walking loads, they predicted the bone density and its elastic constants with an acceptable approximation. The object of the present study is to extend the above analysis to obtain the distribution of these same parameters in the case of the human mandible applying normal mastication loads. The density and anisotropy distributions obtained have been validated with data found in the literature (Arendts and Sigolotto, 1989; Schwartz-Dabney et al., 1991). 2 2 MATERIALS AND METHODS 2 Materials and Methods 2.1 FE model: geometry and materials The position of a set of points in the surface of the human mandible was obtained using a coordinate measuring machine. For the sake of simplicity, measures were limited to the right half of the jaw. Later, an operation of symmetry with respect to the median plane in the symphyseal region was applied. Once the outer surface of the mandible was obtained, the internal volume was meshed with linear 8-noded hexahedral elements (type C3D8 of the elements library of ABAQUS R ). Measures was also limited to the basal bone. The teeth geometry was approximated based on the few teeth that were still present and the alveolar process, altered by the individual’s edentulism, was extrapolated from the basal bone. A layer of elements with a thickness of 0.2mm was used to simulate the periodontal ligament, similarly to Korioth et al. (1992). The FE model had a total of 77,490 elements and 88,836 nodes. It is shown in figure 1. The material properties of bone are defined in section 2.3, while non-remodelling materials, the periodontal ligament and teeth, were attributed elastic linear isotropic behaviour (see table 1). Teeth are essentially composed of dentin, surrounded by a layer of enamel. This layer covers a part of the teeth, above the gingiva and has not been considered here. 2.2 FE model: boundary and loading conditions In the FE model, the forces exerted by the masticatory muscles were imposed as external loads, distributed in the insertion area of each muscle. Figure 1 highlights in different colours the various groups of nodes where the different muscles were inserted (Hylander, 1992). The orientation of these forces were taken from a similar model made by Korioth et al., 1992, and the magnitude, from the same source used by these authors, Nelson 1986, adding some other load cases not modelled by Korioth et al. Boundary conditions were imposed on the nodes of the joint surface of the condyles and on the nodes of the teeth corresponding to each type of bite. During canine and incisive clenching, the mouth is closed, or practically closed, depending on the size of the food being cut. When the mouth is closed, the action of the mastication muscles confronts the anterior surface of the condyle with the posterior surface of the articular eminence in the temporal bone. The articular surface of both condyles was fixed in the canine and incisive clenching (see figures 1 and 2). The mastication forces are the result of the pressure in the teeth-food contact. In the present 3 2.2 FE model: boundary and loading conditions 2 MATERIALS AND METHODS model, displacements were simply restrained at the nodes of the surface of the teeth that come in contact with the food. This way, the reactions in those nodes represent mastication forces. Koolstra et al. (1988) and Haraldson et al. (1988) show that these mastication forces have a vertical component and a small component transverse to the axis of the jaw. Displacements in these two directions were restrained in canine and incisive clenching, in the canine cuspids and in the incisal borders of the incisors (figure 2). Mastication with molars were modelled differently. Mastication produces cyclic movements of opening and closure of the mouth with a small lateral deviation (Hylander, 1992), called chewing cycles. The instant of maximal bite force practically coincides with the centric occlusion (Graf, 1975; Hylander, 1992): the mouth is closed and the condyles at their back position in contact with the articular eminence of the temporal. It will be assumed that the food thickness prevents the ipsilateral condyle from contact. Consequently, when mastication is carried out with the right side, for example, the articular surface of the left condyle was fixed (figure 2) and the right condyle were assumed free to move. Mastication forces in the transverse and axial directions are a consequence of the resistance that food offers to be crushed, but very small. The vertical force is the component of highest magnitude, being responsible for the grinding of the food. Vertical displacements were restrained in the occlusal face of the corresponding molars in order to simulate these forces as reactions. According to Carter et al. (1987), bone remodelling depends on the maximal stresses that the bone bears throughout its load history. Assuming that mastication is a pseudostatic process, these maximum values can be obtained by solving a static problem in which the forces developed by the masticatory muscles at the moment of centric occlusion are applied, together with the commented displacement boundary conditions, at the teeth and condyles. The load history of the mandible was simplified by assuming a mastication pattern referred to as “alternating bilateral”: a succession of mastication with the right molars (RM load step) and with the left molars (LM load step). Manns and D´ ıaz (1988) established that 75% of the population follows this pattern, as opposed to the 10% who presents a simultaneous bilateral pattern (food is located among the molars on both sides), and the other 15% with either left or right unilateral mastication. The food is first cut by the incisors and then a symmetrical incisive bite is applied (I load step), involving the four incisors. Following comes a canine bite with the right side (RC load step), where food is cut with the second incisor, the canine, and the first premolar of that side. After this, a left canine bite is applied (LC load step), symmetrical to the previous one. Finally, unilateral mastication 4 2 MATERIALS AND METHODS 2.3 Internal bone remodelling model is alternated with the first and second molar on either side, up to a total of 15 load cases. The complete sequence is: I-RC-LC-RM1-LM1-RM2-LM2-RM1-LM1-RM2-LM2-RM1-LM1-RM2-LM2, where RM1, for example, is a mastication with the first right molar. No distinction was made between the mastication forces with the first and second molars due to the lack of data. It must be pointed out that in the long-term, the order of application of the load cases has not much influence in the results [17] and only the number of cycles of each load affects those results. Therefore, these sequences may simulate any random sequence, within some limits, if they have the same loads in the same number and in the same proportion. The remodelling response of the bone is not significantly affected by the order in which loads are applied (Beaupr´ e et al., 1990; Jacobs, 1994). It is, however, influenced by the number of daily cycles corresponding to each load case. It has been supposed that the daily number of cycles is n=500, distributed among the different activities in the same proportion as seen in the previous sequence. Yet, instead of superimposing all the activities in a day, it has been assumed that only one activity is developed each day, with the previous sequence being a sequence of days. Jacobs (1994) proved that, on a long term level, grouping the load cases this way, does not affect the results significantly, provided that the grouping time (one day in this case) is short enough (Jacobs, 1994; Doblar´ e and Garc´ ıa, 2002). 2.3 Internal bone remodelling model The remodelling model based on the theory of internal variables, developed by Doblar´ e and Garc´ ıa (Garc´ ıa, 1999; Doblar´ e and Garc´ ıa, 2002) has been used in this work. The mechanical properties of bone depends on the porosity and the fabric tensor, ˆ H, (Cowin and Hegedus, 1976). Doblar´ e and Garc´ ıa defined a remodelling tensor Hthat includes both the amount of material and the anisotropy. The eigenvectors of ˆ Hare parallel to the orthotropy directions and the influence of porosity (or equivalently the apparent density ρ , or the bone volume fraction, vb) was given by Beaupr´ e et al. (1990) and Hernandez et al. (2001). Beaupr´ e et al. E=   2014 ρ 2.5si ρ ≤1.2g/cm3 1763 ρ 3.2si ρ >1.2g/cm3,(1) Hernandez et al. E=84370v2.58 b α 2.74 (2) where α represents the ash fraction and varies due to the mineralization of bone tissue, a process not considered here. Considering an average value of α =0.6 (Garcia-Aznar et al.), equation (2) 5 3 RESULTS becomes E=3388 ρ 2.58, which gives larger values for E than the Beaupr´ e correlation (1). This equation was used in a previous study (Martinez et al. 2003) and those results will be compared with the results obtained using equation (2). The mechanical stimulus, Y, is defined in this model as the thermodynamic variable associated to the remodelling tensor: Y= ∂ Ψ(H, ε ) ∂ H(3) with Ψthe free energy function. Another tensor J, was defined in Garc´ ıa-Aznar (1999) to differently weigh the deviatoric and octaedric parts of the stimulus, by means of a parameter ω ∈[0,1]. J=1− ω 3trY I+ ω devY,(4) Remodelling criteria are scalar combinations of Jand define the resorption, formation and equilibrium ranges. Two functions, grand gfwere defined for the remodelling criteria, such as figure 3: gf(J,Ψ∗ t,w)≤0gr(J,Ψ∗ t,w)>0 resorption gr(J,Ψ∗ t,w)≤0gf(J,Ψ∗ t,w)>0 formation gr(J,Ψ∗ t,w)≤0 and gf(J,Ψ∗ t,w)≤0 dead zone (5) where Ψ∗ tand ware respectively the “reference stimulus”, or “attractor state”, and the “dead zone width” (Huiskes et al., 1987). Beaupr´ e et al. (1990) provided remodelling curves for femur and cranium, the latter had a lower reference stimulus and a very small slope in resorption. This agrees with the observations of Turner (1999): the reference stimulus experiences a long-term adaptation to the external stimulus. So, weight-bearing long bones, like the femur, has higher values of the reference stimulus than protective flat bones, like the cranium. The evolution of the remodelling tensor, H, can be found in Garc´ ıa-Aznar, 1999, and in Doblar´ e and Garc´ ıa, 2002, ˙ H=f(J,Ψ∗ t,w)(6) which is integrated using a forward Euler explicit integration scheme. 3 Results All the simulations start from an unrealistic situation: isotropic and homogeneous density distribution ρ =0.5g/cm3. After applying mastication loads during a certain time, the anisotropy and density distribution changes reaching a remodelling equilibrium situation. 6 3 RESULTS The mastication habits of the individual influences that equilibrium situation. In this sense, an important factor to be analyzed is the relation between the number of bites (incisive and canine) and the number of mastications. The sequence described above, from now on called S1, was used in a previous study (Mart´ ınez et al. 2003), whose results are compared with another sequence, S2, which only includes mastications with the molars: RM1-LM1-RM2-LM2... The use of the two mentioned correlations between Eand ρ (equations (1) and (2)) are also compared. Lastly, a sensitivity analysis was made, dealing with the influence of parameters ω and Ψ∗ t. The first one, ω , measures the significance of the anisotropy of the stimulus in the remodelling response. Doblar´ e and Garc´ ıa (2002) used ω =0.1 in their analysis of the proximal femur. The values compared here are slightly higher ω =0.18,0.24,0.3 and 0.33. As commented above, the reference stimulus, Ψ∗ t, is different from one bone to another, being the high loaded bones those with greater reference stimulus. The mandible does not bear loads as high as the femur and consequently its reference stimulus should be lower. Doblar´ e and Garc´ ıa used a value of Ψ∗ t=50MPa/day for the femur, together with a dead zone width of 2w=25MPa/day. Here a reference value of Ψ∗ t=10MPa/day has been adopted as initial reference, keeping the dead zone width as half of the reference stimulus. Other values analyzed were Ψ∗ t=25 and 50MPa/day. The remodelling equilibrium may be characterized by a negligible variation of bone mass. This can be interpreted as a convergence criterium. In order to check that convergence the variable conv was defined. conv =Rv˙ ρ dV Rv ρ dV (7) The evolution of this variable is shown in figure 4 for the simulation that includes the S2 sequence, uses the correlation of Hernandez (2), ω =0.3 and Ψ∗ t=10MPa/day. Most of the results that follow correspond to this case, which will be called the reference simulation from now on. From figure 4, it can be stated that the global remodelling equilibrium, i.e. convergence, has been reached after 368 load steps. Figure 5 shows the final bone density distribution in the reference simulation. The results of density practically coincide in all simulations. If a density limit is established at 1.92g/cm3to distinguish between cortical and trabecular bone (Beaupr´ e et al., 1990), practically all of the mandible’s surface resulted cortical bone. This coincides with what actually happens. However, two different areas of cortical bone with low density can be distinguished: the coronoid process and the pterygo-masseteric tuberosity. In this zone the cortical layer is thinner (Arendts and Sigolotto, 1989 and 1990) and it is also less rigid, therefore with lower density. Figure 6 shows some sections of the mandible with their corresponding distribution of bone 7 3 RESULTS density. Some computer tomographies were taken from the actual mandible at the same sections and are shown below. The similarity between the numerical results and the CTs is quite noticeable, except, perhaps, for the upper third of the sections. It must be remembered that the geometry of this region, the alveolar process and the teeth, could not be measured and it was only approximated. The molar CT section shows a thin layer of dense bone covering the hollow, left by the tooth. This tooth was lost previously to death and external bone remodelling changed that region. These changes have not been taken into account, since the tooth was simulated to be present. In all sections, a central area of trabecular bone surrounded by a layer of cortical bone can be distinguished. This tubular structure is usual in the diaphysis of long bones. Bone structure is optimum from a strength point of view (Currey, 1984) and nature puts the bone tissue to maximize its stiffness with the least weight. A tubular section is without a doubt the best option for resisting the bending and torsion that mastication loads produce. Global equilibrium is determined by the conv variable. It does not, however, evidence the areas where this remodelling equilibrium has not been locally reached. In order to evaluate this convergence, the evolution of density and elastic constants was analyzed in 60 control points placed throughout the mandible. After 368 days of loads in the reference simulation, all points reached remodelling equilibrium. The rest of the simulations required a similar number of loading days to reach convergence. Figure 7 shows the evolution of density and the elasticity moduli Ea,Etand Er at point P, highlighted in figure 6 (at the first right molar, in the labial side). Ea,Etand Ercorrespond to the elasticity moduli in the following directions, respectively: axial (parallel to the axis that runs through the corpus of the mandible), tangential (contained in the section perpendicular to the previous axis and tangent to the profile of that section), and radial (perpendicular to the previous two directions). Finally, table 3 shows the average values of the elasticity moduli in the cortical bone layer of the symphyseal region: at the incisors and at the first molar. These values are compared with the experimental ones obtained by Schwartz-Dabney et al. (1991) for the symphyseal and mentonian region and with those ones obtained by Arendts and Sigolotto (1990). The latter ones were averaged through the entire mandible, which almost completely distorts them. Table 3 shows that using the Beaupr´ e’s correlation (1), the elasticity moduli are significantly lower than using the correlation of Hernandez (2), and less similar to those obtained experimentally. Regarding with the anisotropy, the largest stiffness resulted in the axial direction, followed by the tangential and the radial directions, in all the simulations and for all the points of the mandibular 8 3 RESULTS corpus. This result is logical since the axial direction is the one requiring larger stiffness to resist the bending stresses produced by mastication. The two simulated mastication sequences lead to very similar results, almost identical in the density distribution, which consequently was not shown. This was expected, given that they differ in only a few load cases. Despite that slight difference, it can be affirmed that results of sequence S2 are more similar to the experimental ones (see table 3) than those of sequence S1, apart from being more reasonable. It is convenient to remember that in S1 there is one bite for every five mastications, which seems clearly excessive. This proportion between load cases depends on the eating habits of the individual and type of food. However it has a scarce influence on the results, except for the case of very particular mastication patterns. The influence of parameter ω on the results of the model is quite notable. This parameter weighs the influence of the deviatoric part of the stimulus, i.e. of the load, on the remodelling response. This deviatoric part is larger in the mandible than in the femur, because bending, the main load that the femur resist, must be added to the torsion produced by mastication loads. Thus, in the mandible, it is more recommendable to use larger values of ω than the 0.1 used by Doblar´ e and Garc´ ıa (2002) for the femur (see table 3). The larger the value of ω , the larger the degree of anisotropy in the cortical bone layer: an increase in the axial stiffness and a decrease in the radial one, is obtained, the transverse being only vaguely affected. The choice of ω =0.3 seems the most reasonable as it leads to results more similar to the experimental ones, especially in sequence S2. The influence of parameter Ψ∗ twas also analyzed. The sequence S2, the correlation of Hernandez and ω =0.3, with three different values of Ψ∗ t: 10, 25 and 50MPa/day were chosen for the sensitivity analysis. A comparison of the density distribution at the region of the first molar is given in figure 8. Great differences can be seen from one to another, the one corresponding to Ψ∗ t=10MPa/day giving a thicker cortical layer than the others. This result is more similar to reality, as can be seen by comparing figures 6 and 8. A higher reference stimulus makes net formation only to appear in zones with very high stress level, thus resulting in a very thin layer of cortical bone. This result shows that the reference stimulus used by Doblar´ e and Garcia (2002) for the femur, Ψ∗ t=50MPa/day, is not proper for the mandible, that bears not so high loads. The elastic properties are very similar in those points with the same density, that is, the cortical bone layer with maximal density (see table 4). 9 REFERENCES REFERENCES Ψt Ψ∗ t ww ˙r Resorption zone Dead zone Formation zone Figure 3: 0 100 200 300 400 0,00 0,02 0,04 0,06 0,08 conv days Figure 4: 16 REFERENCES REFERENCES ρ (g/cm3) Figure 5: 17 REFERENCES REFERENCES abc •P Figure 6: 18 REFERENCES REFERENCES                    ρ                   days E(GPa) ρ (g/cm3) Figure 7: 19 REFERENCES REFERENCES abc Figure 8: E(MPa) ν Dentin (Craig and Peyton, 1958) 17600 0.25 Periodontal ligament (Widera et al., 1976; Ralph, 1982) 3 0.45 Table 1: Muscle Orientation of the forces Magnitudes of the forces (N) X Y Z Incisive Canine Molar R L R L R L Superficial masseter +0.419 +0.207 -0.885 76.2 76.2 87.6 110.4 106.6 38.1 Deep masseter -0.358 +0.546 -0.758 21.2 21.2 37.5 47.3 45.7 16.3 Anterior temporalis +0.044 +0.149 -0.988 12.6 12.6 85.3 22.1 102.7 80.6 Middle temporalis -0.500 +0.221 -0.837 5.7 5.7 45.9 19.1 57.4 50.7 Posterior temporalis -0.855 +0.208 -0.474 3.0 3.0 31.8 19.7 40.8 40.8 Medial pterygoid +0.372 -0.486 -0.791 136.3 136.3 96.1 82.2 169.6 82.2 Lateral pterygoid +0.757 -0.630 +0.174 61.9 61.9 28.7 62.1 33.5 23.9 Table 2: 20 REFERENCES REFERENCES S1 S2 1st molar incisive 1st molar incisive Hernandez ω =0.18 Ea17.4 17.8 17.4 17.9 Et16.6 16.0 16.7 16.0 Ψ∗ t=10MPa/day Er14.7 14.4 14.6 14.3 Hernandez ω =0.24 Ea18.3 19.6 19.0 20.0 Et17.0 16.6 17.5 16.1 Ψ∗ t=10MPa/day Er13.1 12.9 12.7 12.8 Hernandez ω =0.30 Ea20.8 24.8 22.1 23.8 Et19.5 16.5 19.0 16.6 Ψ∗ t=10MPa/day Er10.4 10.3 10.1 10.6 Hernandez ω =0.33 Ea23.7 27.3 22.5 27.6 Et19.7 16.2 21.7 17.7 Ψ∗ t=10MPa/day Er8.9 9.4 8.6 9.0 Beaupr´ e ω =0.30 1Ea15.5 21.5 - - Et11.0 14.9 - - Ψ∗ t=10MPa/day Er10.4 9.3 - - Schwartz-Dabney et al. (1991) Ea23 - 23 Et15 - 15 Er10 - 10 Arendts and Sigolotto (1990) Ea17.3 Et8.2 Er6.9 Table 3: First molar Incisors Ψ∗ t=50MPa day Ψ∗ t=25MPa day Ψ∗ t=10MPa day Ψ∗ t=50MPa day Ψ∗ t=25MPa day Ψ∗ t=10MPa day Ea23.9 23.9 23.8 22.1 22.1 22.1 Et16.7 16.7 16.6 17.0 16.9 19.0 Er10.2 10.2 10.6 11.5 11.6 10.1 Table 4: 21 View publication statsView publication stats