A new biomechanical approach to cranial suture function: the role of contact elements in linear and nonlinear models

Abstract

Understanding cranial sutures and how they relieve and dissipate stress is essential to assess their role in cranial biomechanics and to develop highly accurate predictive models. This involves examining how ontogeny affects cranial sutures, as well as their morphology and function, and how these changes through time may impact essential biomechanical loadings such as chewing or direct biting. In this work, we study the cranial sutures of Crocodylus niloticus in detail using contact elements under finite element analysis. Contact elements permit the creation of a physical relationship between two bones that are in contact and even the configuration of these relationships, for example, in terms of movement or flexibility. The definition of bone contacts may require linear and/or nonlinear computational solutions to attain higher accuracy. Herein, skull geometry is tested to determine how bones may be altered by different types of contacts under various conditions. As predicted, the absence of sutures or cranial kinesis leads to a reduction in stress distribution across the skull, whereas sutures and cranial kinesis help the skull relieve stress and prevent certain bones from sustaining high stress levels. The type of contact used in individual sutures has a significant effect on model outcomes. Additionally, feeding behaviors significantly impact cranial biomechanics, reflecting the influence of other variables that may be applied to the models. As highlighted by the results, in order to obtain accurate results when analyzing fossil taxa, the nature of the cranial sutures should be taken into account. Therefore, developing predictive models based on living taxa is invaluable, because it facilitates the study of extinct taxa for which there is a lack of information on the fibrous joints due to poor or no preservation in the fossil record.

Non-technical Summary

This study proposes a new computational method using a 3D digital model to understand how the zone of contact between cranial bones (namely sutures) works mechanically. In cranial mechanics, only movements in which the bones remain in constant contact are allowed. Therefore, we reproduced all these movements in four cranial models using different contact properties to validate the new methodology. The results obtained suggest that first, when the skull allows less movement between bones, stress is concentrated on certain points of the skull; and second, when the skull allows more movement between cranial bones, the stress tends to be dispersed in other bones, protecting the skull's physical integrity. Finally, we used these computational models to reproduce different types of predatory feeding behavior observed in animals such as crocodiles and alligators. We conclude that the new method to model the contact between bones can be applied to fossils of extinct animals that do not preserve soft tissue as cranial sutures.

Introduction

The role of sutures in cranial biomechanics of extinct and extant vertebrates has been a long-debated subject among functional and evolutionary biologists and paleontologists. Cranial sutures are boundary areas between bones, comprising a soft tissue component and the contacting bone edges (Curtis et al. 2013). Sutures are, however, not only sites of bone deposition during growth, but also determine the biomechanics of the skull (Bailleul et al. 2016). Essentially, they serve as critical elements in the dissipation and reduction of stress (Curtis et al. 2013).

The morphology of cranial sutures strongly influences the forces experienced during different passive (e.g., resisting impacts) or active (e.g., feeding) activities of the vertebrate skull (Curtis et al. 2013). Nevertheless, cranial sutures are not equal across skulls and could widely vary externally and internally from one skull to another (Kathe 1995, 1999). The different roles of these structures through ontogeny and how they evolved across the fossil record remain unclear (Brochu 1996).

From younger to older skulls, sutures are progressively replaced by bone, facilitating alterations in length, width, and shape of the head (White et al. 2021). The ossification of cranial sutures also modifies their biomechanical properties, influencing the forces experienced and transmitted by the craniofacial complex during growth (Sharp et al. 2023). Consequently, variations in the characteristics of sutures may result in disparate biomechanical properties and outputs in cranial models. In the fossil record, the tissue filling the sutures is no longer present, and instead, the morphology of the union between bones is the only evidence remaining. Nevertheless, recent developments in digital imaging and modeling techniques have enabled the creation of versatile tools for the reconstruction of various soft-tissue structures (Lautenschlager 2016), but many times it is not possible to recognize suture soft tissue (or it is not preserved) in the computed tomography (CT) data or the restoration in deformed and disarticulated specimens makes it impossible to apply these techniques.

In this work, we introduce an alternative, novel methodology for suture modeling based on the morphology of the junction between bones, as this is more likely to be preserved, that can be applied in instances where existing techniques cannot.

Finite element models have shown that cranial sutures can influence biomechanical results such as stress and strain in the skull of mammals, reptiles, and birds (Kupczik et al. 2007; Bright 2012; Curtis et al. 2013; Moazen et al. 2013; Cuff et al. 2015). Based on current understanding sutures exhibit nonlinear material properties that cause them to behave differently in tension and compression (Popowics and Herring 2007). This is not addressed in current works that model sutures by means of computational methods whereby linear and static approximations are employed due to their computational efficiency and the associated errors are often tolerable. However, in biomechanical models involving soft tissues, the incorporation of nonlinearities needs to be considered, although the inclusion of this would imply an increase in the computational cost because of the increased complexity and the need to use iterative solvers (Marcé-Nogué 2022). Consequently, the incorrect representation of sutures as soft-tissue structures could lead to different results by over- or underestimating stresses in the studied specimens.

Furthermore, studies of extinct species that want to include sutures must face the problem of the lack of information about the sutures due to nonpreservation. Therefore, a new methodology is needed to enable the inclusion of sutures without including them. Finite element analysis (FEA) provides the opportunity to add contact elements. This permits the establishment of a physical relationship between two surfaces that are in close proximity, as well as the configuration of these relationships in terms of separation, stiffness, or damping, among other mathematical formulations (Wriggers 2002). In other words, to apply contact elements in cranial structures, only the behavior between adjacent bones needs to be defined, instead of having to create new bodies that represent suture soft tissues. Therefore, the proposed approach is especially appropriate for models lacking information on suture soft tissue, as is common in the fossil record, where only the suture morphology is preserved.

This work aims to demonstrate that contact elements in FEA models can improve the modeling of the cranial sutures, resulting in higher accuracy results. However, it is very difficult to elucidate how the stresses acted in fossil taxa; thus, it is important to first strengthen the studies in living taxa for their paleobiological implications. Suture function(s) are generally studied in in vivo specimens, as this is the most direct evidence to understand the forces related to sutures (Rayfield 2007). In this work, the living taxon Crocodylus niloticus is the case study to test different configurations of contact definitions. Each configuration is discussed in detail and reveals ontogenetic implications for cranial suture closure. It is a prior required step to later reproduce the methodology in extinct taxa.

Materials and Methods

Digital Model

The skull of a well preserved Crocodylus niloticus was digitally scanned on an industrial CT scanner (YXLON Y.TU450.D09) at the Institut Català de Paleontologia Miquel Crusafont (Sabadell, Spain). The original specimen, corresponding to an adult male, is housed at the Museu de Ciències Naturals de Barcelona, under the label MZB 2003-1423. As a result of the CT scanning, 500 slices were obtained with a pixel size of 567 μm. Raw CT data and derived 3D models of each cranial bone are available via Morphosource for research purposes (https://www.morphosource.org/projects/000600962?locale=en).

The tomographic raw data were imported into VSG Avizo software (Thermo Fisher Scientific) to obtain two different 3D digital models: (1) a complete 3D model of the whole skull, with all the bones merged to generate a unique geometry; and (2) a 3D model of the skull comprising the different bones (individually segmented and separated from the rest) (Fig. 1). Irregularities in the surface that appeared due to the generation of the models following reconstruction and segmentation were repaired using refinement and smoothing tools from Geomagic Wrap 2021 (3D Systems) and converted to a CAD format (Marcé-Nogué et al. 2011).

Figure 1. Types of sutures observed on Crocodylus niloticus (MZB 2003-1423) based on computed tomography (CT) scan. Left, Dorsal, ventral, and lateral views of the analyzed specimen with each bone digitally segmented. Right, Suture classification based on Kathe (1999). A, Vertical walls based on frontal–postorbital suture. B, Shelves based on maxilla–nasal suture. C, Interdigitated based on lacrimal–maxilla suture.

Models with Different Contact Definition

In a finite element model, contacts are defined according to the relationship of movement that is allowed when two separated surfaces of different bodies touch each other (Marcé-Nogué 2022). This relationship of movement is defined in two directions according to how bones can move perpendicularly (separate) and how bones can move in the tangential plane (slide) with respect to the other body. When bones in contact are allowed to separate, the solving process requires an iterative nonlinear process to reach the solution. To summarize, a linear contact can be defined as a bonded contact when separation and sliding are not allowed and as a no-separation contact when separation is not allowed but sliding in the tangential plane is allowed. Otherwise, nonlinear contacts include a frictionless contact, when separation and sliding are allowed, and rough contact, when separation is allowed but sliding in the tangential plane is not allowed.

As the loads transmit differently depending on the type of contact, four models were created with different contact types between bones to study their performance.

Model A: The presence of contact elements was omitted by using a unique solid model without bone separation. This model constitutes a comparative base model against which all other models are compared to test the influence of cranial sutures.

Model B: A linear model with all the bones of the skull “stuck” together using bonded contacts that allow neither normal nor tangential movements between contacting bones.

Model C: A linear model in which the relationship between all the bones in contact was defined as “no separation” using contacts that allow tangential movement between bones in contact but do not enable normal movement.

Model D: The relationship between all the adjacent contacting bones was defined as a function of the type of cranial sutures. It represents a model with well-defined sutures, as it combines a mix of the preceding linear contacts as well as nonlinear contacts. As a starting point, we would expect this model to represent the closest to reality for each cranial suture, as each cranial suture is referred to as one contact type depending on its inner suture morphology examined using CT scanner–generated images. In the case of contacts defined as rough and frictionless, the problem is nonlinear, implying a certain number of iterations to reach the solution.

In model D, nonlinear contacts can be adjusted via the stabilization damping factor and the normal stiffness factor. The stabilization damping factor provides resistance to damp the motion between the contacting surfaces and prevents rigid body motion. We used 0.1 in model D to ensure the convergence of the solution. The normal stiffness factor controls the amount of penetration between surfaces in contact. A value of 0.1 is usually appropriate when bending dominates.

Sutures of Model D

Cranial sutures represent different external and internal morphologies with potentially different functions and biomechanical implications. Herein, regarding model D, cranial sutures have been categorized following Kathe (1995, 1999). Up to eight different cranial suture morphologies have been recognized in vertebrates (Kathe 1995), but all types are not present in all taxa. In this case, the types of recognized sutures were divided into vertical walls, shelves, and interdigitated, depending on the kind of contacts present between adjacent bones (Fig. 1).

For these reasons, the raw CT data (i.e., the tomographic slices) of the analyzed C. niloticus specimen were carefully examined by LG-E using VSG Avizo to identify the different types of sutures previously described. Each bone contact was studied individually, and depending on the morphology of each structure, the type of suture corresponding to each bone was determined (Fig. 1, Table 1). It should be remarked that the same bone could present different types of sutures, as it may be in contact with multiple different bones. In the computational model, each kind of cranial suture has been modeled using a different type of contact (Table 1). The criteria used to correlate each suture type to the corresponding type of contact are as follows: Sutures that have an interdigitated morphology that showed possibility of sliding, hence, no tangential contact, were proposed as a no-separation contact. Sutures presenting the shelves morphology that showed uneven surfaces that did not allow for sliding movements, opposite movements to the no-separation contact, were proposed as rough sutures. Finally, the vertical walls suture type that does not restrain any type of movement was proposed as a frictionless contact. It should be noted that no suture morphology resembling either bonded contacts or frictional contacts with partial friction was recognized in the analyzed taxon, and therefore none of them were used in any bone contact in Model D.

Table 1. Classification of all the sutures observed in Crocodylus niloticus. Contacts between bones (connection) are divided into three different types of sutures (ST): interdigitated, shelves, and vertical walls; contacts (C) are denoted as no separation, rough, and frictionless in the finite element analysis (FEA) model.

¹These cells should be in orange and the text in white - using the same colors as in the uppermost cell (with the text Interdigitated (ST)–no separation (C) ) as they represent different types of contacts.

Finite Element Model

A structural static analysis was performed using the finite element package ANSYS 2021 R2 on a HP Z6 G4 Workstation with 96 GB (8 cores × 12 GB) and 1.90 GHz. Isotropic and homogeneous elastic properties were assumed for the bones that formed the skull in all the models (E = 6.65 GPa and v = 0.35) (Fortuny et al. 2016). Although the specific mechanical properties of bone for the investigated individuals are unknown, the use of general values for bone does not alter results in a relative comparison between models (Gil et al. 2015), and models with heterogeneous mechanical properties of bone closely matched models assuming homogeneous properties (Strait et al. 2005).

Three loading scenarios were defined using an extrinsic approach that has been applied previously using the three behavioral categories employed by crocodilians to catch, subdue, and process prey (McHenry et al. 2006), namely: a bilateral bite, an axial twist, and a lateral bite. The bilateral case was created using two parallel forces of 100 N among the principal tooth-row into the −y-axis direction. The jaw joint was fixed in the y-axis, whereas the condyle was fixed in all directions to prevent free body rotation. The axial twist was created using two parallel forces of 100 N in the same position as for bilateral biting, but one of the forces is facing in the upward direction. The jaw joint on the lateral of the up direction is suppressed, and the jaw joint on the lateral with a force in the direction of the −y-axis is fixed on the y-axis. The condyle is fixed to prevent free movement. Finally, the lateral case was defined using two forces of 100 N facing in the direction toward the −z-axis. The condyle was fixed, and both jaw joints were suppressed to allow lateral displacement (Fig. 2).

Figure 2. Boundary conditions applied to models A–D depending on the feeding behavior studied in each case. In the bilateral case (A), two parallel forces of 100 N were applied between the principal tooth-row and the pal. The jaw joint was fixed in the −y-axis, and the condyle was fixed in all directions. In the lateral bite (B), two forces of 100 N were applied in the direction of the −z-axis, with the jaw joints suppressed and the condyle fixed. The axial twist case (C) was created using two parallel forces of 100 N, with one of them facing in the upward direction. Only one jaw joint was suppressed, and the other one was fixed as the condyle.

Models were meshed with an adaptive mesh of hexahedral elements. The mesh of the model consisted of about 2 million nodes and 1.4 million elements for model A and about 9 million nodes and 6 million elements for models B, C, and D. The last three models are more complex, because all the bones are present. In all cases, the size of the mesh was tested in front of result variations.

Finally, an alternative scenario is presented (see Supplementary Appendix S1) that reproduces, in a simple manner, the behavior of different intrinsic adductor forces during biting, achieved through a bending movement. The objective of this scenario is to validate the FEA models in comparison with published data from in vivo specimens (Metzger et al. 2005).

Results

The application of different types of bone contacts implies changes in the results obtained in each model. We evaluated Von Mises stress and displacements of the skull (Fig. 3). Von Mises stress is a good predictor of the strength of the skull and favors easier and more comprehensible comparisons between models due to the combination of all the stress tensor Cartesian components into a single value. Although it is usually applied to ductile materials, it remains an accurate criterion for predicting fracture location and stresses in bone tissue when the material is modeled using isotropic properties (Doblaré et al. 2004). To better understand how the changes in the stress develop in the bones of the skull, we computed the average value of stress in each bone of the skull (Fig. 4, Supplementary Table S1). We used the mesh-weighted arithmetic mean (MWAM), which gives the mean value of stress in non-uniform meshes (Marcé-Nogué et al. 2016). This was not applied to the nonsutured model because it was built as a single bone.

Figure 3. Von Mises stress and displacement distribution of models A–D of a Crocodylus niloticus skull under different feeding behaviors (bilateral, lateral, twist). Maximum values are colored in red. Due to the significantly lower values in the case of the axial twist feeding behavior, two different scales had to be used for a more accurate visualization of the displacement data.

Figure 4. Mesh-weighted arithmetic mean (MWAM) for each bone of models B–D under different feeding behaviors. As the nonsutured model (model A) was constructed as a single bone, MWAM cannot be applied. Anatomical Abbreviations: BO, basioccipital; EC, ectopterygoid; EO, exoccipital; F, frontal; J, jugal; L, lacrimal; M, maxilla; N, nasal; P, parietal; PL, palatine; PM, premaxilla; PO, postorbital; PRF, prefrontal; PT, pterygoid; Q, quadrate; QJ, quadratojugal; SO, supraoccipital; SQ, squamosal.

Feeding Scenario 1: Bilateral Biting

In the nonsutured and sutured linear models (models A–C) under a bilateral feeding behavior (Fig. 3), the high values of stress are mainly placed on the posterior part of the skull, with the basioccipital, the exoccipital, and the quadrate being the most affected bones. Furthermore, the bonded model (model B) is less affected compared with models C and D. Among the models, the mixed model (model D) shows the highest values of stress along the entire skull. Observing the stress reported in each bone (Fig. 4), the bilateral scenario presents higher stresses distributed across the skull, even though, the maximum values for the same bone, the basioccipital, are lower than in the lateral bite scenario. The maximum value of average stress (Supplementary Table S1) is 4.22 MPa present on the mixed linear–nonlinear model (model D), but lower in the nonseparated case (model C: 3.83 MPa), and even lower in the bonded case (model B: 3.26 MPa). Bilateral biting behavior (Fig. 3) deformation has higher values in the anterior part of the skull, with the maximum values of deformation on the premaxilla. The bilateral biting behavior does not result in elevated stress values in any of the lateral skull regions, as the deformation is distributed equally among them.

Feeding Scenario 2: Lateral Bite

The results obtained under a lateral feeding behavior (Figs. 3, 4) show that major stress values are in the posterior part of the skull in all nonsutured or sutured linear models (models A–C). These three models (models A–C) show a similar distribution of stress over the skull. In the case of the mixed linear–nonlinear model (model D), the maximum values, also located on the basioccipital, are higher than in preceding models. In model D, higher values are not only concentrated in the posterior part of the skull, high values of stress are also located in the anterior part, above the nasal and the maxilla. The maximum MWAM value (Supplementary Table S1) on linear models (models B and C) under lateral load case is located, in both cases on the basioccipital, with higher mean stress values on the nonseparated model (model C: maximum value of 5.75 MPa), than on the bonded model (model B: maximum value of 5.04 MPa). In the case of the mixed linear–nonlinear model (model D), the higher mean value is also placed on the basioccipital with a maximum value of 5.80 MPa. The deformation resulting from the lateral feeding behavior (Fig. 3) in all analyzed models is higher in the anterior part of the skull, particularly on the premaxilla, and progressively decreases to the posteriormost part of the skull.

Feeding Scenario 3: Axial Twist

Stress distribution on the nonsutured and sutured linear models (models A–C) under axial twist feeding behavior (Fig. 3) are mainly present in the center of the cranium, particularly on the maxilla, prefrontal, and nasal bones. Even if the maximum stress values are located at the same area, the values obtained for the nonsutured model (model A) are lower along the skull in general terms, being highest on the nonseparated model (model C). Under an axial twist behavior, the mixed linear–nonlinear model (model D) presents higher stress throughout the cranium. Under axial loads, the maximum values are lower than in the other two types of feeding behavior simulations (lateral and bilateral cases). The mixed linear–nonlinear model (model D) presents the highest mean values in this case (Supplementary Table S1), because the maximum value is 1.18 MPa, located in the right palatine bone. The maximum mean values in the other models (Supplementary Table S1) are 0.65 MPa, located in the prefrontal bone (corresponding to the case of the model without separation, model C), and even lower in the united model (model B): 0.58 MPa, located in the right nasal bone. The total deformation values under axial twist load case (Fig. 4) vary according to the direction of the x-axis force, being higher on the lateral edges when positive forces are applied and lower on the lateral edges when negative forces are applied. In all models, the deformation majorly affected the left lateral side of the skull. For the mixed linear–nonlinear case (model D), the left quadratojugal is even more affected on the left.

Discussion

FEA is a commonly used tool in functional morphology, because it allows the easy digitization of biological structures of both living and extinct taxa to solve biomechanical problems in highly complex biological geometries using computational analyses (Rayfield 2007). Although these types of analyses have become prevalent within comparative biomechanics fields, there are still important limitations that need further investigation. It cannot be ignored that FEA creates models that can be very useful, provided one keeps in mind that all models make assumptions and that they are, by definition, not literal representations (Anderson et al. 2012). Therefore, continuing to study this methodology is important to achieve more realistic approaches with fewer assumptions each time, as preceding works already demonstrated that some variables used during the FEA can impact the resultant outcome (Rayfield 2019). Furthermore, not all variables affect the results to the same extent, clouding the results, as it is not known a priori if these variables could cause potential modeling limitations (Walmsley et al. 2013). In this work, we propose a new methodology that enables the inclusion of the behavior of the sutures, without including the sutures themselves as a body in the FEA model, via contact elements. For this reason, the methodology presented herein has great potential to be used in finite element models of fossil species where soft tissues are not preserved. Moreover, the contact formulation available in FEA software also allows the user to configure and modify the parameters to adapt the contact to a linear or a nonlinear behavior of the soft tissue desired (Marcé-Nogué 2022), as done, for example, when modifying the value of the stiffness and the damping of the contact (see “Materials and Methods”; Wriggers 2002).

Discussion of Sutures

This work reveals how under the same analyzed feeding behaviors and identical variables the results obtained in our models are different when cranial sutures are defined differently. As the results suggest, the cranial sutures of Crocodylus niloticus have an important role in stress distribution and total deformation. In this regard, the use of extant taxa represents a valuable case study for testing the efficacy of new computational approaches to properly place and deeply analyze all sutures. Our results clearly demonstrate that depending on the type of conditions applied among the C. niloticus cranial suture network, different results are obtained, thus confirming that sutures directly influence skull biomechanics. The results of stress obtained reveal the capacity of dissipation in the different models, with the nonseparated case (model C) being the one that, apparently in a more optimal way, facilitates the skull not concentrating stresses over a specific bone by relieving stresses over the whole skull.

This phenomenon contrasts with the models that do not allow movement between bones (models A and B). This response was also noted by Curtis et al. (2013), who concluded that patent sutures (comparable to model C) distributed the stresses and strains better throughout the skull than fused ones (comparable to models A and B). In the mixed linear and nonlinear cases (model D), stresses are concentrated, with exceptionally high values, over particular bones that have nonlinear sutures, although it should be keep in mind that no-separation contacts are dominant in model D (whereas in model C all sutures are considered to be no-separation ones). However, stresses are in general also greater than in the other models, probably because the no-separation sutures dissipated the stress over the whole skull. Based on these results, model D shows higher values over the skull (produced by linear contacts) and concentrated focuses of stress among certain bones (produced by nonlinear contacts).

All in all, the mixed nonlinear and linear model D is probably the closest to the ontogenetic stage (adult) of the analyzed individual, but is maybe not the most optimal one for feeding behaviors, as skulls represent a compromise between optimal (adult) function for these loads and other constructional constraints that explain why sutures vary throughout the skull. On the other hand, even with frictionless sutures, the nonseparated model (model C) would represent the case where these types of sutures ossify, creating more contact between bones over the whole skull; the bonded and no-suture models (models B and A), respectively, represent the last stages, wherein sutures start to fuse and lose movement. Overall, this new methodology also has interesting implications for understanding how ontogenetic stages (based on cranial suture closure) cause cranial biomechanical differences. The models developed in FEA can be used as hypotheses regarding the mechanical behavior of a specimen (Porro et al. 2013). The suitability of the models is dependent on the objective or scenario under study. Herein, for the specimen analyzed in this study, a mixed linear–nonlinear model is the most accurate, as it depicts the development in suture morphology over the skull and through the age of the specimen.

Feeding Behavior

Feeding behaviors affect the biomechanics of crocodile skulls, as well as the suture morphology. This conclusion is in accordance with Walmsley et al. (2013), who reported that feeding behaviors greatly influenced their results. According to our results, the bilateral biting scenario is the feeding behavior that shows especially high values on the posterior part of the skull.

Work on living taxa benefits from a posteriori validation of the results on the analyzed feeding behaviors. As morphology is the result of a summary of genetics, development, and environment (Rayfield 2007), particular mechanisms and functions are the consequence of the adaptation of each taxon to a concrete behavior, biological role, and performance, all of which are relevant parts of their ecology. Hence, the results obtained allow the extrapolation of ecological data from the examination of the response of its cranial suture morphology (Rayfield 2007). Both morphological and ecological outcomes, in combination, help us understand the ecomorphology of the taxon. Ecomorphology links the interrelationship between morphological variation with the ecology of the anatomical structure by studying the material's composition, its arrangement, and the physical properties that could affect all its levels of organization (Bock 1994). Diet in crocodilians undergoes considerable changes across size and age, and even with habitat, so different individuals of the same taxon would present different outcomes when the effects of feeding behaviors are studied (Fergusson 2010). One of the reasons for this variability is that sutures display ontogenetic changes, as they tend to fuse through time (Curtis et al. 2013). The analyzed specimen was an adult male specimen of C. niloticus, one of the largest species of all living crocodilians. Commonly, large crocodilians feed by taking off parts of their prey by rotating along their longitudinal axis (axial twist scenario) or kill their prey by suffocating them under the water instead of killing by crushing and biting (bilateral scenario) (Cleuren and de Vree 2000). The results obtained here are a demonstration of how the skull adapted to the most common feeding behavior of the animal, the axial twist, and less to the more uncommon, but not precluded, bilateral feeding behavior.

Future Direction and Implications for Paleobiological Studies

The present study is an example of how it is achievable to examine the morphology to understand how cranial sutures work in different scenarios and how this approach could be used for extinct taxa (Rayfield 2007). Further investigations with this methodology could be applied to get important new paleobiological data (Rayfield 2007) that are impossible or difficult to get from other sources, such as data related to soft tissues, which are rarely preserved in the fossil record. Moreover, this work achieves the aim of using nonlinearities in functional morphology models. In fact, nonlinearities are eminently suitable when modeling soft tissues as contacts between separated bones (Marcé-Nogué 2022), and this work highlights the usefulness of nonlinearities in FEA models.

Learning about the ecology of extinct animals is a complex process, because there is usually scarce evidence or remains that could help to disentangle how the animal lived and how its ecological niche was (Witmer 1995). The present model could serve as reference for future studies of cranial morphology on extinct taxa with no evidence of suture remains of soft tissue and could provide character correlation between morphological and evolutionary hypotheses.

Conclusions

The present work reveals that the type of contact used in each suture clearly affects the outcome of the FEA, demonstrating the importance of considering cranial sutures in any biomechanical analyses. Our results show how the lack of sutures or the presence of sutures with restrictions of movement caused reduction on the stress distribution over the skull, while sutures with movement helped the skull to relieve stress and prevented bones from enduring high stress values.

Comparison with different feeding behaviors is useful to further understand the implications of external variables, revealing that feeding behaviors highly influence skull biomechanics. Bilateral feeding behavior caused the most major stress on the bones, in contrast to axial twist, which caused the least stress. In addition, our results demonstrate that FEA is a useful method to comprehend a suture's function, especially in fossil taxa, and biomechanics when contact elements are considered, provide new clues to understand their biomechanical role. All the contact types defined in this work can be used to test different biomechanical hypotheses, because they provide different movement behaviors between adjacent bones, adapting the requirements of the desired FEA model to the objective of the research. These results, combined with previous work, bring us a step closer to creating more realistic models that could be applied in future paleontological investigations, as they generate new ecological and evolutionary hypotheses to be tested. Especially in fossil taxa, where no or poor preservation of soft tissues occurs, contact elements are a good solution to model cranial sutures.