Ontogeny of bite force in a validated biomechanical model of the American alligator

Developing computational models that accurately estimate biological and biomechanical functions remains difficult for researchers interested in complex systems such as the feeding apparatus of vertebrates. Even modeling superficially simple systems such as the akinetic, two-part skull of crocodylians poses challenges given the complex, 3D nature of the bones and the jaw muscles that actuate them. Nonetheless, accurately characterizing the biomechanical performance of the feeding apparatus is critical to our understanding of the evolution of the skull. Many modeling techniques have been advanced in recent years, such as finite element analysis (FEA), multi-body dynamic analysis and BoneLoad (Davis et al., 2010), but the predictions generated by these techniques often remain to be validated by comparisons with in vivo data. Validated modeling may be used to carry out in silico experimental investigations into the performance of the feeding apparatus in extant taxa. Moreover, the use of validated models is critical to studies of form and function in fossil taxa.

The feeding apparatus of extant crocodylians produces the highest measured bite forces among extant tetrapods (Erickson et al., 2003). Counter-intuitively, crocodylians do not display the typical morphology of other hard-biting tetrapods such as tegu lizards, hyenas or Tyrannosaurus, in which the skull is dorsally heightened, expanding the attachment area of temporal muscles and resisting dorsoventral bending of the rostrum (Molnar, 1998; Metzger and Herrel, 2005; Tseng and Stynder, 2011; Schaerlaeken et al., 2012). Instead, crocodylians evolved a dorsoventrally flattened skull, which is hypothesized to be an adaptation for aquatic ambush predation (Iordansky, 1973). This configuration is associated with a suite of biomechanical modifications of the feeding apparatus. The consequences of this flattened skull include a shift in temporal muscles toward more mediolateral orientations as well as increased resistance to mediolateral bending and axial torsion of the cranium (Busbey, 1995). The discordance between the morphology of the crocodylian skull and that of most animals with hard bites suggests that study of the former group can yield valuable insights into the anatomical determinants of feeding performance.

The morphology and performance of the crocodylian feeding apparatus have been the focus of numerous investigations over the past century. The early comparative and functional morphological studies of Iordansky (1964) and Langston (1973) set the stage for more quantitative biomechanical investigations of the crocodylian feeding apparatus such as the functional morphological investigations of Sinclair and Alexander (1987), Busbey (1989, 1995) and Cleuren et al. (1995). These early works used 2D static equilibrium lever analyses to estimate the magnitude and orientation of bite force and joint forces in the parasagittal plane. More recently, advances in computational power have facilitated in silico quantitative studies of feeding biomechanics. FEA provides tools to computationally investigate stress and strain (Daniel and McHenry, 2001; Metzger et al., 2005) and reaction forces (Porro et al., 2011, 2013) under a wide range of loading conditions. Researchers have applied this technique to investigate how feeding behavior loads the rostrum in comparative samples of extant crocodylians (McHenry et al., 2006; Pierce et al., 2008), comparing among archosaurs (Rayfield et al., 2007; Rayfield and Milner, 2008), and extinct crocodylomorphs such as thalattosuchians (Pierce et al., 2009). Using multi-body dynamics, Bates and Falkingham (2012) compared estimates of bite force of Tyrannosaurus with estimates from various tetrapods including a large individual of Alligator mississippiensis. These authors found that peak (impact) bite forces in A. mississippiensis matched previously reported maximum in vivo data.

Numerous in vivo studies have provided crucial data and advanced our understanding of the performance of the feeding apparatus in crocodylians. Researchers using electromyography (EMG) and X-ray cineradiography found that all major groups of adductor muscles were active during crushing bites but showed that the pterygoideus muscles were inactive during holding bites (Van Drongelen and Dullemeijer, 1982; Busbey, 1989; Cleuren and de Vree, 1992; Cleuren et al., 1995). Metzger et al. (2005) measured in vivo strain in the cranium during biting and found that these data broadly agreed with computational predictions. Most recently, Porro et al. (2013) found that in vivo strain magnitudes in the mandibles generally surpass predictions from FEA.

In vivo bite force recordings using force transducers are an invaluable source of data for further understanding the ontogenetic scaling and comparative biomechanics of the crocodylian feeding apparatus. Erickson and colleagues conducted a series of studies that measured bite force in a growth series of American alligator (Erickson et al., 2003; Gignac and Erickson, 2015) and extant crocodylian species (Erickson et al., 2012, 2014). Across a size range of 31.5 to 405.5 cm total length, these studies found a positive allometric relationship (b) between maximum bite force and a variety of body size proxies, such as total length (b=2.62; isometry=2), snout–vent length (b=2.59; isometry=2), mass (b=0.79; isometry=0.667) and head length (b=2.75; isometry=2). Erickson et al. (2003) hypothesized that this positively allometric increase in bite force may be responsible for intraspecific niche partitioning in A. mississippiensis. Furthermore, with the exception of Gavialis gangeticus, all extant crocodylians of a given size have equivalent bite forces. These studies have provided a solid in vivo basis against which in silico predictions may be tested.

Gignac and Erickson (2016) compared estimates of A. mississippiensis bite force derived from static bite force modeling with their previously published in vivo findings. The modeling techniques employed by the authors reliably calculate bite force in A. mississippiensis. These authors used photographs of dissections to precisely measure muscular physiological cross-sectional area (PCSA). Although this method is successful in calculating bite force in A. mississippiensis, it relies on access to cadaveric specimens, and, as such, is not applicable to the fossil record.

Although most previous modeling studies investigated the effects of muscle force on cranial forces in the skulls of crocodylians, few relied on anatomically detailed muscular attachment geometry from an ontogenetic series of animals. Of course, some studies of extant crocodylian biomechanics have benefited from dissection of cadaveric specimens to inform muscle force reconstructions (Porro et al., 2011; Gignac and Erickson, 2016), but these methods are not immediately applicable to questions of cranial function in extinct crocodyliforms. As such, accurate methods for muscle modeling in fossil taxa should be pursued to better understand the evolution of this system. Thus, despite the long history of functional analyses of crocodylian skulls, our understanding of the biomechanical environment of the feeding apparatus requires improvement.

Here, we tested the hypothesis that bite force predictions using digital models of skull morphology, myology and 3D computational modeling will be consistent with in vivo data. Specifically, we hypothesized that our model will generate bite forces consistent with the positively allometric relationship previously reported in in vivo bite force data from an ontogenetic series of A. mississippiensis (Erickson et al., 2003). Bite force was calculated using 3D lever mechanics and FEA. Although finite element software was used to create models, map muscle attachment sites and interface with BoneLoad in order to distribute muscle forces, we used FEA to calculate bite force; we did not investigate stress or strain distributions or deformation in the present study. Barring biologically unrealistic deformation, forces calculated with FEA should converge with results obtained with 3D lever mechanics.

Five frozen, unpreserved specimens of A. mississippiensis Daudin 1802 were obtained from Rockefeller Wildlife Refuge (Grand Chenier, LA, USA), and a single dry skull was obtained from a private collector. Skull lengths ranged from 4.8 to 45.4 cm (Fig. 1, Table 1), corresponding to total body lengths of approximately 38–326 cm (Woodward et al., 1995). The smallest individual was CT scanned at the University of Missouri Biomolecular Imaging Center (Siemens Inveon MicroCT, Siemens Medical Solutions USA Inc., Malvern, PA, USA); the largest individual was scanned at the University of Missouri School of Medicine Department of Radiology (Siemens Somatom Definition Scanner, Siemens Medical Solutions USA Inc.). All other animals were scanned at the University of Missouri School of Veterinary Medicine (GE LightSpeed VCT CT scanner, GE Medical, Milwaukee, WI, USA).

Stacked images were manually segmented in Avizo 9 (Visualization Sciences Group, SAS, Merignac, France; Fig. 2A), and three-dimensional models of skeletal anatomy were created (Fig. 2B). Using Geomagic Studio 13 (Geomagic Inc., Research Triangle Park, NC, USA), models were aligned to world axes and cleaned to remove features that unnecessarily increased computational time. Meshes were constructed with four-noded tetrahedra in Strand7 (G1D Computing Pty Ltd, Sydney, Australia; Fig. 2C). Four-noded tetrahedral bricks were used to construct finite element models. Finally, muscle attachments were ‘mapped’ onto the models using Strand7 following the methods of Grosse et al. (2007; Fig. 2D). All models had at least 500,000 elements to ensure models behaved convergently (McCurry et al., 2015; Table S1). The dimensions of the models are: x is positive in the left lateral direction, y is positive in the dorsal direction and z is positive in the rostral direction. All models were tested at 5 deg of gape. Finite element models are available online (https://osf.io/jmpck/). We modeled unilateral, left-sided, static, crushing bites at the most caudally located maxillary tooth, where bite force is theoretically highest and where Erickson et al. (2003) measured bite force. Muscles were assumed to contract maximally, as reported by previous EMG studies (Busbey, 1989; Cleuren et al., 1995).

All muscle terminology follows Holliday and Witmer (2007). In the present study, muscles modeled were M. adductor mandibulae externus superficialis (mAMES), M. adductor mandibulae externus medialis (mAMEM), M. adductor mandibulae externus profundus (mAMEP), M. adductor mandibulae posterior (mAMP), M. pseudotemporalis superficialis (mPSTs), M. pseudotemporalis profundus (mPSTp), M. pterygoideus dorsalis (mPTd), M. pterygoideus ventralis (mPTv) and M. depressor mandibulae (mDM). In extant crocodylians, most cranial muscles have substantial mediolateral components (mAMEM, mAMEP, mPSTp), rostrocaudal components (mPTd), or both (mAMES, mAMP, mPSTs, mPTv, mDM; Fig. 3).

In this study, muscle force was distributed over the surface area of attachment of the origin rather than modeled as a single vector. Each face of a tetrahedral element belonging to a muscle origin bore a portion of the total force, directed at the centroid of the muscle insertion. The computational toolkit BoneLoad version 7 (Davis et al., 2010) was used for distributing muscle forces across attachment sites and to calculate moments about axes. BoneLoad was originally used in modeling bite forces in phyllostomid bats, and its predictions are well supported by in vivo measurements (Davis et al., 2010; Santana et al., 2010). BoneLoad uses the geometry of muscle attachments (Fig. 3) and magnitude of muscle forces to automate the calculation of moments about an axis of rotation and distribute these muscular forces (Fig. 2E). These muscular force distributions were then used in both 3D lever analysis (Fig. 2F) and FEA (Fig. 2G).

Bite force was also calculated using FEA. In FEA, the user constrains the degrees of freedom at specified nodes. The FEA software then calculates the force required at each constraint to maintain equilibrium. Methods described by Strait et al. (2005) were followed to assign constraints. A single node at the tip of the caudal tooth and a single node in the middle of the articular surface of each quadrate bone were constrained in all three translational and all three rotational degrees of freedom. Material properties of alligator mandibular cortical bone were assigned to all elements of the FEM following Zapata et al. (2010). The scope of this project prohibited the inclusion of cranial sutures in the models; further, the material properties of sutures in A. mississippiensis (or indeed, the material properties of cranial sutures in any reptile) are unknown. Porro et al. (2011) found that although including sutures in finite element models affects stress and strain distributions in the alligator mandible, reaction forces including bite force were not dramatically affected. Although we were not investigating stress and strain distributions in the skull in the present study, we would expect artificial concentrations of stress and strain near bite points and muscle attachments (Curtis et al., 2013). For these reasons, cranial sutures were not included here. Because FEA provides force orientations, the component of force in each dimension is reported in addition to overall magnitudes (mediolateral, F_x; dorsoventral, F_y; rostrocaudal, F_z; total, F_sum).

Biomechanical models are useful only insofar as they produce consistent results that are at least broadly comparable with in vivo data. To validate this method, bite forces calculated using both 3D LM and FEA were compared with in vivo bite force data reported by Erickson et al. (2003). Erickson and colleagues measured maximum in vivo bite force in an ontogenetic series of A. mississippiensis using force transducers. To assess how bite forces calculated in this study correspond to in vivo data, ordinary least squares regression was conducted on bite force calculated with both FEA and 3D LM against skull length using R (http://www.R-project.org/). Because Erickson et al. (2003) did not report skull lengths, we used the published relationship of skull length against snout–vent length (Woodward et al., 1995) to calculate skull lengths for individuals in the study of Erickson et al. (2003). Ordinary least squares regression is justified over standardized major axis regression, because although skull lengths were presumably not measured without error, the error is likely to be low and the ratio of this error to the error in either directly measured or estimated bite force is also low. To compare slopes of regressions of log-transformed bite force on log-transformed skull length between Erickson et al.’s (2003) data and our results, we used a linear model with data source (in vivo, FEA and 3D LM), skull length and the interaction term. This analysis of covariance model allows each source of bite force data to have a separate slope while allowing comparison between slopes.

In the case of significantly different slopes between modeled and measured bite force data, we used the Johnson–Neyman technique to determine the region in which there is no significant difference in slope (Johnson and Neyman, 1936; White, 2003). The Johnson–Neyman technique compares two regressions and provides upper and lower values of the independent variable between which slopes do not significantly differ. We used the Johnson–Neyman technique to compare in vivo data with both FEA and 3D LM forces. All code for analysis is available online (https://osf.io/jmpck/).

Bite force typically refers only to the compressive (i.e. dorsoventral) components of force acting on a food item, as this is the component of force that will do work on food. FEA calculates forces acting in all dimensions, but lever mechanics calculates only forces acting perpendicular to the plane containing the axis of rotation and the point of force application. Similarly, the bite force transducers used by Erickson et al. (2003) measured only forces acting perpendicular to the long axis of the cantilever. In both cases, forces are dorsoventrally oriented. The results of FEA presented below therefore only include the magnitude of dorsoventral force, and statistical analyses were performed on only the dorsoventral component of bite force, which was typically ∼90% of total bite force.

Regression of muscle length, volume and force against skull length showed that these parameters scaled isometrically in all muscles with the exception of mAMEP and mPSTp, in which muscle volume and force scaled with slight negative allometry (slope estimates of 1.59 and 1.52, respectively; see Table 2 and Table S1). Fig. 5 shows reconstructed muscle force and the proportion each muscle contributes to total muscle force. Note that mAMP and mPTd together account for approximately two-thirds of muscle force in our model. However, our methods likely underestimate the force of mPTv; see Discussion.

Bite force estimates ranged from 49.3 N in the smallest individual (both methods) to 3460 N in the largest individual (3D LM). Bite force estimation with 3D LM and FEA yielded nearly identical results (Table 3). Magnitudes of total bite forces calculated with FEA and 3D LM differed by <6%. Whereas bite forces calculated in intermediately sized individuals matched in vivo data well, bite force in larger and smaller individuals diverged from in vivo data, with lower force estimates in larger individuals and higher force estimates in smaller individuals relative to in vivo data. The percentage contribution of a muscle to bite force is not necessarily the same as its percentage contribution to total muscle force (Fig. 5; Tables 4 and 5), because muscles vary in attachment site geometry in the crocodylian adductor chamber. For comparisons with in vivo bite force data, we only considered the dorsoventral component of bite force. However, the conditions of static equilibrium demand that forces be balanced in all three dimensions. Therefore, our FEA also calculated rostrocaudal and mediolateral components of bite force (Table 6). Bite points experienced medially and rostrally acting forces in addition to dorsoventral force.

Regression of bite force against skull length showed that both 3D LM- and FEA-calculated bite forces do not significantly differ from isometry (3D LM 95% confidence interval 1.64–2.06; FEA 95% confidence interval 1.63–2.05). By contrast, in vivo data from Erickson et al. (2003) showed positive allometry (95% confidence interval 2.51–2.61). Application of the Johnson–Neyman technique on both sources of calculated bite force data against in vivo results of Erickson et al. (2003) revealed that both samples had a region of non-significant difference of slopes. For bite force calculated with 3D LM, median values of lower and upper skull length in the region where slopes were not significantly different were 9.8 to 17.9 cm, respectively. Between these sizes, 3D LM predicts a slope that does not significantly differ from in vivo data. For bite force calculated with FEA, median values of skull length were 9.5–17.4 cm (Fig. 6). Between these sizes, FEA predicts a slope that does not significantly differ from in vivo data.

Biomechanical modeling offers researchers powerful tools with which to test hypotheses of feeding performance of extant and extinct taxa. In vivo bite force data of wild crocodylians are challenging to obtain, and in vivo measurements are obviously not possible in extinct taxa, making computational modeling necessary to explore patterns of form and function in the group. Accurate computational methods can model biting under varying conditions of tooth contact, gape and muscle recruitment, and thus modern computational methods are an excellent option for investigating the relationship between morphology, bite force and resulting cranial forces.

The two biomechanical modeling techniques used in this paper produce results consistent with each other. Like other validated models (Davis et al., 2010; Santana et al., 2010), our alligator models faithfully perform within an acceptable range of in vivo bite force. The good performance of the model further demonstrates the utility of high-fidelity muscle inputs. However, bite force in most individuals differed somewhat from predicted in vivo bite forces from similarly sized animals. These differences between in silico and in vivo techniques shed light on the challenges of modeling complex feeding function and cranial biomechanics. Causes for this mismatch may be divided into (1) differences between modeled bites and in vivo bites and (2) submaximal model performance.

First, the mechanism of bite force production in our models differs from the direct measurements of Erickson et al. (2003). Whereas we modeled static, crushing bites, Erickson et al. (2003) included some unknown amount of momentum: the maximum force during ‘aggressive, snapping’ bites. Therefore, the peak force reported by these authors likely included some degree of impact force resulting from rapid mandibular deceleration. Daniel and McHenry (2001) suggested that ‘dynamic loading due to rapid deceleration’ likely plays a role in maximal forces experienced by the skull. Because the present study modeled static, crushing bites, which have no momentum contribution from impact forces, estimated maximum bite forces are presumably below peak forces experienced by the cranium.

Second, the modeling techniques employed here underestimate the force of mPTv because models do not adequately capture aponeurotic muscle attachments. Crocodylians have a complex tendinous skeleton among the adductor mandibulae and pterygoideus muscle bellies (Iordansky, 1964, 2000; Schumacher, 1973; Busbey, 1989) that is challenging to model. Traditional tomography techniques fail to image these tendons with enough reliability to create digital models. Because the tendinous attachments of mPTv are missing (‘U-tendon’ of Iordansky, 1964; ‘pterygoideus-tendon aponeuroses’ of Schumacher, 1973; ‘posterior pterygoid tendon’ of Busbey, 1989), we underestimate the total surface area of muscle attachment and thus the total estimated force. mPTv accounts for approximately one-third of A. mississippiensis jaw muscle mass (Busbey, 1989; Cleuren et al., 1995) but only about 3% of the total muscle force in our model (Table 2). Moreover, the methods employed in this study place force vectors directly between attachment sites. mPTv originates on the edge of the pterygoid flange, courses caudoventrolaterally, then passes around the ventral border of the mandible to attach on the lateral surface of the angular bone. Force vectors oriented from muscle origin to insertion will therefore pass through the body of the mandible, resulting in erroneously collinear forces between the cranial and mandibular attachments of this muscle. Although BoneLoad accounts for muscle fibers that pull on a curved surface of bone (‘muscle wrapping’ sensu Grosse et al., 2007), it does not account for subsequent changes to the course of a muscle fibers, such as occur when a muscle passes around a bone or other structure. We suggest that the former phenomenon (muscle forces distributed across a curved attachment surface) be called ‘muscle traction’ and the latter anatomical phenomenon be called ‘muscle wrapping’. To better reconstruct muscle function accurately, future neontological studies could use contrast-enhanced CT imaging of soft tissues (Gignac et al., 2016) to facilitate the inclusion of the tendinous skeleton and associated musculature as well as properly orient force vectors in wrapping (Moazen et al., 2008; Gröning et al., 2013). However, tendons are rarely preserved in the fossil record, leaving studies of muscle function in extinct taxa to rely on inferential methods. As this method was developed to apply to fossil crocodylomorphs and other vertebrates, tendinous attachments were not included in muscle attachments.

Because mPTv makes up a sizeable proportion of jaw muscle mass (Busbey, 1989; Cleuren et al., 1995), errors in modeling this muscle may be particularly deleterious to model fidelity. To explore the effects of altering the magnitude and orientation of mPTv force on model performance, we used our largest specimen to calculate bite force under three additional scenarios. To orient the force vector of mPTv more correctly, we used DiceCT-based scans to determine the angles of insertion of mPTv on the lateral surface of the articular bone. We then oriented muscle force along this adjusted vector rather than towards the muscle's cranial attachment site (Fig. 7). To account for the underestimated PCSA of mPTv, we scaled muscle force magnitude by the ratio of muscle force calculated for mPTv by Gignac and Erickson (2016) over that of our own mPTv for equivalently sized animals (approximately 5.25 times). Gignac and Erickson (2016) dissected cadaveric specimens to calculate PCSA. This is inapplicable to fossil taxa, but presumably yields more accurate PCSA data. We then combined these orientation and magnitude variables. Bite force for the 45.4 cm-long specimen was originally 3200 N. Reorienting mPTv force raised bite force to 3377 N, while scaling its force by 5.25 raised bite force to 4160 N. Combined, these effects resulted in a bite force of 4515 N. Both models in which mPTv force was scaled up resulted in substantially higher bite force. Of course, both sources of correctional data would be unavailable without access to fresh cadaveric specimens. Because a primary goal of this study was to develop a method with applicability to the fossil record, we present these results but do not incorporate them into the broader workflow or statistical analysis in this paper.

An integrative understanding of the feeding apparatus requires an accurate, comprehensive characterization of muscular inputs and their concomitant impacts on joint and cranial function. 3D computational analysis of musculoskeletal behavior is time intensive; however, these techniques allow researchers to investigate performance in conditions and numbers that cannot be replicated in vivo. Researchers are able to visualize muscle resultants and cranial forces in three dimensions, correlate muscle inputs with cranial performance, and uncover loading environments of key cranial structures such as jaw joints and sutures. This modeling workflow can assess each muscle's contribution to various cranial forces, characterize joint reaction force magnitude and orientation to better understand the gross anatomical and microanatomical adaptations joints have to loading environment, and investigate intracranial joints, secondary craniomandibular joints and dual joint systems. The jaw joint minimally resists all dorsoventrally oriented input forces that do not contribute to bite force. As such, it plays a key role in modulating and dispersing forces in the feeding apparatus. We would expect that the jaw joint will have a tissue composition that is well suited to its loading regime. Therefore, biomechanical studies and histological investigations can provide reciprocal illuminations into how joint morphology and composition correlate with loading environment.

The 3D modeling approaches used in this study are well suited for application to the fossil record. Because these techniques use the area of muscle attachment site rather than the weight of dissected muscles to estimate PCSA, they can be used to make accurate estimates of muscle force, moments about axes and reaction forces in extinct animals. Effective use of osteological correlates (Holliday, 2009) and the extant phylogenetic bracket (Witmer, 1995; Holliday and Witmer, 2007) can constrain reconstructions of muscle attachment location, size and shape. Extant relatives can also inform the reconstruction of myological parameters such as pennation angle and fiber length.

The evolution of the modern crocodylian skull involved substantial changes to the skull (Langston, 1973; Busbey, 1995). In contrast to the platyrostral skulls of crocodylians, the earliest members of crocodylian-line archosaurs had oreinirostral skulls, such as the rauisuchian Postosuchus (Chatterjee, 1985) and the sphenosuchian Sphenosuchus (Walker, 1990). In these animals, the skull is dorsoventrally deep and mediolaterally narrow. The quadrates were dorsoventrally oriented (Walker, 1990) and, at least in sphenosuchians, were still not rigidly sutured to the braincase (Langston, 1973). A mobile suspensorium alongside an open palatobasal joint (Langston, 1973; Busbey, 1995) has been interpreted as evidence that these early ancestors of crocodylians were kinetic (Walker, 1990). The evolution of the crown group therefore involved substantial reorientation of adductor muscles along with a reduction in kinetic potential. The transformation from oreinirostry to platyrostry presumably required crocodylians to evolve higher mass or pennation of adductor muscles to achieve equivalent bite forces to their fossil ancestors.

The biomechanical origins of the pterygoid buttress are also poorly understood. In crocodylians, the pterygoid buttress articulates with the medial surface of the mandible in what has been described as an ‘open’ or ‘sliding’ joint (Schumacher, 1973). Some researchers have hypothesized that it braces the mandible against ‘reverse-wishboning’ (Iordansky, 1964; Schumacher, 1973; Busbey, 1995). Porro et al. (2011) included the pterygoid buttress as a constraint of the mandible. Recently, the pterygoid buttress has been suggested to represent a key innovation underlying the crocodylian feeding apparatus (Holliday et al., 2015). Although mediolaterally acting muscle forces will cancel out when symmetrically recruited, they will load structures on which they act. In crocodylians and other taxa with substantial mediolateral components to muscle force, these forces are likely to be a primary source of stress in the skull.

The loss of cranial kinesis along with the elaboration of the pterygoid flange into a novel craniomandibular articulation represent key features of crocodylian evolution (Holliday et al., 2015). Indeed, many of the hallmark features of the crocodylian skull, including an expanded retroarticular process, laterally attaching pterygoideus muscles, a bony secondary palate and broad scarf joints followed the evolution of the pterygoid buttress. The methods used in this paper will be applied to assess the kinetic status of fossil relatives of crocodylians and to investigate the role of the pterygoid buttress in the feeding apparatus of these derived archosaurs.

This study is one of the first to use both 3D LM and FEA to investigate the production of bite force in an ontogenetic series of A. mississippiensis. The use of anatomically accurate muscle attachments is key to the success of the models, and the good agreement between the two methods lends support to these techniques. The modeling techniques in this study can be used to assess the effect of changing muscle size and orientation during the evolution of the modern crocodylian skull. Key features of the crocodylian skull may have permitted novel muscular morphologies. These methods will allow researchers to test hypotheses linking bony features such as the loss of kinesis, secondary palate, scarf joints and the pterygoid buttress with muscular innovations, such as generally enlarged adductor mass, laterally inserting mPTv and the extensive cranial tendinous skeleton.

We thank Ruth Elsey and the staff of Rockefeller Wildlife Refuge for providing specimens. We thank the University of Missouri Biomolecular Imaging Center, the University of Missouri School of Medicine Department of Radiology and the University of Missouri School of Veterinary Medicine for scanning specimens. We thank Betsy Dumont, Larry Witmer, Laura Porro and Kent Vliet for helpful discussions. We thank two anonymous reviewers, whose comments greatly enhanced the quality and clarity of the manuscript.