Rethinking the use of finite element simulations in comparative biomechanics research

Introduction

Structure-function relationships underlie many hypotheses about patterns of morphological disparity in organisms (Lauder, 1981; Lauder & Thomason, 1995). Among the major tools in estimating the functional performance of vertebrate structures in particular is the use of biomechanical simulations to predict traits like bite force, structural stiffness, mechanical efficiency, etc. (Richmond et al., 2005; Ross, 2005). In the past 15 years, the adaptation of finite element (FE) modeling, an engineering method for solving load-deformation scenarios using principles of continuum mechanics, to biological questions has been greatly expanded across the broad spectrum of organismal study systems. An important feature of FE approaches is the ability to quantitatively test functional morphological hypotheses, contrasting it with largely qualitative conclusions and inferences in classic comparative and functional anatomical approaches A search in the “Web of Science” database returned 768 publications between 2005 and 2020 on using FE models in the study of evolution and biology (Fig. 1A). Research using FE analysis of the vertebrate skeleton covers topics such as inferring locomotory and masticatory performance in the vertebrate fossil record (Rayfield, 2007), morphofunctional evolution in entire clades (Pierce, Angielczyk & Rayfield, 2008), to evolutionary optimization of functional morphological attributes (Polly et al., 2016), and others.

Studies using the FE method to test hypotheses about vertebrate structure-function fall into two major categories: 2D and 3D analyses. 2D FE models are typically derived from photographs of specimens, whereas 3D FE models are typically derived from computed tomography (CT) or surface scans. The main trade-off between 2D and 3D approaches is model sample size vs time investment in building each model (Morales-García et al., 2019). 2D models are quicker to build and allow for larger taxonomic sample sizes, but the extent of model simplification restricts its application to structures whose function can be reasonably approximated in two dimensions. 3D models can provide a fuller characterization of the morphology at hand but are much more time-consuming to build. As such, most studies using 3D FE models examine fewer than 10 taxa (Fig. 1B).

Given the steady increase in research studies employing FE methods to address comparative biomechanical and evolutionary questions, it is critical for both practitioners of FE analyses and researchers considering using the FE toolkit to test biomechanical hypotheses to be able to design studies efficiently using this time-intensive method. To highlight the important issue of sample size in 3D FE studies of structure-function relationships, here I ask the question: are time-bottlenecked small-sample 3D FE datasets adequate in reproducing performance-ecology correlations observed in broader taxonomic datasets? Given the already rich literature in using 3D FE analyses in comparative skull biomechanics research, I take a meta-analysis approach to addressing this question using mainly published studies.

The study system I use to demonstrate the effect of taxonomic sample size on performance-ecology relationships is the skull. Although the FE method is applicable to any morphological system whose geometry and biophysical boundary conditions can be digitized and parameterized, most studies in vertebrates employing FE modeling in a comparative context have done so to study skull biomechanical performance (Ross, 2005). Nevertheless, the effects of taxonomic sample size are expected to be shared in part by other study systems. Therefore, the findings from these analyses are expected to be relevant to researchers in comparative biology, paleobiology, bioengineering, and biomedical engineering fields that use multi-taxon comparative FE datasets to test structure-function hypotheses.

Survey Methodology

I tallied the number of taxa included in published FE studies in general (Fig. 1A) and specifically of the vertebrate skull from 2005 to 2020 for 3D-model based analyses (Fig. 1B). Total number of peer-reviewed publications using FE methods were extracted from the Web of Science database (accessed 22 December 2020) using the key words “finite element” + “evolution” + bio*.

A second survey of the literature was conducted on FE studies that specifically address vertebrate skull biomechanics in a comparative context by searches in both Web of Science and Google Scholar (both accessed 15 December 2020) using key words “finite element” + “skull” + phylo*. The year of publication was constrained to between 2005 and 2020. The total number of unique species studied in each surveyed publication was counted. The surveyed publications were further vetted by removing all studies that used 2D FE models (Fig. 1B). Out of the 28 studies obtained from the skull FE survey (Table S1), the FE model construction methodology used in each survey publication was noted and a composite FE model dataset was constructed (see next section).

Results

Baker’s Gamma correlation coefficient between FE-based taxon clusters and feeding ecology-based taxon clusters in the full data set is 0.0498 (bootstrapped 95% CI [0.0019–0.1260]), indicating weak to no association between taxon groupings generated by FE traits vs ecological traits (Fig. 2A). Of the bootstrapped subsamples, datasets with 12 or more taxa returned correlation coefficients within the range calculated for the full dataset. Datasets with 11 taxa or fewer exhibit increasingly large correlation coefficients as sample size decreased. All taxonomic sample sizes at n = 25 or smaller contain at least 1 replicate with correlation coefficient above 0.25, and all taxonomic sample sizes at 13 or fewer contain replicates with correlation coefficients of 0.75 or more.

Baker’s Gamma correlation coefficient between FE-based clusters and phylogenetic grouping in the full dataset is 0.0171 (bootstrapped 95% CI [0.0033–0.2451494]), indicating weak to no association between FE traits and phylogenetic structure (Fig. 2B). Resampled taxonomic datasets with 7 or more taxa returned correlation coefficient values within the range observed in the full dataset. Datasets with 6 or fewer taxa exhibited increasingly large correlation coefficients. Resampled datasets with 22 or fewer taxa contained at least one replicate with correlation coefficient above 0.25; datasets with 10 or fewer taxa contained at least one replicate with correlation coefficient value above 0.75.

Baker’s Gamma correlation coefficient between feeding ecological clusters and phylogenetic grouping in the full dataset is 0.0498 (bootstrapped 95% CI [0.0015–0.1037]), indicating weak to no association between ecological groupings and phylogenetic groupings (Fig. 2C). Resampled datasets with 12 or more taxa returned correlation coefficient values within the range observed in the full dataset. Datasets with 17 or fewer taxa contain at least one replicate with correlation coefficient larger than 0.25; datasets with 8 taxa or fewer contain at least one replicate with correlation coefficient larger than 0.75.

Cophenetic correlation coefficients that consider branch length information from the dendrograms returned broadly similar results to Baker’s Gamma correlation coefficient analyses (Figs. 2D–2F). The main differences are found in (1) the FE to phylogenetic structure correlation, which showed that resampled datasets with 13 taxa and above returned correlations not significantly different from those estimated in the full dataset (compared to 7 taxa and above using Baker’s Gamma coefficient), and (2) feeding ecological to phylogenetic structure correlation, which showed that resampled dataset with 10 taxa and above returned correlations similar to those in the full dataset (compared to 12 taxa and above using Baker’s Gamma coefficient).

Discussion

The use of small (<10 taxa) comparative FE model datasets in the majority of published 3D FEA studies of vertebrate skulls likely results in overestimates of the feeding ecological association of simulated biomechanical traits. Relative to a weak or no correlation “full” dataset of 27 taxa, resampled datasets containing 11 or fewer taxa exhibit significantly elevated correlation coefficient estimates for the relationship between FE value and ecological groupings (Fig. 2). Small taxonomic datasets also exhibit stronger association with phylogenetic groupings in both ecological and FE-based clusters, albeit to different extents. Ecological groupings up to ~10-taxon datasets show significantly higher correlation to phylogenetic groupings, and FE-value based groupings show significantly higher correlation to phylogenetic groupings at between 5 and 12 taxa or fewer (Fig. 2). Therefore, at small taxonomic sample sizes of 3–5 taxa, there is relatively high association to phylogenetic structure in the FE data groupings compared to larger taxonomic samples.

The elevated correlation coefficients in smaller taxonomic samples are in part driven by the prevalence of outliers in the bootstrap replicates. The smaller the taxonomic sample, the higher the quantity of high correlation coefficient replicates. This observation suggests taxon sampling choice could have a substantial effect on the resulting presence/absence of significant FE data to ecological grouping correlations. At the low end of the sampling spectrum, all three-sample datasets are expected to return correlation coefficients of 0.5 or higher, with a coefficient of 1.0 defining the upper quartile. This suggests that the practice of using a small number of taxa to perform FE simulations to interpret the overall performance-ecology association of a larger taxonomic clade runs the risk of finding spurious high correlation results when the underlying full dataset exhibits weak or no correlation. In other words, the high number of outliers and the elevated mean correlation coefficients at taxonomic sample sizes smaller than ~11–12 taxa translate to higher incidences of false positives relative to the full dataset.

Despite the inability for small taxonomic datasets to replicate performance-ecology correlations and phylogenetic structure of larger datasets, these findings do not render small taxonomic sample FE studies obsolete. However, the results do highlight the importance of “calibrating” the research question at hand to the appropriate taxonomic sample to lessen potential biases from false positives. Results from small taxonomic sample FE studies could remain useful if the research question is focused on taxon-specific performance-ecology comparisons, rather than extrapolation to broader clade level correlations. Nevertheless, the simultaneous elevated correlation in performance-ecology, performance-phylogeny, and ecology-phylogeny relationships invites caution in interpreting functional correlation when a phylogenetic one is equally plausible based on available data.

To further assess whether the significant correlation of FE and feeding ecological trait data with phylogenetic structure represents elevated phylogenetic signal at smaller taxonomic samples, I conducted post-hoc analyses to estimate phylogenetic signal in the FE and feeding ecology data. I used a multivariate implementation of Blomberg’s K (Adams, 2014) to calculate phylogenetic signal in FE values in the full dataset as well as a similar bootstrap series of smaller resampled datasets (Fig. 3A). Results indicate that there is a concomitant increase in K with decreasing taxonomic sample size, as has been previously observed for simulated datasets (Münkemüller et al., 2012). Increasing K and increasing FE data-feeding ecology correlation with smaller taxon samples suggest the presence of confounding ecological and phylogenetic factors.

I also estimated phylogenetic signal in the categorical ecological data using the delta statistic (Borges et al., 2019). In contrast to the elevated phylogenetic signal in FE data at small taxon sample sizes, no elevations in phylogenetic signal are observed in either diet breadth or trophic level categorical data (Figs. 3B and 3C). However, the delta statistic is known to exhibit low sensitivity in detecting phylogenetic signal at small taxon sample sizes (<20) (Borges et al., 2019), so it is unclear whether the feeding ecology-phylogenetic structure correlation (Figs. 2C and 2F) is unrelated to phylogenetic signal or whether phylogenetic signal is undetected by current methods. Comparative methods have been shown to be affected by high variance and low power in estimating phylogenetic signal and other parameters at small taxon samples (Boettiger, Coop & Ralph, 2012); the elevated FE data, ecology, and phylogeny correlations in the bootstrapped samples of the current study (Fig. 2) appear to be similarly affected by uncertainties at small taxon sample sizes.

Although the bootstrap analysis involves random resampling and it is not possible to pinpoint which specific taxa or phylogenetic factors produce outsized effects on FE-ecology correlations at small taxonomic sample sizes, some general recommendations for small-sample comparative FE studies can still be made. The categorical groupings in the ecological trait data were converted into a continuous distance framework in the cluster analyses via Gower distance measures (see “Methods” section). The result is a relatively clumped dendrogram with polytomies representing ecologically ‘identical’ taxa in the context of the input variables (diet breadth and trophic level). The overall dataset (n = 27) contains multiple samples in each ecological grouping or clump, providing more evenly represented samples in each group than is possible in smaller randomized resampled datasets. As such, a recommended sampling strategy for smaller taxonomic datasets might be to focus on maximizing the even sampling of taxa representing unique ecological trait combinations and to avoid asymmetric sampling of some ecologically similar taxa over others. Such a sampling strategy should reduce spurious high-correlation outcomes when using distance-based comparison methods to study FE-ecology relationships in smaller datasets.

For comparative analyses using FE data to study broader clade-level associations between biomechanical and ecological traits, the time-consuming nature of constructing 3D FE models remains a bottleneck to achieve comparable scope to other data sources such as sequence or geometric morphometric shape data. Given the sensitivity of the correlation coefficient between FE values and ecological categories on choice in taxonomic sampling and sample sizes, I suggest some rethinking in future comparative FE-based biomechanics research design. Rather than focusing on collection of comparative FE data in increasingly large samples of taxa, I posit that most comparative biomechanics research using the FE method would be better served with a theoretical morphology approach (Polly et al., 2016).

The fusion of morphospace analysis (using methods such as morphometric morphometrics or other multivariate trait data) and FE analysis has already been demonstrated to be a fruitful approach to test hypotheses about evolutionary optimality and the relative importance of multiple selective forces in explaining morphological disparity (Stayton, 2009; Tseng, 2013; Dumont et al., 2014; Polly et al., 2016). Theoretical morphological models constructed at extremes and/or regular intervals over a given morphospace reduces the subjective nature of taxon selection during FE analyses by summarizing the range of morphological variation that underlies the subsequent FE simulation outcomes, rather than relying on taxonomic-specific interpretations that may be more sensitive to outlier and sampling effects. In this morphospace-driven context, the biomechanical performance context of the taxonomic dataset at hand is only indirectly dependent on the choice of taxon sampling, assuming morphospace sampling is representative of the morphological disparity in the clade studied. A shift to comparative FE analyses based on a morphospace framework leverages the comparative power of the method, especially when the practice of experimental model validation has yet to become standard practice to permit the evaluation of absolute magnitudes of FE simulation outcomes (e.g., Strait et al., 2005; Bright & Rayfield, 2011).

Conclusion

Time as a limiting factor in applications of 3D FE simulations in comparative biomechanics research has a direct effect on limiting the taxonomic breadth of biomechanics research using comparative FE analysis. A consequence of this limitation is the presence of significant biases in performance-ecology correlation coefficients driven by small taxonomic sample sizes and outliers. Future advances in comparative biomechanics research using FE modeling may depend on a bifurcation in the application of this method. On the one hand, small taxonomic sample studies can remain useful for interpreting taxon-specific biomechanical adaptations, if carefully designed with consideration of phylogenetic structure, ecological trait representation among taxa, and preferably integrated with model validation. On the other hand, research efforts in quantifying clade-level form-function associations could be better served through theoretical morphological approaches of representing morphological disparity, rather than building increasingly large FE datasets of taxon-specific values that are more vulnerable to sampling outlier effects even at larger sample sizes. Continued improvements in model construction efficiency and accuracy are key to solving the time bottleneck issue in using FE methods in broad comparative studies. Finite element analysis is a once chic biomechanical modeling method in comparative biology that has come of age, and continued methodological and application development should help to build its analytical rigor to be on par with any other comparative methodology.

Supplemental Information