Importance of the variable periodontal ligament geometry for whole tooth mechanical function: A validated numerical study

Journal of the mechanical behavior of biomedical materials 67 (2017) 61–73

Contents lists available at ScienceDirect

Journal of the Mechanical Behavior of Biomedical Materials journal homepage: www.elsevier.com/locate/jmbbm

Importance of the variable periodontal ligament geometry for whole tooth mechanical function: A validated numerical study

MARK

⁎

Anneke Nikolausa, John D. Curreyc, Tom Lindtnerb, Claudia Flecka, , Paul Zaslanskyb,⁎⁎ a b c

Department of Materials Engineering, Institute of Technology Berlin, 10623 Berlin, Germany Department for Restorative and Preventive Dentistry, Charité – Universitaetsmedizin Berlin, 13353 Berlin, Germany Department of Biology, University of York, York YO10 5DD, United Kingdom

A R T I C L E I N F O

A BS T RAC T

Keywords: Tooth displacements PDL strain Alveolar bone stress Finite element analysis Mastication response Three-rooted tooth-PDL-bone complex

When mammalian teeth breakdown food, several juxtaposed dental tissues work mechanically together, while balancing requirements of food comminution and avoiding damage to the oral tissues. One important way to achieve this is by channeling mastication forces into the surrounding jaw bone through a thin and compliant soft tissue, the periodontal ligament (PDL). As a result, during a typical chewing stroke, each tooth moves quite substantially in its anchor-site. Here we report a series of experiments, where we study the reaction of threerooted teeth to a single chewing event by finite element (FE) modelling. The nonlinear behaviour of the PDL is simulated by a hyperelastic material model and the in silico results are validated by our own in vitro experiments. We examine the displacement response of the complete tooth-PDL-bone complex to increasing chewing loads. We observe that small spatially-varying geometric adjustments to the thickness of the PDL lead to strong changes in observed tooth reaction movement, as well as PDL strain and bone stress. When reproducing the regionally varying thickness of the PDL observed in vivo, FE simulations reveal subtle but significant tooth motion that leads to an even distribution of the stresses in the jaw bone, and to lower strains in the PDL. Our in silico experiments also reproduce the results of experiments performed by others on different animal models and are therefore useful for overcoming the difficulties of obtaining tooth-PDL-bone loading estimates in vivo. This data thus enhances our understanding of the role the variable PDL geometry plays in the tooth-PDL-bone complex during mastication.

1. Introduction Teeth in mammals are intensively used working tools, made to reduce and comminute food so that it can be further processed by the gut. They usually have a crown consisting of one or several cusps that protrude into the oral cavity, and one or multiple roots that anchor the teeth into the jawbone, through a specialised ligament attachment system (Nanci, 2008). The outer regions of the crown consist of a hard and wear-resistant highly-mineralised enamel layer covering a bonelike dentine core (Kinney et al., 2003). The roots are also made of dentine and are lined with cementum which is a low density bony tissue (Ho et al., 2009). Teeth are bound to the jawbone via cementum and through a thin 'membrane', the periodontal ligament (PDL), made of collagen fibres surrounded by a gel-like matrix of proteoglycans and cells (Berkovitz, 1990). The compliant nature of the PDL allows some mobility of teeth in their sockets. When muscles apply forces for the purpose of chewing, the lower

⁎

jaw moves and brings the teeth in both upper and lower arches into contact with the chewed substance. The forces crush the food (Mohl, 1988, Lucas, 2004) and any excess strain-energy that is not consumed by comminution, is distributed across the crown and into the roots. This leads to deformation of the PDL and, consequently, dissipation into the supporting tissues. This supporting and cushioning role of the PDL may also circumvent irreversible damage to the cyclically-loaded tooth materials, and may thus contribute to extending functionality of the teeth. Usually additional mastication cycles follow the first chewing stroke, such that reduction of the food particles takes place, until the food is ready to be swallowed. Teeth vary in size and shape (Lucas, 2004) and consequently the load distributions in enamel, dentine or cementum, as well as the overall response of the entire tooth-PDL-bone complex to mastication loading (Bakke, 2006; Chen, 2009), are not simple to understand. This is further complicated by the existence of multiple roots, particularly in the distal more highly loaded regions of the dental arch. The usual

Correspondence to: Institute of Technology Berlin, Materials Engineering, Straße des 17. Juni 135 - Sekr. EB-13, 10623 Berlin, Germany. Corresponding author. E-mail address: claudia.fl[email protected] (C. Fleck).

⁎⁎

http://dx.doi.org/10.1016/j.jmbbm.2016.11.020 Received 11 March 2016; Received in revised form 1 November 2016; Accepted 24 November 2016 Available online 25 November 2016 1751-6161/ © 2016 Elsevier Ltd. All rights reserved.

Journal of the mechanical behavior of biomedical materials 67 (2017) 61–73

A. Nikolaus et al.

physiological forces acting on teeth arise due to the action of muscles, and are applied approximately along the tooth long axis. Reported maximal chewing loads for common foods measured on human molars are somewhere between 70 N and 150 N, based on direct measurements in teeth of human volunteers (Anderson, 1956; Boldt et al., 2012; Mioche and Peyron, 1995), but maximum voluntary biting forces are known to exceed 500 N (e.g. Raadsheer et al., 1999; Bakke, 2006; Waltimo and Könönen, 1995). Average chewing loads in pigs are somewhat higher than for humans: biting forces of 240–510 N have been measured on porcine third premolars (Bousdras et al., 2006). Chewing loads lead to displacements of teeth in their sockets (e.g. Boldt et al., 2012; Gathercole, 1987; Naveh et al., 2012a; Parfitt, 1960) with an overall non-linear response of the tooth-PDL-bone complex. The latter has no simple analytic description but has been addressed by 3D imaging and numerical simulations under variable simplifying assumptions. Naveh et al. (2012a) reported measurements obtained from a loading system situated within a micro computer tomography (μCT) scanner. Following the motion of a rat molar in the bone socket in 3D, they identified specific contact areas between the root surface and the jaw bone, appearing at increased loads. This resulted in a deflection of the tooth as it moves into the socket, described as a 'seesaw' tooth motion. Interestingly, tooth tilting following loading was also observed in compression on a bucco-lingual tooth slice from a pig premolar (Qian et al., 2009) both in laboratory experiments and in computer simulations. Ho et al. (2013) report seemingly matching root-bone surfaces in human molars, areas that come into contact under compressive load. Finite element (FE) methods have been in engineering use for some time now, and have also found extensive application in studying the outcome of loads on teeth, reported as early as 1973 (Farah et al., 1973). A PubMed search of the key words 'finite element' and 'tooth' results in more than 1000 published papers, with the number ever increasing. The FE method is often used to study tooth-biomaterial interactions, e.g. the influence of restorations and implants on load distribution in and around treated teeth. Notably, however, there are differences between different FE models (e.g. 2D versus 3D), and despite important insights arising from many studies, only few address normal function of whole, intact teeth under mastication-relevant conditions. Furthermore, most of these studies addressed single-rooted teeth, and very little is known about three rooted teeth, such as upper human molars. Contemporary technical developments provide means to reproduce, in great detail, all the components of the tooth-PDL-bone complex, and to study the reactions of these components to load. Finer details make it possible to assess previously hidden aspects of the load transfer dynamics as well as reveal deformation and response-to-load of every zone in each of the tooth-PDL-bone components. Essentially all of the dental FE models that have analysed intact tooth responses to load originated from the field of orthodontics, as listed in Table 1. Orthodontic forces (Ren et al., 2003) are however significantly lower than mastication forces, and they are static, acting over much longer timescales than mastication forces. Further, they are often applied in non-physiological trajectories and they induce permanent repositioning of teeth within the jawbone. Such low magnitude long-acting loads allow the PDL materials to flow over time, thereby decreasing the acting forces, and they result in bone remodelling (resorption and formation) in the jaw. Thus, tooth displacement in the socket due to application of orthodontic forces (e.g. Cattaneo et al., 2005; Chen et al., 2014; Field et al., 2009; Huang et al., 2016; Jones et al., 2001; Lombardo et al., 2014; Middleton et al., 1996) is significantly different from the reversible tooth displacement occurring during mastication, although all these tooth movements always lead to substantial deformation of the PDL (Hohmann et al., 2011; Pietrzak et al., 2002). To describe the complex non-linear behaviour of the PDL, specific material models were developed (e.g. Favino et al., 2013; Huang et al., 2016; Natali et al., 2007; Poppe et al., 2002; Toms and Eberhardt, 2003; Wang et al., 2012). Two of these models were used to

Table 1 PDL thickness and material models used in published FE models of teeth. Where linear elastic material models were used, the Young's modulus (E) and Poisson's ratio (ν) are listed. If several FE models with different PDL thicknesses have been analysed they are separated by bullet points. Ref.

PDL thickness

FE models of orthodontic loading Cattaneo et al., 2005 Non uniform

Dorow and Sander, 2005 Hohmann et al., 2011

Huang et al., 2016 Jeon et al., 1999 McCormack et al., 2014

0.2 mm – 0.05 mm – 0.1 mm – 0.2 mm – 0.3 mm – Non uniform 0.2 mm 0.2 mm 0.2 mm

Natali et al., 2007

Non uniform

Poppe et al., 2002 Provatidis, 2000

Non uniform 0.229 mm

Toms and Eberhardt, 2003

– Non uniform – 0.25 mm

Vollmer et al., 1999

0.2 mm

Wang et al., 2012

0.25 mm

FE models of mastication loading Keilig et al., 2016 0.2 mm Menicucci et al., 0.25 mm 2002 Natali et al., 2004 0.2 mm – 0.2 mm Ona and – 0.4 mm Wakabayashi, 2006 Pietrzak et al., 2002 0.2 mm Qian et al., 2009 Rees, 2001 Ren et al., 2010 Schrock et al., 2013

Non uniform 0.3 mm 0.2 mm Non uniform

Tuna et al., 2014

0.18 mm

Material model PDL

Non-linear Linear elastic (E=0.044 MPa, or 0.17 MPa, ν=0.45) Linear elastic (E=0.01 MPa, 0.1, or 1, ν=0.45) Linear elastic (E=0.1 MPa, ν=0.45)

Hyperelastic Linear elastic (E=0.667 MPa, ν=0.49) Fibre matrix model (Ematrix=1 MPa, νmatrix=0.45, Efibre=1000 MPa, νfibre=0.35) Anisotropic hyperelastic constitutive model Bilinear Linear-elastic (E=0.68 MPa, ν=0.49) and different fibrous models Linear elastic (E=0.303 MPa, 0.208 MPa, 0.143 MPa, 0.179 MPa, or 0.25 MPa, ν=0.45) Non-linear elastic Bilinear elastic Deviatoric viscoelastic Volumetric viscoelastic Tension-compression volumetric viscoelastic Several bilinear elastic models 3D non-linear visco-elastic spring elements Non-linear constitutive law One model consisting of a non-linear elastic and linear elastic phases Large strain non-linear elastic isotropic 'split' law Non-linear, viscoelastic Linear elastic (E=50 MPa, ν=0.49) Linear elastic (E=68.9 MPa, ν=0.45) Linear elastic (E=1 MPa, 3 MPa, 5 MPa, 7 MPa, 9 MPa, or 11 MPa, ν=0.45) Non-linear contact model

examine low mastication loads (up to 3 N) on a canine (Pietrzak et al., 2002; Natali et al., 2004), and although different approaches were employed, both studies showed results in good agreement with values determined experimentally (Parfitt, 1960). None of the mastication FE studies, known to the authors, provide information about the influence of the PDL geometry due to realistically high mastication forces ( > 50 N), which is of high relevance for multi-rooted teeth. Due to its complex architecture and non-linear deformation behaviour, modelling of the tooth-PDL-bone complex in multi-rooted teeth requires paying attention to substantial geometric and technical details. Once created and validated, however, such models are able to numerically reproduce in silico a whole range of observations related to responses to loads, providing insights that are not available even from a very large range of real-world experiments. μCT scans of whole toothPDL-bone segments make it possible to model the tooth-PDL-bone geometry in great detail. This includes not only the major enamel and dentine tissues, but also the cementum or trabecular regions of the bone. The PDL exhibits non-linear elastic behaviour, with large 62

Journal of the mechanical behavior of biomedical materials 67 (2017) 61–73

A. Nikolaus et al.

deformations, the latter violating the basic assumptions of simulations using linear-elastic theory. Consequently, other material models such as viscoelastic (Fung, 1993; Zhang et al., 2007) or hyperelastic (Martins et al., 2006; Weiss et al., 1996) models must be used. Viscoelastic models specifically provide interesting insights into time- and velocitydependent aspects of deformation, but they require corresponding experiments for validation, which is not necessary when investigating and modelling single mastication events at one single chewing speed. For the latter, hyperelastic models provide excellent workarounds. From all of the above, it appears that during loading, and due to the PDL suspension mechanism in the bone, teeth tilt while moving along defined paths, guided by contact points between root and bone (Naveh et al., 2012a; Qian et al., 2009; Ho et al., 2013). We therefore hypothesize that in multi-rooted teeth, a predefined functional design exists, directing tooth motion in the sockets. We assume that the PDL geometry has an important role in this functional design, the inner mechanisms of which are largely not understood. Consequently, in this work, we report on the effects of the PDL geometry for the mechanical response of a highly loaded three-rooted tooth, studying real and simulated fresh porcine tooth-PDL-bone complexes. We center on understanding the response to initial mastication loading, commencing at the onset of food comminution with the early phase of the first chewing stroke. We examine how the geometry of the PDL influences the stress and strain distribution into the most highly-loaded (furcation) bone tissue. A series of FE models is used to assess the impact of subtle geometry changes in the PDL architecture. The most prominent findings arise from comparisons of a simple PDL geometry of uniform thickness with geometries of an intricate PDL microarchitecture of variable thickness, as is observed in nature. We validate our computer simulation results by comparison with laboratory measurements of deformation of a series of fresh mechanically-loaded pig teeth. Consequently we are able to numerically reproduce tooth displacement patterns and tooth-bone contact patterns typical of the first cycle of mastication, which helps us understand design principles and attributes of this important component of the oral cavity. Fig. 1. Generation of our tooth-PDL-bone FE models: (a) 3D view and typical sections through reconstructed tomographic data of a μCT scan of a three-rooted pig tooth; (b) Exploded view of the different entities identified in the structure and forming the main components of the FE mesh of the tooth-bone complex; (c) 3D rendering of the integrated FE model of the whole tooth in the bone. Outermost elements of the mesh are shown on the left-hand side and are given partially transparent shading on the righthand side, to help visualise the components of the system in 3D. The x-, y- and zdirections of the Cartesian coordinate system applied for the FEA correspond with the mesio-distal, bucco-lingual and occlusal-apical directions of the tooth, respectively, as indicated.

2. Materials and methods 2.1. Loading of entire tooth-PDL-bone complexes Four lower hemi-mandibles of adult sows containing permanent teeth were obtained from a nearby abattoir shortly after butchering, and were mechanically tested while still fresh (all preparation and testing was performed within 24 h after the animals were slaughtered, and the jaws and teeth were kept wet throughout storage, preparation and testing). After removing all external soft tissue (gingivae) to facilitate gripping, each jaw was embedded in quick-set dental gypsum and mounted upright with teeth oriented to replicate the in vivo occlusal plane. This ensured reproducible tooth loading, with a main force component applied orthogonal to the occlusal plane. Six premolars were mechanically tested in different jaws. Compression was applied with a vertical steel rod under displacement control, with a rate of 50 mm/min (see Appendix A - Determination of loading velocity) in a UTS universal test machine (type 009.00, UTS Testsysteme GmbH, Ulm, Germany, crosshead-displacement measured with a resolution of 0.001%, span 1220 mm, force determined using a 2 kN sensor with a resolution of 0.01%).

an effective pixel resolution of 12 µm. Data were reconstructed using Nrecon (v. 1.6.4, Skyscan, Kontich, Belgium), cropped, downscaled to obtain isotropic 120×120×120 μm3 voxels and segmented (separated and labeled) manually using ImageJ/Fiji (1.47 n, Rasband, 19972015). 2.2.2. FE mesh creation from µCT data The tooth and bone components of the tooth-PDL-bone complex were identified based on characteristic mineral densities and geometries observed in the µCT scan (see Fig. 1a and Fig. 2b and c). The lowdensity interphase layer between enamel and dentine, which we name 'sub-DEJ dentine' (Zaslansky et al., 2006) is poorly reproduced in the 3D tomographic data due to limitations of the density resolution of the polychromatic µCT X-ray source (Salmon and Liu, 2014). Thus, 3D outlines of enamel and dentine in the crown were used to define the outer bounds of the lower density zone between the two higher-density materials. The segmented data of enamel, dentine and bone was smoothed and meshed with tetrahedral elements using Amira (v. 5.4.3, FEI, Hillsboro, Oregon, USA) whereas the sub-DEJ dentine geometry was generated by filling the gap to match the inner surface of enamel and the outer surface of crown dentine, using HyperMesh (v.

2.2. Finite element model creation 2.2.1. Sample preparation and 3D imaging A three-rooted, average sized premolar with parts of the adjacent teeth and the surrounding bone of another jaw were cut out of the fresh mandible, and stored in a 70% ethanol-water solution for preservation during imaging. The sample was scanned in a lab-µCT (Skyscan 1172 micro-CT, BrukerCT, Kontich, Belgium) with an energy of 100 keV and 63

Journal of the mechanical behavior of biomedical materials 67 (2017) 61–73

A. Nikolaus et al.

Fig. 2. Images of the 3D μCT data used for model creation and resulting FE model. (a) Typical 3D renderings of the FE model tooth observed from a distal and from a buccal viewpoint. (b and c) Slices in the μCT volume reconstruction data of the three-rooted pig tooth and corresponding 2D sections through the FE model exhibiting the 'variable PDL': (b) bucco-lingual section, (c) mesio-distal section. The dashed lines indicate where the slices (b) and (c) intersect. The assigned material properties for the FE model are identified and listed (E - Young's modulus, ν - Poisson's ratio, μ1, α1, and D1 - Ogden coefficients (for hyperelastic modeling based on Ogden (1972)). The references used are listed in Table S2 and marked with an asterisk (*). Measurements used to determine the width of the PDL are indicated, the values are given in mm.

12, Altair Engineering Inc., Troy, Michigan, USA). Cementum, which essentially has the same density as dentine, was defined for simplicity to span the outer 0.2 mm surface layer of each of the three roots. The soft low-density PDL is transparent in the µCT data obtained at such energies. Two different PDL mesh geometries were created between the roots and the bone, as described in the following section. To represent the microstructure of the cortical and cancellous bone in our FE models, the cortical bone area and struts of the trabecular bone were modelled ‘as bone’ according to their geometric appearance (see Fig. 1 c half transparent right-hand side of the bone) based on the µCT data. Therefore we did not homogenise the material properties of the cancellous bone volume. Parts of the neighbouring teeth on either side of the studied premolar were reproduced in our FE models to be able to examine effects of bone deformation beyond the zone immediately surrounding the simulated premolar. The enamel caps of these neighbouring teeth were removed to reduce the number of elements in the models, because of the large gaps existing between adjacent premolar pig teeth and because we found that this does not influence any of our results. Fig. 3. Results of experimental loading of six freshly-mounted three-rooted pig teeth in vitro: each tooth was loaded occlusally similar to in vivo loading conditions, with the tooth fixed within the bone, as depicted schematically by the inset graphic. The loading rate was 50 mm/min. Load displacement responses of left and right second and fourth premolars were used to establish three archetypical load-displacement curves. These correspond to a regular, intermediate (fit mean) response, a second, more compliant case (fit soft) and a relatively resilient response (fit hard). The three idealised tooth-PDL-bone mechanical loading scenarios (dashed lines) were used as a basis for FE simulations. The fits are described by the given exponential growth functions with the listed constants.

2.2.3. Uniform versus variable PDL geometries The width of the PDL, measured in the original µCT data, varies from region to region of the root by an order of magnitude, spanning ~50 to ~500 µm (see Figs. 3a and b) with a median thickness of 210 µm. The median was determined by measuring the PDL in the µCT data in 21 slices evenly distributed over the root height at a total of 613 locations. Two model prototypes were created in order to test the significance of reproducing regional differences in the thickness of this 64

Journal of the mechanical behavior of biomedical materials 67 (2017) 61–73

A. Nikolaus et al.

dentine obtained from the buccal sides of teeth is E=3.5 GPa, whereas the average for lingual sides is E=9.7 GPa (Zaslansky et al., 2006). Consequently, an average of E=6 GPa was used for the less-mineralised sub-DEJ dentine. The mechanical property differences in dentine occur mainly along the longitudinal axis of the tooth (see Table S2). To partially account for the significant local differences known to exist in the Young's modulus of dentine yet to avoid making the model too complex, we chose to differentiate between crown (E=20 GPa), cervical (E=15 GPa) and root (E=10 GPa) zones in dentine. The trabecular structure of the bone was modelled representing the true geometric allocation of the trabecular struts (see transparent half of Fig. 1c). Both cortical bone and trabecular bone within the struts have reasonably similar elastic properties (Turner et al., 1999). Therefore, we assigned an average value of E=13 GPa (Ashman and Rho, 1988). The Poisson's ratios used are listed in Fig. 2 (Benazzi et al., 2012, Dorow and Sander, 2005).

soft membrane. In the 'uniform PDL' FE models, a PDL with a uniform width of the median thickness of our sample (210 µm) was generated surrounding the outer rim of cementum, using HyperMesh. Because 200 µm is a PDL thickness frequently used in FE tooth studies in the literature (Table 1), results of a model with a 200 µm thick uniform PDL are also provided. The mesh of the bone was then locally adjusted to seamlessly connect with the PDL component of the model. For the second model, termed 'variable PDL' FE model, a significant amount of geometric detail was reproduced, and the PDL was assumed to fill the uneven gap observed by µCTbetween the three roots and the bone including the furcation zone beneath the crown. Zones where the bone model was eroded due to the automatic mesh creation steps were corrected manually, to match the original tomographically-imaged dimensions. To ensure fidelity, the final model was scrutinised and compared with the original width of the PDL, as measured at many locations in the reconstructed µCT volume. Two slices exemplifying the goodness of fit of this FE model and the scan are shown in Fig. 2.

2.2.6. Mechanical property assignment to the PDL Following mechanical loading, teeth displace largely due to deformation of the PDL. Much evidence exists that the PDL experiences large deformations when teeth are loaded (Naveh et al., 2012b; Yoshida et al., 2001; Dorow et al., 2003). Hyperelastic material models (Ali et al., 2010) are useful to simulate large non-linear strains (Weiss et al., 1996; Martins et al., 2006) under constant strain rates, and consequently such models were used to study the initial PDL response to load in our in silico experiments. A hyperelastic 'First Ogden' model (Ogden, 1972) was used, following comparison against several other FE model approaches, including linear-elastic (Tanne et al., 1987), viscoelastic (Su et al., 2013) and various hyperelastic models ('Marlow' (Marlow, 2003), 'Arruda-Boyce' (Arruda and Boyce, 1993), 'MooneyRivlin', 'Polynomial'). For additional simulation details see Appendix A - Hyperelastic material model for the PDL. The PDL properties were validated based on our mechanical tests. Through iterative adjustments, we found the material parameters μ1=0.4, α1= 6.4, and D1=0.1 to best reproduce tooth movement for the 'variable PDL' model. Using these parameters, we analysed the difference between simulation results and results obtained experimentally (see next section).

2.2.4. Simulation conditions All simulations were solved using Abaqus 6.14-4 (SIMULIA, Providence, Rhode Island, USA). The calculations were performed using the North-German Supercomputing Alliance hardware (HLRN, Zuse Institute, Berlin, Germany). The 'variable PDL' and the 'uniform PDL' FE models contain 6,276,898 and 4,760,240 elements, respectively. For all FE simulations, symmetry boundary conditions were imposed in the x-direction on the mesial and distal sides of the model. This allows free tooth movement in all directions, thus reproducing the mechanical environment of the tooth-bone-complex as part of a jaw segment in the best achievable way. The lower rim of the tooth-jaw complex (x-y-plane) was defined as having fixed boundary conditions. After assignment of mechanical properties (see following section), simulations were performed using a displacement velocity of 50 mm/ min (see Appendix A - Determination of loading velocity) in z-direction (occluso-apical). This load was delivered to the tooth over a small wear facet extending over 1150 elements on the uppermost surface of the crown (see Fig. S1). During the simulations, increasing loads were applied until a reaction force of 400 N was reached. The convergence of the mesh (Erdemir et al., 2012) was successfully tested for all models. Further technical details of the simulation settings are given in Appendix A.

2.2.7. Validation by the results of the mechanical tests The force displacement responses of six teeth, loaded up to 400 N load, are shown in Fig. 3. Some teeth appear to have a 'softer' response, while others have a 'harder' response, exhibiting a more bowed or steeper graph, respectively. These data were used to establish the load response of 'typical' tooth-PDL-bone complexes of different stiffnesses. A mean value was calculated from all measurements and fitted with the exponential growth function F =A1 +A2 e−z / φ where F is the loading force, A1, A2 and φ are constants and z represents the tooth displacement. To account for the range of stiffnesses observed in different tooth-PDLbone complexes, two additional exponential growth fits were derived, termed 'fit hard' and 'fit soft', fitted to the hardest and the softest measurements. The resulting fits and the corresponding fit parameters are also shown in Fig. 3. For additional validation considerations see Appendix C. Comparing the in vitro and in silico forces at the same applied displacements, we find a small difference of 4.3 N or less (mean difference=2.5 N). Fig. 4 shows that as expected the FE model 'variable PDL' excellently reproduces the 'fit mean' in vitro load displacement reaction.

2.2.5. Mechanical property assignment for tooth and bone The material properties used in our FE models are based on a range of values found in the literature (see Table S2) and we used intermediate values (marked by *) listed in the graphic representation in Fig. 2. For simplification, all materials were modelled as isotropic and homogeneous and all the mineralised tissues (tooth and bone) were modelled as linear-elastic materials, a reasonable approximation allowing us to focus this work on the PDL-geometry effects and on the dynamics during the first cycle of mastication. While the elastic material properties of the tooth and bone components are relatively well agreed-on, much less certainty exists regarding the best way to represent the mechanical behaviour of the PDL (Fill et al., 2011; Fill et al., 2012; Minch, 2013) and how to model it. Some technical aspects are presented in Appendix A and a survey of considerations related to the choice of tooth-bone material properties is given in Appendix B and summarised in Supplementary Table S2. Interestingly, even though the properties in Table S2 are from a range of different mammals, they are similar, supporting the notion that teeth differ by shape rather than by changes to the material properties (Constantino et al., 2012). Preliminary tests using models employing different literature values revealed that most materials properties variations other than for the PDL had little effect on our displacement simulation results. We assumed an average Young's modulus (E) of 80 GPa for enamel and of 5 GPa for cementum. The average modulus for human sub-DEJ

3. Results 3.1. Tooth movement under masticatory load The 'variable PDL' and 'uniform PDL' models exhibit increasing, albeit different tooth displacements with loading. An overall stiffer load displacement response is seen in the case of the 'uniform PDL', 65

Journal of the mechanical behavior of biomedical materials 67 (2017) 61–73

A. Nikolaus et al.

Fig. 4. Force-displacement results for the uniform (dotted lines) and variable (solid line) PDL thickness models obtained in silico. The mean in vitro fit (dashed line) is given for comparison with the simulations and reveals the excellent match between our simulation and the results obtained in the lab. The solid (red) curve is the load displacement response of the FE model with the 'variable PDL' and the significantly stiffer (dotted light grey and dark purple) lines are for the models with the’uniform PDL'. Theshaded region shows the range between the in vitro 'fit hard' and 'fit soft' (see Fig. 3).

although we used the same materials properties in all components. This is best demonstrated by direct comparison of the displacements calculated with the 'uniform PDL' models versus the 'variable PDL' model at the same forces (Fig. 4). The reaction with a 210 µm thick PDL is only about 4% softer than with a 200 µm thick PDL. The 10 µm difference in PDL thickness in the ‘uniform PDL’ models has therefore only minor effects on the results. The displacements along the three Cartesian directions (X, Y, Z), approximately coinciding with the mesio-distal, bucco-lingual, and occluso-apical tooth orientations, are shown for the 'variable PDL' model at a load (reaction force) of 150 N in Fig. 5a–c and for the 'uniform PDL' models in Fig. 5d–f. The latter reveal overall smaller displacements than those observed in the ‘variable PDL’ model (results from the ‘200 µm uniform PDL’ model are very similar to results from the ‘210 µm uniform PDL’ model). As expected, due to the loading set-up, the largest displacements occur approximately along the long tooth axis, corresponding to the Zdirection (Fig. 5c and f). In both models, the tooth moves very little in the mesio-distal (X) direction (Fig. 5 a and d) despite lacking contacts with neighbouring teeth. We observed no more than 3% of the displacement magnitude seen for the Z-direction. Even for higher loads exceeding 400 N, this movement is far below the motion needed to contact the neighbouring teeth (note gaps to neighbouring teeth seen in Fig. 2c). The removal of neighbouring enamel from the FE model, therefore, does not influence our results in pig premolars and we conclude that the healthy three-rooted tooth does not significantly affect its neighbours when loaded occlusally. Importantly, however, the ‘variable PDL’ differs from the ‘uniform PDL’ models in the amount of bucco-lingual tooth rotation observed. Notably, in the 'variable PDL' model, in addition to the Z-direction, a significant amount of tooth motion is also observed along the Ydirection, far greater than the movement observed in the 'uniform PDL' models (Fig. 5b and e). This difference between the models is also observed in the large gradient of displacements observed along the Zdirection of the ‘variable PDL’ as compared to the ‘uniform PDL’ models. Together this suggests that a significant inwards (lingual) rotation takes place. Tooth positioning in the jaw certainly plays a role

Fig. 5. Graphic representation of the different tooth displacements simulated using the 'variable PDL' (left column, panels (a–c) and the 200 μm 'uniform PDL' (right column, panels (d–f) FE models. Displacements in the three coordinate directions Ux on the top, Uy in the middle, and Uz on the bottom, are shown for a load of 150 N. The dotted area marked in (a) corresponds to the position of the images and zones of analysis shown in Fig. 6a.

(see Fig. S1). Note that in all cases, the three-rooted pig tooth is loaded as similar as possible to the natural setting in the mouth during chewing. Despite the identical loading axis and tooth orientation in the uniform and variable models, significant differences in the Y-displacements are seen. Thus, the differences in the geometry of the two PDLs variable versus uniform - lead to marked differences in the initial mastication load-induced displacements, mainly appearing in the

66

Journal of the mechanical behavior of biomedical materials 67 (2017) 61–73

A. Nikolaus et al.

Fig. 6. Stress and strain response in the bone and PDL, in the variable and uniform PDL FE models: (a) Identification of the regions ‘furcation’ and ‘lateral’ used for calculation of mean stress and mean strain, in a magnified 2D view of the area indicated by the dashed line in Fig. 5a. Note the subtle PDL geometry differences in the 'variable PDL' FE model (a - left-hand side) and the 'uniform PDL' model (a - right-hand side). The PDL in the 'uniform PDL' model has a constant thickness around the tooth roots but exhibits different thicknesses along the slice shown due to interaction between the slice and the complex geometry of the root. Importantly, the 'variable PDL' model exhibits a relatively thick PDL in the 'furcation' observed between the lingual and buccal roots and a comparatively thin PDL at the buccal side of the lingual root ('lateral'). (b,c) Mean von Mises stresses in the furcation and lateral bone regions, with increasing load on the tooth. (d,e) Mean effective strains evolving in the PDL with increasing loads. The stress/strain values are means of > 1000 elements in 3D in the volumes 'furcation' and 'lateral' and are only marked in the 2D slices in (a) for graphical clarity. The shaded grey region in the plot shown in panel (b) marks a range of stresses where we note a gradual shift of load transfer into the lateral bone area, as indicated by the change in curvature observed for both bone zones.

tissue surrounding the bone tissue between the lingual and buccal roots, and the volume between the inner side of the lingual root and the bone residing between the roots. These volumes are marked as 'furcation' and 'lateral' in the schematic Fig. 6a. Here, important mechanical interactions appear to take place between the different components of the tooth-PDL-bone complex. We find that the differences in local width of the PDL between the two FE models strongly affect the local strain/stress distributions in the PDL tissue and also in the supporting bone. Fig. 6b and c show the mean von Mises stresses in the bone elements ( > 1000 elements in

lateral, Y- (bucco-lingual) direction. 3.2. PDL strain and bone stress under masticatory load To better understand the behaviour of the tooth-PDL-bone complex, we identified regions where the most significant stresses and strains were observed. These appear to concentrate beneath the tooth crown, near the branching zone of the roots, highlighted in Fig. 6a. The importance of the PDL thickness is demonstrated by comparing the strains in PDL volumes in the most affected regions: the volume of PDL 67

Journal of the mechanical behavior of biomedical materials 67 (2017) 61–73

A. Nikolaus et al.

models containing adjacent teeth were compared with additional FE models where the remnants of the neighbouring teeth were excluded. The displacements of the tooth in its socket as well as the strain and stress distributions in the tooth-PDL-bone complex were not affected by the removal of the neighbouring teeth.

3D), in the 'lateral' and 'furcation' regions (marked in Fig. 6a) under increasing loads. In the following and for simplicity, we refer to this mean von Mises stress as 'bone stress'. For loads on the tooth up to 250 N, all models result in similar bone stress values and responses in the 'furcation' area (Fig. 6b and c dark green and dotted, 'furcation bone' marked curves). For higher loads, the two PDL geometries lead to a significantly different 'furcation' bone stress: in case of the 'uniform PDL', the 'furcation' bone stress levels off, that is, it increases with a decreasing monotonic rate (Fig. 6c dark green and dotted, 'furcation bone' marked curves). In contrast, in the case of the 'variable PDL', the graph exhibits an overall slope that appears rather linear (Fig. 6b dark green curve). Closer inspection reveals a slight 's' shape, where the initial rate of stress increase becomes somewhat milder (marked by the shaded area in Fig. 6b) at about 100 N up until about 250 N, followed at higher loads by an increase in the rate of stress build-up in the bone. The difference between both PDL geometries becomes very significant in the 'lateral' bone areas (Fig. 6b and c light green and dotted, 'lateral bone' marked curves). In the 'uniform PDL' models, the ‘furcation’ and ‘lateral’ graphs have the same shape but different values: bone stress in the 'lateral' area is markedly lower than the bone stress determined at the 'furcation', but both increase with decreasing rates (Fig. 6c). In the 'lateral' region in the 'uniform PDL' case, we observe stress levels of 35% to 40% of the bone stress levels observed in the 'furcation' area (Fig. 6c). Conversely, the 'variable PDL' model shows a very balanced force-distribution picture (Fig. 6b), with an initially moderate slope in the 'lateral' region bone stress, observed for loads below 100 N, increasing twofold in the shaded graph area. In this 'variable PDL' model, bone stress in the 'lateral' and 'furcation' areas reaches very similar values at high mastication loads. Note that overall, the stresses sustained by the near-PDL bone are higher in the 'variable PDL' model (Fig. 6b) as compared with the 'uniform PDL' models (Fig. 6c), but they do not exceed 15 MPa, i.e. they are well below known values of bone strength. In the ‘uniform PDL’ the thickness difference between the 200 µm and 210 µm thick PDL has no influence on the 'lateral' bone stress (Fig. 6c) and has only a small influence on the 'furcation' bone stress, decreasing with a thicker uniform PDL (Fig. 6c). Calculated in a similar way to von Mises stresses, the mean effective strains (see Appendix D) in the PDL volume abutting the same 'lateral' and 'furcation' bone regions (marked in Fig. 6a) are shown in the diagrams in Fig. 6d and e. Referring to this as 'PDL strain' for simplicity, it is used as an indicator for the higher bounds on deformations taking place in the PDL during the first mastication stroke that we model. In all cases, the PDL strain increases under load in an exponentially decaying form (Fig. 6d and e). Curiously however, and different to what is seen for the stresses in bone, the 'furcation' PDL strain is much lower than the 'lateral' PDL strain in the 'variable PDL' model (Fig. 6d). The opposite is seen in the 'uniform PDL' models where ‘furcation’ PDL strains are higher than ‘lateral’ PDL strains (Fig. 6e). For loads of up to 170 N the PDL strain response predicted by both models is similar in the 'lateral' area (Fig. 6d and e dark purple and dotted, 'lateral PDL' marked curves), despite the markedly different PDL geometry and the tooth motion that we observed. For loads higher than 170 N the ‘lateral’ PDL strain in the ‘uniform’ models is higher than in the ‘variable’ model and in the ‘furcation’ area between the roots, PDL strain in the ‘uniform’ models is higher than in the ‘variable’ model (Fig. 6d and e). Note how, in the 'furcation' area, the PDL in the 'variable PDL' model (Fig. 6d light purple curve) appears to be more protected, as it deforms less than 50% of what is observed in the 'uniform PDL' models (Fig. 6e light purple and dotted, 'furcation PDL' marked curves).

4. Discussion Damping and strain distribution into the supporting bones can be considered as main mechanical functions of the PDL, in both singleand in multi-rooted teeth (Mohl, 1988). The soft-tissue microstructure and the macroscopic geometry emerged via natural selection through a long evolutionary history (Gaengler, 2007; Agrawal and Lucas, 2003) suggesting that the PDL architecture confers an important survival advantage. Our results show the significance of the variable PDL geometry for the stress and strain distributions in highly loaded three-rooted teeth which we analysed during the very first loading phase at the onset of mastication. By using a hyperelastic PDL material model with a highly detailed FE mesh, we observe details of how the tooth-PDL-bone complex forms an efficient low maintenance and longlasting load delivery apparatus. We perceive this complex as a compound tool, optimised for coping with the external, occasionally large muscle forces, capable of breaking down food of variable hardness and composition, while protecting the non healing teeth against damage due to overload or fatigue. To fulfil its function properly in the more heavily-loaded post-canine regions, the tooth-PDL-bone complex sustains the forces that teeth apply to comminute food, while ensuring that excess strain energy is reliably dissipated into the supporting jaw bone tissues. This is partially achieved through tooth micro-mobility, where the roots and the PDL define routes for tooth motion into the sockets without any damage to any of the tooth-PDLbone complex components. We limit our discussion here to observations during the first stroke in a series of chewing cycles for which our FE model is validated by in vitro experiments. It is well known, that with repeated loading, regardless of whether chewing forces decline or not (Gibbs et al., 1981; Schindler et al., 1998) viscoelastic effects become increasingly important (Dorow et al., 2002, Sanctuary et al., 2005). Although time effects are beyond the scope of the present work, our results provide valuable insights into the general way that teeth move about during mastication, because we examine a range of PDL stiffnesses (see Appendix C) and a range of applied forces at a loading rate that is typical for mastication. Understanding the first stroke of the complete chewing series is a step towards comprehension of the structurefunction relations and dynamics of the tooth-PDL-bone complex during mastication. Our main findings are that (i) guiding tooth movement is an important component of the healthy tooth-PDL-bone complex function in three-rooted teeth, and is facilitated by local changes in PDL thickness. This leads to (ii) an even and seeminglyoptimised distribution of bone stresses and PDL strains, and consequently (iii) increases the endurance of the system. 4.1. PDL geometry guides tooth movement Tooth motion following the application of load differed significantly between the two simulated PDL morphology models. Although we used the same set of material parameters and almost identical geometries for all in silico experiments, the 'variable PDL' model load-displacement reaction leads to a combination of substantial axial (on-axis) and lateral (bucco-lingual) motion of the tooth, with the PDL response being far softer than for the 'uniform PDL' models (shown in Fig. 4). These results do not originate from the PDL material parameters, but are rather a true consequence of the slightly different PDL geometries. We observe that the subtle architectural differences are important, as may be seen near the root furcation (Fig. 6a). The 'variable PDL' geometry is thicker at the 'furcation' than the assumed ~200 µm that

3.3. Importance of simulating neighbouring teeth To estimate the possible influences of neighbouring tooth crowns and roots in the bone on the movement and deformation behaviour of the tooth-PDL-bone complex in our three-rooted pig premolar, the FE 68

Journal of the mechanical behavior of biomedical materials 67 (2017) 61–73

A. Nikolaus et al.

between the buccal side of the lingual root and the surrounding bone under load. We certainly expect differences in both the geometry involved (rat versus pig, molar versus premolar) and in the loading trajectories used. Nevertheless similar root-bone contact areas dominate tooth movement in both species. While the magnitude and direction of tooth displacement differ between these ex vivo and in silico evaluations, the fundamental similarity of tooth movement between the rat compression and the pig model experiments is striking. Despite the different size, habitat and diet of different mammals, and for teeth that even differ in the number of roots (four roots in the rat molar versus three roots in the pig premolar), the architectural similarity of the tooth-PDL-bone complex dominates tooth function, suggesting that guided tooth motion is an evolved functional principle. We conclude that in multi-rooted teeth, the establishment of tooth contact with the inter-radicular bone is part of the design and of the 'normal' tooth-PDL-bone function. Therefore, reproducing this movement is key to understand the structure-function relations of the normal tooth-PDL-bone complex.

we used in the 'uniform PDL' models. This appears to allow a larger extent of tooth movement both laterally and along the loading direction for the 'variable PDL' model, when compared to the 'uniform PDL' geometry. The situation is reversed when moving away from the furcation: between the inner side of the lingual root and the bone (marked as 'lateral' in Fig. 6a) the PDL becomes significantly thin: in the 'variable PDL' model the thickness here can decrease threefold, down to 70 µm, which is thinner than the 'uniform PDL'. Lateral movement for the tooth in the 'variable PDL' model is therefore restricted by this thin 'lateral' layer which presumably guides the root down the socket. Interestingly, Jang et al. (2015) reported experiments that highlight the importance of the PDL thickness in the furcation area. They showed that young rats raised on a soft diet develop a thinner PDL at the furcation which restricted tooth movement. Their results nicely fit our pig simulation findings - although they are from a very different animal. The main function of the PDL is to transfer the load from the tooth into the jaw bone and to dissipate strain energy, thereby protecting the tooth. It also includes rich innervation that forms an important part of the oral neurological feedback loop, necessary for proper function of the mastication system as a whole (Mohl, 1988; Hildebrand et al., 1995). It is likely that tooth motion activates the PDL receptors leading to neurological activity which has been shown to be proportional to tooth bucco-lingual rotation (Trulsson, 2006). This may enable sensation of the hardness of the food and the amount of force needed (Mioche and Peyron, 1995). The tooth-PDL-bone complex therefore also includes design considerations that facilitate a sensory function in the response to load. The differences in microarchitecture between the variable and uniform PDL geometries lead to markedly different contact areas between tooth and bone. In the 'variable PDL' model a contact region forms in the 'lateral' area under increasing load, while in the 'uniform PDL' models the contact region appears much higher up along the root, near the 'furcation' area. Thus, if root motion activates mechanoreceptors, the tooth-PDL-bone complex architecture may include ‘evolutionarily optimised’ locations for these sensors. Note, that in our FE model, no real contact is possible between root and bone, because the PDL material cannot vanish in the simulations, whereas in reality, the gel-like PDL shrinks substantially with time. While the tooth clearly exhibits a rocking movement in the 'variable PDL' FE model, this kind of 'seesaw motion', as previously named (Naveh et al., 2012a), is much smaller in the 'uniform PDL' model. The PDL architecture in the 'variable PDL' ensures a sliding movement of the tooth along the inner (buccal) side of the lingual root, coupled with some tilting in the bucco-lingual direction. In the 'uniform PDL' model, displacements appear to be much more restricted in all directions: displacements in the occluso-apical direction are much smaller (compare Fig. 5c and f), and in spite of a greater lateral PDL thickness, there is less movement in the bucco-lingual direction. This is significant, and is caused by formation of a ‘contact point' between root and bone in the furcation area, before the full deformation potential of the lateral PDL is attained. The comparable thin uniform PDL in the furcation area thus restricts and dominates tooth motion and is of greater significance to the tooth-PDL-bone function in the 'uniform PDL’ models than the relatively thick uniform PDL in the lateral area. We conclude that despite the overall membrane-like, quasi-uniform geometry of the real PDL, subtle architectural details, specifically thickness, are needed to ensure orderly distribution of loads, an activity that appears to be coupled with well defined tooth motion. While additional insight might also arise from simulation of the collagen fibrils as defined entities in the PDL (e.g. McCormack et al., 2014), our findings highlight the need to reproduce 'variable PDL' geometries during numerical simulations of the tooth-PDL-bone complex. Our numerical simulations not only reproduce the motion of the pig premolar, they show surprising cross-species similarity with experimental findings of the load response in the rat tooth-PDL-bone complex. Naveh et al. (2012a) observed the formation of contact areas

4.2. Bone stress and PDL strain One of our more striking results relates to the stress state in the bone surrounded by and supporting the roots (inter-radicular bone). Stresses in this region appear to vary markedly between the FE models, which is easier to observe at higher loads, and is most notable in the lateral bone areas. While the bone stress rises for the 'variable PDL' model it appears to level off in the 'uniform PDL' models (Fig. 6b and c) suggesting that in such FE simulations, strain energy must be shifting to other parts of the tooth-PDL-bone complex. In other words, the 'variable PDL' geometry is extremely important for evenly distributing stresses into bone. In this model, the stress response of the bone coincides with the bucco-lingual tooth movement. Following an initially pronounced stress increase in furcation bone, it increases with a decreasing slope (Fig. 6b grey area) while the stress-increase with load shifts to the lateral bone area (Fig. 6b), coinciding with the observed movement of the tooth in the lingual direction. This is coupled with a reduction in PDL width between the bone and the root. After a lateral ‘contact region’ is formed between the root and bone, the bone stresses in both the furcation and the lateral areas continue to increase similarly with increasing external forces on the tooth. The significance of the lateral 'contact region' is even more clear when we examine FE models using a softer PDL, as the PDL is deformed more easily (see Fig. S2) and the ‘contact region’ appears sooner. We found that in the 'variable PDL' model the PDL is not compressed as much as in the ‘uniform PDL’ models (Fig. 6d,e). We thus put forward the notion that the geometry of the variable PDL thickness in the tooth-PDL-bone complex protects the PDL from overload damage. 4.3. PDL geometry is central to the endurance of the tooth-PDL-bone complex Our findings, directly derived from the FE model, suggest that the variable geometry of the PDL has an important role in circumventing damage in both the hard and soft components of the complex. Protection to teeth arises from energy dissipation into/via the PDL and the PDL geometry plays a decisive role in the manner by which stress is distributed into the bone. We observe that with increasing loads, the bone stress continues to rise while the PDL strain levels off. Thus, damage to the PDL is prevented while the bone is increasingly loaded, though with stresses still within its load-bearing capacity: reported compressive strength values of typical bone span 130– 270 MPa (Currey, 2002). This indirectly suggests that the 'uniform PDL' models poorly reproduce the natural loading state of the bone by teeth. Compression damage to the bone would probably induce pain which would likely trigger a jaw-opening reflex leading to immediate 69

Journal of the mechanical behavior of biomedical materials 67 (2017) 61–73

A. Nikolaus et al.

of minor importance in the pig, but that the PDL needs to be modelled as accurately as possible, specifically regarding its variable geometry. In other species, e.g. humans, where tooth contacts are tighter, these probably have an additional effect on tooth motion paths and magnitudes. It is however interesting to consider that in multi-rooted teeth, the roots may nevertheless have an important guiding role in tooth motion independent of adjacent teeth.

load release. Furthermore, bone is capable of undergoing healing, unlike the tooth tissues. We conclude that the tooth-PDL-bone complex is designed to prevent irreversible damage by channeling overloads into the bone, protecting the tooth tissues from irreparable damage, such as rupture, cracking or crushing. The PDL itself is also safeguarded from ripping or crushing due to overload damage in the 'variable PDL' geometry, where a thicker PDL is found above the furcation and a thinner PDL exists at the lateral root contact areas. It is well-known in clinical dentistry that PDL regeneration does not take place without medical intervention, and remains a significant challenge (Andreasen, 2012). Thus, the PDL appears to be designed such that, ideally, damage does not occur to the soft tissue under natural, non-pathological conditions. In this manner, the design of multi-rooted teeth protects the PDL against overload and directly contributes to the endurance of the toothPDL-bone complex.

5. Conclusions This work reveals that the micro-geometry of the PDL plays a central role in the structure-function relations of the tooth-PDL-bone complex under masticatory load, specifically for three-rooted teeth. The maximum difference in the PDL geometry between the 'uniform PDL' and 'variable PDL' FE models is only a fraction of a millimetre in thickness, in some areas, and much lower in most regions. Nevertheless, both geometries lead to substantially different tooth displacements, and markedly differing bone stress and PDL strain results. Our results highlight the importance of PDL geometry, having a thicker PDL layer in the region of the 'furcation' between the roots, and a relatively thin layer in the 'lateral' area of the buccal root. The incorrect PDL thickness above the 'furcation' between the lingual and buccal roots in the 'uniform PDL' FE models does not only reduce the possible occluso-apical displacement of the tooth, but - surprisingly reduces movements in the lateral (bucco-lingual) direction and thus results in a significantly different tooth motion pattern, an underestimation of bone stresses and an overestimation of the PDL strains. We therefore conclude that the PDL geometry, specifically the thickness, is a key factor to guiding tooth movement in a reproducible and stable way in multi-rooted teeth. This suggests that the PDL function goes well beyond acting as an anchorage system similar to other entheses in the body: rather, it appears that the PDL also has a function to guide the tooth motion in the socket, acting as a joint with cartilagelike lubricant properties, facilitating establishment of contact and transfer of loads in different parts of the tooth socket. Further, the tilting movement reproduced in the 'variable PDL' model and in agreement with experimental results is a design concept of threerooted teeth, dissipating mastication strain energy by channeling overloads into the bone, that is the main component in the system able to heal, if mechanically induced overload damage should occur.

4.4. Validation and limitations Similar to other FE and experimental models, our simulations have limits. Even though we assume homogeneous, linear-elastic material properties for bone and tooth constituents (besides PDL), we learn much about the design principles of the system. By use of a hyperelastic model, we provide high spatial resolution which is important for the thin-layered architecture of the PDL. The simplifications relating to microstructural and property details (e.g. fibre orientations in the PDL, anisotropy, viscoelasticity) are of minor importance for the characterisation of tooth displacement in its socket under normal function, the focus of this work. Future studies, however, will have to address the effects of these details. The von Mises stress and effective strain parameters are used as proxies that ease understanding of the overall stress/strain states in isotropic materials. While of arguable use for biological materials, these parameters are easier to compare, and they provide a quantitative measure of the differences between the two models. More elaborate models are needed to obtain information on the significance of anisotropy on the analysis of the stress and strain state. The PDL is a membrane suspending the tooth in the bone with viscoelastic - non-linear, time and rate dependent - properties, as it is composed of collagen fibres surrounded by a gel-like matrix. While our hyperelastic PDL material model is useful to model the non-linear elastic characteristics of the PDL, it neither includes cyclic and velocity effects of viscoelasticity nor the different flow behaviour of the PDL gellike matrix and fibres. As the focus of our work is on tooth motion and single mastication strokes and not the cyclic dynamics of mastication, we assume that the time depending effects are unimportant for our research question. Furthermore our PDL material parameters are validated to in vitro experiments, which were all performed with the same, chewing-relevant loading rate of 50 mm/min in freshly excised jaws all loaded in a reproducible way. Our results are therefore valid for this rate, whereas other loading rates remain to be investigated. Additional simulations (Supplementary material) show that, even with softer or harder PDL properties - as to be expected due to viscoelastic effects during further loading cycles - we reach results that lead to similar conclusions (see Appendix C). To investigate tooth movement using an FE model and simulating mastication loads, we found that modelling of the neighbouring teeth is

Conflict of interest The authors declare no conflict of interests. Acknowledgements The authors thank Martin Burke and Matthieu Chappaz for the pigtooth compression tests. We gratefully acknowledge the North-German Supercomputing Alliance (HLRN) for its computing resources and Dr.Ing. Baumann for his advice. We thank the Extrusion Research and Development Center TU Berlin for use of HyperMesh. AN is indebted to the TU for financial support. PZ and CF are grateful for financial support of the German Research Foundation (DFG) through priority programme SPP 1420: 'Biomimetic Materials Research: Functionality by Hierarchical Structuring of Materials'. The BCRT is acknowledged for the use of the lab-μCT.

Appendix A. Supplementary experimental details Determination of loading velocity The chosen loading velocity is based on measurements of the chewing forces on a porcine third premolar (Bousdras et al., 2006) and a scaling

70

Journal of the mechanical behavior of biomedical materials 67 (2017) 61–73

A. Nikolaus et al.

law between the time for one chewing cycle and the body mass, given by Druzinsky (1993). We used a mean force of 375 N. This corresponds to a tooth displacement of l=0.22 mm (see Fig. 3). The time for one chewing cycle scales with body mass by TChew=125.89M 0.128 with TChew the chewing time (in ms) and M the mass (in g) of the animal (Druzinsky 1993). Assuming an average mass of 120 kg for an adult sow, this gives TChew = 563 ms. Taking an estimated loading time on the tooth of 30% (Menegaz et al., 2015) up to 50% of each chewing cycle, results in t=169–283 ms. This gives a range of loading speeds between 78 mm/min and 47 mm/min. Note, that these are only approximate estimates, and we therefore chose a loading speed of 50 mm/min for our in vitro and in silico experiments. Simulation settings The element types used were linear tetrahedra (C3D4) for all materials, except for the two 'uniform PDL' and the cementum meshes. The latter have different elements, because for the uniform thicknesses of 200 µm and 210 µm, a simple layer offset can be used to produce wedge elements (C3D6). Based on its Poisson's ratio of 0.49, the PDL was defined to be an 'almost incompressible material' (Andersen et al., 1991) for which hybrid elements were used. Typically, modelling of materials with a high Poisson's ratio ( > 0.45) leads to numerical instabilities. These are resolved in Abaqus by use of hybrid elements (Simulia, 2014) that give stable solutions for incompressible materials. The calculation of the state of deformation/rotation for each load level is either based on comparison with the undeformed initial geometry yielding a so-called 'linear geometry simulation', or it is based on the already deformed geometry under preceding, lower loads, resulting in so-called 'non-linear geometry simulation'. For simulations including hyperelastic materials such as the PDL in our models, it is necessary to use 'non-linear geometry simulations'. Hyperelastic material model for the PDL The PDL experiences large deformations, violating linear elastic theory which requires small deformations. Hyperelastic material models circumvent this limitation and were therefore used here to model the PDL. There are several forms of strain energy potentials available and implemented in the solver Abaqus 6.14, to model isotropic hyperelastic materials, e.g. Ogden, 1972; Arruda and Boyce, 1993; Marlow, 2003; Yeoh, 1993. The strain energy potential in the first Ogden form describes our PDL:

U=

2μ1 α12

( λ1α1 + λ 2α1 + λ 3α1 − 3) +

1 ( J el − 1)2 . D1

U is the strain energy per unit of reference volume and μ1, α1, and D1 are material parameters. λi = J −1/3 λ i are the deviatoric principal stretches, with λ i the principal stretches and J the total volume ratio. J el is the elastic volume ratio. Using hyperelastic material models already implemented in the FE solver makes the non-linear elastic modelling of the PDL easily reproducible. The PDL exhibits time and rate dependent behaviour and is therefore a viscoelastic material. Using a hyperelastic material model offers the advantages of static FE simulations, in terms of the amount of detail that can be calculated within a reasonable timeframe. We are aware of the limitations of such simulations for dynamic processes such as mastication. However, as the focus of this study is on the first biting stroke and the related cascade of events, and as exemplified by the validation experiments/results, the hyperelastic approach is an excellent, practical, first approximation experiment. Appendix B. Material properties of the tooth-bone complex The tooth-PDL-bone complex is composed of various materials spanning a wide range of stiffnesses. The hardest and stiffest material, enamel, is found on the outer chewing surfaces of teeth and has rather anisotropic mechanical properties, as it is composed of parallel mineral crystal bundles that form rods radiating outwards to the occlusal surface and that incorporate and are surrounded by organic sheaths (He and Swain, 2008). Some measurements of the Young's modulus of enamel are listed in the supplementary Table S2. Young's modulus varies with orientation and position in the enamel cap (Cuy et al., 2002). Dentine is a bone-like material, a nano-composite of collagen fibrils and nanometer-sized mineral particles. It is neither homogeneous nor isotropic: at the micrometre length scale, dentine contains elongated tubules that traverse the tissue, and at the millimetre length scale, dentine contains regions that are more mineralised than others, notably the bulk of the crown, which is encased by softer dentine externally towards the enamel cap and in the roots. Young's modulus values of dentine are listed in Table S2. The dentine roots are encased by a variable-thickness layer of cementum. It is a bone-like material and typical Young's moduli are also given in Table S2. Zaslansky et al. (2006) characterised a 200–300 µm soft zone below the DEJ in human teeth, the sub-DEJ dentine, and showed that this zone is structurally different from bulk dentine. The stiffness estimates for the sub-DEJ dentine fall in the range of 3–12 GPa, which is well below the stiffness reported for bulk crown dentine. The PDL is a fibrous connective tissue (Natali, 2003), which secures the tooth to the alveolar bone. It displays a non-linear stress-strain behaviour (Dorow et al., 2003) and various works report a large range of Young's modulus values spanning six orders of magnitude (Fill et al., 2011; Rees and Jacobsen, 1997). Several reviews exist on the material properties of the PDL (Fill et al., 2011, 2012; Minch, 2013), a survey of models and properties is listed in Table 1. Bone is an anisotropic material (Currey, 2002) comprising cortical bone and trabecular bone. Cortical bone is the hard outer shell-like region of long bones. Trabecular bone forms a porous, foam-like continuum (Natali, 2003). Measurements of the Young's modulus of single trabeculae are often performed by nanoindentation and therefore the specimens are usually dehydrated, which results in higher measured Young's moduli as compared to the wet state. In Table S2, the results of different measurements of the Young's moduli of trabecular bone are listed.

71

Journal of the mechanical behavior of biomedical materials 67 (2017) 61–73

A. Nikolaus et al.

Appendix C. Model validation The load displacement responses of different teeth vary somewhat. This is demonstrated by two of the force displacement curves of teeth in opposite sides of the lower jaw (Fig. 3). They are from premolars harvested from the same sow, tested on the same day under the same wet conditions, both with the same loading speed of 50 mm/min. The only difference is that one (marked by stars) is from the left and the other (marked by squares) from the right half of the jaw. Nevertheless, one (stars) shows a much stiffer reaction than the other one (squares). To verify that our FE results incorporate the full range of our lab measurements, we produced three fitting curves representing a 'fit hard', 'fit mean', and 'fit soft' case (Fig. 3). Three FE models with the 'variable PDL' were validated to the three cases by adjusting the hyperelastic properties of the PDL, listed in Table S1. With these parameters the load displacement curves of each FE model coincide nicely with the corresponding measurement, see Fig. S3. The simulation 'mean' is the model to which we refer in the main text as 'variable PDL' FE model. Fig. S2a and b show the von Mises bone stress as mean over several elements in the regions 'lateral' and 'furcation' for the three FE models. Especially for the first 80–120 N of load, the bone stress reaction is similar for the three models. Also note that both the decreases in stress gradient (i.e. rate of stress change) in the furcation area and the increases in the lateral bone stress gradient are reproduced by both the 'soft' and 'hard' model: in case of the 'soft' model the changes appear 'earlier', at lower loads, and for the 'hard' model 'later', at higher loads. Therefore, we conclude that changing the PDL stiffness to adapt the load displacement behaviour only has a minor influence on the resulting bone stress, and, furthermore, subtle differences in microanatomy between the different teeth are of minor importance as compared to the PDL geometry. The effective strain in the adjacent PDL is more influenced by its stiffness, see Fig. S2c and d. A softer PDL results in higher tooth displacements for the same load and therefore in higher PDL strains. This effect is reversed for the 'hard' PDL. Nevertheless, the strain in the 'lateral' region is always higher than in the 'furcation' region and the exponential decay form of the curves remains. Therefore, our conclusions for the 'mean' model are transferable to the 'soft' and 'hard' models. Consequently, we have reason to believe that our results and conclusions for the mean 'variable PDL' are valid within the range of load displacement responses measured for several pig teeth, shown in Fig. 3. Appendix D. Effective strain The effective strain εeff was calculated, in a manner similar to the von Mises stress, by:

εeff = 2/3 ε112 + ε22 2 +ε332 − ε11 ε22 − ε11 ε33 − ε22 ε33 +

3 (ε12 2 + ε132 +ε232 ) 4

where ε11, ε22 , ε33, ε12 , ε13 and ε23 represent the components of the strain tensor. Appendix E. Supporting information Supplementary data associated with this article can be found in the online version at http://dx.doi.org/10.1016/j. jmbbm.2016.11.020.

properties of the periodontal ligament. J Orofac. Orthop. 63 (2), 94–104. Dorow, C., Krstin, N., Sander, F.-G., 2003. Determination of the mechanical properties of the periodontal ligament in a uniaxial tensional experiment. J. Orofac. Orthop./ Fortschritte Kieferorthopädie 64 (2), 100–107. Dorow, C., Sander, F.-G., 2005. Development of a model for the simulation of orthodontic load on lower first premolars using the finite element method. J Orofac. Orthop. 66 (3), 208–218. Druzinsky, R.E., 1993. The time allometry of mammalian chewing movements: chewing frequency scales with body mass in mammals. J. Theor. Biol. 160 (4), 427–440. Erdemir, A., Guess, T.M., Halloran, J., Tadepalli, S.C., Morrison, T.M., 2012. Considerations for reporting finite element analysis studies in biomechanics. J. Biomech. 45 (4), 625–633. Farah, J.W., Craig, R.G., Sikarskie, D.L., 1973. Photoelastic and finite element stress analysis of a restored axisymmetric first molar. J. Biomech. 6 (5), 511–520. Favino, M., Gross, C., Drolshagen, M., Keilig, L., Deschner, J., Bourauel, C., Krause, R., 2013. Validation of a heterogeneous elastic-biphasic model for the numerical simulation of the PDL. Comput. Methods Biomech. Biomed. Eng. 16 (5), 544–553. Field, C., Ichim, I., Swain, M.V., Chan, E., Darendeliler, M.A., Li, W., Li, Q., 2009. Mechanical responses to orthodontic loading: a 3-dimensional finite element multitooth model. Am. J. Orthod. Dentofac. Orthop. 135 (2), 174–181. Fill, T.S., Carey, J.P., Toogood, R.W., Major, P.W., 2011. Experimentally determined mechanical properties of, and models for, the periodontal ligament: critical review of current literature. J. Dent. Biomech. 2 (1), 312980. Fill, T.S., Toogood, R.W., Major, P.W., Carey, J.P., 2012. Analytically determined mechanical properties of, and models for the periodontal ligament: critical review of literature. J. Biomech. 45 (1), 9–16. Fung, Y., 1993. Biomechanics: mechanical Properties of Living Tissues. Springer-Verlag, New York, USA. Gaengler, P., 2000. Evolution of tooth attachment in lower vertebrates to tetrapods. In: Teaford, Mark F., Smith, Moya Meredith, Ferguson, Mark W.J. (Eds.), Development, Function and Evolution of Teeth. Cambridge University Press, United Kingdom, 173–185. Gathercole, L.J., 1987. In-vitro mechanics of intrusive loading in porcine cheek teeth with intact and perforated root apices. Arch. Oral. Biol. 32 (4), 249–255. Gibbs, C.H., Mahan, P.E., Lundeen, H.C., Brehnan, K., Walsh, E.K., Sinkewiz, S.L., Ginsberg, S.B., 1981. Occlusal forces during chewing—influences of biting strength and food consistency. J. Prosthet. Dent. 46 (5), 561–567. Hildebrand, C., Fried, K., Tuisku, F., Johansson, C.S., 1995. Teeth and tooth nerves. Progress. Neurobiol. 45 (3), 165–222. He, L.H., Swain, M.V., 2008. Understanding the mechanical behaviour of human enamel

References Agrawal, K.R., Lucas, P.W., 2003. The mechanics of the first bite. Proc. R. Soc. Lond. B: Biol. Sci. 270 (1521), 1277–1282. Ali, A., Hosseini, M., Sahari, B.B., 2010. A review of constitutive models for rubber-like materials. Am. J. Eng. Appl. Sci. 3 (1), 232. Anderson, D.J., 1956. Measurement of stress in mastication, Part II. J. Dent. Res. 35 (5), 671–673. Andreasen, J.O., 2012. Pulp and periodontal tissue repair-regeneration or tissue metaplasia after dental trauma. A review. Dent. Traumatol. 28 (1), 19–24. Arruda, E.M., Boyce, M.C., 1993. A three-dimensional constitutive model for the large stretch behavior of rubber elastic materials. J. Mech. Phys. Solids 41 (2), 389–412. Ashman, R.B., Rho, J.Y., 1988. Elastic modulus of trabecular bone material. J. Biomech. 21 (3), 177–181. Bakke, M., 2006. Bite force and occlusion. Semin. Orthod. 12 (2), 120–126. http:// dx.doi.org/10.1053/j.sodo.2006.01.005. Benazzi, S., Kullmer, O., Grosse, I.R., Weber, G.W., 2012. Brief communication: comparing loading scenarios in lower first molar supporting bone structure using 3D finite element analysis. Am. J. Phys. Anthropol. 147 (1), 128–134. Berkovitz, B.K.B., 1990. The structure of the periodontal ligament: an update. Eur. J. Orthod. 12 (1), 51–76. Boldt, J., Knapp, W., Proff, P., Rottner, K., Richter, E.-J., 2012. Measurement of tooth and implant mobility under physiological loading conditions. Ann. Anat. - Anat. Anz. 194 (2), 185–189. Bousdras, V.A., Cunningham, J.L., Ferguson-Pell, M., Bamber, M.A., Sindet-Pedersen, S., Blunn, G., Goodship, A.E., 2006. A novel approach to bite force measurements in a porcine model in vivo. Int. J. Oral. Maxillofac. Surg. 35 (7), 663–667. Cattaneo, P.M., Dalstra, M., Melsen, B., 2005. The finite element method: a tool to study orthodontic tooth movement. J. Dent. Res. 84 (5), 428–433. Chen, J., 2009. Food oral processing—a review. Food Hydrocoll. 23 (1), 1–25. Chen, J., Li, W., Swain, M.V., Ali Darendeliler, M., Li, Q., 2014. A periodontal ligament driven remodeling algorithm for orthodontic tooth movement. J. Biomech. 47 (7), 1689–1695. Constantino, P.J., Lee, J.J., Gerbig, Y., Hartstone‐Rose, A., Talebi, M., Lawn, B.R., Lucas, P.W., 2012. The role of tooth enamel mechanical properties in primate dietary adaptation. Am. J. Phys. Anthropol. 148 (2), 171–177. Currey, J.D., 2002. Bones: Structure and Mechanics. Princeton University Press, New Jersey, USA. Dorow, C., Krstin, N., Sander, F.-G., 2002. Experiments to determine the material

72

Journal of the mechanical behavior of biomedical materials 67 (2017) 61–73

A. Nikolaus et al. from its structural and compositional characteristics. J. Mech. Behav. Biomed. Mater. 1 (1), 18–29. Ho, S.P., Yu, B., Yun, W., Marshall, G.W., Ryder, M.I., Marshall, S.J., 2009. Structure, chemical composition and mechanical properties of human and rat cementum and its interface with root dentin. Acta Biomater. 5 (2), 707–718. Ho, S.P., Kurylo, M.P., Grandfield, K., Hurng, J., Herber, R.-P., Ryder, M.I., Altoe, V., Aloni, S., Feng, J.Q., Webb, S., Marshall, G.W., Curtis, D., Andrews, J.C., Pianetta, P., 2013. The plastic nature of the human bone–periodontal ligament–tooth fibrous joint. Bone 57 (2), 455–467. Hohmann, A., Kober, C., Young, P., Dorow, C., Geiger, M., Boryor, A., Sander, F.M., Sander, C., Sander, F.G., 2011. Influence of different modeling strategies for the periodontal ligament on finite element simulation results. Am. J. Orthod. Dentofac. Orthop. 139 (6), 775–783. Huang, H., Tang, W., Yan, B., Wu, B., Cao, D., 2016. Mechanical responses of the periodontal ligament based on an exponential hyperelastic model: a combined experimental and finite element method. Comput. Methods Biomech. Biomed. Eng. 19 (2), 188–198. Jang, A.T., Merkle, A.P., Fahey, K.P., Gansky, S.A., Ho, S.P., 2015. Multiscale biomechanical responses of adapted bone–periodontal ligament–tooth fibrous joints. Bone 81, 196–207. Jeon, P.D., Turley, P.K., Moon, H.B., Ting, K., 1999. Analysis of stress in the periodontium of the maxillary first molar with a three-dimensional finite element model. Am. J. Orthod. Dentofac. Orthop. 115 (3), 267–274. Jones, M.L., Hickman, J., Middleton, J., Knox, J., Volp, C., 2001. A validated finite element method study of orthodontic tooth movement in the human subject. J. Orthod. 28 (1), 29–38. Keilig, L., Drolshagen, M., Tran, K.L., Hasan, I., Reimann, S., Deschner, J., Brinkmann, K.T., Krause, R., Favino, M., Bourauel, C., 2016. In vivo measurements and numerical analysis of the biomechanical characteristics of the human periodontal ligament. Ann. Anat. - Anat. Anz. 206, 80–88. Kinney, J.H., Marshall, S.J., Marshall, G.W., 2003. The mechanical properties of human dentin: a critical review and re-evaluation of the dental literature. Crit. Rev. Oral. Biol. Med. 14 (1), 13–29. Lombardo, L., Scuzzo, G., Arreghini, A., Gorgun, Ö., Ortan, Y.Ö., Siciliani, G., 2014. 3D FEM comparison of lingual and labial orthodontics in en masse retraction. Progress. Orthod. 15 (1), 1–12. Lucas, P.W., 2004. Dental Functional Morphology: How Teeth Work. Cambridge University Press. Marlow, R.S., 2003. A general first-invariant hyperelastic constitutive model. Constitutive Models for Rubber III, Eds. Busfield, J., Muhr, A., A.A.Balkema, Netherlands, 157–160. Martins, P., Natal Jorge, R.M., Ferreira, A.J., 2006. A comparative study of several material models for prediction of hyperelastic properties: application to silicone? Rubber and soft tissues. Strain 42 (3), 135–147. McCormack, S.W., Witzel, U., Watson, P.J., Fagan, M.J., Gröning, F., Agarwal, S., 2014. The biomechanical function of periodontal ligament fibres in orthodontic tooth movement. PLoS ONE 9 (7), e102387. Menegaz, R.A., Baier, D.B., Metzger, K.A., Herring, S.W., Brainerd, E.L., 2015. XROMM analysis of tooth occlusion and temporomandibular joint kinematics during feeding in juvenile miniature pigs. J. Exp. Biol. 218 (16), 2573–2584. Menicucci, G., Mossolov, A., Mozzati, M., Lorenzetti, M., Preti, G., 2002. Tooth–implant connection: some biomechanical aspects based on finite element analyses. Clin. Oral. Implant. Res. 13 (2), 334–341. Middleton, J., Jones, M., Wilson, A., 1996. The role of the periodontal ligament in bone modeling: the initial development of a time-dependent finite element model. Am. J Orthod. Dentofac. Orthop. 109 (2), 155–162. Minch, L., 2013. Material properties of periodontal ligaments. Postȩpy Hig. Med. Dośw. (Online) 67, 1261–1264. Mioche, L., Peyron, M.A., 1995. Bite force displayed during assessment of hardness in various texture contexts. Arch. Oral. Biol. 40 (5), 415–423. Mohl, N.D., Zarb, G.A., Carlson, G.E., Rugh, J.D., 1988. A Textbook of Occlusion. Quintessence Pub Co., Illinois. Nanci, A., 2008. Ten Cate's Oral Histology: Development, Structure, and Function. Elsevier Health Sciences, Pennsylvania, USA. Natali, A.N., Pavan, P.G., Scarpa, C., 2004. Numerical analysis of tooth mobility: formulation of a non-linear constitutive law for the periodontal ligament. Dent. Mater. 20 (7), 623–629. Natali, A.N., Carniel, E.L., Pavan, P.G., Bourauel, C., Ziegler, A., Keilig, L., 2007. Experimental–numerical analysis of minipig's multi-rooted teeth. J. Biomech. 40 (8), 1701–1708. Naveh, G.R., Shahar, R., Brumfeld, V., Weiner, S., 2012a. Tooth movements are guided by specific contact areas between the tooth root and the jaw bone: a dynamic 3D microCT study of the rat molar. J. Struct. Biol. 177 (2), 477–483. Naveh, G.R., Lev-Tov Chattah, N., Zaslansky, P., Shahar, R., Weiner, S., 2012b. Tooth– PDL–bone complex: response to compressive loads encountered during mastication – a review. Arch. Oral. Biol. 57 (12), 1575–1584. Ogden, R.W., 1972. Large deformation isotropic elasticity-on the correlation of theory and experiment for incompressible rubberlike solids. Proc. R. Soc. Lond. A. Math. Phys. Sci. 326 (1567), 565–584.

Ona, M., Wakabayashi, N., 2006. Influence of alveolar support on stress in periodontal structures. J. Dent. Res., 1087–1091. Parfitt, G.J., 1960. Measurement of the physiological mobility of individual teeth in an axial direction. J. Dent. Res. 39 (3), 608–618. Pietrzak, G., Curnier, A., Botsis, J., Scherrer, S., Wiskott, A., Belser, U., 2002. A nonlinear elastic model of the periodontal ligament and its numerical calibration for the study of tooth mobility. Comput. Methods Biomech. Biomed. Eng. 5 (2), 91–100. Poppe, M., Bourauel, C., Jäger, A., 2002. Determination of the elasticity parameters of the human periodontal ligament and the location of the center of resistance of singlerooted teeth. J. Orofac. Orthop./Fortschritte Kieferorthopädie 63 (5), 358–370. Provatidis, C.G., 2000. A comparative FEM-study of tooth mobility using isotropic and anisotropic models of the periodontal ligament. Med. Eng. Phys. 22 (5), 359–370. Qian, L., Todo, M., Morita, Y., Matsushita, Y., Koyano, K., 2009. Deformation analysis of the periodontium considering the viscoelasticity of the periodontal ligament. Dent. Mater. 25 (10), 1285–1292. Raadsheer, M.C., van Eijden, T., van Ginkel, F.C., Prahl-Andersen, B., 1999. Contribution of jaw muscle size and craniofacial morphology to human bite force magnitude. J. Dent. Res. 78 (1), 31–42. Rasband, W.S., 1997–2015. ImageJ. U. S. National Institutes of Health, Bethesda, Maryland, USA. 〈http://imagej.nih.gov/ij/〉. Rees, J.S., 2001. An investigation into the importance of the periodontal ligament and alveolar bone as supporting structures in finite element studies. J. Oral. Rehabil. 28 (5), 425–432. Ren, L.-M., Wang, W.-X., Takao, Y., Chen, Z.-X., 2010. Evaluation of the biomechanical characteristic of tooth supporting structure under occlusal load. J. Funtai Oyobi Fummatsu Yakin/Jpn. Soc. Powder Powder Metall. 57 (5), 298–305. Ren, Y., Maltha, J.C., Kuijpers-Jagtman, A.M., 2003. Optimum force magnitude for orthodontic tooth movement: a systematic literature review. Angle Orthod. 73 (1), 86–92. Salmon, P.L., Liu, X., 2014. MicroCT bone densitometry: context sensitivity, beam hardening correction and the effect of surrounding media. Open Access J. Sci. Technol. 2. Sanctuary, C.S., Wiskott, H.A., Justiz, J., Botsis, J., Belser, U.C., 2005. In vitro timedependent response of periodontal ligament to mechanical loading. J. Appl. Physiol. 99 (6), 2369–2378. Schindler, H.J., Stengel, E., Spiess, W.E.L., 1998. Feedback control during mastication of solid food textures—a clinical-experimental study. J. Prosthet. Dent. 80 (3), 330–336. Schrock, P., Lüpke, M., Seifert, H., Borchers, L., Staszyk, C., 2013. Finite element analysis of equine incisor teeth. Part 1: determination of the material parameters of the periodontal ligament. Veter-. J. 198 (3), 583–589. Su, M.-Z., Chang, H.-H., Chiang, Y.-C., Cheng, J.-H., Fuh, L.-J., Wang, C.-Y., Lin, C.-P., 2013. Modeling viscoelastic behavior of periodontal ligament with nonlinear finite element analysis. J. Dent. Sci. 8 (2), 121–128. Tanne, K., Sakuda, M., Burstone, C.J., 1987. Three-dimensional finite element analysis for stress in the periodontal tissue by orthodontic forces. Am. J Orthod. Dentofac. Orthop. 92 (6), 499–505. Toms, S.R., Eberhardt, A.W., 2003. A nonlinear finite element analysis of the periodontal ligament under orthodontic tooth loading. Am. J. Orthod. Dentofac. Orthop. 123 (6), 657–665. Trulsson, M., 2006. Sensory‐motor function of human periodontal mechanoreceptors. J. Oral. Rehabil. 33 (4), 262–273. Tuna, M., Sunbuloglu, E., Bozdag, E., 2014. Finite element simulation of the behavior of the periodontal ligament: a validated nonlinear contact model. J. Biomech. 47 (12), 2883–2890. Turner, C.H., Rho, J., Takano, Y., Tsui, T.Y., Pharr, G.M., 1999. The elastic properties of trabecular and cortical bone tissues are similar: results from two microscopic measurement techniques. J. Biomech. 32 (4), 437–441. Vollmer, D., Bourauel, C., Maier, K., Jäger, A., 1999. Determination of the centre of resistance in an upper human canine and idealized tooth model. Eur. J. Orthod. 21 (6), 633–648. Waltimo, A., Könönen, M., 1995. Maximal bite force and its association with signs and symptoms of craniomandibular disorders in young Finnish non-patients. Acta Odontol. Scand. 53 (4), 254–258. Wang, C.-Y., Su, M.-Z., Chang, H.-H., Chiang, Y.-C., Tao, S.-H., Cheng, J.-H., Fuh, L.-J., Lin, C.-P., 2012. Tension-compression viscoelastic behaviors of the periodontal ligament. J. Formos. Med. Assoc. 111 (9), 471–481. Weiss, J.A., Maker, B.N., Govindjee, S., 1996. Finite element implementation of incompressible, transversely isotropic hyperelasticity. Comput. Methods Appl. Mech. Eng. 135 (1), 107–128. Yoshida, N., Koga, Y., Peng, C.-L., Tanaka, E., Kobayashi, K., 2001. In vivo measurement of the elastic modulus of the human periodontal ligament. Med. Eng. Phys. 23 (8), 567–572. Zaslansky, P., Friesem, A.A., Weiner, S., 2006. Structure and mechanical properties of the soft zone separating bulk dentin and enamel in crowns of human teeth: insight into tooth function. J. Struct. Biol. 153 (2), 188–199. Zhang, W., Chen, H.Y., Kassab, G.S., 2007. A rate-insensitive linear viscoelastic model for soft tissues. Biomaterials 28 (24), 3579–3586.

73