Department of Biomedical Engineering, University of Florida
P.O. Box , Gainesville, FL 32611, USA firstname.lastname@example.org Abstract
This study explored the asymmetrical distribution of cortical bone in the Macaca fascicularis mandible. The morphology of the mandible is very intriguing and it has attracted much attention due to its complexity. The objectives of the present study were to use Finite element analysis (FEA) to test the hypothesis that the cortical asymmetry is related to strain energy density (SED) and to make predictions about the remodeling activity in the mandible.
FEA was used to create a realistic mandibular model, to mimic in vivo physiologic loading and to analyze the model. The most common loading environment was simulated: combined loads (bending, torsion and direct shear) with a tilted biting load. Linear regression analysis was used to correlate the thickness and SED values. SED criterion was applied to predict the remodeling activity.
The regression results suggest that SED and thickness are related and a strong positive correlation exist between them. The SED values found at midcorpus always exceed the values found at the mandibular base. The SED values found at the lateral midcorpus were higher than the values at the medial midcorpus. The values were not distributed uniformly throughout the mandibular bone or within the lazy zone interval.
The study supports the functional relationship theory between form and function of the mandible and based on the SED profile, it concluded that the mechanostat theory might not be valid for the mandible. The current work is intended as a contribution to the continued research on the biology of bone and is aiming to enhance our understanding of primate masticatory biomechanics.
Keywords: mandible, finite element method, strain energy density, thickness, remodeling
In the last 25 years, extensive research on macaque mastication explored mandibular anatomy, mandibular movements during mastication, investigated biting and reaction forces occurring during mastication, portrayed the stress-strain behavior of the mandibular bone and overall, expanded our understanding of primate masticatory biomechanics. The macaque model is an excellent model for studying mastication, not only because of abundant available data, but also because it is a primate model. Studies on the primate skull are regularly used as reference for studying human masticatory biomechanics.
The mandible is characterized by a very odd and fascinating geometry: the cortical bone is distributed asymmetrically throughout the entire mandible. Despite of extensive research on the mastication system, the biomechanical justification for this unique, asymmetrical distribution of cortical bone is still ambiguous. A direct relationship between mandible form, function, and mechanical load history, although crucial from a biomechanical point of view, was often assumed but has never been established. Despite an abundant record of biomechanical studies on mandibular morphology and profiles of strain (Hylander 1979a, Daegling and Green, 1991, Daegling 1993, 2002, 2004, Dechow and Hylander, 2000) nothing is known about the relationship between the local differences in mandibular bone mass and SED. Understanding the functional morphology of the mandible is critical for uncovering the evolutionary transformations in facial bones form and expanding our knowledge of primate origin.
The present study concentrates on the relationship between strain, SED and local differences in bone mass and thickness. An improved mandible model, the mandible with masticatory muscles, is used to simulate the physiologic loading conditions which occur during mastication. Masticatory muscles (left temporalis muscle, left masseter-pterygoid sling, right temporalis muscle and right masseter-pterygoid sling) are simulated as individual vectors. A FE analysis is performed in which the mandible is subjected to combined loading: torsion, direct shear and parasagittal bending. Principal strain values and SED data are recorded. SED values are used to evaluate the remodeling process and to predict if the remodeling activity will be present or absent.
Regional variation in cortical bone
The asymmetrical pattern of cortical bone distribution in the mandible is unique and fascinating. Even more intriguing is that this cortical asymmetry is stereotypical among anthropoid primates regardless of variations in mandible dimensions or dietary preferences (Daegling 2002, Daegling and Hotzman 2003). Considerable differences in cortical bone can be observed between the basal or alveolar regions, symphysis or molar region, and medial or lateral aspects of the mandible. The mandibular thickness varies significantly throughout the mandible (Daegling 1993, Futterling et al, 1998). The mandible is vertically deep and transversely thick. In the molar region, the lingual aspect of the mandibular corpus is thinner than the lateral aspect. The distribution of cortical bone changes from the molars toward the symphysis, such that under the premolars the thin lingual bone is much less apparent. The base of the mandibular corpus in the molar region is the thickest part. At midcorpus, the mandibular corpus is thicker on the lateral aspect than on the medial aspect (Daegling 1993). In addition, experimental studies showed that not only the geometrical properties but also the mechanical properties differ significantly throughout the mandible. The mandible is very stiff in the longitudinal direction and usually stiffer on the medial aspect than on the lateral aspect (Dechow and Hylander, 2000).
Both functional and non-functional explanations have been presented over the years but currently there is no consensus regarding the unusual cortical distribution in the mandible. One of the non-functional theories which tried to explain the cortical asymmetry is that the mandibular corpus is deep and thick to accommodate large teeth, more specifically their long roots (Hylander 1988). However, this theory was not accepted as the roots do not extend all the way down to the mandibular base. Some studies show that there is actually no relationship between the mandibular corpus dimensions and teeth size (Daegling and Grine, 1991).
Numerous studies tried to offer a functional justification and to relate the mandibular corpus size to the functional demands of the mandible during mastication. These studies made assumptions of how the cortical bone should be distributed in the mandible for the mandible to resist the masticatory loads and stresses. Hylander proposed several possible explanations as to why the cortical bone is distributed the way it is in the mandibular corpus: the corpus should be vertically deep to counter parasagittal bending stress and transversely thick to counter torsional stress and wishboning (Hylander, 1979 a,b, 1988).
The mandible has a very intricate, asymmetrical shape. It is materially anisotropic, with a huge variation in material properties throughout the entire structure. These studies emphasized the difficulties in analyzing the mandible due to considerable regional variation in thickness, cortical area, size, shape and mechanical properties throughout the bone.
Loading patterns and strain gradients
Numerous studies explored a functional relationship between the form and the function of the mandible (Hylander 1979a, b, 1984, Demes 1984, Russell 1985, Hylander et al., 1987, 1998, Lahr and Wright, 1996, Ross and Hylander, 1996). A large body of research explored the relationship between the stress and strain patterns and the mandible morphology (Hylander 1979a, Daegling and Green, 1991, Daegling 1993, 2002, 2004, Dechow and Hylander, 2000). Although extensive research exists, a functional correlation between the mandibular morphology and the stress and strain patterns has never been established and it is still one of the most controversial issues in physical anthropology.
Different regions of the mandibular corpus are loaded differently during mastication. In vivo experiments brought evidence that the macaque mandible is subjected to a combination of bending and torsion during mastication (Hylander, 1979b, 1981, 1984,; Hylander and Crompton, 1986, Hylander et al. 1987, Hylander and Johnson, 1997). The distribution of stresses and strains during mastication were inferred from theoretical and experimental studies (Knoell, 1977, Hylander, 1984, Bouvier and Hylander 1996, Daegling and Hylander, 1997, 1998, Dechow and Hylander, 2000). Specifically, during the mastication, the mandible is primarily twisted about its long axis and sheared perpendicularly to its long axis. In addition, the mandible is subjected to parasagittal and transverse bending (Hylander 1979b). The simultaneous application of twisting, bending and direct shear during mastication is a possible explanation for the unusual asymmetrical distribution of cortical bone in the mandibular corpus (Demes et al., 1984, Daegling 1993). According to theoretical studies, the shear stresses resulted from torsion and direct shear add up on the lateral aspect and are subtracted on the medial aspect of the mandibular corpus.
The strain history for the facial bones of Macaca Fascicularis is well-documented nowadays. Strain magnitudes are available for various skull regions: mandibular symphysis (Hylander 1984), zygomatic arch (Bouvier and Hylander 1996, Hylander and Johnson 1997), supraorbital bar (Hylander et al., 1991, Bouvier and Hylander 1996) and mandibular corpus (Hylander 1979 a, b, Hylander 1986, Hylander and Crompton, 1986, Hylander et al., 1998, Bouvier and Hylander 1996, Dechow and Hylander, 2000, Daegling and Hotzman, 2003). Usually the strains recorded are in the 250-1000 range which is considered the functional interval (Fritton et al., 2000, Wood and Lieberman, 2001, Daegling 2004). In conclusion, non-uniform and steep strain gradients were found for the Macaca facial bones: high strain have been found in the anterior zygomatic arch and in the mandibular corpus while low strains have been found in the posterior portion of the zygomatic arch and supraorbital bar. Particularly, in the mandible’s case, experiments show that very low strains (below 200) as well as very high strains (2000) are present. It seems that for the mandible, by evaluating the experimental strain gradients, the form does not always follow the function, or at least not all parts of the mandible are designed as to maximize strength and minimize bone tissue (Daegling and Hylander 1997, Daegling 1993).
Strain density energy and functional adaptation
Since more than 100 years ago, it was proposed and accepted that the bone is an optimized structure (Roux, 1881, Koch, 1917, Lanyon, 1973, Rubin, 1984). An optimized structure exhibits maximum strength with minimum amount of material. Many researchers proposed that especially the appendicular skeleton (bones of the limbs) acts an optimal force-resisting structure (Lanyon 1973, Pauwels 1980, Rubin 1984, Frost 1986, Rubin et al. 1994). Consequently, other bones, such as the bones of the facial skeleton, could be considered optimized structures.
Frost proposed first the mechanostat theory according to which bones adapt to mechanical loads in order to sustain those loads without hurting or breaking (Frost 1998, Schoenau and Frost 2002). The mechanostat is a combination of non-mechanical factors (hormones, calcium, vitamins, etc), mechanical factors (loads, strains, etc), modeling and remodeling mechanisms, thresholds and possibly other mechanisms. Frost proposed that the mechanostat model is applicable “in all amphibians, birds, mammals, and reptiles of any size, age and sex” (Frost 1998). Four mechanical usage windows or strain ranges are usually defined: below 50me (disuse characterized by bone loss), between 50-1500me (the adapted window or lazy zone, normal load), 1500-3000me (mild overload characterized by bone gain) and above 4000me (irreversible bone damage) (Frost 1994, Mellal et al, 2004). Most of the values are expected to be generally situated in the adapted window range or in the lazy zone interval and therefore bone homeostasis is predicted. Homeostasis means that bone resorption and bone formation are in equilibrium and therefore the values should be near uniform throughout the bone.
Following the idea that bones should be optimized structures, it has been proposed that the facial bones could be optimized as well for countering and dissipating mastication forces and consequently, there is a functional correlation between the morphology and function of the mandible (Hylander 1979a). In 1984, Demes proposed a theory which supports the hypothesis that form follows function in the mandible (Demes, 1984). He proposed a very interesting theoretical explication as to why the mandible is vertically deep and transversely thick in the molar region by using superposing torsional and occlusal loads. Demes used shear and bending moment diagrams to prove his theory. The mandible is vertically deep to counter the bending stress and transversely thick to counter the combined effects of torsion and direct shear. Moreover, shearing and torsional stresses add up on the lateral side and are subtracted on the lingual aspect of the mandible which correlates with mandibular corpus being thicker on the lateral aspect and thinner on the lingual aspect. Daegling and Hotzman performed several in vitro experimental strain analyses on human mandibles by superposing torsional and occlusal loads to test Demes’ theory (Daegling and Hotzman, 2003). The study partially supported Demes’ theory and showed that the lingual strains are indeed diminished and the lateral basal corpus strains are increased when the mandible is subjected to combined loading. However, the authors obtained different results for the midcorpus and alveolar aspects of the mandible. Various other researchers supported the hypothesis according to which the facial bones are especially optimized for countering and dissipating mastication forces. In 1985, Russell proposed a novel theory for that time regarding the morphology of the facial bones (Russell, 1985). The author postulated that the stress obtained from chewing hard food leads to developing more pronounced supraorbital region.
On the other hand, many researchers challenge the functional correlation theory and usually use experimental stress and strain data to prove their point. If bone is an optimized load bearing structure, there should be near uniform stress and strain levels throughout the bone. A large body of experimental work proves that the facial bones and mandibular bone in particular, exhibit a totally different behavior than the behavior expected for an optimized structure (Hylander 1979b, 1984, Daegling 1993, Daegling and Hotzman 2003, Hylander and Johnson, 1997, Futterling et al, 1998, Dechow and Hylander, 2000). It is recognized that facial bones cannot behave as perfectly optimized structures. Uniform stress and strain levels can be obtained for structures with uniform cross-section and material properties. The mandible presents a huge variation in geometrical and mechanical properties and therefore a variation in stress and strain values throughout the mandible is expected. Moreover, to obtain uniform stress and strain levels, a structure has to be loaded axially in tension or compression or in pure shear. The mandible is subjected during mastication to torsion, bending, direct shear or any combinations of these. Consequently, it is expected that stress and strains gradients to vary in a certain range throughout the mandible during mastication (Ross and Metzger, 2004).
Torsional and occlusal loads which occur during mastication will produce a certain stress and strain environment in the mandible. Changes in stress and strain patterns will impact cell remodeling activity. Activation of cell remodeling will result in bone resorption or bone formation. Cell remodeling will therefore have an impact on the bone mass. According with this proposed relationship, variation in bone mass and cell remodeling should be directly related to the stress and strain environment and consequently to the applied load.
The functional adaptation of the mandible is triggered by mechanical or non mechanical stimuli. Today it is accepted that mechanical stimuli govern bone adaptation (Cowin, 2001). The most common mechanical stimuli are: strain, stress, strain energy, SED, strain rate and fatigue microdamage. Strain energy density has been considered by many researchers a valid stimulus for bone remodeling (Huiskes et al. 1987, Katona et al. 1995, Cowin, 2001). SED is the rate of variation in bone density (Mellal et al. 2004). Brown and his colleagues investigated twenty-four mechanical parameters that are related to functional adaptation in bone (Brown et al, 1990). The results of the study reveal that only four parameters are directly related to remodeling: strain energy density, shear stress and tensile principal stress and strain.
In 1986, Fyhrie and Carter proposed a novel mathematical theory which described the relationships between cancellous bone apparent density and stress (Fyhrie and Carter, 1986). They assumed that bone is a “self-optimizing” material. Huiskes and his colleagues were among the first to consider strain energy density the stimulus for remodeling instead of strain (Huiskes et al. 1987). They developed an adaptive-remodeling model and used strain energy density to predict the shape or bone density adaptations. Fyhrie and Carter developed later another theory using strain energy density as remodeling stimulus. Their study showed that strain energy density can successfully predict remodeling activity in the femur (Fyhrie and Carter, 1990).
Since then, the strain energy density is successfully used to investigate remodeling patterns in bones (Katona et al. 1995, Turner and Pidaparti, 1997, Barbier et al. 1998, Cowin, 2001, Mellal et al. 2004). A strain energy density criterion was developed in which SED is the main stimulus for bone remodeling. The rate of change of apparent density at a particular location in the mandible is described by the following formula:
where is the apparent density, tis the time, B and k constants that quantify bone gain or loss, and u is the strain energy density. The area in which no net change of bone density occurs, the zone between bone densification and bone resorption, is called the lazy zone (bone homeostasis). A lazy zone can be expressed by using the following formulas, where s is expressed in percents and represents the extent of the lazy zone:
Theoretical and experimental studies on the mandible seem to convey conflicting conclusions regarding a possible relationship between the strain field, remodeling and the distribution of bone mass. One of the principal reasons for this ongoing discrepancy arises from the complexity of the mandible. To create a biomechanical model of the mandible is a very challenging task. Theoretical models (solid, hollow, asymmetrical ellipse models) can offer insight into the mandible biomechanics but because they do not account for the intricate materials properties of the mandible, they cannot be used for complex biomechanical analyses. The discrepancy in results might arise also from using different methods to test the mandible. In vivo or in vitro experimental methods are considered limited field methods. The mandible is usually analyzed only in certain regions “of interest”. Other limitations contribute to the problem: the load magnitude and the loading environment cannot be controlled in an in vivo experiment. In vitro methods can control the load magnitude but have difficulties is recreating the in vivo environment. The boundary conditions applied in vitro (loading environments, how the mandible was supported, how it was constrained, etc) are not quite similar with the in vivo boundary conditions. The results will be greatly affected by all these changes. Moreover, experimental studies can only estimate from a couple sites “of interests” the strain gradients for the entire mandible.
The best solution to combine many of the advantages of in vivo and in vitro methods is the finite element method. The finite element model of the mandible, even if it is a simplified representation, can be successfully used to estimate the real mandible behavior. The geometry and realistic material properties can be simulated in the model. The load magnitude and the loading environment can be controlled during the analysis. The stress and strain results can be obtained throughout the entire mandible not just in some regions of interest. Most important, an FEA can predict regions in the model with maximum stress and/or maximum strain values. Strain gradients throughout the model can be easily calculated, plotted and displayed. However, the success of FEA strongly depends on several significant factors: accurate geometry of the FE model, realistic material properties and replicating physiologic boundary conditions (restraints and applied loads).
Several hypotheses may be tested by examining the asymmetrical cortical bone distribution in the mandibular corpus. (1) Does the transverse cortical thickness of the mandibular corpus at various locations have a predictable relationship to SED values? (2) If a relationship exists between cortical thickness and SED, the remodeling activity should be present or absent? (3) Is the Frost mechanostat theory applicable to all bones, including the mandible? In other words, is the mandibular bone in homeostasis?
The current work aims to use FEA to answer questions related to functional morphology that the limited field methods cannot resolve. The objectives of this study are to attempt to explain patterns of cortical asymmetry in the mandible using FEA, to test the hypothesis that mandibular thickness is related to strain and SED and to make predictions about the remodeling activity using SED stimulus. The null hypothesis was that there is no relationship between cortical thickness, strain environment and SED and at all locations near uniform SED values will be measured.
Materials and methods
Strain energy density criterion
The lazy zone interval was calculated using the parameters published by various authors, k=0.004 J/g, B= 1 (g/cm3)2/ MPa time units and s = 10% (Weinans et al. 1992, Turner and Pidaparti, 1997, Mellal et al. 2004) (Figure 1). The lazy zone interval is:
0.0036 0.0044 MPa/(gcm-3)
The density () was estimated using the relationship between the elastic modulus and density (Carter and Hayes, 1977, Huiskes et al. 1987, Weinans et al. 1992):
E = 3790 x 3
The elastic modulus of the model was chosen as E = 17 GPa therefore the density is =1.65 g/cm3. Consequently, the calculated interval for the strain energy density is:
Two mandible models were obtained initially through 3D reconstruction from 90 CT scans: a dentate and an edentulous model. Both models were imported into the MSC Patran finite element analysis package (MSC Software Corporation, Santa Ana, CA). As demonstrated in the previous work (Marinescu et al, 2004), the edentulous model performed always better than the dentate one in recreating the experimental strains. The dentate model was too stiff and consequently, the strain values obtained from it were extremely low when compared with the experimental strain values. Because of the limits of spatial resolution in conventional CT, the periodontal ligament was not visualized and consequently impossible to model. As a result, no interface existed between the teeth and alveoli. The edentulous model yields more accurate results and eliminates the challenging task of modeling the periodontal ligament.
The FE analyses were performed using the edentulous mandible model. The most complex case of heterogeneity and directional dependence studied previously was chosen for the analysis. The analyzed model was heterogeneous and transversely isotropic. Material properties were assigned using a local coordinate system: the 3-axis was defined as following the length of the mandible, with the 1-axis mediolaterally oriented and 2-axis superoinferiorly oriented. The sets of material properties were assigned for cortical and trabecular bone. Cortical bone was assigned material properties in the frontal plane and in the longitudinal direction (E1,2cortical = 13 GPa, G12 = 5 GPa, E3cortical = 17 GPa, G23 = 6.91 GPa, n12 = 0.3 and n23 = 0.229). The trabecular region was modeled as isotropic (Etrabecular = 1.5 GPa and trabecular = 0.3). Three regions, right (posterior corpus and ramus), anterior corpus, and left (posterior corpus and ramus) were defined. One local coordinate system was defined for each region. The finite element analysis package was used to assign material properties for cortical and trabecular bone to each region, according with their local coordinate system.
The validation of a FE model is a very challenging task and according with some researchers, it is virtually impossible to find a perfect match (Dalstra et al, 1995). To validate a model usually FE strains are plotted against experimental strains and their correlation is determined (Keyak et al, 1993, Dalstra et al, 1995, Lengsfeld et al, 1998, Stolk et al. 2002, Barker et al. 2005, Anderson et al, 2005). Linear regression analyses are used in FE studies to correlate the measured experimental strains with the FE predicted strains. Maximum predicted principal strain is paired with maximum experimental principal strain and minimum predicted principal strain is paired with minimum experimental principal strain. The correlation coefficient and the p-value are determined.
The comparison between experimental and finite element principal strain data was presented in Marinescu et al., 2004. The maximum and minimum principal experimental strains (lateral strain gauge: 755 and - 221; medial strain gauge: 273 and -108) and their ratio (lateral strain gauge: 3.41; medial strain gauge: 2.52) were recorded during the strain gage experiment. The maximum and minimum principal finite element strains (lateral strain gauge: 528 and - 174; medial strain gauge: 154 and -84) and their ratio (lateral strain gauge: 3.03; medial strain gauge: 1.83) recorded for the edentulous homogeneous isotropic model, fully constrained bilaterally at condyles and angles were usually lower than the experimental strain values. During the strain experiment, the mandible was supposed to be restrained by a symmetrical steel fixture. However, we believe the specimen deflected during the experiment because of the mandible asymmetrical geometry and movement occurred in the transverse direction at the constraint locations. More congruent finite element results were obtained by simulating the altered boundary conditions (relaxing the degrees of freedom in the transverse direction): maximum and minimum principal finite element strains (lateral strain gauge: 768 and -243; medial strain gauge: 233 and -292) and their ratio (lateral strain gauge: 3.16; medial strain gauge: 0.80).
A linear regression analysis was performed to correlate the FE strains with the experimental strains obtained in the current study. The results of the linear regression analysis show very significant correlation between measured and predicted strain for the edentulous homogeneous isotropic model. The regression results strongly support the validity of the mandible FE model (r = 0.99, P-value=0.0001).
A new parameter, masticatory muscles, was introduced in an attempt to improve previous modeling efforts and simulate more realistically the physiologic loading environment (Marinescu et al., 2004). The human masticatory system is well described and analyzed in the literature (Pruim et al., 1978, 1980, Koolstra and van Eijden, 1992, 1997a, 1997b, 1999, Koolstra 2002, 2003 Korioth and Hannam, 1994a, 1994b, Korioth and Johann, 1999). The masticatory system can be modeled as a constrained lever model and it is based on the assumption that the forces in the mandible are related through a triangle of support (Greaves, 2000, Spencer, 1998). The unconstrained lever model uses equilibrium equations to determine the relationship between the forces in the mandible (Spencer, 1998, 1999). The muscle resultant force (M), the bite force (B) and the joint reaction forces (JW and JB) are determined through spatial relationship (Figure 2). The model assumes that the working and the balancing side muscle forces are equal and the muscles resultant force lies in the midline.
The information available about the Macaca masticatory system is limited compared with the information available about the human masticatory system. The documentation available concerning Macaca masticatory system consists usually of physiologic data or in vivo electromyographic data (EGM). The internal architecture, fiber length and cross sectional area of macaque masseter muscle and pterygoid muscle were described in detail by Anton (Anton 1999, 2000). There are many studies that offer electromyographic data for the Macaca masticatory system (Hylander and Johnson, 1994, Hylander et al., 2000).
There are multiple problems in simulating the masticatory system. The extremely complex system consists of: the mandible, which is moved in respect with the skull, two intricate temporomandibular joints and masticatory muscles. The mandible has an extremely complicated irregular geometry and it is characterized by a large variation in mechanical properties. The temporomandibular joint is the articulation between the condyle and the temporal bone. The condyle and the temporal bone are separated by a cartilaginous structure called the meniscus. The masticatory system consists of many muscles with different shapes and sizes. It is impossible to determine the contribution of each muscle from in vivo studies. To simulate the masticatory muscles, the muscles forces or the resultant masticatory force have to be estimated. Both the magnitude and the direction of the force vectors have to be known. Because the masticatory system is such a complex system, the experiments performed are very challenging and the experimental data very difficult to collect and validate. The mastication models are usually based on many assumptions and simplifications (Pruim et al. 1980, Meyer et al. 2000, Erdman et al 2002, Wagner et al. 2002, Feller et al. 2003). The most used assumptions in simulating the mastication muscles are:
The forces of the mastication muscles can be simulated as vectors.
The forces of the masseter and medial pterygoid muscles can be added and simulated through a single vector (the masseter-pterygoid sling).
The main mastication forces simulated are the forces exerted by the masseter-pterygoid sling and the temporalis muscle.
In the current study, the masseter-pterygoid sling and the temporalis muscles were simulated by applying two loads at condyles and angles, on both sides of the mandible (Figure 3). Several researchers attempted to determine the maximum possible mastication force using mathematical models or experimental models, by using electronic strain gauges (Howell and Manly, 1948). A few of the magnitudes reported for the maximum possible mastication force are: 250 N (Koolstra and van Eijden, 1992) and 300 N (Van Ruijven et al. 2002) for human mandibles and 154-258 N for monkey mandibles (Reitzik et al. 1978). The magnitude of the maximum resultant mastication force was therefore considered equal to 300N. The current model was inspired by the unconstrained lever model, in which the muscle resultant applies a force that is resisted at the biting point. The working and the balancing side muscle forces are considered equal. Therefore four loads of 75N each were applied to the mandible model to account for the left temporalis muscle, the left masseter-pterygoid sling, the right temporalis muscle and the right masseter-pterygoid sling.
The orientation for the masseter-pterygoid sling was chosen according to the data available for the Macaque masseter muscle. The fibers of the Macaque masseter muscle are oriented between 65 and 85 to the occlusal plane (Anton 1999). The sine of the masseter muscle angle for monkeys is 0.93 which gives a 68 angle from the horizontal, occlusal plane (Dechow and Carlson, 1990). Only the angles between the muscles and the occlusal plane were considered in the Dechow and Carlson study. The angle of inclination of the masseter-pterygoid sling in the frontal plane was estimated at 45. The estimation was necessary because no information was available regarding the orientation of the masseter-pterygoid sling in the frontal plane. Therefore, in the current study, the masseter-pterygoid sling vector was inclined 45 in the frontal plane and 68 from the occlusal plane. The sine of the temporalis muscle angle for monkeys is 0.87 which gives a 60 angle of inclination from the occlusal plane (Dechow and Carlson, 1990). The temporalis was considered to act only in the sagittal plane. The insertions of the masticatory muscles were simulated by points on the mandibular angles and the coronoid process of the mandible, on each side.
The model was constrained at condyles (Hart et al. 1992, Futterling et al. 1998). An equal number of nodes (three) were totally constrained at condyles on each side. The mandible was subjected to a combined loading regime (Figure 4). The mastication force is tilted by 15 in the frontal plane, toward the right side of the mandible (Daegling and Hotzman, 2003). Previous studies considered more realistic to simulate a tilted bite force. Experimental work showed that a vertical mastication force is highly questionable. It is believed that the lateral component of the bite force contributes to increase twisting of the mandible during mastication (Daegling and Hotzman, 2003). The current FEA is simulating simultaneous application of bending, torsion and direct shear, a combined loading pattern which is believed to occur most often during mastication. The tilted mastication force (100 N magnitude) was applied to the left second molar.
The experimental data available in the literature (Daegling 1993) for the cortical thickness for Macaca fascicularis was used in this study for making the correlation assessment between cortical thickness and SED values. Six mandibular sections are used in experimental studies: M3 (third molar), M2 (second molar), M1 (first molar), P4 (fourth premolar), P3 (third premolar), C/I2 (canine/incisor) (Figure 5). The variation in transverse cortical thickness in the mandible is presented in Table 1.
The same mandibular sections (CI2, P3, P4, M1, M2, M3) and the same mandibular regions (lateral midcorpus, basal and medial midcorpus) presented in the experimental work were identified in the FE mandible model. The principal strains and the SED values were automatically calculated for the entire mandible model using the finite element analysis package. The SED values were manually extracted at each section, for each region.
For each mandibular section, SED values were extracted for the lateral, basal and medial region. First, the six sections were identified in the model by using the sections described in the experimental study. Second, for each section, three regions were identified: lateral midcorpus, basal and medial midcorpus. One common node and its neighboring elements were identified for each region. The SED values for the five or six neighboring elements situated at the each location were found and their values averaged. Consequently one average SED value was obtained for each region. Eighteen averaged SED values were obtained in the end: three values (lateral, basal and medial) for each of the six sections.
Usually the experimental strains recorded on the mandible’s lateral aspect, below the left second molar are in the functional interval of 200-1000 (Hylander 1979b, Dechow and Hylander 2000, Daegling 2004). Hylander reported the highest strain values from in vivo strain gage experiments performed during mastication, strain values of up to 2000 (Hylander 1984). Finite element principal strain results obtained on the lateral aspect, below the left second molar during mastication were very similar with the published experimental strain results: maximum principal strain 1178, minimum principal strain - 1114. The values were situated in the middle of the 200 - 2000 published strain range.
According with the published experimental data, the cortical thickness varies considerably within each mandibular region and between mandibular sections. In this study, using FEA, SED values were obtained at the same locations and were correlated with available thickness data. Generally, the SED profiles and the thickness profiles vary in a similar manner within the chosen mandibular sections. The SED values were recorded for the following mandibular regions: the lateral midcorpus (Table 2), base (Table 3) and medial midcorpus (Table 4).
For the lateral region, the SED values tend to increase slightly from the anterior corpus toward the molar region (Figure 6a). The maximum values are found in the molar region. The SED values recorded for the basal region also tent to increase slightly from the anterior corpus toward the molar region, with small values in the anterior corpus and higher values toward the molar region (Figure 6b). The overall trend for the medial region is the decrease of the SED values from the anterior corpus toward the molar region (Figure 6c). The highest values are encountered in the anterior corpus. The thickness values on the medial aspect decrease correspondingly, from the anterior corpus toward the molar region. An arbitrary path was created in the finite element geometry along the midcorpus and the overall SED profile was plotted along the path (Figure 7).
The most interesting finding is that the values compared have a similar trend and the pattern is consistent across various mandibular sections. The values found at the lateral or the medial midcorpus are always higher than along the mandibular base (Figure 8 a, b). Also, the values found at the lateral midcorpus generally exceed the values found at the medial midcorpus (Figure 9). The only exception is the value for the most anterior region considered. Only in this region (C/I2), the SED value on the medial aspect is higher than the SED value recorded for the lateral aspect.
The SED values obtained from all the sections and regions were plotted on the same graph and the lazy zone interval was determined. The SED values obtained in the current study range from 0.00283 U 0.01213 (MPa). As it can be seen from Figure 10, only a few values are within the lazy zone interval. The few values situated within the lazy zone interval were collected from the medial midcorpus and basal area. None of the SED values collected from the lateral midcorpus fall within the lazy zone interval.
The mandibular thickness varies significantly throughout the mandible. The base of the mandibular corpus in the molar region is the thickest part. At the symphysis, the mandibular thickness is greatest along the medial aspect. However, in the molar region, the lingual aspect of the corpus is thinner than the lateral aspect. Especially at midcorpus, the mandibular corpus is thicker on the lateral aspect than on the medial aspect (Daegling 1993).
The thickness varies slightly on the lateral midcorpus, increasing from the anterior corpus toward the molar region. A similar trend is observed for the SED values calculated for this region, except for the most anterior section. The SED values increase slightly from the anterior corpus toward the molar region. The maximum values are found in the molar region for both variables. If the value for the most anterior section is considered an outlier and is excluded, the correlation is positive strong (r=0.673, p-value=0.12).
The base of the mandibular corpus is very thick with the highest values in the molar area. The SED values recorded for the basal region tent to increase from the anterior corpus toward the molar region. The higher values, calculated for the base of the mandibular corpus, are found also in the molar region. The correlation between the SED values and thickness is positive strong (r=0.737, p-value=0.05).
On the medial aspect, the cortical thickness decreases from the anterior section toward the molar section, having the largest value at the anterior corpus. The overall trend for the medial region is the decrease of the SED values from the anterior corpus toward the molar region. The highest SED values are also encountered at the anterior corpus. The correlation between the SED values and thickness for the medial midcorpus is positive strong (r=0.634, p-value=0.09).
As it can be seen, perfect linear relationship between thickness values and calculated SED values were not found, but positive strong correlations were found for all the mandibular regions. The results suggest that the two variables, SED and thickness, are directly related and a strong relationship exist between them. More important, the thickness and the SED profiles have very similar trends for each region considered. Given the relatively small number of sites, the consistency of the results within sections is very significant. The values found at midcorpus always exceed the values found at the mandibular base. The findings are in agreement with the results published by other researchers (Daegling and Hylander, 1988, 2000). This result also confirms the experimental observation according to which the values at midcorpus are usually very different than the values obtained at the basal corpus (Daegling and Hylander, 1988). Vertical shear stress due to parasagittal bending will vary as a parabola from zero at the bottom to a maximum at the centroid (midcorpus) and zero at the top and bottom. Therefore, using stress distribution diagram, higher values were expected at midcorpus than at the base. Positive strong correlations between thickness values and calculated SED values were found for all the mandibular regions. Therefore the current model supports the functional relationship theory between the form and the function of the mandible.
Another significant result is related to lateral versus medial aspect of the mandible. Under the molars, the cortical bone is thicker on the lateral aspect than along the medial aspect. Generally, the SED values found at the lateral midcorpus are higher than the SED values for the medial midcorpus. At the symphysis and at the anterior corpus, the cortical bone is distributed exactly the opposite, being thicker on the medial aspect than along the lateral aspect. This variation is also observed from the SED values calculated for the most anterior section (C/I2). Only for this section, the SED value for the medial aspect exceeded the value for the lateral midcorpus. The mandible model was subjected to a combined loading: superposition of bending, torsion and direct shear. The shear stresses resulted from torsion and direct shear theoretically will add up on the lateral aspect and subtract on the medial aspect of the mandibular corpus (Demes 1984). According to this theory, high stress, strain and consequently SED values are expected on the lateral aspect of the mandibular corpus and lower stress, strain and SED values are expected on the medial aspect.
The strain energy density values calculated by FEA are in the same range with SED values obtained in similar studies (Mellal et al, 2004). The calculated lazy zone interval was: 0.00594 U 0.00726 (MPa). The SED values obtained in the current study range from 0.00283 U 0.01213 (MPa). Only a few of the calculated SED values are within the lazy zone interval. The values were collected from the medial midcorpus and basal area. None of the SED values collected from the lateral midcorpus fall within the lazy zone interval. The majority of the SED values are higher than the values comprises in the lazy zone interval, especially for the lateral and medial midcorpus. Higher SED values usually predict bone gain. The majority of the SED values collected from the base are lower than the values comprises in the lazy zone interval. Lower SED values usually predict bone loss. Therefore, the current model predicts that the remodeling activity is present in all the mandibular regions.
According to this model, the mandibular bone is not in an equilibrium state or in bone homeostasis. This model does not support the theory according to which the mandible is an optimized structure. The values were not uniform throughout the bone or very close to the values within the lazy zone interval. However, non-uniform SED values may not be an indication of a high bone turnover; it could also signify that the mechanostat theory is not valid for the mandible.
On the other hand, the calculated lazy zone interval and consequently the calculated strain energy density interval depend strongly on the parameters chosen, the relationship between the elastic modulus and the apparent density, and the values selected for the elastic modulus. The calculated lazy zone interval could be too narrow. Unfortunately the predictions for the bone gain or bone loss are strictly made in relation to the lazy zone interval.
The current research project explored using finite element methods a very important and controversial topic: the unique morphology of the mandible. A complex FE model was developed and analyzed by considering the geometrical complexity of the mandible, heterogeneous material properties and mastication muscles. Thickness and SED values were compared in order to determine their relationship. Positive strong correlations were found for SED values and thickness for all the mandibular regions. Therefore, the model supports the functional relationship theory between mandibular form and function. According with the SED profile, the mandibular bone is not in homeostasis and consequently mandibular remodeling is predicted. Moreover, the mechanostat theory might not be applicable to the mandible. The current project aims to investigate the controversial relationship between form and function and extend our knowledge in mandible biomechanics.
Table 1 Regional cortical thickness for macaque jaws. Units in mm. (from Daegling DJ 1993. The relationship of in vivo bone strain to mandibular corpus morphology in Macaca fascicularis. J Hum Evol 25: 247–269).
Table 2 SED values and cortical thickness data for the lateral midcorpus region.
Table 3 SED values and cortical thickness data for the basal region.
Table 4 SED values and cortical thickness data for the medial midcorpus region.
Figure 1 Calculated lazy zone interval. The adapted window or the lazy zone is the functional strain or strain energy density interval for which bone homeostasis is predicted. Homeostasis means that bone resorption and bone formation are in equilibrium. Values higher than the values within the lazy zone predict bone gain and values lower than the values within the interval predict bone loss.
Figure 2 The unconstrained lever model. The unconstrained lever model uses equilibrium equations to determine the relationship between the forces in the mandible. The muscle resultant force (M), the bite force (B) and the joint reaction forces (JW and JB) are determined through spatial relationship. The model assumes that the working and the balancing side muscle forces are equal and the muscles resultant force lies in the midline.
Figure 3 The masseter-pterygoid sling and the temporalis muscles simulation. Top and lateral view of the FE mandible. The insertions of the masticatory muscles were simulated by points on the mandibular angles and the coronoid process of the mandible, on each side. Four loads of 75N each were applied to the mandible model to account for the left temporalis muscle, the left masseter-pterygoid sling, the right temporalis muscle and the right masseter-pterygoid sling. The masseter-pterygoid sling vector was inclined 45 in the frontal plane and 68 from the occlusal plane. The temporalis vector was inclined 60 from the occlusal plane.
Figure 4 Mandible model subjected to combined loading. The model was constrained at condyles. An equal number of nodes (three nodes) were totally constrained at condyles on each side. The mastication force is tilted by 15 in the frontal plane, toward the right side of the mandible. The model was subjected to combined loading: simultaneous application of bending, torsion and direct shear.
Figure 5 Mandibular sections: M3 (third molar), M2 (second molar), M1 (first molar), P4 (fourth premolar), P3 (third premolar), C/I2 (canine/incisor) (from Daegling DJ 1993. The relationship of in vivo bone strain to mandibular corpus morphology in Macaca fascicularis. J Hum Evol 25: 247–269).
Figure 6 SED and thickness data for various mandibular regions: a) lateral midcorpus region. For the lateral region, the SED values tend to increase slightly from the anterior corpus toward the molar region; b) basal region. The SED values recorded for the basal region also tent to increase slightly from the anterior corpus toward the molar region, with small values in the anterior corpus and higher values toward the molar region; c) medial midcorpus region. The overall trend for the medial region is the decrease of the SED values from the anterior corpus toward the molar region.
Figure 7 SED profile along midcorpus. An arbitrary path was created with the finite element analysis package along the midcorpus and the overall SED profile was plotted along the path. a) Mandible model and an arbitrary path. b) SED values were plotted along an arbitrary path created on the mandibular corpus.
Figure 8 Regional SED values: SED values are consistently higher at a) the lateral midcorpus than basally and b) at the medial midcorpus than basally. Vertical shear stress due to parasagittal bending will vary as a parabola from zero at the bottom to a maximum at the centroid (midcorpus) and zero at the top. Therefore, using stress distribution diagram, higher values were expected at midcorpus than at the base. Mandibular sections: M3 (third molar), M2 (second molar), M1 (first molar), P4 (fourth premolar), P3 (third premolar), C/I2 (canine/incisor).
Figure 9 Regional SED values: SED values are generally higher at the lateral midcorpus than at the medial midcorpus. The only exception is at symphysis (C/I2) where the medial SED values are higher than the lateral SED values. The same variation is observed in cortical bone distribution: under the molars, the cortical bone is thicker on the lateral aspect than along the medial aspect. At the symphysis and at the anterior corpus, the cortical bone is distributed exactly the opposite, being thicker on the medial aspect than along the lateral aspect. Mandibular sections: M3 (third molar), M2 (second molar), M1 (first molar), P4 (fourth premolar), P3 (third premolar), C/I2 (canine/incisor).
Figure 10 Regional SED values and the calculated lazy zone interval. Only a few of the calculated SED values are within the lazy zone interval. The values were collected from the medial midcorpus and basal area. None of the SED values collected from the lateral midcorpus are within the lazy zone interval. The majority of the SED values are higher than the values comprises in the lazy zone interval, especially for the lateral and medial midcorpus. The majority of the SED values collected from the base are lower than the values comprises in the lazy zone interval.