Beware the black box: investigating the sensitivity of FEA simulations to modelling factors in comparative biomechanics

Aims

Here we investigate the sensitivity of models in a broad scale comparative Finite Element (FE) dataset to different values of several modelling factors, to determine the extent by which the pattern of results is changed by the choice of different input values. The specific focus is on factors associated with material properties, scaling, linear load cases, and bite position.

Our approach is to make use of a previously compiled comparative dataset, which drew conclusions relating to form and function in extant crocodilians (Walmsley et al., 2013). As in the previous study, we simulate biting, shaking, and twisting feeding behaviours, which are typically used by crocodilians to process prey items. We explore many of the modelling factors inherent in the growing body of comparative biomechanical studies, and explicitly test the extent to which these factors influence, or change the pattern of results.

Factors affecting FE analysis

Finite Element Analysis (FEA) is a computational technique commonly used in engineering disciplines, whereby complex structures are discretised in order to approximate their mechanical response (behaviour) to applied loads. In recent years FEA has become increasingly prevalent in the fields of comparative biomechanics (McHenry et al., 2006; Oldfield et al., 2012; Walmsley et al., 2013), paleontology (Degrange et al., 2010; McHenry, 2009; McHenry et al., 2007; Tseng & Wang, 2010; Wroe, 2008; Wroe et al., 2013), biology (Dumont, Piccirillo & Grosse, 2005; Wroe et al., 2008), medicine (Chen et al., 2012; Omasta et al., 2012), and anthropology (Wroe et al., 2010), as improvements in computational capabilities mean lower entry level costs for researchers. In the context of comparative biomechanics FEA offers a number of advantages:

1. Biological structures included in comparative analyses often differ in size and shape, whilst structure-function questions typically focus on the role of shape. In Finite Element (FE) models, differences between specimens can be easily standardized through scaling (Dumont, Grosse & Slater, 2009; McHenry, 2009; Snively, Anderson & Ryan, 2010; Tseng, 2008; Walmsley et al., 2013).

2. Experiments can be quickly changed to test new hypotheses, simply by changing boundary conditions and loading.

3. Experimental time is greatly reduced, and since simulations are digital many more tests can be performed on a single specimen than would be feasible working with live animals or ex vivo specimens.

4. Hypotheses on form and function implications for extinct taxa can be tested (McHenry et al., 2007; Oldfield et al., 2012; Plotnick & Baumiller, 2000; Snively & Theodor, 2011; Tseng, 2008; Wroe et al., 2007a).

5. When combined with mesh deformation/warping (O’Higgins et al., 2011; Parr et al., 2012), purely theoretical morphotypes can be generated to help tease out the important features of shape that effect the structural response.

Despite the many advantages of FEA, there are limitations to the conclusions that can be drawn from the results (Rayfield, 2007). Finite element models are complex and informative simulations require deliberate choices for multiple factors (listed below, and here termed modelling factors). In many instances, biologically relevant empirical data for each factor are lacking and researchers necessarily make assumptions about realistic/plausible values to use as input variables for these modelling factors (McHenry et al., 2006).

The goal of many comparative analyses is to discover the pattern of differences in biomechanical performance between different models; that is, the relative order of the models’ performance under specific loads (e.g., strongest to weakest) and the degree by which they vary. While sufficient accuracy is critical in mechanical and biomedical engineering, for many comparative biomechanical studies the accuracy of absolute results is not required as long as the interspecific pattern of results resembles the actual biological pattern (McHenry et al., 2007; Oldfield et al., 2012; Parr et al., 2012; Rayfield, 2005; Tseng, 2008; Walmsley et al., 2013; Wroe et al., 2007a; Wroe et al., 2010). Whether the FEA results do reflect reality can only be examined if the results of the analysis can be compared with empirical data, a process termed validation. Although validation data has been used in a number of biological FE analyses (Bright, 2012; Bright & Rayfield, 2011a; Bright & Rayfield, 2011b; Gröning et al., 2009; Kupczik et al., 2007; Liu et al., 2012; Metzger, Daniel & Ross, 2005; Panagiotopoulou et al., 2010; Panagiotopoulou et al., 2012; Rayfield, 2011; Ross et al., 2005; Strait et al., 2005; Tsafnat & Wroe, 2011), the data required to validate models are difficult to obtain. Many comparative biomechanical analyses are thus run without validation (McHenry, 2009; McHenry et al., 2006; McHenry et al., 2007; Oldfield et al., 2012; Walmsley et al., 2013; Wroe et al., 2007a; Wroe et al., 2008). Combined with the lack of data on realistic values for modelling factors, this creates a degree of uncertainty about the validity of the results. In many cases, researchers assume (either explicitly or implicitly) that the precise value of input variables for modelling factors will not alter the pattern of results, as these values are kept constant across the models in the analysis, and the results obtained will be a valid reflection of the pattern. Whilst this is a logically plausible approach, it is seldom tested.

In the absence of the required empirical data to validate FE models, the sensitivity of results to the choice of input values for modelling factors can be explored through sensitivity analysis. In such an analysis, the input value for one factor is varied across models, while all other values are held constant; thus the effect upon the pattern of results can be quantified. If the pattern of results does not change markedly for different values, then the analysis is deemed relatively insensitive to the precise values chosen for that modelling factor. Where this is the case, the assumption – that the results of a comparative analysis can be informative, even in the absence of empirical data on modelling factors and absence of model validation – remains logically plausible (although still untested). If, however, the pattern of results is strongly affected by the precise values used for modelling factors, then the analysis is sensitive to input parameters and its results are only informative if input parameters are founded upon empirical data.

Investigations into the sensitivity of FEA simulations have been performed in relation to a number of different modelling factors. These investigations have targeted input values associated with scaling (Dumont, Grosse & Slater, 2009), material properties (Bright & Rayfield, 2011b; Cox et al., 2011; Gröning, Fagan & O’Higgins, 2012; Kupczik et al., 2007; Panagiotopoulou et al., 2010; Porro et al., 2011; Reed et al., 2011; Strait et al., 2005; Tseng et al., 2011; Wroe et al., 2008), muscle activation (Fitton et al., 2012; Ross et al., 2005; Tseng et al., 2011), sutures (Bright, 2012; Kupczik et al., 2007; Porro et al., 2011; Reed et al., 2011; Wang et al., 2010), bite position (Cox et al., 2011; Fitton et al., 2012; Porro et al., 2011; Wang et al., 2010), ligaments (Gröning, Fagan & O’Higgins, 2012; Gröning, Fagan & O’Higgins, 2011; Wood et al., 2011), mesh density (Bright & Rayfield, 2011a; Tseng et al., 2011), mesh warping (O’Higgins et al., 2011), jaw joint constraint (Gröning, Fagan & O’Higgins, 2012; Tseng et al., 2011), orientation of muscle force (Bright & Rayfield, 2011b; Cox et al., 2011; Gröning, Fagan & O’Higgins, 2012; Grosse et al., 2007), muscle loading application (Grosse et al., 2007; Kupczik et al., 2007; Wroe et al., 2008), FEM element type (Bright & Rayfield, 2011a; Dumont, Piccirillo & Grosse, 2005), and subcortical geometries (Panagiotopoulou et al., 2010). This growing body of literature has helped to identify those modelling factors which most affect simulation results (information which is invaluable to comparative studies where validation is unfeasible or impossible). However, the majority of sensitivity analyses to date involve a single specimen; there are limited instances of sensitivity analyses involving either multiple specimens of one species (Kupczik et al., 2007), or multiple species (Cox et al., 2011). As an important goal of many comparative analyses is to ascertain the pattern of relative biomechanical performance between taxa, multi-factorial, multi-species sensitivity analyses allow us to assess how suitable FEA is for comparative studies in the absence of validation.

Modelling factors

Modelling factors in comparative FEA are specific aspects of model set-up that can influence results in comparative simulations. Common modelling factors include, but are not limited to: scaling, material properties, simulated feeding behaviour, linear versus non-linear load cases, sutures, bite position, muscle activation schemes, muscle proportions, number of muscles, constraints, and how muscles are modelled in the FE simulation. Each of these factors can often be implemented in a number of different ways; for example, muscles are typically modelled either as beams or as point loads, and both of these implementations contain subsets of options, such as beam geometry and directionality of the point load. In addition to the numerous combinations in which they can be sensibly assembled, the cascade of assumptions within each modelling factor leads to a very large parameter space, which can potentially produce appreciably different results if sensitivity to these is high.

Exploring this parameter space is logistically complex. In the absence of empirical data that can be used to select realistic input values for the various factors, exploration of parameter space provides a sensitivity analysis but does not necessarily improve model accuracy. In this paper we present a comprehensive sensitivity analysis of the following five modelling factors, which are each specifically relevant to questions about skull optimisation (minimized stress/strain) for given feeding scenarios in crocodilians.

FEA models

For this analysis many methodological aspects (particularly those relating to data acquisition, CT processing, and surface/solid mesh generation) are identical to the previous study (Walmsley et al., 2013). In summary: CT data was processed in MIMICS v11 (MATERIALISE, Belgium), surface meshes of the cranium and mandible were optimised before forming the foundation to generate suitable solid meshes using Harpoon (SHARC), and FE simulations were performed using Strand7 (www.strand7.com).

High-resolution finite element model construction was based on previously published protocols (Bourke et al., 2008; Clausen et al., 2008; McHenry et al., 2007; Moreno et al., 2008; Wroe et al., 2007a); see Walmsley et al. (2013) for specific details. In the present study, we varied the parameters (described below in detail) of several modelling factors: Material properties: models are simulated with either isotropic-heterogeneous or -homogeneous material properties; Scaling: models are scaled to a consistent volume, surface area, or length; Feeding behaviour: models are loaded to simulate biting, shaking, and twisting feeding behaviours; Linear load case combinations: loads are scaled to 2 metrics per feeding behaviour (see below for details); Bite position: simulations are performed at front, mid, and back bite positions. Each of these variables is altered whilst holding all others constant, allowing all possible combinations of variables to be investigated across the seven species simulated. Specific feeding behaviours are strictly functional variables and are not considered to be a target of this sensitivity analysis – it is nonsensical to investigate strength under twisting as an indicator of strength under biting – but these load cases do increase the parameter space investigated in this study by a factor of three, and including them gives some insight into how sensitivity to a modelling factor is affected by functionally different loading conditions. Whilst the breadth of modelling factors explored here is extensive, it is by no means complete; for example, our simulations do not account for the influence of structures such as sutures, which have an important effect on biomechanics (Bright, 2012; Kupczik et al., 2007; Porro et al., 2011; Reed et al., 2011; Wang et al., 2010).

Data collection and presentation

We collected, assessed and here present data in multiple formats to ascertain the degree that multi-dimensional variation of input parameters (for common modelling factors) has on the results and their interpretation. In brief (outlined in detail below), the presented formats are:

• Signal: qualitative visual comparison between pairs or triplets of sets within one species.

• Rank: qualitative comparison between pairs of conditions between multiple species.

• Percentage Difference and Mean Percentage Difference: quantitative comparison between pairs of conditions within one species.

• Pattern and Standardised Pattern: quantitative comparison between all conditions, for multiple species.

• Standardised Pattern Difference (SPD): quantitative comparison between sets, for multiple species, and allowing comparison of qualitatively- vs quantitatively-ordered data.

• Shape correlations: quantitative comparison between shape differences and set differences, for multiple species.

As in Walmsley et al. (2013) the results assessed here are the 95% von Mises strain values; this 95% von Mises strain constitutes the largest elemental (individual brick) value of strain in the model if the highest 5% of all elemental values are ignored. It is important to note that measuring strain in this way only accounts for the magnitude of strain; it compares the magnitude of upper strain values but does not include any information on the type of strain (i.e., compressive or tensile), or upon the location of that strain within the anatomical structure. Each individual result is a combination of specific values for each parameter; we use condition to describe that combination of parameters. In total 108 unique conditions exist for each species model, each describing the type of scaling (volume, surface, or length), bite position (front, mid, or back), feeding type (bite, shake, or twist), material properties (HET or HOM), and specific LLCs (one of a possible two per feeding type) used in an individual simulation. A set is used to describe an arbitrarily ordered group of conditions ( $X=\left\{x_1,x_2,\dots ,x_n\right\}$ ) with a common parameter. When comparing between two sets, the number of conditions in each set depends upon the number of values (parameter options) for the modelling factor; where there are two parameter options (i.e., for material properties and LLC) the number of conditions in a set is 54, and where there are three parameter options (i.e., for scaling, feeding type, and bite position) there are 36 conditions per set. Thus, a ‘volume’ set would include all 36 conditions where the model is scaled to volume, and a ‘HOM’ set includes all 54 conditions where the model has isotropic homogeneous material properties.

Signal

Signal is plotted as the microstrain value for each condition, with the set of conditions for each value plotted along the x axis. For signal, all sets for a modelling factor can be plotted simultaneously (e.g., Figs. 1A and 1B). This gives a chart that superficially resembles a signal trace. By treating the strain response of each specimen as a signal to a set of unique conditions, differences between the variables for each modelling factor can be compared at a holistic level, and the closer one signal follows (or tracks) to another, the less influence that modelling factor has upon the results (Figs. 1A and 1B).

Material properties (isotropic HET vs isotropic HOM)

Results for HET and HOM models closely match each other across all conditions, both qualitatively and quantitatively. The signal of HOM models tracks closely to HET (Fig. 2) and percentage differences between each pair of conditions were consistently low, averaging <6% for all species excluding C. intermedius and C. johnstoni, which each averaged ≈10% (Table 2 and Fig. S1). Between conditions, the greatest differences were for those that involved twisting, but mean percentage difference values remained below 10% (Table 2). Between species models the largest differences were for C. johnstoni and C. intermedius, while the smallest were for M. cataphractus and C. novaeguineae (Fig. S1).

Consistency in ranking (Table 3) was very high, with 24 of the 54 condition pairs predicting identical rankings, and a further 22 pairs differing in the rank of 2 models only. Of the remaining condition pairs, 5 differed in rankings by 3 or 4 species models, and there was 2 instances of ‘5 out’(for a detailed account of how well each condition pair predicted rank see Fig. S2). Charts of pattern (Fig. 3) and standard pattern (Fig. 4) likewise show that for each HET-HOM condition pair, strain values are very similar (horizontal parts of the trace). Standardised pattern difference (SPD) showed high consistency between conditions when ordered either by rank or mean standardised pattern difference (Fig. 5A), with low values of mean SPD (<0.1ε_Mc) throughout. Mean percentage difference showed no correlation with shape as measured by ΔPC1 and ΔPC2 (Fig. 6B).

Scaling

Strain in volume-scaled models closely matched that of surface-scaled models across all condition pairs. The signal of models track closely (Fig. 7), and consistency in rankings is high (22 identically ranked condition pairs and 9 pairs that differ by two models, out of a total of 36 condition pairs in the set). Similarly, pattern (Fig. 3), standard pattern (Fig. 4), and standard pattern difference (Figs. 5C–5E), all show very small differences between volume- and surface-scaling.

For each of these qualitative and quantitative measures, length-scaled models exhibited quite different strain values from both volume- and surface-scaled models across all condition pairs. Rankings (Table 3) for length- vs volume-scaled models had 14 identical predictions, and 6 conditions that differed in the order of 2 species, out of a total of 36 conditions; for length- vs surface-scaled the equivalent numbers were 11 and 3 respectively. Plots of signal (Fig. 7) indicate that, while scaling to length follows the same broad trend as volume- and surface-scaling, for individual conditions it consistently shows the greatest deviation (akin to noise) in its signal.

This pattern of results is also evident in percentage difference values; volume- vs surface-scaled models show the smallest differences in strain response overall, with the average for individual species ranging from ≈1% for T. schlegelii, through to ≈8% for C. johnstoni (Table 2 and Fig. S3). Volume- vs length-scaling shows much larger averages for individual species, from ≈2% for C. intermedius, through to ≈32% for C. moreletii (Table 2 and Fig.S4). Similarly surface- and length-scaling show large differences, from <1% for C. intermedius, through to ≈34% for O. tetraspis (Table 2 and Fig. S5). Crocodylus intermedius shows very little difference between all three scaling parameters, with a maximum average difference of ≈2%, and absolute max difference of ≈5% for volume- and length-scaled simulations (Table 2 and Fig. S6). Between species models the largest and smallest differences were identical for volume- and length-scaling, and surface- and length-scaling, with C. intermedius and C. novaeguineae displaying the smallest differences, and O. tetraspis and C. moreletii displaying the largest (Table 2, Figs. S4 and S5). Between volume- and surface-scaling C. intermedius and T. schlegelii show the smallest differences, while C. johnstoni and O. tetraspis show the largest (Table 2 and Fig. S3).

Between conditions, for all species models and each scaling comparison, the greatest differences were for those that involved twisting, although this was less pronounced in C. intermedius between surface- and length-scaling (Table 2, Figs. S3–S5). These large differences are also apparent in pattern (Fig. 3), and standard pattern (Fig. 4), which both show larger variation between scaling parameters for twist when compared to either bite or shake. Additionally, the largest differences for SPD are overwhelmingly dominated by twist feeding behaviours, which all fall in the worst half of SPD, and those conditions also involving MeM linear load cases consistently perform worst of all (Figs. S6B, S7B and S8B).

Qualitative and quantitative measures of sets gave inconsistent results for comparisons between length- and either volume- or surface-scaled models, in that those conditions that predict identical rank show some of the largest differences in SPD (Figs. 5D–5E, Figs. S7 and S8). Comparing between volume- and surface-scaled models shows much higher consistency between qualitative and quantitative measures; conditions with the smallest variation in SPD were predominantly identical or near predictions of rank (Fig. 5C and Fig. S6).

Mean percentage differences between volume- and surface-scaled models show no correlation with shape, as measured by ΔPC1 and ΔPC2; however, the larger variation in results between length- and both volume- and surface-scaled models showed some correlation with shape (Figs. 6A, 6D and 6E). In both cases mean percentage difference correlated well with ΔPC2, and poorly with ΔPC1, with surface-length comparisons showing r² values of 0.67 and 0.48, and volume-length 0.85 and 0.39, for ΔPC2 and ΔPC1 respectively.

Linear load cases

Results between LLC models show large variation across conditions, both qualitatively and quantitatively. For some combinations of species model and conditions, TeT/MeM results correlate well with NoLLC/ELA, whilst for others the agreement is low. Conditions involving both volume-scaling and twisting show good correlations across all species, while conditions involving length scaling and biting show consistently poorer correlations, ranging from an average percentage difference of 13% for C. novaeguineae through to 58% for O. tetraspis (Table 2 and Fig. S9); signal also shows large differences for those conditions (Fig. 8). The largest deviations in SPD were always biting conditions, with the very worst also involving length scaling (Fig. S9B).

With respect to species models, C. novaeguineae shows good correlations, while O. tetraspis shows poor correlations overall, but even this is inconsistent; the signal for TeT/MeM models tracks NoLLC/ELA closely for volume-scaled twist conditions, but tracks poorly for length-scaled biting (Fig. 8A). Percentage difference between each pair of conditions shows considerable variation (Fig. S9), ranging from a minimum of 2% through to maximum of 59% for O. tetraspis. The largest mean percentage differences were for O. tetraspis (avg. ≈21%), and C. johnstoni (avg. ≈17%), with C. novaeguineae showing the smallest ≈6% (note that the M. cataphractus models have zero differences since LLCs are equal for this species model) (Table 2, Fig. S9).

Pattern (Fig. 3) and standardised pattern (Fig. 4) show small and large differences between LLCs; for example, in the C. intermedius models, shake and twisting conditions show small differences between LLCs, while large differences are seen for biting conditions. This variation in quantitative pattern is illustrated by the SPD (Fig. 5B), where mean standard pattern difference varies from almost indistinguishable through to ≈ 0.5ε_Mc. Additionally SPD for individual species models ranged from almost indistinguishable from M. cataphractus, to more than 0.8 of that benchmark (Fig. S10).

Consistency in ranking (Table 3) was low, with only 16 of the 54 condition pairs predicting identical rankings, and a further 6 pairs differing in the rank of 2 models only. Of the remaining conditions, 10 were out by 3 or 4, and 21 reported substantially different rankings (out by 6 or 7). With respect to consistency between qualitative and quantitative results, predictive rank displayed an appreciable spread when ordered by SPD, but the smallest differences in pattern are still dominated by identical or near (‘2 out’) predictions of ranked order (Fig. 5B and Fig. S10); conditions that were qualitatively consistent were also quantitatively similar.

High variation in LLC results showed some correlation with shape; mean percentage difference showed good correlation with ΔPC2 (r² = 0.77), but poor correlation with ΔPC1 (r² = 0.21). This suggests that sensitivity of models to LLC is related to shape (Fig. 6C), particularly those aspects of shape captured within PC2 – inter-rami angle, followed closely by symphyseal length and mandibular width; see Figure 19 in Walmsley et al. (2013). Plots of signal support this observation, in that where differences are large, TeT/MeM over-predicts compared to NoLLC/ELA in species models with long and narrow rostra (Figs. 8D, 8E and 8G), and under predicts for those with more robust, broad and short rostra (Figs. 8A and 8B).

Bite position

Results between front, mid, and back bite positions exhibit appreciably large differences across all conditions, both qualitatively and quantitatively. All three bite positions show poor correlation in signal (Fig. 9), poor predictive rank (Table 3), large percentage differences (Table 2, Figs. S11–S13), as well as large differences in pattern (Fig. 3), standard pattern (Fig. 4), and standard pattern difference (Figs. 5F–5H and Figs. S14–S16).

The overall waveform of signal for front, mid, and back bite positions remains reasonably consistent across all conditions for each species, mainly varying in amplitude (Fig. 9); however, this variation is large (e.g., M. cataphractus) and not uniform throughout conditions (e.g., O. tetraspis shows smaller variation in bite conditions than in shake).

Pattern (Fig. 3) shows large differences between all three bite positions, where response decreases and compresses across all species models when moving from front to back positions. Although somewhat less noticeable, standard pattern (Fig. 4) also shows reasonable differences for all species models. SPD also shows large differences, with individual species models extending beyond 0.4 of M. cataphractus for most conditions (Figs. S14–S16), and averaging >0.1 of M. cataphractus across most conditions for front and mid, and front and back condition pairs (Figs. 5F–5G).

While percentage differences are typically large between all bite positions, mid and back show the smallest overall, with the average ranging from 29% for C. moreletii to 39% for C. novaeguineae (Table 2 and Fig. S13). Front and back show the largest difference, ranging from 45% for C. moreletii, through to 66% for T. schlegelii (Table 2 and Fig. S12), and similarly for front and mid, the average ranges from 26% for C. moreletii, through to 48% for T. schlegelii (Table 2 and Fig. S11).

Ranked order of specimen is highly sensitive to bite position, with a large proportion of simulation conditions resulting in substantially different predictions (Table 3). Of the 36 possibilities 15 were out by 5 or more between front and mid condition pairs, 21 for front and back, and 11 for mid and back. While identical predictions were only observed 8 times for front and mid bite position pairs, once for front and back, and twice for mid and back, slight differences were somewhat more frequent, particularly between front and back, and mid and back condition pairs (Table 3).

Comparisons between either front and back or mid and back bite condition pairs show a low consistency between qualitative and quantitative results, in that those conditions that predict identical rank show large differences in SPD, and the smallest differences in SPD are consistently very poor predictions of rank (Figs. 5G–5H and Figs. S15–S16). Between front and mid condition pairs, good predictors of rank spread appreciably when ordering conditions by SPD, although the best half of conditions ordered by SPD predominately consists of good predictors of rank (Fig. 5F and Fig. S14).

Between conditions, for all species models and each bite position comparison, the largest differences were for those that involved shaking – with the exception of C. novaeguineae between front and mid positions, whose largest differences were for twist conditions (Table 2, Figs. S11–S13). These large differences are also apparent in pattern (Fig. 3), and standard pattern (Fig. 4), where much larger variation is apparent between bite positions for shake compared to either bite or twist. The smallest SPD is also dominated by conditions involving biting, although somewhat less pronounced between mid and back positions; while the largest are dominated by conditions involving twisting, specifically those also involving HET material properties, which is less pronounced for front and mid positions (Figs. S14–S16).

Mean percentage differences show no correlation with shape, as measured by ΔPC1 and ΔPC2, for all three bite position comparisons (Figs. 6F–6H).

Material properties (Isotropic HET vs Isotropic HOM)

Qualitatively and quantitatively the selection of either HET or HOM material properties (as we calculated these) made little difference in the interpretation of results. This is evident from the small differences in signal (Fig. 2), percentage difference (Table 2 and Fig. S1), pattern (Fig. 3), standard pattern (Fig. 4), and standard pattern difference (Fig. 5A and Fig. S2), as well as the large proportion (46 of 54) of conditions that predict identical or near (2 out) specimen rankings (Table 3). Interestingly these differences are small despite HOM material properties for all species models being calculated from the average of M. cataphractus, and not from their own HET average (Table 4).

The fact that conditions involving twisting displayed the greatest sensitivity to the selection of material properties may relate to differences in material stiffness at the outer surface of HET models compared with HOM models; during elastic torsional loading material furthest from the axis of rotation carries a higher proportion of the load (Spotts & Shoup, 2004). This result suggests that torsional loads may be at least as important as bending loads in determining the distribution of cortical bone within beam-shaped skeletal elements.

Scaling

Qualitatively and quantitatively, scaling to either surface or volume made little practical difference upon the results or their interpretation. This is evident from the small differences in signal (Fig. 7), percentage difference (Table 2 and Fig. S3), pattern (Fig. 3), standard pattern (Fig. 4), and standard pattern difference (Fig. 5C and Fig. S6), as well as the large proportion (31 of 36) of conditions that predict identical or near (2 out) specimen rankings (Table 3).

Comparing length- to either volume- or surface-scaling made a bigger difference in the results, displaying large differences in signal (Fig. 7), percentage difference (Figs. S4 and S5), all measures of pattern (Figs. 3, 4, 5D, 5E and Figs. S7–S8), and a large proportion of inconsistent rank predictions (Table 3). The higher sensitivity to length-scaling is related to the spectrum of skull shape in crocodilians, ranging from longirostrine through to brevirostrine taxa (Busbey, 1995; Langston, 1973; McHenry et al., 2006). Scaling to length is arguably appropriate for exploring the consequences of different head length morphologies and symphyseal morphologies, however this needs to be used very carefully; a brevirostrine animal with the same head length as a longirostrine would be a much larger animal with a much stronger skull. Differences between length- and either volume- or surface-scaling appear to be a function of shape, where the largest differences are seen in both relatively shorter and broader (O. tetraspis and C. moreletii) or longer and narrower (T. schlegelii), skulls than M. cataphractus (Table 2 and Fig. 7); additionally this is supported by the strong correlations with ΔPC2 scores (Fig. 6).

The differences in results between all three scaling parameters are a function of the proportional difference between the linear scaling factors (LSF) used to scale models to volume, surface, and length (Fig. 10). Larger proportional differences between LSFs directly translate to larger differences in the response of models after scaling to one parameter or another. This explains why length-scaled models have such different results to both volume- and surface-scaled models; the difference between the LSFs of length- compared to both volume- and surface-scaling is proportionally larger than between volume- and surface-scaling.

Similar to material properties, conditions involving twisting display the greatest sensitivity to the selection of scaling parameters, consistently showing the largest absolute percentage difference across all species (Table 2, Figs. S3–S5), and dominating the largest standard pattern difference (Figs. S6B, S7B and S8B), particularly those conditions also involving MeM Linear Load Cases. In this regard conclusions relating to twisting feeding behaviours should be considered carefully, since the selection of one scaling parameter over another has a substantial influence over how the results would be interpreted.

Results show a high sensitivity both qualitatively and quantitatively to simulations where models are scaled to length as opposed to either surface or volume, and while this doesn’t speak to the appropriateness of one scaling parameter over another, it does suggest that the selection of length as a scaling technique should be well justified, since it is likely to dramatically change the pattern of results and their interpretation.

Linear load cases

Selection of appropriate LLC is important, and the interpretation of results would be largely dependent on which were used in the simulations. This is evident from the large difference in signal for most species across a number of simulation conditions (Fig. 8), the high proportion (26 of 54) of conditions that badly (4 or more out) predict rankings (Table 3), and large SPD – averaging > 0.1ε_Mc for most simulation conditions (Fig. 5B). Qualitatively conditions that involve biting show the greatest sensitivity to the selection of LLCs, showing very poor predictions of ranked order in addition to accounting for all of the largest differences in SPD (Fig. S10B).

In this analysis most conditions that show good predictions of rank also show the smallest variation in SPD (Fig. 5B and Fig. S10B). This means that those conditions that show good predictions of rank are also quite similar in regards to their pattern of results, and thus selection between the LLCs presented here becomes somewhat arbitrary since each yield similar results. However, this information could only be acquired through an extensive sensitivity analysis such as this, and is unlikely to remain true for other comparative datasets.

Quantitatively, absolute percentage differences (Fig. S9) vary considerably with respect to scaling parameter, feeding behaviour, and the model species, displaying both very large and very small difference for different combinations of these parameters. The only distinctive trend is that the largest difference for all species occurring under conditions combining biting and length-scaling. This range of difference in the results suggests that LLCs are far more sensitive to combinations of factors than to any one factor, particularly those combinations relating to the variation in skull shape which changes the values used for each LLC. In the example of shaking and the two LLCs used here, one simulates an identical lateral force across all species, and the other simulates a constant ratio of outlever-length to lateral force – i.e., shaking identical mass at an identical frequency. For the applied force for each of these simulations to be identical (and thus the microstrain results), scaling must be such that outlever-length is identical for each model. In this way the small differences between LLCs for shaking manifest as a result of outlever-length being very close to that of M. cataphractus at the rescaled size, and is most evident in conditions involving length-scaled shaking where differences in out-lever length are smaller for all species (Fig. S9).

In many comparative analyses an arbitrarily selected (normally equal) load is simulated on all specimens, with the prevailing logic that after size is accounted for, all that remains to influence biomechanical response is shape. Importantly, by simulating identical forces across all species, information about the functional aspect of that feeding behaviour is lost. In shaking, simulation of an equal lateral force results in each animal NOT shaking a prey item of the same mass at an identical frequency. Conversely by simulating identical mass and frequency, the selection of an appropriate scaling parameter becomes much more important since the simulated force is calculated by outlever-length, which is an aspect of shape determined by the scaling parameter. Similarly the forces calculated for twisting and biting would also be influenced by aspects of shape that can be over- or under-stated as a result of scaling parameter selection.

Bite position

Qualitatively, selection of either front, mid, or back bite positions in simulations is important, and interpretation of results would be largely dependent on which were used. This is evident from the large differences in signal between all three bite positions (Fig. 9), the small proportion (8, 1, and 2 of 36 for front-mid, front–back, and mid-back comparisons respectively) of conditions that predict identical rankings (Table 3), and the large differences in SPD, averaging >0.1ε_Mc for most simulation conditions (Figs. 5F–5H). From a quantitative point of view, those conditions involving bite show the least sensitivity to the selection of bite position across all species, consistently showing the smallest absolute percentage differences (Figs. S11–S13), and additionally dominating the smallest differences in SPD for all comparisons (Figs. S14–S16). However, absolute percentage difference and SPD for bite conditions is much larger than that seen for either material properties (Figs. S1 and S2) of volume- vs surface-scaling (Figs. S3 and S6).

The combination of large differences in pattern (Fig. 3), standard pattern (Fig. 4), standard pattern difference (Fig. 5), the small number of identical rank predictions (Table 3), the large absolute percentage differences (Table 2), and the large differences in signal (Fig. 9) between all three bite positions, together illustrates a very important point for comparative biologists. Broad assertions (or interpretations) about skull optimisation for a specific feeding type cannot be inferred from simulations at a single bite position, since the pattern of results between specimen changes dramatically depending on the selection of bite position. For example some skulls may be better optimised for back, or mid biting, than they are for front biting, so simulations of front biting should only be used to make interpretations about that specific behaviour and not extended to biting in general. For a more comprehensive understanding of skull optimization for a specific feeding behaviour when multiple bite positions are feasible, each bite position must be analysed separately and conclusions drawn from the aggregation of all data. Further to this point, where observational data relating to feeding behaviours is available, it should be incorporated into the simulations so that comparisons remain logical in the context of their biological reality. If species A, B, and C are all know to engage in shake feeding behaviours but species B tends to grip prey at its mid bite position, while species A, and C tend to grip prey at a front bite position, the most logical comparison is not simulating all 3 shaking at a front bite position, but A and C at front, and B at mid. Simulations performed in this way, guided heavily by accurate observational data, are likely to better reflect biological reality, and additionally increase confidence in the results they provide.

Feeding behaviour

For comparative simulations conclusions can only be drawn for the specific feeding behaviour being compared, and blanket conclusions relating to performance cannot be inferred from one behaviour to another. For instance if the results from a simulation relating to biting suggests one specimen performs better than another, this does not mean that for a different feeding behaviour (i.e., twisting or shaking) the same relationship exists. While simulations relating to different feeding behaviour are used here, we do not compare predictions about overall skull performance between different feeding behaviours, since they are functionally incomparable.

Overall patterns

Model sensitivity varied between modelling factors; the highest sensitivity was for bite position with an average percentage difference >30% for all bite position comparisons (Fig. 11A). All other modelling factor comparisons averaged <20%, with volume- vs surface-scaling, and material property selection showing the smallest differences, both averaging ≈5%. Individual feeding behaviours show varied degrees of sensitivity to modelling factors, with bite being the least sensitive, averaging <30% for all modelling factors (Fig. 11B). Linear Load Case comparisons are not directly comparable between feeding behaviours (since each scale loads differently — see methods); however, bite load cases were highly sensitive to LLC, far more so than either shake or twist. This is likely due to the functional difference between the TeT and NoLLC conditions; TeT simulates a standardised bite force across all species models, so that all but M. cataphractus are biting either above or below their calculated maximal bite force, while NoLLC simulates maximal muscle recruitment for each species model. Shake shows the highest sensitivity to all bite position comparisons (average >40%) compared to either bite or twist, and very low sensitivity (average <10%) to all other modelling factors (Fig. 11C); the high sensitivity to bite position is likely a result of the applied loads being a function of outlever-length, which changes dramatically between bite positions. While less sensitive to bite position than shake, twist also shows high sensitivity to scaling, specifically to length- vs either volume- or surface-scaling (Fig. 11D).

Recent comparative analyses have directed a lot of attention to the importance of scaling to either volume or surface (Dumont, Grosse & Slater, 2009) in addition to the need for accurate material properties (Wroe et al., 2007b). However, our results show that our models are not nearly as sensitive to these factors as they are to bite position or linear load cases (at least, for the way we have modelled these). Interestingly, these are generally accepted in the literature without question. In particular, bite position has a large influence on the pattern of results, and the high sensitivity of models to bite position emphasises the importance of using empirical data on behaviour as input variables for any comparative modelling analyses; specifically, that it’s much more important than the validation of material properties of bone or the selection of either volume- or surface-scaling.