Human feeding biomechanics: performance, variation, and functional constraints

Construction of solid models

The seven specimens chosen for analysis were CT-scanned at Penn State’s Center for Quantitative Imaging (pixel size = 0.16 mm) and the 2D digital image stacks were used to create seven solid meshes (Fig. 3) using Mimics v 14.0 (Materialise, Ann Arbor, MI, USA), following the methods outlined by Smith et al. (2015a) and Smith et al. (2015b). Mandibles corresponding to the seven crania (except for BERG and KSAN2, which lacked mandibles; see below) were also scanned so that they could be used to direct muscle force vectors in the loading simulations described below. The crania were solid-meshed at similar densities using tet4 elements (element count: GRGL = 2,118,350; BERG = 1,928,931; KSAN1 = 1,620,112; KSAN2 = 1,392,417; MALP = 1,323,093; TIGA = 2,059,433; WAFR = 1,831,053). Solid meshes were then imported as Nastran (NAS) files into Strand7 (Strand7 Pty Ltd, NSW, Sydney, Australia) FEA software.

We created two sets of human FEMs that differed in their assigned muscle force and bone properties. One set of human FEMs (“ALL-HUM” models) was assigned human properties for bone tissue and masticatory muscle force, whereas chimpanzee properties were applied to the second set (“CHIMPED” models). The ALL-HUM models provide the most realistic assessment of human cranial mechanics, in terms of the predicted strains and bite forces. These models also allow for a more thorough examination of intraspecific variation in humans. In contrast, the CHIMPED models permit direct comparisons between our humans FEMs and our previously analyzed FEMs of chimpanzees and fossil hominins (Smith et al., 2015a; Smith et al., 2015b). These comparisons focus on shape-related differences in mechanical performance that are free of the effects of differences in cranial size and bone material properties. Therefore, the comparisons between the CHIMPED human models and the chimpanzee data from Smith et al. (2015a) and Smith et al. (2015b) most directly address our mechanical hypothesis described above because the hypotheses relate specifically to the mechanical consequences of shape differences.

Material properties of tissues

Human cortical bone material properties assigned to the ALL-HUM models were collected from various locations across the craniofacial skeletons of two fresh-frozen human cadavers (female, aged 22; male, aged 42) by measuring their resistance to ultrasonic wave propagation (Ashman et al., 1984; Peterson & Dechow, 2002; Schwartz-Dabney & Dechow, 2002; Wang & Dechow, 2006; Wang, Strait & Dechow, 2006; see Supplemental Information). Previous studies show that freezing has only a very minimal effect on ultrasonic measurements and elasticity of cortical bone (Zioupos, Smith & Yuehuei, 2000). For each location sampled, the elastic (Young’s) modulus in the axis of maximum stiffness (E3) was averaged between the human donors and used to distribute spatially heterogeneous isotropic material properties throughout the seven human FEMs using a method (Davis et al., 2011) analogous to the diffusion of heat through a highly conductive material. To achieve this, values at each of the sampled locations, which ranged from 17.92–25.52 GPa (mean = 20.61 GPa, SD = 1.92), were converted to temperatures and distributed throughout the cortical volume of the FEM. The elastic modulus of cortical bone was then set to vary with temperature during the subsequent loading analysis, with any thermally-induced strains removed from the analysis. For Poisson’s ratios, models were each assigned the average of the sampled locations (v23 = 0.293). The same procedure was used to diffuse chimpanzee material properties to the CHIMPED model variants using data collected from a cadaveric female chimpanzee at 14 craniofacial regions (Smith et al., 2015a; Smith et al., 2015b). In both the ALL-HUM and CHIMPED sets of model variants, homogeneous isotropic properties were used to model both trabecular bone (E3 = 637 MPa; v23 = 0.28) and enamel (E3 = 80,000 MPa; v23 = 0.28), following Smith et al. (2015a) and Smith et al. (2015b).

Muscle forces and constraints

Jaw adductor muscle forces were applied to both sets of FEMs for the anterior temporalis, superficial masseter, deep masseter, and medial pterygoid under the assumption that the chewing muscles were acting at peak activity levels on both sides of the cranium. These loads allow an estimate of the maximum bite force produced by each individual. In the ALL-HUM variants, muscle forces were applied based on muscle physiological cross-sectional area (PCSA) data reported by van Eijden, Korfage & Brugman (1997), with forces corrected to account for pennation and differences in gape during fixation using formulae from Taylor & Vinyard (2013). Corrected PCSAs were then used to calculate forces in Newtons (N) such that each cm² of muscle was equivalent to 30 N (Murphy, 1998). These unscaled forces were applied to the “average” specimen (GRGL), while the six “extreme” variants were applied forces that were either scaled up or down based on differences in model size (Table 5), with size represented by model volume (i.e., the summed volume of all tet4 elements in mm³) to the two-thirds power. This muscle force scaling procedure removes the effects of differences in model size on stress, strain, and strain energy density (SED) from the mechanical results (Dumont, Grosse & Slater, 2009; Strait et al., 2010). The CHIMPED model variants were also assigned forces that were scaled dependent on their size using PCSA data from an adult female chimpanzee (Strait et al., 2009; Smith et al., 2015a; Smith et al., 2015b). However, rather than scaling the FEMs around the “average” specimen (GRGL), we scaled the forces applied to the CHIMPED models (see Table 5) from the baseline chimpanzee model used for scaling purposes (PC1+) in the analysis by Smith et al. (2015b), permitting size-free comparisons between humans and chimps. For both sets of muscle loadings, plate elements modeled as 3D membrane were “zipped” at their nodes to the surface faces of tet4 elements representing each muscle’s origin. The scaled muscle forces for each set of analyses were applied using Boneload (Grosse et al., 2007) to the normal surfaces of the plate elements as tractions directed toward their respective insertions on the mandible, with the mandible slightly depressed and the condyles translated onto the articular eminences (Dumont et al., 2010). Mandibles were only used here to direct these vectors. In the case of the BERG specimen, which was lacking its mandible, a scaled version of the GRGL mandible was used to define the orientation of muscle force vectors. Similarly, a scaled version of the KSAN1 mandible was used to replace the missing mandible in KSAN2.

For both sets of biting simulations, each of the seven FEMs was oriented such that one of three axes (i.e., X, Y, or Z) was parallel to the occlusal plane. Each model was constrained at a single node against translation in all axes at the working-side TMJ, while the balancing-side TMJ was constrained only in the superoinferior and anteroposterior directions (Strait et al., 2009; Smith et al., 2015a; Smith et al., 2015b), thus creating an axis of rotation around the TMJs. Models were subjected to simulations of left premolar (P³) and left molar (M²) biting by constraining a node in the center of occlusal surface in each tooth, respectively, in the superoinferior direction. These constraints generated strains in the craniofacial skeleton, as well as reaction forces at the TMJs and bite point, upon the application of muscle forces.

Analysis of model output parameters

Following Smith et al. (2015a) and Smith et al. (2015b), we displayed global strain patterns using strain maps. These maps are analogous to histograms in that they illustrate strain magnitudes at thousands of nodes simultaneously, but have the added advantage of preserving spatial information. In addition, we collected strain data generated by each FEM from surface elements at 14 locations across the craniofacial skeleton (Fig. 4). These locations correspond to those included in previous in vitro and in silico (e.g., FEA) studies on primate feeding biomechanics (e.g., Hylander, Johnson & Picq, 1991; Hylander & Johnson, 1997; Ross et al., 2011; Smith et al., 2015a; Smith et al., 2015b). At each location, we examined several strain metrics from each of the seven FEMs in order to understand patterns of deformation. These included maximum principal strain (tension), minimum principal strain (compression), maximum shear strain (maximum principal strain–minimum principal strain), von Mises strain (distortional strain or non-isometric strain), and strain energy density (SED, the strain energy stored at a given point). Additionally, strain mode, the absolute value of maximum principal strain divided by minimum principal strain, was recorded for each location. This measure indicates whether tension or compression is dominant at a given location.

Data on the reaction forces generated at constrained nodes (i.e., the bite point and two TMJs) were recorded in Newtons (N). Reaction forces at the P³ and M² were recorded relative to the occlusal plane, while reaction forces at the left and right TMJs were recorded and compared relative to a user-defined “triangle of support” Cartesian coordinate system, with one of three axes perpendicular to a reference plane defined by the “triangle of support” formed by the constrained nodes at the bite point and two articular eminences (Smith et al., 2015a; Smith et al., 2015b). The efficiency of bite force production at a given bite point in each model was also compared using the mechanical advantage (MA), a measure of masticatory muscle efficiency or leverage, calculated as the ratio of bite force output to muscle force input.

In the evaluation of our mechanical hypothesis, we first inspected data collected from the ALL-HUM models for large levels of intraspecific variation that could potentially invalidate the functional significance of our results. Strain magnitudes and SED at each of the 14 sampled locations were examined for large differences between individuals, in addition to a comparison of coefficients of variation (CVs) at specific locations. Differences in the spatial patterning of strain magnitudes between the ALL-HUM models were also compared using strain maps, in addition to variation in biting efficiency (i.e., MA). Lastly, we also calculated CVs for von Mises strain and MA in the CHIMPED model variants for direct comparison with the chimpanzee CVs reported by Smith et al. (2015b) using the Fligner-Killeen test for equal CVs.

To analyze relative mechanical performance in our human FEMs, we focused on comparisons between the CHIMPED humans and our previously analyzed FEMs of chimpanzee crania (Smith et al., 2015b). Specifically, we compared the magnitudes of von Mises strain, considered to be a key metric in assessing regional bone strength (Keyak & Rossi, 2000), at the 14 sampled locations, as well as differences in biting efficiency, between humans and chimps. We tested for significant differences between species using the Mann-Whitney U test.

Variation in strain magnitude and spatial patterning

Box-plots of strain and SED distributions recorded from the ALL-HUM models at the 14 sampled locations during premolar (P³) and molar (M²) biting are shown in Fig. 5 (see also Tables S1 and S2). Despite notable differences in craniofacial morphology between the models, comparisons of strain magnitudes reveal strong similarities. For P³ biting, the highest strain magnitudes were experienced at the working nasal margin (Location 12), although on average higher tensile strain magnitudes were generated at the working and balancing postorbital bars (Locations 4 and 5). During M² biting, the working zygomatic root (Location 8) was subjected to the highest strain magnitudes, except that tension was greatest at the balancing postorbital bar. During both bites, low strain magnitudes were generated along the supraorbital torus (Locations 1–3), the balancing zygomatic root (Location 9), balancing infraorbital (Location 11), and the zygomatic bodies (Locations 13 and 14). All FEMs of human crania were found to exhibit this general pattern.

Some regions of the face did exhibit large differences among individuals. In particular, the FEMs were found to differ in von Mises strain magnitude by as much as 210% at the nasal margin, which also has the highest CVs for all forms of strain during both P³ and M² biting (Table 7), with the exception of minimum principal strain at the working dorsal orbital (Location 2) and balancing infraorbital (Location 11) during P³ biting, SED at the working dorsal orbital (Location 2) during P³ biting, and the balancing zygomatic body (Location 14) for both bites.

Strain mode was nearly always compressive or tensile at a given location across the seven ALL-HUM models (Fig. 6), with a few exceptions. During premolar biting, only 3 locations varied with respect to strain mode (Locations 1, 10, 11), with only one FEM differing from the other models in each case. These three locations also differed in strain mode during molar biting, with Locations 1 and 10 exhibiting slightly higher levels of variation, in addition to variation in strain mode at Location 4.

By comparison with CHIMPED FEMs, humans were found to exhibit lower levels of shape-related variation in von Mises strain magnitude and lower CVs than chimpanzees at the 14 sampled locations (Table 8). However, results of the Fligner-Killeen tests reveal that only 3 of the 14 “gage sites” exhibit significant differences in CV values. Specifically, humans were found to exhibit a significantly lower CV at the zygomatic arches during both P³ and M² biting, and at the working infraorbital during P³ biting.

Variation in the spatial patterning of strain concentrations

Despite some large differences in strain magnitude, the spatial patterning of strain distributions was similar across the ALL-HUM models. The color maps during P³ biting (Fig. 7) reveal two predominant deformation regimes that are common across the seven FEMs: (1) superior displacement of the anterior maxilla in proximity to the loaded P³, which creates highly tensile and compressive (hence highly shearing) strains surrounding the root of the nasal margin, compression along the nasal margin, and compression at the working zygomatic root; and (2) frontal bending of the zygomae under the inferiorly directed pulling action of the masticatory muscles, which generates tension at the zygomatic body and near the zygomaticomaxillary junction, particularly at the working-side, and deforms the orbit such that it is tensed along an inferolaterally-oriented axis and compressed along a superolaterally-oriented axis.

The color maps of strain patterning during M² biting were also generally similar across the ALL-HUM models (Fig. 8). As expected, all models exhibited lower strain magnitudes in the lower maxillary region during molar biting compared to premolar biting, but higher concentrations of compressive strain at the working zygomatic root. Molar biting was also associated with the same type of frontal bending, zygomatic torsion, and orbital deformation that was observed for premolar biting, with relatively large concentrations of strain at the postorbital bars, orbital margins, and medial infraorbital.

In their study of chimpanzee biomechanical variation, Smith et al. (2015b) compared color maps of principal strain magnitudes in their six models with the scales normalized to an average of 10 landmarks (Locations 1–5, 8–12). They suggest that, by illuminating similarities and differences between individuals in the concentrations of relatively high and low strain concentrations through this normalization step, such “relative strain” maps strain may be particularly informative in comparative analyses of craniofacial function. When viewed in this manner (Fig. 9), the CHIMPED human models more clearly reveal a shared pattern of facial deformation that differs from that of chimpanzees under identical loading conditions, which was predominantly characterized by torsion of the zygoma and resulting orbital deformation under the inferiorly-directed masseteric muscle force.

Variation in bite force production and efficiency

The ALL-HUM models exhibit moderate differences in bite force production and efficiency (mechanical advantage, MA) at P³ and M² bite points (Table 9). With respect to bite force production, humans generated premolar bite forces that ranged from 333–507 N when loaded with scaled masticatory muscle forces. The MA range for premolar biting was 0.34–0.43 with all but one individual (WAFR) occupying a narrower range of 0.39–0.43. Molar bite forces ranged from 496–756 N. In terms of leverage, most FEMs exhibited molar MAs of 0.57–0.64, but with the WAFR model again being considerably less efficient (0.53).

When compared to the chimpanzee data in Smith et al. (2015a), the CHIMPED human models analyzed here were found to exhibit somewhat lower ranges of variation in biting MA. However, results of the Fligner-Killeen tests reveal no significant differences in CV values between the species at either the P³ (chimp = 8.67, human = 5.65; p = 0.18) or M² (chimp = 8.11, human = 6.67; p = 0.13) bite point.

Variation in reaction forces generated at the TMJs

During premolar biting, all seven of the ALL-HUM models generated strongly compressive reaction forces at both TMJs (see Table 9), similar to the results for chimpanzees (Smith et al., 2015b). However, unlike in chimpanzees, M² biting generated distractive (tensile) reaction forces at the working-side TMJ that would have “pulled” the mandibular condyle away from the articular eminence in five of the seven models. In order to remove distractive forces, these models required reductions in the muscle force applied to the balancing-side, which ranged from 5 to 15% (see Table 9). Interestingly, when loaded with chimpanzee muscle forces, all seven of the CHIMPED human models exhibit distractive forces in the working TMJ during M² biting, with larger muscle force reductions required to eliminate the distraction (see below).

Structural stiffness of the human craniofacial skeleton

Direct comparisons of shape-related mechanical performance between our human FEMs and our previously analyzed chimpanzee FEMs (Smith et al., 2015a; Smith et al., 2015b) were permitted by the CHIMPED models. These comparisons reveal that the human craniofacial skeleton is less stiff and experiences von Mises strains that are elevated relative to those experienced by chimpanzees when subjected to identical loading conditions (Fig. 10). Several of the sampled locations were found to experience significantly higher magnitudes in humans during both P³ and M² biting following the results of Holm-Bonferroni-corrected Mann-Whitney U tests (Table 10). These included the working nasal margin (Location 12), postorbital bars (Locations 4 and 5), working zygomatic root (Location 8), and the working dorsal orbital (Location 2). However, strains at the mid-zygomatic arches in humans were within the range observed for chimpanzees (which are extremely variable). Additionally, human zygomatic bodies were found to be structurally stiff, with significantly lower von Mises strain magnitudes than chimpanzees.

Human bite force production and mechanical efficiency

Analysis of our CHIMPED human FEMs reveals that human crania are capable of generating bite forces with higher mechanical efficiency than chimpanzees (Fig. 11). Pairwise comparisons using the Mann-Whitey U test demonstrate that these differences are significant at both P³ (U = 1.5, z = −2.73, exact p = 0.003) and M² (U = 1, z = −2.79, exact p = 0.002) bite points. However, unlike chimpanzees, all seven of the CHIMPED human models generated highly distractive (tensile) reaction forces at the working-side TMJ during molar biting. Therefore, molar biting in humans increases the risk of having the muscle resultant vector fall outside the triangle of support. To bring the joint back into compression, a reduction in balancing side muscle force of 15–30% was required (Table 11).

Craniofacial strength: Is the human face weak?

Our results suggest that the modern human craniofacial skeleton is structurally less strong, in terms of resistance to masticatory stress, than that of chimpanzees when subjected to identical loading conditions (i.e., same properties and constraints, muscle forces scaled to model size). In the CHIMPED variants of our human FEMs, most of the locations analyzed experienced von Mises strain magnitudes that were elevated relative to chimpanzees, in particular the working nasal margin, the postorbital bars, the working zygomatic root, and the working dorsal orbital region. Exceptions to this pattern include the zygomatic arches, where strains were bracketed by the range of values seen in chimp FEMs, and the prominence of the zygomatic body (i.e., the “cheek bone”), which is apparently strong in modern humans.

During unilateral P³ biting, the nasal margin of modern humans experienced von Mises strains that were on average more than 350% greater than chimpanzees. Similarly, previous investigations identify the “root” of the nasal margin to be an area of high stress and strain during masticatory loading in humans (Endo, 1965; Endo, 1966; Arbel, Hershkovitz & Gross, 2000; Szwedowski, Fialkov & Whyne, 2011; Maloul et al., 2012). This region is often described as a pillar-like structure (Benninghoff, 1925; Bluntschli, 1926), or section of a frame-like structure (Görke, 1904; Endo, 1965; Endo, 1966), that resists mainly compression during anterior tooth biting. The results of our analysis are in general agreement with these findings, except that tension at the nasal margin was also found to be high in magnitude, indicating intense bending and shearing of the lower maxillary region during anterior tooth biting (see Figs. 7 and 9).

In addition to the nasal margin, the postorbital bars of the human FEMs were also found to experience highly elevated von Mises strain magnitudes compared to chimpanzees. However, adjacent regions, including the zygoma/zygomatic body (“cheek bone”) region and zygomatic arch, were found to be similar in strength to the lower end of the chimpanzee range. Mechanical analyses of Paranthropus boisei and Australopithecus africanus (Smith et al., 2015a) show a similar pattern of relatively low strains in the zygomatic body. Smith et al. (2015a) suggest that the structural strength of the zygomatic body in australopiths could be adaptively significant, offering as one possibility that it serves to reduce strains in the nearby zygomatico-maxillary suture. In pigs, it has been demonstrated that unfused sutures can fail at relatively modest stress levels (e.g., Popowics & Herring, 2007), so some bony facial regions may serve to shield nearby sutures from masticatory stresses rather than bone itself (Wang et al., 2012). Among smaller-faced modern human crania, the zygomatico-maxillary suture may be especially prone to experiencing relatively large masticatory stresses. In our FEMs, the largest strains in this region of the mid-face were generated medial to the zygomatico-maxillary suture. The location of these elevated strain magnitudes corresponds roughly to the location of facial fractures experienced commonly during physical altercations (Ellis, El-Attar & Moos, 1985). Facial fractures are also common at the postorbital bar, as opposed to the zygomatic body or zygomatico-maxillary suture, when the zygomatic body is exposed to traumatic blows (Ellis, 2012; Pollock, 2012). Therefore, it is possible that the strength of the human zygomatic body, and perhaps the relative weakness of the postorbital bar, is related to diverting stress from sutures that might otherwise fail under relatively lower stress magnitudes.

Like the zygomatic body (“cheek bone”) region, humans were found to exhibit lower average von Mises strains and markedly lower peak strains than chimpanzees at the mid-zygomatic arch, although human values were bracketed by the range of chimp values. This potentially reflects differences in arch length. Specifically, the size of the temporalis muscle, which is correlated with the area of the infratemporal fossa (Weijs & Hillen, 1984), is significantly reduced in humans compared to that of chimpanzees (Taylor & Vinyard, 2013). Demes & Creel (1988) show that the area of the infratemporal fossa is nearly half that of chimpanzees, meaning that the total length of the zygomatic arch is also reduced. Bone strain analyses demonstrate that the arch is subjected to sagittal bending, as well as torsion along its long axis (e.g., Hylander, Johnson & Picq, 1991; Hylander & Johnson, 1997; Ross, 2001; Ross et al., 2011). Predictions based on beam theory therefore suggest that a decrease in the length of the arch will lessen these bending and torsional moments, whereas a reduction in the height and/or breadth of the arch will weaken it under bending and shear, respectively.

Functional interpretations based on the morphology of the zygomatic arch are complicated by the fact that the temporalis fascia has been hypothesized to stabilize it from the inferiorly-directed pulling action of the masseter muscle (Eisenberg & Brodie, 1965). Curtis et al. (2011) tested this hypothesis using FEA and found that models that do not include the temporalis fascia will overestimate strains in the arch and surrounding regions, including the postorbital bar and infraorbital. However, they also found that their models lacking a fascia generated strains more similar in magnitude to those collected during in vivo experiments (Hylander, Johnson & Picq, 1991; Hylander & Johnson, 1997; Ross, 2001; Ross et al., 2011). Similarly, previous FEA studies on primate crania that have not included a modeled fascia (e.g., Ross et al., 2005; Ross et al., 2011; Strait et al., 2005) also find broad agreement with in vivo data. Therefore, we did not feel that it was necessary to include this structure in our FEMs. Importantly, Curtis et al. (2011) did not actually model the temporalis fascia, rather, they applied external forces along the margin of the attachment of the fascia. This procedure assumes that the load transferred to bone by the fascia is evenly distributed around its perimeter. However, the fascia is subjected to load by the inferiorly directed force produced by those temporalis fibers that arise off of the deep surface of the fascia. This force should elevate tension in the fascia along its superior margin (i.e., where it arises off of the superior temporal line) while reducing tension along its inferior margin (i.e., along the arch). This factor may mitigate the role of the fascia in resisting the contraction of the masseter muscle.

Although the brow ridges are not thought to play an important role in masticatory stress resistance (e.g., Picq & Hylander, 1989; Hylander, Johnson & Picq, 1991; Ravosa, 1991a; Ravosa, 1991b; Ravosa et al., 2000) it is interesting to note that our human FEMs experienced higher von Mises strain magnitudes than chimpanzees at all three of the supraorbital sites examined, particularly during premolar biting. Between the human and chimpanzee samples, differences were found to be greatest at the working and balancing dorsal orbitals, not the dorsal interorbital, supporting the idea that the brow ridge cannot be modeled as a bent beam (Picq & Hylander, 1989; see also Chalk et al., 2011). The fact that the smaller brows of humans experienced elevated strain magnitudes during biting could be interpreted as meaning that large brow ridges are an adaptation to resist masticatory loads. However, a wealth of experimental data on humans and non-human primate species has shown (e.g., Hylander, Johnson & Picq, 1991; Ravosa et al., 2000; Szwedowski, Fialkov & Whyne, 2011; Ross et al., 2011; Maloul et al., 2012) that strains along the supraorbital margin are relatively low during biting and chewing, which is supported by the results presented here. Therefore, it is more reasonable to interpret differences in supraorbital morphology between humans and chimpanzees as being related to some non-dietary function, and that the resulting increases in brow ridge strain among humans are experienced as a secondary byproduct. For example, Moss & Young (1960) suggest that a large separation is formed posterior to the orbits when brain size is small, forming a supraorbital ridge. When brain size is large, the frontal bone is more steeply inclined posterior to the orbits, forming a vertical forehead rather than a large torus. A byproduct of this missing bar of bone above the orbits among modern humans could be that strain magnitudes are mildly elevated in that region.

Overall, our findings show that the human craniofacial skeleton is weaker than that of chimpanzees when subjected to feeding loads. These findings support the hypothesis that dietary changes involving a shift to softer and/or more processed foods along the modern human lineage has led to masticatory gracilization and reduced structural strength of the bony facial skeleton (e.g., Lieberman et al., 2004). However, in their biomechanical analysis, Wroe et al. (2010) recently found that although the human cranium is less robust, it experiences low peak strains and an even distribution of facial strain magnitudes compared to extant apes and fossil australopith species. Differences between our results and those of Wroe et al. (2010) could reflect differences in the way muscle loads were applied to the models in each analysis and/or the manner in which models were constrained. For example, we applied both normal and tangential tractions over entire muscle areas using Boneload (Grosse et al., 2007), whereas Wroe et al. (2010) loaded their models with muscles modeled as straight pre-tensioned beam elements. However, we conducted a sensitivity analysis to explore this possibility further (see Supplemental Information) and found that these differences in methodology only resulted in small differences in strain magnitude at most locations across the craniofacial skeleton.

Another possible explanation for the differences between our study and the study by Wroe et al. (2010) relates to the magnitudes of the applied muscle forces. Wroe et al. (2010) subjected their FEMs to three sets of simulated biting on various teeth. In their first simulation of the three, FEMs were assigned a set of species-specific muscle forces (or muscle force estimates) from the literature. In a second simulation, models were scaled to the surface area of their chimpanzee model and re-loaded using chimpanzee muscle forces. Lastly, in the third simulation, models were scaled to the surface area of their chimpanzee model and loaded with muscle loads required to generate an equivalent bite force. In this third simulation, the high biting leverage offered by the retracted human face meant that the muscle forces required to generate a bite compared to the other hominoids examined were relatively low. Therefore, Wroe et al. (2010) concluded that the human facial skeleton may in fact be well-adapted to resist masticatory stresses generated during high magnitude biting. Importantly, however, mean element von Mises stresses were found to be relatively high in their human FEM during the second simulation, where FEMs were scaled to the same surface area and loaded with equivalent muscle forces. This is the most similar of their three scaling procedures to the scaling performed here (scaling muscle forces to model volume^2/3), which we believe is the best means for removing the effects of size on comparisons of mechanical performance (e.g., Dumont, Grosse & Slater, 2009; Strait et al., 2010).

Bite force production and efficiency: are humans suited to produce large biting forces?

When analyzed using human bone and muscle properties (i.e., ALL-HUM models), our human FEMs produced bite forces of 333–507 N at the premolar (P³) and 496–756 N at the molar (M²). These results are similar to, but lower than, previous estimates of human bite force production using both 2D and 3D modeling techniques (e.g., Wroe et al., 2010; Eng et al., 2013). For example, using skeletal measurements and data on muscle cross-section, Eng et al. (2013) recently estimated that humans are capable of producing approximately 660–1106 N of M² bite force, while Wroe et al. (2010) estimated a maximum unilateral M² bite force of 1109–1317 N using FEA. However, our M² bite force results are bracketed by bite force transducer data collected from various western populations, which range from approximately 368 N (Sinn, de Assis & Throckmorton, 1996) to around 911 N (Waltimo, Nystram & Kananen, 1994), although Inuit males have been shown to produce an average of 1277 N in M² bite force (Waugh, 1937). Therefore, our results for bite force production lie within and do not exceed the known range of in vivo variation exhibited by recent human populations.

Because chimpanzees have absolutely and relatively larger jaw adductor muscles than humans (e.g., Taylor & Vinyard, 2013), it is no surprise that the chimp FEMs were capable of producing more forceful bites than our human FEMs when loaded with species-specific muscle forces (compare data in Table 9 to Smith et al. (2015b), Table 4). However, when loaded with muscle forces scaled to remove differences in size (as in the CHIMPED model variants), we found that humans are more efficient producers of bite forces, in terms of biting leverage, consistent with the findings of Wroe et al. (2010). Specifically, the MA for P³ biting in humans ranged 0.39–0.47, compared to 0.32–0.42 in chimpanzees (Smith et al., 2015b), with only two chimps overlapping the human range. Humans were found to exhibit even more elevated leverage during M² biting (0.60–0.71), with only one individual overlapping the chimpanzee range (0.49–0.61). When comparing these data using the Mann-Whitey U test, humans were found to be significantly more efficient at producing bite forces at both mesial and distal bite points. The CHIMPED humans were even found to exhibit a biting efficiency similar to that observed in australopiths (Smith et al., 2015a). In fact, P³ MA in P. boisei (0.40) and A. africanus (0.41) were near the lower end observed in humans. The FEM of A. africanus also generated M² bites with similar efficiency (0.62) to humans, whereas P. boisei produced more mechanically efficient (0.75) molar bites (Smith et al., 2015a).

Our data on bite force efficiency in humans support previous findings that have demonstrated the mechanical advantage of modern human bony facial architecture compared to both non-modern humans and non-human primate species (e.g., Spencer & Demes, 1993; O’Connor, Franciscus & Holton, 2005; Lieberman, 2008; Lieberman, 2011; Wroe et al., 2010; Eng et al., 2013). Using estimates of muscle leverage from 2D measurements (Lieberman, 2008; Lieberman, 2011), humans have been shown to achieve high biting leverage through a marked degree of facial retraction (orthognathism), which reorients the muscles of mastication relative to the tooth rows. As noted above, we found that our human FEMs produced bite forces with leverage ratios similar to those observed in A. africanus and P. boisei (Smith et al., 2015a). However, australopiths achieve high biting leverage through an anterior positioning of the chewing muscles relative to the tooth rows (Rak, 1983; Strait et al., 2009; Strait et al., 2010; Smith et al., 2015a). In humans, the midfacial region is “tucked” beneath the anterior cranial fossa (Lieberman, McBratney & Krovitz, 2002; Lieberman et al., 2004; Lieberman, 2008; Lieberman, 2011), which similarly places bite points in a position that offers higher mechanical advantage to the jaw adductors.

Although the human cranium can theoretically produce mechanically efficient bite forces, the production of unilateral molar (M²) bite force is limited by the risk of TMJ distraction, as predicted by the constrained lever model (Greaves, 1978; Spencer, 1998; Spencer, 1999). Specifically, we found that all seven of the CHIMPED human FEMs experienced a highly distractive (tensile) reaction force at the working-side joint during molar biting. These forces have the effect of “pulling” the mandibular condyle from the jaw joint, increasing the risk of joint dislocation (Spencer, 1998; Spencer, 1999). As noted in the introduction, the soft tissues of the mammalian jaw joint are well suited to resist compressive joint reaction forces, but are poorly configured to resist distractive joint forces that “pull” the mandibular condyle from the cranial base (Greaves, 1978; Spencer, 1998; Spencer, 1999). In contrast, only one of the six chimpanzee FEMs analyzed by Smith et al. (2015a) generated a tensile force at the working TMJ, and this reaction was only very weakly tensile (12.7 N). Similarly, Smith et al. (2015b) found that their FEMs of P. boisei and A. africanus lacked working-side distraction and were able to produce “stable” bites on both the premolars and molars, offering these species the ability to produce maximally forceful molar bites with limited risk of causing pain and/or damage to the TMJ capsule.

Interestingly, when loaded with human muscle forces (i.e., ALL-HUM), two of the human FEMs (TIGA and WAFR) were capable of maintaining weakly compressive reaction forces at both TMJs during molar biting. Additionally, balancing side force reductions required to eliminate distraction in the remaining models were proportionately less (5–15%) than when applying chimpanzee forces (15–30%). Comparisons of the muscle loads applied to the models and their force ratios in the ALL-HUM and CHIMPED models (see Tables 9 and 11) reveal that chimpanzees devote a higher proportion of muscle strength to anteriorly-positioned muscle compartments (superficial masseter and anterior temporalis) compared to more posteriorly-positioned ones (deep masseter and medial pterygoid). Therefore, it is tempting to suggest that changes in human jaw muscle force ratios may have coincided with the retraction of the lower face during human evolution in order to reduce the risk of TMJ distraction. Likewise, if the repositioning of cranial elements for reasons other than food processing (Lieberman, 2008; Zink & Lieberman, 2016) led to an increase in biting efficiency but the generation of working side joint distraction during molar biting, the overall reduction of chewing muscle size in Homo could also be viewed as a result of positive selection rather than relaxed selection so as to lessen these distractive forces.

Our findings that humans are limited in their ability to produce forceful unilateral molar bites are supported by data on bite force and muscle activity in humans. Spencer (1995) and Spencer (1998) tested some predictions of the constrained lever model and found that humans produced bite forces that increased as the bite point moved from the incisors to the first molar. However, moving from M¹ to M³, bite forces were found to decrease as a result of the decreasing balancing force muscle recruitment required to avoid joint distraction. Spencer (1995) also notes that most of the participants (8 of 10) in his analysis reported pain near the working-side TMJ when biting forcefully using the back molars. In addition to this study, Hylander (1977) suggests that specialized anterior tooth biting and increased masticatory muscle leverage may be related to the high incidence of third molar reduction and agenesis among modern Inuit due to the increased risk of distraction when biting on these teeth, although the results of our single pre-historic Arctic FEM (TIGA) provide no support for this hypothesis. Similarly, Spencer (2003) demonstrates that seed predating New World primates with adaptations for increased anterior bite force have relatively small third molar roots.

As discussed above, Wroe et al. (2010) analyzed human feeding biomechanics within a comparative context. One of the principal findings of their analysis, supported by the data presented here, is that humans are capable of generating bite forces with higher mechanical efficiency than chimpanzees. Wroe et al. (2010) use this as evidence to argue that human craniofacial evolution may have been influenced by selection for powerful biting behaviors. However, the results of this study showing the comparative weakness of the human cranium combined with the increased risk of jaw joint distraction during molar biting leads us to interpret the increased biting leverage exhibited by humans, which is particularly high among recent populations (Spencer & Demes, 1993; O’Connor, Franciscus & Holton, 2005), to be a byproduct of human facial orthognathism, which may be at least partly related to facial size reduction. Human facial flatness may also have been acquired through selection for some non-dietary function. For example, Lieberman (2008) and Lieberman (2011) suggests that the marked degree of facial retraction exhibited by modern human crania could be related to changes in brain size and cranial base flexion. However, Ross (2013) shows that basicranial flexion cannot produce significant facial retraction on its own. Alternatively, Holton et al. (2010) propose that dietary shifts leading to reduced facial strain magnitudes among early human species may have led to reduced facial growth and earlier fusion of the maxillary sutures, and thus smaller and more retracted facial skeletons.

Although the majority of the morphological and mechanical evidence is not consistent with the hypothesis that the human masticatory apparatus has experienced recent selection for high magnitude biting, the results of our analysis cannot reject the hypothesis that, in addition to changes in diet and tool use, increases in muscle force efficiency during human evolution could have led to relaxed selection for large chewing muscle size and reductions in facial size (Wroe et al., 2010) or that humans benefited from increased biting leverage when using submaximal forces by exerting less energy per bite. Our results for premolar biting leverage also do not conflict directly with the hypothesis that anterior tooth biting could have been selectively important in humans. However, the reduced size of the premolar teeth in humans increases the risk of tooth crown fracture (Constantino et al., 2010). Therefore, studies on premolar size and strength are not consistent with the hypothesis that humans are particularly well adapted for forcefully loading their anterior teeth, but such studies have yet to be conducted on incisors or canines, which are the more likely to be used during paramasticatory activities. For example, Hylander (1977) identifies features of the modern Inuit craniofacial skeleton that he argues to be adaptations for powerful biting behaviors using the incisors, although our single pre-historic Arctic FEM (TIGA) was not found to be exceptional in this regard. Additionally, Spencer & Ungar (2000) show that incisor bite force leverage varies in relation to the intensity of incisor tooth use among some Native American populations. Similarly, it is possible that differences in anterior tooth use among “archaic” members of the genus Homo are reflected in mechanical differences between the species. In particular, the Neanderthals (H. neanderthalensis) exhibit a number of derived characteristics hypothesized to be adaptations for forceful incisor biting (e.g., Brace, 1962; Smith, 1983; Trinkaus, 1983; Trinkaus, 1987; Rak, 1986; Demes, 1987). Notably, Spencer & Demes (1993) show that Neanderthals exhibit high incisor bite force leverage relative to H. heidelbergensis (but not modern H. sapiens). In order to maintain functional use of the posterior dentition (i.e., avoid TMJ distraction), Spencer & Demes (1993) further show that the molar tooth row in Neanderthals was anteriorly shifted, resulting in the characteristic retromolar gap.

Data on enamel thickness seemingly contrasts with the hypothesis that humans have experienced relaxed selection for powerful biting behaviors. Specifically, a number of studies find that recent human populations exhibit thick molar enamel (e.g., Martin, 1983; Martin, 1985; Olejniczak et al., 2008; Smith et al., 2006; Vogel et al., 2008), which has been interpreted as a primitive retention. However, notwithstanding disagreements over the significance of enamel thickness (Grine, 2005), Smith et al. (2012) recently show that “thick” molar enamel in humans is primarily the result of small coronal dentine areas. They found that enamel area in humans is reduced, but there was a disproportionately large reduction in dentine to enamel as human teeth were evolving smaller size, resulting in a relatively “thick” enamel cap. Thus, Smith et al. (2012) argue that the dichotomy between thick and thin enamel is an oversimplification.