The intervals method: a new approach to analyse finite element outputs using multivariate statistics

FEA models

Plane models of 11 mandibles of Cingulata (Table 1), each one corresponding to a different species, were created according to the methodology summarized by Fortuny et al. (2012). The models were created using the ANSYS FEA Package (Ansys Inc.) v.15 for Windows 7 (64-bit system) to obtain the von Mises stress distribution.

Two main masticatory muscles (i.e., temporalis and masseter) were included in the model as a vector between the centroid of the muscular attachment in the mandible and the centroid of the equivalent muscle attachment in the skull following the modelling approach used in Serrano-Fochs et al. (2015). To compare the models, a scaling of the values of the forces was applied according to a quasi-homothetic transformation in the FEA models (Marcé-Nogué et al., 2013) using the plane model of Chaetophractus villosus as a reference. This method corresponds to an adaptation of the scaling methods proposed by Wroe, McHenry & Thomason (2005) and Dumont, Grosse & Slater (2009) for plane models. This procedure was performed to apply the appropriate force in each model, thus allowing the comparison of the stress results when the specimens differ in size.

The information for each analysed species regarding the area of the mandible, insertion areas, forces (musculature applied force per unit area (N/mm²)), thickness and the scale factor in the quasi-homothetic transformation can be found in Table 1.

The boundary conditions were defined and placed to represent the loads, displacements, and constraining anchors that the structure (i.e., mandible) experiences during its function. The mandible was constrained in the x and y direction at the most anterior part and fixed in the x and y directions on the condyle at the level of the mandibular notch (Fig. 1) following the procedures described in Serrano-Fochs et al. (2015) and Marcé-Nogué et al. (2016).

Isotropic and linear elastic properties were assumed for the bone. In the absence of data for Cingulata or any other closer relative, as well as lacking data for any mammalian clade with a similarly shaped jaw, we decided to apply the mandibular material properties of Macaca rhesus: E (Elasticity Modulus) = 21,000 MPa and v (Poisson coefficient) = 0.45 (Dechow & Hylander, 2000). We chose the available properties of Macaca rhesus because it has a wide range of habitats and diet which resembles omnivorous or generalist insectivorous armadillos (Richard, Goldstein & Dewar, 1989). In addition, it has been shown that in a comparative analysis these values are not crucial (See Gil, Marcé-Nogué & Sánchez, 2015 for discussion).

As primary data, we obtained the von Mises stress distribution of each one of the analysed species. Von Mises stress is an isotropic criterion used to predict the yielding of ductile materials determining an equivalent state of stress (Reddy, 2008). Considering bone as a ductile material (Dumont, Grosse & Slater, 2009) and according to Doblaré, García & Gómez (2004) when isotropic material properties are defined in cortical bone, the von Mises criterion is the most adequate for comparing stress states.

Quantitative stress data

Von Mises stress values were obtained for all the elements of each finite element model. Values of stress (and strain) can be obtained for each element because FEA mathematically solves stress (and strain) inside each element, while displacements are computed at the nodes (Zienkiewicz, Taylor & Zhu, 2013). ANSYS provides a complete data file with these results that can be easily manipulated in spreadsheets and other software. This means that the number of stress data points for each mandible is the same as the number of elements of the mesh.

The new methodology proposed here divides the values of stress into N equal intervals (Fig. 1). Therefore, each interval will contain all the elements within a certain threshold, each one being defined by a lower threshold T_lower and an upper threshold T_upper. The range of the intervals is constant across the sample, the distance from T_lower to T_upper being the same in all the intervals. For an interval Φ_i the lower threshold coincides with the upper threshold of the interval Φ_i−1 and the upper threshold coincides with the lower threshold of the interval Φ_i+1.

Once all the stress values of a single specimen were obtained they were subdivided into different intervals. When all the elements of this specimen were allocated into an interval, the total number of elements in each interval EΦ_i (Eq. (1)) was computed and the total area of each interval AΦ_i (Eq. (2)) was calculated from the individual area of each element. These values must fulfil the Eqs. (1) and (2). Then, the percentage area of each interval was computed in relation with the total area of the model for each specimen, following the Eq. (3). (1)${\displaystyle \sum _{i=1}^{N}E{\Phi }_i=E=\text{Total number of Elements in the mesh}}$ (2)${\displaystyle \sum _{i=1}^{N}A{\Phi }_i=A=\text{Total Area of the Model}}$ (3)${\displaystyle A\Phi \left[\text{\%}\right]=\frac{A{\Phi }_i}{\sum _i^{N}A{\Phi }_i}=\frac{A{\Phi }_i}{A}.}$ After carrying out this procedure, a new set of variables was generated (AΦ [%]), each one representing a different interval of stress values. Each variable represents the amount of area (as a percentage) of the original model having this range of stress values. If this procedure is performed in several models, these new variables can be used in later analyses comparing the different models.

A Fixed Upper Threshold FT_upper has to be defined in the interval Φ_N−1 to allow including in the interval Φ_N all the highest values of stress. Some of these high values of stress represent artificial noise, which are numerical singularities produced in the results of FEA at the points where the displacement boundary conditions are applied (Marcé-Nogué et al., 2015b). The presence of this artificial noise in a numerical model is due to the mesh, since as the elements get smaller, the artificial high stress values become even higher (Dumont, Grosse & Slater, 2009). To select the FT_upper we chose not to consider the maximum value of 2% of the higher values of the model based in the suggestions of Walmsley et al. (2013) and Marcé-Nogué et al. (2016).

Convergence of the data for the case study

To enable the comparison between different species it is necessary to specify the number of intervals that will be applied to the specimens under analysis, since this number will define the number of variables (A[%]). A large number of intervals will generate an excessive number of variables to work with, which can be difficult to interpret in later analyses. On the other hand, a limited number of intervals may yield incorrect results as they could oversimplify the stress pattern. Therefore, we decided to carry out multiple Principal Components Analyses (PCAs) as a guide to choose the minimum number of the intervals.

We performed these PCAs based on a different number of intervals to establish the threshold in which the data converged (i.e., when adding more intervals yielded similar patterns in the PCAs). We named this threshold value N_PCA.

For the 11 mandible models of our case study (Fig. 2), the PCAs were performed using nine different sets of intervals: N = 3, 5, 10, 15, 20, 25, 50, 75 and 100. In this case the FT_upper was fixed to 0.1 MPa, which represents the appropriate value where less than the 2% of the higher values for all the eleven models are in its interval Φ_N. Once all the PCAs were carried out we performed Major Axis Linear regressions between equivalent PC scores of one set of intervals and the next one (e.g., regression between the PC1 scores of the PCA of 3 intervals with PC1 scores of the PCA of 5 intervals). The coefficient of determination (R²) was then used to assess the convergence of the obtained results (Document S5). For each PC (1, 2 and 3 only, as they represent more than 80% of the variance for all the cases) we computed the coefficient of determination for successive pairs of PCs. The number of intervals at which R² reached a plateau (i.e., it does not further increase) was considered the N_PCA.

When the variables are all in the same units, it is usually preferred to carry out the PCA on the variance–covariance matrix. Since our variables were in the same units, we performed the aforementioned PCAs using the variance–covariance matrix. Nevertheless, if some variables, despite being in the same units, have a noticeably larger variance than others they will have a higher weight on the PCA, which might obscure more subtle relationships between those variables. In these kind of situations PCAs using the correlation matrix are more adequate (i.e., this is equivalent to perform the PCA using standardized variables). In the context of the present study this is relevant, since areas with high stress are usually small in all specimens. This implies that the variance for this specific interval (which represents one of the variables used to perform the PCA) will be small compared to the variance of an interval representing lower stress values, which occupy very large areas. Since the interval representing higher stresses could be still informative in comparative terms between species, it is important to account for this situation. By performing the PCA on the correlation matrix, a variable that has a smaller variance in absolute terms can be informative in relative terms. For this reason, in the case study we carried out PCAs using both kinds of variable matrices (i.e., variance–covariance and correlation).

Validation of the data for statistical analyses

Document S4 show the percentage of area for each interval AΦ[%] using different N intervals, as well as the variance loadings for the PCAs.

Based on the patterns observed in the plots displaying the first two PCs (Fig. 3) for each set of intervals it is clear that the distribution in stress space became more or less the same after N = 50. Additionally, the results of the regressions of the PCs indicate that the scores are almost completely correlated for N = 50 (Fig. 4 and Document S5), therefore we chose to compare the different species with the set of 50 intervals (i.e., N_PCA).

Case study

Using a N_PCA = 50, the PCA developed on the covariance matrix showed more than 90% of the variance concentrated in the first PC, while using the correlation matrix the first PC explained only 75.4% of the variance.

The bi-plot of stress space defined for the first two PCs of the PCA performed on the variance–covariance matrix allow to identify different areas (Fig. 5A). The right side of this bi-plot (Fig. 5A) is occupied by species that can be described as having a very large area of a specific range of values (i.e., low to very low stress), whilst downwards on the left of this stress space are specimens characterized by moderate stress values. Finally, the upper left part of this stress space is exclusively occupied by individuals showing values intermediate between the former ones. Based on the loadings of the variables, is evident that the proportionally larger areas with low stress dominate this PCA (Fig. 5B).

The results obtained from this PCA make it possible to establish certain patterns. A proportionally larger area of intermediate stresses characterizes the insectivores, while proportionally larger areas with very low stress values characterize the omnivores (Fig. 5A). However, within this trend two species showed a distinct pattern in PC2: T. matacus and Z. pichiy. Both are located in the upper part of the plot (Fig. 5A), characterized mainly by one specific value of low stress area. Nonetheless, the second PC only explains 16.8% of the variance.

According to a hypothesis suggested in Serrano-Fochs et al. (2015), high levels of stress along the mandible would represent a fragile mandible with a reduced capability to chew or process hard items. Supporting this hypothesis, in the present case the insectivores are in the left part of the stress space, defined by intermediate to low stress levels. This can be expected considering that these species feed mainly on ants without chewing them; thus their mandibles are expected to show a higher stress level. On the other hand, very low stress values would be indicative of mandible with a higher ability to bear higher loads. This would be in agreement with the location of the omnivorous taxa at the right side of the stress space characterized by intermediate stress values (E. sexcintus, Ch. villosus, Ch. vellerosus and Z. pichiy, Fig. 5A).

Nevertheless, this PCA did not detect any difference in the stress pattern between specialist and generalist insectivores. All the loadings for the first two PCs belong to low or intermediate stress values, whilst the loadings of the variables representing intermediate to high stress have almost negligible values. This arises because the percentage of areas with high values is small in all analysed mandibles (the variables representing intermediate and high stresses have very small loadings for PC1 and PC2). Since we carried out the PCA using the variance–covariance matrix, the larger absolute values of the lower stress areas implied larger variance for those variables, thus completely masking the variability that might exists in the variables representing high stress values.

In the case of the PCA using the correlation matrix the distribution of the loadings of the variables (Fig. 5B) is noticeably more homogeneous (more variables have an important contribution to the PCs analysed) without being so strongly influenced by just a few specific variables as was the case of the PCA performed using the variance–covariance matrix (where variables representing low stress values overcome the rest of the variance). The PCA of the correlation matrix successfully distinguishes the main three diets, with omnivore-carnivore species on lower-left area of the plot. PC2 separates specialist insectivores from the rest of species, while Chlamyphorus truncatus (i.e., a generalist insectivorous species exhibiting a very particular diet due to its completely fossorial life style) is located near the omnivore/carnivore species in the negative part of PC1 and PC2.

Within this PCA, omnivore/carnivore species were characterized by very low stress values (lower-left quadrant of the plot), whilst generalist insectivores showed a proportionally larger area of intermediate stress values than the rest of the species. Specialist insectivores have proportionally larger areas of high stress.

As was the case when performing the PCA using the variance–covariance matrix, the fossorial generalist insectivore—Chlamyphorus truncatus was located with the omnivores. In spite of the fact that insects correspond to large part of its diet, this species also feeds on worms, snails, and a small proportion of roots (Borghi et al., 2011). This consumption of relatively harder objects that need more oral processing than just tiny insects, might explain the biomechanical results obtained here.

As similarly found by the variance–covariance PCA, Z. pichiy corresponds to the omnivorous species with proportionally less area of low stress. This species is considered to show a preference for soft-bodied species (e.g., larvae, tarantula spiders, pods) (Redford, 1985), thus being the least carnivorous taxon among all the omnivore/carnivore armadillos (Superina, Pagnutti & Abba, 2014). Indeed, microwear studies classified Z. pichiy together with the insectivorous species (Green, 2009). In concordance with this information, our results show that Z. pichiy is located in the areas of the stress space between carnivore and insectivore species. Finally, it has been suggested that T. matacus has a seasonal diet based on fruits and pods representing more than 50% of its consumed items during the dry season (Bolkovi, Caziani & Protomastro, 1995), being the the generalist insectivores species showing the lowest stress.

Influence of the mesh in the results

The methodology presented in this work is based on obtaining the percentage areas containing certain values of stress. To validate the results, the influence of coarse meshes and non-homogeneous meshes on the computed areas was also tested in the Supplementary Information, where methods and results are discussed (Document S1).

According to the results obtained (Figs. S5A and S5B and Tables S1 and S2, our proposed approach is not affected by the size of the elements (i.e., we obtained extremely low values of the relative error ($REr{r}_i^j$)). In addition, when the homogeneity of the mesh was tested, the absolute error (AErr^j) and its variance showed very low values as well Figs. S5C and S5D and Tables S3 and S4, hence confirming that the homogeneity of the mesh is not affecting the newly generated variables.

Overall, this means that under our proposed methodology the size and the homogeneity of the mesh do not affect the obtained results, thus providing a robust method that is effective irrespective of these factors. Therefore, it is not necessary to apply mesh correction techniques such as generating a Quasi-Ideal Mesh (QIM) as proposed by Marcé-Nogué et al. (2016).

Convergence of the data in the case study; which N is adequate?

There were important differences between the PCAs depending on the N value. When this value was low (i.e., N = 3 and 5) the obtained results were noticeably different, whereas from N = 15 and above the PCA results were very similar, as shown by the dispersion plots of the first and second PCs (Fig. 4).

In the case study, the PCA convergence was N_PCA = 50. This means that at least for structures such as armadillo jaws, 50 intervals are enough to describe the sample variability in terms of distribution of stress values per area. In fact, for the purpose of this analysis (i.e., to distinguish between dietary groups) around 15 intervals would have been enough to describe at least the main eco-morphological categories.

Whether N = 50 is an adequate number of intervals for other cases is something that needs to be addressed in the future and that might depend on the structure and sample under analysis. It is likely that for more complex geometries or if species from other taxonomic groups are included, the N_PCA will differ. Based on our present results, we suggest that the N_PCA should be determined on a case-by-case basis.

Comparative biomechanical performance of extant armadillos

Using a N_PCA = 50, the PCA of the variance–covariance matrix showed differences between omnivores and both groups of insectivores (Fig. 5A). These observed differences are concordant with previously published results (Serrano-Fochs et al., 2015; Marcé-Nogué et al., 2016), although the present results provide a better distinction of the groups when compared to similar study that only analysed the stress values collected at specific points (Serrano-Fochs et al., 2015) (Fig. 5A). The mandibular ecomorphology of Cingulata has been previously studied using morphometrics, microwear and biomechanics, among other techniques (Fariña, 1985; Fariña, 1988; Vizcaíno & Fariña, 1997; Vizcaíno & Bargo, 1998; Vizcaíno, De Iuliis & Bargo, 1998; De Iuliis, Bargo & Vizcaíno, 2000; Fariña & Vizcaíno, 2001; Vizcaíno et al., 2004; Vizcaíno, 2009; De Esteban-Trivigno, 2011). Nevertheless, the results obtained applying our proposed approach are more accurate when distinguishing between dietary groups. It is noteworthy that none of the previously proposed methods to compare mandibular performance using FEA applied to different species of armadillos (Serrano-Fochs et al., 2015; Marcé-Nogué et al., 2016) were able to effectively distinguish the specialist from the generalist insectivores. Consequently, the approach presented here, especially when using the correlation matrix in the PCA, seems to be able to detect subtler differences, thus representing a useful way to characterize and understand the biomechanical behaviour of the mandibles in relation to diet.

The bi-plot of the first two PCs of the PCA computed using the correlation matrix displayed T. matacus and Z. pichiy in a location that was slightly different than the rest of species of their corresponding dietary group. Classifying diets is always a complicated matter, as it is reducing the existing dietary variability to only a few categories. As explained in the results section, these two species have a diet that differs from the general trend of their groups. In a similar fashion, Chl. truncatus occupies a position along the omnivore/carnivore group instead of showing stress values similar to other insectivores. As previously mentioned, even though this species has been usually grouped with insectivores, Chl. truncatus feeds on worms, snails and roots; all items that need more processing if compared to simply swallowing small insects. It is also relevant that the PCAs of the stress intervals were able to detect these differences in diet, including some subtle differences that have been elusive using other techniques. Even though further studies are required, the proposed procedure seems to be a promising approach to be tested in future ecomorphological studies considering other taxonomic groups.

The PCA computed using the correlation matrix provides information on the variability of all the variables, giving them all the same importance through standardization. This explains why it was possible to distinguish between the two insectivore’s diets even if their differences were rather subtle (Fig. 5B). Nonetheless, interpreting the loadings from the correlation-based PCA can be misleading if it is not properly done. The loadings of the PCA performed using the variance–covariance matrix can be easily interpreted; in the present study, a variable (AΦ[%]) with a higher loading for a given PC represents a larger actual area for that interval of stress values in the mandibles of the specimens located in the direction of that loading. However, in the PCA carried out using the correlation matrix (i.e., where the variables have been standardized), a variable with a higher loading indicates that this variable has a proportionally larger variability in that area of the stress space (but not necessarily in absolute terms). For example, in the Fig. 5B C. unicintus is located in region of the stress space that corresponds to higher loads for high and very high stress variables. This does not mean that the species in question has an area with high stresses larger than its area with small stresses in absolute terms. Rather, it indicates that these high stress variables are proportionally larger as compared with the rest of species. It is important to bear this in mind when using correlation matrix-based PCAs to interpret the results obtained from the intervals method proposed here. As discussed above, overall both PCAs (i.e., using the variance–covariance or correlations matrices) have shown to be very useful for identifying ecomorphological patterns based on stress values obtained from FEA.

We have successfully reduced thousands of data points to just 50 variables that are still biologically meaningful. In addition, we managed to use the output data from FEA, which has a different number of data for each specimen (i.e., due to the differences in the number of elements in each model) and represent it with the same 50 variables in each specimen under analysis making them comparable in a statistical framework. We have shown that at least for our case study this approach works better than previously published examples which are based on either the visual inspection of the stress maps, collecting stress data on specific points, or comparing the complete stress value distributions (Serrano-Fochs et al., 2015; Marcé-Nogué et al., 2016). We consider that for comparative analysis where the main aim is to test the same functional hypothesis under equivalent loads in different species, this method can be generalized to other taxonomic groups and structures. Additionally, it is promising that in addition to the PCAs other techniques might be applied using this newly generated stress variables (intervals) depending on the sample and the biological hypotheses being tested (a matter that should be explored in future studies).

Although the proposed approach implies losing the spatial information of the stress, the case study analysed here strongly supports that the new interval variables are still biologically meaningful and easily interpretable, thus being very useful to describe and analyse the data. Our proposed approach is especially appropriate in cases where the amount of stress might vary but the stress spatial distribution is highly similar. In these situations, it is very difficult to interpret the obtained results, as well as to answer the tested biological hypothesis just by visually inspecting the coloured stress maps. It is not clear yet whether this will work or not in cases where the analysed structures are very different (e.g., when comparing dissimilar morphologies at a higher taxonomic level where usually anatomical differences tend to be prominent). For this reason, we recommend that the newly proposed intervals method be used together with other traditional approaches, such as the visual inspection of the stress and/or deformation maps.

Further research is required to test the advantages and limitations of the intervals method when analysing other anatomical structures and/or taxonomic groups, as well as testing it using three-dimensional models (where the same method can be applied just by simply using volume instead of area in the equations presented here).