Finite element modeling of occlusal variation in durophagous tooth systems

Teeth are intimately involved in food acquisition and processing; therefore, preventing tooth fracture, and the resultant loss of tooth function, is important. Work on fracture failure during food processing has focused primarily on bunodont teeth, such as the molars of humans or other primates. The dentistry community has been a strong motivator for this, with a focus on better understanding the mode of fracture in bunodont teeth under various loads. For example, radial cracks are the most prevalent mode of failure in simple bi-layered spheres when loaded with hard food items, but marginal cracks and semilunar chipping dominate when loading with softer foods (Qasim et al., 2005, 2007). Other work has focused on the role diet has played in human evolution. Primates that have to process large hard prey items prevent tooth fracture by thickening the enamel caps of their teeth (Lucas et al., 2008). In fact, enamel thickness is one of the aspects of tooth morphology that determines the mode of tooth failure in generalized bunodont teeth (Lawn and Lee, 2009). In addition to resisting tooth failure, hominid tooth morphology was also influenced by the ability to break food items (Berthaume et al., 2011). A mix of different cusp morphologies on the same tooth, as opposed to all sharp or all blunt for instance, takes these opposing evolutionary pressures into account to optimize bunodont tooth function: creating high stress in prey objects while minimizing stresses in the tooth enamel (Berthaume et al., 2013, 2014).

While the bulk of tooth fracture literature focuses on bunodont dentition, the function of other tooth morphologies has also been addressed; for example, carnivores, especially hypercarnivores, have modified teeth to pierce and cut soft tissues. Work on piercing teeth has focused on the length and bladed aspects of puncturing tooth morphology, approaching tooth failure as a functional trade-off with puncturing ability (Freeman and Lemen, 2007; Van Valkenburg and Ruff, 1987). At the same time, other work has looked at more detailed aspects of morphology, such as tooth composition, and found that the smoothed tips of conical puncturing teeth reduce the likelihood of chipped teeth, which allows for thinner enamel (Lawn et al., 2013). However, when puncturing soft tissues, stress concentrations are shifted to the margins of tooth crowns. These stress concentrations can lead to failure, but may be mitigated by the addition of a cingulum, a reinforced ledge of enamel that wraps around the base of many mammalian teeth (Anderson et al., 2011). Cutting teeth, like shark teeth or the carnivore carnassials, are often notched. These notches reduce the work needed to process malleable prey (Anderson and LaBarbera, 2008; Anderson, 2009; Anderson and Rayfield, 2012) though they can also concentrate stresses in the tooth, thus making tooth failure more likely (Whitenack et al., 2011).

Bunodont dentition shears and pulps food, caniniform teeth and carnassials pierce and cut through flesh, but crushing teeth have the straightforward job of transmitting the compressive force required to shatter a prey item. Across vertebrate taxa, crushing teeth are characterized with subjective, and ultimately uninformative, stereotypical terms, ‘flattened’ or ‘molariform’ (Mara et al., 2010; Mehta, 2009; Summers, 2000; Wilga and Motta, 2000), that serve to obscure the great diversity of tooth forms associated with durophagous diets. Teeth associated with hard-prey crushing diets can vary in occlusal convexity, and some even have cusps (Fig. 1). However, this diversity of tooth forms is not entirely expected – some tooth shapes are better able to crush hard prey than others, which should lead to convergence on this design (Crofts and Summers, 2014). Of course, the ability of teeth to crush prey items is only one selective pressure affecting tooth shape; the ability to withstand high forces without breaking must also be shaping teeth. The arms race between durophagous predators and their prey has had a profound impact on ecological structures throughout time (Vermeij, 1977), but the question still stands: what evolutionary pressures have shaped hard-prey crushing teeth through time?

The goals of this paper were threefold: to determine the effect of occlusal concavity/convexity on strain in a crushing tooth, to quantify the role of a centrally located stress concentrator on strain in the tooth, and to determine whether the ability of different tooth shapes to resist strain is sensitive to variations in prey size. To do this, I analyzed three series of canonical tooth shapes as finite element (FE) models, with each series varying by a single aspect of tooth morphology (Fig. 2). In this way I determined how changes in tooth morphology would affect strain distribution in teeth subjected to occlusal loads simulating either small or large hard-prey items.

For the concave–convex series of models under the small prey item loading regime, total strain energy increases as occlusal surfaces shift from concave to convex (Fig. 3A, Table 1). Maximum first principal strain follows a similar pattern, with concave shapes having low maximum first principal strain and convex shapes having high maximum first principal strain. Fig. 4A shows the changes in first principal strain distribution for the convex–concave model series under the small prey item loading. While the magnitude of strain changes between model morphologies, there is little change in overall strain distribution. For all models in this series, strain is primarily concentrated in the enameloid around the area being loaded, and dissipates through both the enameloid layer and the dentine body.

There is a marked difference between the patterns of total strain energy and maximum first principal strain in the concave morphologies when comparing the large prey item loading regime with the small prey item loading regime. Under the large prey item loading regime, as models go from concave to flat there is a very slight increase in total strain energy and maximum first principal strain; the difference between the most concave model and the flat model is much less than that observed under the small prey item load (Fig. 3B, Table 1). However, as expected given the similarity in loading areas, the overall pattern of increase in total strain energy and maximum first principal strain for the convex morphologies under the large prey item loading regime is similar to the increases seen under the small prey item load regime in these same shapes (Fig. 3A, Table 1). Fig. 5A shows the strain distributions for the convex–concave models under the large prey item loading regime. For convex shapes, strain remains concentrated around the area being loaded, as in the small prey item loading regime. As under the small prey item loading regime, strain is distributed in both the dentine and enameloid of these models. For the concave models, high strain is also concentrated under the site of applied load. However, in these shapes, load is applied to the raised peripheral edges of the occlusal surface, and strain is more concentrated in the dentine than in the overlying enameloid layer. Additionally, rings of high strain develop around the base of both the concave and flat models.

When models with central stress concentrators were loaded under the small prey item loading regime, total strain energy increased as the cusp height increased. In contrast, the maximum first principal strain decreased as the stress concentrator height increased (Fig. 3C, Table 1). In the models themselves, first principal strain is concentrated in the dentine around the area under load and is distributed through both the enameloid and dentine in the immediate area. This pattern remains constant for all stress concentrator heights (Fig. 4B). For the large prey item loading regime, both total strain energy and maximum first principal strain increase with the height of the stress concentrator (Fig. 3D, Table 1). Because the area being loaded changes slightly, the distribution of first principal strain also varies slightly, but is always concentrated around the area being loaded (Fig. 5B). Strain is distributed through both the dentine and enameloid, but for morphologies with lower cusps there is more strain in the dentine layer, but the magnitude remains lower than the strain in the enameloid layer.

For the models with stress concentrators with bases that vary from wide to narrow, there is a similar pattern for both the large and small prey item loads (Fig. 3E,F, Table 1). For both loads, there is little change in the magnitude of total strain energy or maximum first principal stain for most models. But there is a rapid increase in both metrics for the last three morphologies, those with the narrowest stress concentrators. For both loading regimes, first principal strain is concentrated around the area loaded, though this area varies slightly under the large prey item loading regime (Fig. 4C, Fig. 5C). For most tooth models, strain is distributed between the dentine and enamel layers, similar to the pattern seen in the other tooth model series. For the three morphologies with the narrowest stress concentrators, strain begins to form rings around the stress concentrator and is concentrated in the enameloid layer of the model under both loading regimes.

These results demonstrate the evolutionary pressure on hard-prey crushing teeth to resist failure, but this is not the only factor influencing tooth morphology. Comparing these results with previous work on the crushing ability of teeth, there appears to be a trade-off in performance. Testing the force required to break brittle, morphologically identical, 3D printed shells by physical models of the same three series of canonical tooth models tested here tells a different story of tooth optimization (Crofts and Summers, 2014). For large and small prey items, greater strain in convex teeth means a greater likelihood of crack formation than for flat or concave teeth, but convex teeth required less force to break a prey item than flat or concave ones (Fig. 3A,B). Similar relationships are seen in the cusped tooth models: taller cusped teeth are better able to break prey but have higher strain values (Fig. 3C,D), and the narrowest cusps show much higher strain values but are more effective at inducing prey failure (Fig. 3E,F). Given the trade-offs between tooth morphologies that can effectively fracture prey items and morphologies that will resist tooth failure, we might expect that intermediate, ‘ideal’ tooth morphologies would be the de facto tooth shape for hard-prey consumers. This, however, is not the case as there is a wide range of durophagous tooth morphologies.

One explanation for this variation in tooth morphology is the effect of prey shape, size and material properties on tooth failure. Compared with larger prey items, smaller prey generate loads that show a more pronounced change in both maximum strain energy and first principal strain as tooth models go from concave to flat. FE models of 4-cusped bunodont teeth loaded by brittle spheres of varying sizes showed a similar pattern of size-mediated variation in strain (Berthaume et al., 2014). Similarly, hemispherical tooth models worn flat were able to achieve a higher load to critical failure when loaded by flat surfaces than by rounded surfaces (Keown et al., 2012). The interaction between the flat indenter and the increasingly flat occlusal surface serves to spread the applied load (Ford et al., 2009). Given that the concave tooth morphologies of the present study behave similarly to the flat tooth morphology under the large prey loading regime, spreading the applied load will also increase the load to critical tooth failure. In this way, the gradation of concave to flat tooth morphologies represents a greater range of tooth morphospace open to animals that consume large prey items than may be available to those that consume smaller prey items. As concave tooth morphologies also spread the applied load, they can reduce the impact of tooth wear, and allow for increased tooth usage when processing larger prey items as well. While flatter surfaces spread loads and increase the load to critical failure, enamel that has been worn thin can flex and lead to subsurface cracks when loaded by a rounded indenter. Additionally, the discontinuity caused by the flat, worn surface can be prone to chipping, even with a flat indenter (Ford et al., 2009).

It should be noted that this is a discussion of the effects of loading and morphology on a single tooth. Loads can also be distributed across multiple teeth, reducing the applied force to any one tooth. Multiple teeth abutting each other can also allow for stresses and strains to be distributed from one tooth to another, redirect forces to teeth or portions of teeth less likely to fail, and improve overall tooth stability (Nobiling, 1977; Ramsay and Wilga, 2007). The present study is a first step in understanding the effects of loading on hard-prey crushing teeth, and further work is needed to understand how the interaction of multiple teeth would affect the patterns of strain distribution and magnitude.

The patterns of strain that we predict from the FE models are reflective of real-world failure regimes in physical models and in real teeth. When crushing large prey, the highest strain was concentrated in rings around the body of the tooth, for flat and concave teeth. A similar ring of concentrated strain can be seen around the base of the narrowest cusps for large and small prey, unlike other cusp morphologies where high strains are concentrated at the tip of the cusp. This pattern of strain distribution is similar to stress distributions seen in bi-layered epoxy models, which were shown to develop into ring cracks (Qasim et al., 2005). Additionally, a shift in strain to the margins of the model, as seen here in the concave and flat tooth morphologies, can lead to failure at the edge of the model (Qasim et al., 2007). This susceptibility to ring cracks and edge failure may be a reason why concave tooth morphologies are not as common as convex morphologies in nature. Finally, in cusped teeth, this strain pattern indicates a discontinuity between the body of these tooth models and the cusp, and is a site of likely failure.

Having teeth that resist failure is less important for animals that replace teeth frequently than for those that only rarely replace teeth, so the frequency of tooth replacement should be tied to tooth morphology. However, in durophagous animals it is also important to maintain functional tooth sites, which should reduce the rate of tooth replacement (Dalrymple, 1979). This strategy of reducing the rate of tooth replacement to prolong individual tooth function should constrain tooth morphology, favoring tooth morphologies that will be less likely to break. Some animals have developed an interesting solution to the competing pressures to replace teeth infrequently, maintaining function, and to increase the rate of tooth replacement, to remove damaged teeth, by developing a pattern of tooth replacement that maintains one functional crushing surface at a time while replacing another (Neenan et al., 2014). In these lineages, this may allow for more flexibility in replacement rate and morphology.

The results of FE modeling are only as good as the data put into the model. Our model relies on stiffness values for teeth from two species of shark that are not durophagous (Whitenack et al., 2010). There is evidence for a difference in the hardness of tooth materials for cutting versus tearing shark teeth (Enax et al., 2012), and fiber orientation in the enamel of crushing teeth does differ from other tooth types, which may allow these crushing teeth to be more resistant to compression (Preuschoft et al., 1974). Changing the material properties of the different layers of the tooth models could affect how strain is transmitted from one layer to the other and change the magnitude and patterns of maximum principal strain distribution in the models. We also only tested a single thickness of enamel across all tooth models. In mammalian bunodont teeth, variation in enamel thickness plays an important role in preventing tooth failure due to wear versus brittle fracture (Lawn and Lee, 2009). Neither the effects of changing both tooth material properties nor those of varying enamel thickness were addressed in this study and should be pursued in future work. The present study serves as a step towards understanding how shape affects strain distribution through hard-prey crushing teeth, and how this may influence the evolution of different specialized tooth morphologies.

For the first series, I varied the overall occlusal morphology from a deep concave surface, to one that was flat, to a highly convex occlusal surface. This was achieved by setting r=0.4, so that the added curve would cover the whole face of the model's occlusal surface, and varying h from −0.5 to 0.5 in increments of 0.1. This resulted in 11 morphologies (Fig. 2A): five of varying concavity, one with a flat occlusal surface and five with various degrees of convexity. In the second series, I added cusps of varying heights to the center of the occlusal surface of the flat tooth morphology. The width of the base of the cusp was constrained by setting r=0.01, and h ranged from 0 (no cusp) to 0.5, increasing in increments of 0.05. This generated another 11 shapes (Fig. 2B) with cusps of increasing height. The final series of models also had central cusps, but varied in the width of the base of the cusp. Cusp height was fixed (h=0.25), and the base of the cusp ran from r=0.4, which covers most of the occlusal surface, to the narrowest cusp where r=0.01. Beginning at r=0.4, the value of r decreased by 0.05 between each progressive model morphology until r=0.1. From that point, r decreased by 0.02 until r=0.02, which was halved (r=0.01) to create the narrowest cusp. This resulted in a series of 12 models (Fig. 2C).

Tooth models were constructed in the axisymmetric work-flow in COMSOL Multiphysics (v 4.3). To better reflect the structure of real teeth, models were constructed to have an outer layer of brittle, enameloid-like material, over a body of more ductile, dentine-like, material. This was accomplished by duplicating and scaling the initial equation (Fig. 6) in the course of model construction. I did not include a pulp cavity in the tooth models, following Anderson et al. (2011), reasoning that dentine is soft enough to allow for all enameloid deformation. There are few data in the literature on the material properties of non-mammalian tooth tissues, so I used averaged values taken from Whitenack et al. (2010). I set the Young's modulus of the enameloid layer to 70.745 GPa and the Young's modulus of the dentine body to 25.465 GPa, and used 0.3, an accepted estimate for most biological materials, as the Poisson’s ratio for both materials. I used COMSOL's built-in mesh feature to mesh the models, and set the mesh fineness such that there were multiple elements across the depth of the enameloid layer (Fig. 6). To mimic tooth attachment, I anchored the base of the model (Fig. 6), allowing no translation or rotation.

Each model was subjected to two loading regimes (Table 1), both centered over the middle of the tooth model's occlusal surface and running perpendicular to the base of the tooth (Fig. 7). The first loading regime was designed to mimic prey items smaller than the tooth, or with a much smaller radius of curvature. This was achieved by defining the area of the occlusal surface being loaded as a circle with a set radius (0.05; Fig. 7A). For models with narrow cusps, this method of loading is unrealistic, because the load spreads down the sides of the cusp (Fig. 7B). To correct for this, loads were constrained to the tip of the cusp down to a fixed height, approximately the same depth as the small load reached in the convex models (approximately −0.012345). The second loading regime mimics a prey item larger than the tooth or with a much larger radius of curvature. For this loading regime, loads were applied to the leading edge of the model to this same fixed height. For most models, those that were convex or possessing a cusp, this led only to a small change in the area being loaded, but not location of the load. For concave models, however, the large prey item loading regime loads only the peripheral edges of the occlusal surface, versus the middle of the occlusal surface. Similarly, while the small prey item loading regime only loads the middle of the occlusal surface for the flat model, the large prey item loading regime spreads the load over the entire occlusal surface.

The brittle failure of teeth directed both the loading of the models and the types of measurements we could take. Because teeth fail as brittle solids, loads were scaled to the volume of the specific model for each test (Dumont et al., 2009) to allow for comparison between morphologies. Because total load depends only on the volume of the model being tested, not the area loaded, total loads were the same for both loading regimes for each individual model.

Similarly, we used only first principal strain, not Von Mises stresses, as this is more appropriate for brittle solids (Dumont et al., 2009). Maximum principal strain describes the magnitude of strain at the most deformed node in the model, and while there are three principal strains, I measured only the first principal strain because it was consistently the largest tensile strain, and therefore most likely to be associated with failure. It should be noted that because maximum principal strains deal with only a single node, there is the potential for these data to be misleading if the node is anomalous. In addition to the maximum first principal strain, I also gathered data on the total strain energy in each model. This is a measurement of the amount of energy that goes into the deformation of each shape. As there is a threshold for any given material, past which it will break, teeth made of the same material should all have the same threshold. This means that teeth with a higher measured total strain energy will be more likely to pass that threshold and, thus, more likely to break. In the course of modeling, I generated heat maps demonstrating the distribution of strain in each model, allowing prediction of the most likely location of failure in each tooth model morphology.

Thanks to Adam Summers, Tom Daniel, Per Reinhall, Annika Erlebe, Betsy Dumont, Ian Grosse and Christian Sidor for their advice and patience.