Estimating bite force in extinct dinosaurs using phylogenetically predicted physiological cross-sectional areas of jaw adductor muscles

Functional muscle groups

Jaw adductor muscles in archosaurs are largely grouped and named based on developmental biology and various topological criteria such as their relative positioning to nerves and blood vessels (Holliday & Witmer, 2007). However, from a functional perspective, the existing groupings are not necessarily congruent with lines of actions in a lever model. For instance, the Musculus (M) pseudotemporalis superficialis (mPSTs) is topographically and functionally similar to the M. adductor mandibularis externus (mAME) but are developmentally linked to the M. pseudotemporalis profundus (mPSTp), the latter of which is often physically connected to (and indistinguishable from) the M. adductor posterior (mAMP). This largely stems from the fact that the mPSTs and mAME both have cranial attachments in the temporal fossa, while the mPSTp and mAMP both attach onto the quadrate (Fig. 1). Thus, the mPSTs and mAME work together as inter-linked functional in-levers while the mPSTp and the mAMP work together as a separate set of inter-linked functional in-levers.

I distinguish three functional adductor groupings, largely following Rayfield et al. (2001) and Sakamoto (2010), and identified muscle body complexes as follows: (1) the temporal muscle group (mTemp), consisting of mAME and mPSTs; (2) the quadrate muscle group (mQuad), consisting of the mPSTp and mAMP; and (3) the pterygoid muscle group (mPt), consisting of the M. pterygoideus (mPT). Practically, these approximate groupings are necessary as adductor muscles in smaller specimens are often difficult, if not impossible, to separate into the classic topological groupings, and as the goal of this study is to predict A_Phys in fossil species where we do not necessarily have detailed topological information. Furthermore, in the context of both biomechanical modelling and predictive modelling, approximation is often important in obtaining accurate predictions, which is not necessarily the most precise model.

Physiological cross-sectional areas in extant species

A_Phys for extant species (N = 39) were calculated from muscle architecture data collected predominantly from the literature but also from dissections (Sakamoto, Ruta & Venditti, 2019) (Struthio camelus, one specimen; Buteo buteo, three specimens; Larus fuscus, two specimens; Branta canadensis, one specimen; Gallus domesticus, two specimens; File S1; specimens were collected by the Bristol Ornithological Club and were donated to the University of Bristol as part of a clinical veterinary anatomy lab practical, c. 2005–2006 (Sakamoto, 2008)). A_Phys were calculated as: (1) $A_{\mathrm{P}\mathrm{h}\mathrm{y}\mathrm{s}}=(M\mathrm{c}\mathrm{o}\mathrm{s}\theta )/(\rho L)$following Sacks & Roy (1982), where M is the wet weight of the muscle body (g), θ is the mean pennation angle, ρ is the specific density (1.056 × 10⁻³ g/mm³ (Burkholder et al., 1994)), and L (mm) is mean fiber length. In the case of parallel fibers θ is 0° and thus cosθ is 1.

Muscle measurements for A_Phys calculations were taken for two specimens of Buteo buteo and one specimen each of Larus fuscus and Struthio camelus. Muscles were weighed prior to sectioning to obtain M. Muscles were carefully sectioned under a microscope using a sharp scalpel. Incisions were made parallel to the length of the muscle fibers as much as possible. L and θ were measured using ImageJ (Rasband, 2012). As sectioning muscles can artificially truncate fibers and ends of fibers can often be obscured and difficult to determine, fiber lengths were taken at multiple locations on one or more sections through each muscle, the mean of which was taken as L.

In some specimens, A_Phys were approximated using the gross cross-sectional area (A_Gross), as simply the cross-section taken perpendicular to the long axis of the muscle body (Rayfield et al., 2001). A_Gross was measured in one specimen each of S. camelus, B. buteo, and Branta canadensis, and two specimens each of L. fuscus, and Gallus gallus. The muscle body was sectioned roughly perpendicular to the major axis of the muscle body at its widest point, and its A_Gross was digitally measured using Image J. The mean value of the left and right sides was taken as the final A_Gross value. Comparisons between A_Gross measurements and A_Phys calculations (following Eq. (1)) taken from different specimens within the same species, reveal that measured A_Gross values are generally congruent with calculated A_Phys values (Sakamoto, 2008; Persons & Currie, 2011; Snively & Russell, 2007).

Physiological cross-sectional areas in extinct species

For extinct archosaurs, cross-sectional areas of the jaw adductor muscles were estimated as A_Gross using a variant of the dry skull method (Thomason, 1991), whereby cranio-mandibular dimensions (namely the areas of the supratemporal, subtemporal, and mandibular fenestrae) were used to bound the A_Gross of individual jaw adductor muscles. I measured A_Gross of the mAME, mPSTs, mPSTp + mAMP and mPT on photographs and diagrams of reconstructed skulls taken at various angles of view. These are conceptually similar to previously published methods to estimate A_Gross in extinct dinosaurs (Rayfield et al., 2001; Mazzetta et al., 2009; Reichel, 2010). I further applied muscle pennation angle θ = 45° for the mTemp group, θ = 0° for the mQuad group, and θ = 30° for the mPt group, based on average pennation angles in my extant archosaur samples. I applied the effects of pennation on to A_Gross through division of A_Gross by sinθ. This approximates A_Phys in fossil archosaurs.

Predictor variable

I used the width of the skull (W_Sk) as the predictor variable in the phylogenetic predictive models (PPMs). W_Sk was chosen here as it has previously been demonstrated to predict bite force (Herrel et al., 2005a; Sakamoto, 2021; Herrel et al., 2005b) and various jaw adductor muscles well (Fig. 2). It is also readily available from the literature and easy to measure for a vast number of species for which muscle data do not exist, both extant and extinct. The utility of its wider applicability makes simple measures like W_Sk an ideal predictor in predictive modelling. W_Sk were mostly measured directly from osteological and fossil specimens where possible but augmented with data taken from photographs and literature.

Phylogeny

I used an informal supertree of saurians based on the Time Tree of Life (TTOL) (Kumar et al., 2017) with fossil tips inserted manually at the appropriate phylogenetic locations (Sakamoto, Ruta & Venditti, 2019) (Fig. 3). Divergence times for fossil branches are based on first appearance dates (FAD) with terminal tips extended to their last appearance dates (LAD) using the paleotree R package (Bapst, 2012). I used the full range of temporal durations to scale the branches, as this allows for the maximum amount of time possible for trait evolution to occur (Sakamoto, Ruta & Venditti, 2019). Zero-length internal branch lengths were resolved by sharing time with neighbouring branches using the “equal” method (Bapst, 2012; Brusatte et al., 2008). While there are alternative methods to scaling branches, e.g., tip-dating using the fossilised birth-death model (Drummond & Stadler, 2016), phylogenetic regression under Brownian motion, which underlies the PPM framework, is extremely robust to uncertainties in branch lengths (Stone, 2011), so the choice of branch scaling makes minimal impact on PPM.

Phylogenetic predictive modelling

I used a Bayesian PPM (Organ et al., 2007; Sakamoto, 2021) to predict A_Phys in extinct archosaurs. Separate PPMs were fitted on each of the three muscle groups as outlined above with the relevant A_Phys as the response variable and W_Sk as the predictor variable.

Model performance, or prediction accuracy, of each PPM was evaluated in a dataset containing only the extant species (N = 39) first, through Leave-One-Out Cross-Validation (LOOCV). LOOCV procedure largely follows that outlined in (Sakamoto, 2021), and is as follows: (1) the PPM was first fit on the dataset leaving one species out (N − 1) using Markov Chain Monte Carlo (MCMC) generating a posterior distribution of predictive models; (2) the posterior predictive models were used to predict A_Phys for the species that was left out from Step 1 based on the W_Sk and phylogenetic position of that species; (3) the posterior distribution of predictions (posterior predictive distribution) was evaluated against the actual A_Phys value recorded for that species. If the observed value fell outside the vast majority of the posterior predictive distribution (i.e., beyond 95% of the distribution; p_MCMC < 0.05), then it is deemed that the actual A_Phys value is significantly different from the posterior predictive distribution, meaning that the prediction has failed in this particular species. I repeated these steps for all species in the data set (N = 39) and calculated the proportion of species for which the model succeeded in accurately predicting A_Phys out of the total sample size N.

I then predicted A_Phys for 53 fossil species of archosaurs (predominantly theropod dinosaurs). I first fitted a PPM on the N = 39 dataset and generated a posterior distribution of predictive models. I then used the predictive models to generate posterior predictive distributions for all 53 fossil species using their W_Sk and phylogenetic positions. This procedure is largely identical to LOOCV but is conducted in one step instead of one species at a time (Sakamoto, 2021). For 20 of the 53 fossil species, A_Phys measured from reconstructed muscles exist, thus allowing for assessment of the match between PPMs-predicted A_Phys and reconstructed A_Phys in fossil species.

Additionally, I evaluated prediction accuracy of PPMs on an expanded training set (N = 59) that includes A_Phys for select fossil species (N = 20) measured from reconstructed muscles (Sakamoto, Ruta & Venditti, 2019) or taken from literature (Lautenschlager, 2013; Gignac & Erickson, 2017; Bates & Falkingham, 2012; Lautenschlager et al., 2016). Prediction accuracy was evaluated through LOOCV as outlined above.

I then predicted A_Phys for the remaining 33 fossil species of dinosaurs (predominantly theropods). I first fitted a PPM on the N = 59 dataset and generated a posterior distribution of predictive models. I then used the predictive models to generate posterior predictive distributions of A_Phys for all 33 fossil species using their W_Sk and phylogenetic positions.

All model fitting was conducted in BayesTraits V3 over three independent MCMC chains each. The chains were run for 35,000,000 iterations, with the first 25,000,000 iterations discarded as burn-in, and sampled every 10,000 iterations after convergence, to produce a posterior sample of 1,000 predictive models and associated parameters.

Bite force estimation

Using the predicted A_Phys I estimated bite force (F_Bite) for 30 of the 33 fossil dinosaur species for which I predicted A_Phys through the PPM approach. I then compared those against F_Bite estimated for the 19 of the 20 fossil archosaurs based on measured A_Phys reconstructions included in the training set for the PPMs.

For each of the 30 fossil species for which I predicted A_Phys, I took the median value of the posterior predictive distribution for each muscle. Muscle force (F_Musc) was then estimated for each muscle as the product of A_Phys and tetanic stress σ at 0.3 N/mm² (or 300 kPa). The resulting F_Musc was then multiplied by the muscle moment arm to yield the torque of that muscle. I measured relevant moment arms for each muscle following the procedures developed in Sakamoto (2010) (Fig. 4). Muscle moments were summed and divided by the distance between the fulcrum (jaw joint) and bite point and multiplied by two to derive a bilateral F_Bite. F_Bite was estimated for the anterior-most and posterior-most positions along the biting edge (tooth row or beak; Fig. 4).

F_Bite in fossil archosaurs that were included in the PPM training set (N = 19) were estimated in the same way as above but using A_Phys measured from reconstructed muscles as outline in Sakamoto, Ruta & Venditti (2019). I compared F_Bite at the posterior-most positions (maximum F_Bite) between the two sets of fossil species.

Prediction accuracies of PPMs

Prediction accuracies of the PPMs in the initial training set consisting of only extant species (N = 39) was at 87% for all three muscle groups. The prediction accuracies of the PPMs in predicting A_Phys for the 20 fossil species, as compared to their measured A_Phys were 25%, 45% and 35%, respectively for the mTemp, mQuad and mPt groups.

Prediction accuracies of the PPMs in the expanded training set including fossil species (N = 59) were at 95%, 93% and 90% for the mTemp, mQuad and mPt groups, respectively. Out of the 20 fossil species included in the training set, in only two species (Plateosaurus engelhardti and Herrerasaurus ischigualastensis) did the PPMs fail to predict the observed A_Phys.

Bite force estimation

F_Bite estimated for the 30 fossil species based on predicted A_Phys are shown in Table 1. Compared to F_Bite estimated from reconstructed A_Phys in the 19 fossil species, these 30 F_Bite values fall along the expected relationship between F_Bite and W_Sk (Fig. 5). Comparisons between closely related species of similar sizes reveal the accuracy in resulting F_Bite values (Table 1): F_Bite for Deinonychus antirhoppus (predicted A_Phys) with W_Sk of 114.5 mm is 706N, while F_Bite for Dromaeosaurus albertensis (reconstructed A_Phys) with W_Sk of 103 mm is 885N; F_Bite for Carnotaurus sastrei (predicted A_Phys) with W_Sk of 300 mm is 7,172N while F_Bite for Majungasaurus crenatissimus (reconstructed A_Phys) with W_Sk of 300 mm is 7,845N.

Bite force estimates in extinct archosaurs

Using the A_Phys predicted from the PPMs, I estimated F_Bite in several extinct archosaurs. These values can be regarded as reasonably reliable estimates of true F_Bite in these extinct animals, given scaling relation with W_Sk and phylogeny. This owes to the fact that they are highly congruent with F_Bite estimates based on A_Phys measured from reconstructed muscles, which themselves have been demonstrated to be reasonably reliable estimates of F_Bite given size and phylogeny (Sakamoto, 2021). There are still some notable discrepancies between F_Bite estimated here with published values, namely in oviraptorosaurs (Meade & Ma, 2022) and ornithomimosaurs (Cuff & Rayfield, 2015), but also in Deinonychus (Gignac et al., 2010). My estimated F_Bite for oviraptorosaurs are under-estimates of published figures (Citipati osmolskae, 202–225N vs 349.3–499.0N (Meade & Ma, 2022); Incisivosaurus gauthieri, 26–41N vs 53–82.5N (Meade & Ma, 2022)) while those for ornithomimosaurs are over-estimates (Garudimimus, 121N at the tip vs 19N (Cuff & Rayfield, 2015); Ornithomimus, 121N vs 22N (Cuff & Rayfield, 2015); Struthiomimus, 108N vs 57.6N (Cuff & Rayfield, 2015)). The under-estimation of oviraptorosaurs is within the same order of magnitude (×10¹) and thus trivial in a comparative context (Sakamoto, 2021). On the other hand, the over-estimation of ornithomimosaurs is one order of magnitude (×10² vs ×10¹). The likely source of this discrepancy lies in the fact that neither oviraptorosaurs nor ornithomimosaurs were included in my N = 59 training set, thus making the predicted A_Phys reflecting values that are more typical of closely related theropod dinosaurs (Fig. 4; Table 1) as well as regressing to the mean of the training set (as is typical for regression). Future work expanding on the taxonomic sampling of the training set will undoubtedly improve prediction accuracy in specific taxa, especially in clades with unique dietary adaptations such as oviraptorosaurs (Sakamoto, 2010; Meade & Ma, 2022) or ornithomimosaurs (Cuff & Rayfield, 2015).

The case of discrepancy in F_Bite estimates for Deinonychus between this study and that by Gignac et al. (2010) warrants some special attention. Gignac et al. (2010) predicted the forces necessary to puncture bone to the depth observed in Tenontosaurus bones (Gignac et al., 2010; Erickson et al., 1996) and estimated the maximum F_Bite at the posterior-most biting position for Deinonychus at 8,200N, which is an order of magnitude higher than that estimated here at 706N (Table 1). My estimated F_Bite of 706N is in line with those of a similarly sized dromaeosaur Dromaeosaurus at 885N (W_Sk = 103 mm), and a similarly sized basal saurischian Herrerasaurus at 1,937N (W_Sk = 108 mm) (Table 1), but also extant diapsids with similar skull widths, Paleosuchus trigonatus at 1,082N (W_Sk = 121 mm), and Alligator sinensis at 1,084N (W_Sk = 122 mm) (Sakamoto, Ruta & Venditti, 2019). Interestingly the same can be said when the comparison is extended to extant carnivores with similar skull widths, Ursus thibetanus at 871N (W_Sk = 111 mm), Neofelis nebulosa at 1,296N (W_Sk = 119 mm), N. diardi at 1,117N (W_Sk = 115 mm), and Sarcophilus harrisii at 682N (W_Sk = 112 mm) (Sakamoto, Ruta & Venditti, 2019). It is important to note that forces necessary to puncture substrate are not confined to muscle-driven biting and may very well be the product of more aggressive kinetic feeding behaviours involving the whole head and neck (D’Amore & Blumenschine, 2009; Snively et al., 2013; Torices et al., 2018). This is supported by the observation that this bite mark in question was matched to a premaxillary tooth (Gignac et al., 2010), meaning that a long-snouted animal would have had to be capable of generating F_Bite of 3,000–4,000N (Gignac et al., 2010) at the tip of its snout through muscle generated biting, which is not congruent with the available data in comparably-sized amniotes (Sakamoto, Ruta & Venditti, 2019). Given its congruence with a similarly sized dromaeosaur Dromaeosaurus as well as other similarly sized amniotes, there is strong evidence to suggest that my F_Bite estimate for Deinonychus is reasonably accurate for an animal of its size.

Crucially, a comparison of F_Bite based on predicted A_Phys (N = 30) and F_Bite based on reconstructed A_Phys (N = 19) across the extinct non-avian theropods (Fig. 5; Table 1) generally show substantial overlap without any signs of sytemic biases (Fig. 5). This demonstrates that F_Bite estimates based on predicted A_Phys are neither systemically over- or under-estimates compared to F_Bite based on reconstructed A_Phys. Thus, PPMs are a useful approach to expand on reliable F_Bite data based on simple metrics and phylogeny to augment those based on reconstructed muscles. Importantly, this approach allows estimation of A_Phys and F_Bite in taxa (both extinct and extant) where muscle reconstruction is not feasible or possible for any number of reasons.

Discussion of F_Bite for individual species of interest are then valid and worth considering. Of note is that the large-bodied carnivorous dinosaurs, Carcharodontosaurus saharicus and Acrocanthosaurus atokensis, both reaching the size range of Tyrannosaurus rex, have F_Bite that are substantially lower than the latter, at 16,984 and 25,449N respectively, compared to 48,505N of T. rex. Carcharodontosaurus is approximately the same size as T. rex but is here shown to have had F_Bite that is approximately half of the latter. Carcharodontosaurus is typical in build and skull proportion for a theropod dinosaur, so the fact that its F_Bite was only half of that of T. rex is more likely a reflection of just how unique T. rex may have been compared to other theropods of similar sizes. Tyrannosaurus had robust conical-shaped teeth and multiple adaptations in the skull that allowed it to withstand immense forces (Cost et al., 2019; Rayfield, 2005; Snively, Henderson & Phillips, 2006). Multiple lines of evidence also point to habitual bone-crushing and -consumption in T. rex (Gignac & Erickson, 2017; Erickson et al., 1996; Erickson & Olson, 1996). These support the hypothesis that T. rex had at least a partial osteophagous diet, an ecology that was likely different from other theropods.

A similarly, large-bodied carnivorous dinosaur, Spinosaurus aegyptiacus, is here tentatively predicted to have had F_Bite at just under 12,000N, roughly in the same range as Sinraptor (10,845N), Gorgosaurus (13,817N), and Daspletosaurus (16,641N), all substantially smaller theropods. With the caveat that W_Sk for Spinosaurus was simply scaled up from the skull-width ratio of Suchomimus (Sereno et al., 1998), I offer additional support for this taxon to have had unique feeding habits for a theropod of its size. Spinosaurus shows adaptations in the craniomandibular morpho-functional complexes that are advantageous for generating relatively faster shutting speeds with less muscle input force (higher displacement advantage) at the expense of F_Bite (lower mechanical advantage) (Sakamoto, 2010). This would be congruent with a feeding mode relying on fast-snapping jaws rather than slow crushing bites, which is commonly observed in species with semi-aquatic feeding habits, including herons and egrets (Hone & Holtz, 2021).

Wider implications

The phylogenetic predictive framework I present here enhances collection of trait data that may be difficult to obtain across a wide taxonomic sample. Various in vivo measurements (e.g., bite force, muscle architecture) are obviously impossible to collect in extinct taxa, but they may also be difficult to obtain in many extant taxa, especially for those that are exotic, enigmatic, endangered, or dangerous. PPM then allows for predictions of trait data in such taxa provided that more readily available predictor variables (e.g., external physical traits) can be measured for them.

While the case study presented here focused on predicting muscle parameters, the PPM approach is not restricted to soft tissue predictions. The response variable Y of interest can be any single continuously varying trait as long as it exhibits significant relationships with a set of predictor variables X_i in taxa for which both Y and X_i can be measured. Examples include (but not precluded to): predicting body mass from skeletal measurements; predicting skeletal structural strength index from skeletal variables; or predicting metabolism from body mass.

Similarly, this framework is not restricted to any taxonomic group or scope. The taxonomic group of interest can be of any nature and breadth as long as a phylogeny exists that includes the taxa for which the Y variable is to be predicted.

A key feature of phylogenetic regression that makes the PPM framework extremely versatile is that it is extremely robust to uncertainties in terms of phylogenetic relations and branch lengths (Stone, 2011). This means that even a phylogeny with high levels of uncertainties (especially branch lengths for fossil trees) can be used effectively in a PPM framework to predict Y variables of interest. Nevertheless, modern approaches to branch scaling such as the tip-dating approach (Drummond & Stadler, 2016) that dates the divergences and tips simultaneously with topology inference (Gavryushkina et al., 2017; Stadler et al., 2018) can return a set of branch lengths that are guided by data and likelihood (i.e., probabilistic estimates), but these should be conducted using phylogenetic/cladistic data matrices (e.g., molecular or morphological characters), appropriate model of evolution, and tree prior, rather than on a fixed topology with no character data as is increasingly being used (Godoy et al., 2021; Gearty & Payne, 2020). The parameters associated with the latter are likely returning the priors as there are no data informing the likelihood, making the set of scaled branches no better than sampling randomly, but more investigation is needed.

It is important to note however that when considering Y and taxonomic sampling, the Y variable of interest should be broadly homologous across the phylogeny used, especially if the phylogenetic coverage is broad, spanning several higher order taxonomic groups, e.g., synapsids and diapsids. For instance, jaw adductor muscles of mammals (e.g., temporal and masseter muscles) are not directly comparable to those of diapsids (e.g., the mTemp, mQuad and mPt groups used here), but they are homologous within Mammalia across mammalian orders. On the other hand, biomechanical performance measures such as bite force is homologous across a very large portion of the tree of life, e.g., across vertebrates. The consideration here then would be the homology between the predictor variables across vertebrate clades with vastly diverse skull anatomy – e.g., skull measurements across various fish clades may not be precisely homologous, which is exemplified when the comparisons extend to amniotes.

There is no multi-response (multivariate regression) implementation as of date, but this is true for phylogenetic regression in general. This means that applications such as predicting morphospace coordinates would need to be conducted on individual shape (e.g., principal components) axis separately.