Validation experiments on finite element models of an ostrich (Struthio camelus) cranium

Ex-vivo experimentation

An ostrich skull, frozen with jaws in occlusion (a farmed sub-adult, aged 6-9 months, obtained after slaughter for meat production), was thawed and the jaw adductor muscles were dissected out (for an overview of avian cranial musculature see Holliday & Witmer, 2007). Ostriches were used due to their large size and the ability to adhere gauges to the surface of the bone, and as they retain patent sutures in some skull regions into adulthood. The M. pseudotemporalis superficialis was isolated, as this muscle provided the largest and most robust attachment site for loading in the experiments. The muscle was then weighed, before immersion in 15% nitric acid, then bisected to measure individual fibre lengths and pennation angle to determine maximum muscle force (after Biewener & Full, 1992).

Due to difficulties in loading an ex-vivo skull with accurate muscle loads and orientations, a number of methods have been used. These include directly loading the original (albeit formaldehyde preserved) muscles (Kupczik et al., 2007), or screwing material into the skull at the site of muscle attachment on which to apply the forces (e.g., Bright & Rayfield, 2011a). This paper attempts a new method to circumvent some issues associated with previous attempts. Two artificial tendons were made using a strand of carbon fibre, resin and fibreglass. The carbon fibre was looped, with the two loose ends being frayed before being sandwiched between two sets of three fibreglass layers (Fig. 1A). The fibreglass layers were all then impregnated with epoxy resin creating artificial ‘tendons’ that were strong, light, and flexible.

The ‘tendons’ were attached to the M. pseudotemporalis superficialis attachment site (delineated by a bony ridge) on the dorsolateral side of the cranium with cyanoacrylate (M-Bond 200; Vishay Micro-Measurements, Basingstoke, UK), before being screwed in place using three 3.5 mm self-tapping screws (Fig. 1B). Therefore, it should be noted that the loadings replicate only one muscle load and are not truly physiological or representative of a bite. This was considered acceptable, given that the aims of the experiment were specimen-specific validation and testing model sensitivity to the effects of the material properties, sutures and rhamphotheca. During dissection, preparation and testing, the skull and rhamphotheca were kept as moist as possible by frequently covering both with a 50:50 mixture of glycerol and water.

Rectangular planar rosette strain gauges (C2A-06-062LR-350; gauge length 1.52 mm, backing 7.04 × 10.41 mm, Vishay Micro-Measurements, Basingstoke, UK) were attached at 13 sites across the skull (Figs. 2 and 3). The gauge sites were prepared by removing any remaining periosteum with pumice powder (GC-5 Pumice powder; Vishay Micro-measurements, Basingstoke, UK) and the site cleaned with 95% ethanol. Gauges were attached to the sites with cyanoacrylate adhesive (M-Bond 200; Vishay Micro-Measurements, Basingstoke, UK). The gauges were then covered with a waterproof silicon rubber coating for protection (3140 RTV Coating; Dow Corning, Midland, MI, USA). Strains from the 13 gauges were connected to an amplifier (5100B; Vishay Micro-Measurements, Basingstoke, UK), recorded, and converted to principal strains and strain orientations using STRAINSMART 4.01 software (Vishay Measurements Group, Basingstoke, UK). Preliminary drift tests were conducted for all gauges before loading and no thermal drift was observed over a 30 s time period.

The cranium was placed upon a specially designed loading rig adapted from a previous study (Bright & Rayfield, 2011a) which supported the cranium bilaterally under the quadrate at the quadrate/articular (QA) joint with two aluminium bars, and at the anteriormost part of the rhamphotheca with a flat block (Fig. 2). Each of the aluminium bars rested between the quadrate articular surfaces and prevented dorsoventral and mediolateral movement of the posterior of the cranium under loading. To prevent anterior motion during loading, a small, flat, thin plate was attached to the end of each bar that contacted the anterior ventral quadrate during the test. Loads were applied by tying 3 mm low-stretch polyester cord (3 mm Magic Speed; LIROS GmbH, Berg, Germany) to the carbon fibre loops using bowline knots at one end, and hanging balances (HCB 50K100; Kern & Sohn GmbH, Balingen, Germany) at the other. The cord was extended over low-friction pulleys (Size 1 upright block, Barton Marine Equipment, Kent, UK) to apply the load approximating the muscle line of action. The hanging balances (precision of 1 N) were attached to rigging screws (6 mm Fork Bottlescrew; Sea Sure, Hampshire, UK) which allowed for manual tightening of the screws to apply tension without twisting.

Due to the limitations of the amplifier, only four gauges could be tested at any time. During the testing, channel 3 on gauge 6 failed, so gauge 6 was discounted from the study, thus recordings were taken from three groups of four gauges (G1-G4, G5-G9, G10-G13). Each gauge was zeroed, and loads of 32 N were applied to each side by tightening the rigging screws simultaneously. These forces were within the limits of an estimated force for the M. pseudotemporalis superficialis calculated from the dissection (muscle body = 4.18 g, average fibre length = 18.4 mm, average pennation angle = 30.2°; using a muscle stress of 0.3 Nmm⁻² gave an average force of 44.7 N for each side) and caused cranial deformation, but not so large that fracture may occur under loading. All gauges were measured twice during two complete trials. After each batch of three or four gauges was loaded, the residual load remaining in the hanging balance after unloading was measured: right tendon average = 6.86 N [S.D. = 1.46], left tendon average = 5.71 N [S.D. = 0.488]. These residual loads probably resulted from friction, knot tightening and elasticity in the cord, and may have resulted in slight variations in the actual loading conditions between each trial. Thus, average residual loads were deducted from the initial load, to give average applied loads of 25.1 N on the right tendon, and 26.3 N on the left tendon. The data from each gauge during a trial is an average of 10 recordings per second taken for 60 s over which a static load was applied (a total of 600 readings per gauge).

FE models

Prior to dissection, the skull was subject to micro-computed tomography (micro-CT) at the University of Hull, UK using a X-Tek HMX 160 scanner. Due to the size of the specimen, it was scanned in two parts (anterior and posterior skull) taken at slightly different orientations (756 slices and 846 slices, with voxel sizes of 0.1594 mm and 0.1425 mm respectively). Both scan sets were rotated and resampled to the same voxel size (0.1594 mm resolution) in Avizo 6.3 (Mercury Computer Systems, Arlington, Virginia, USA) where the cortical bone, trabecular bone, sutures and rhamphotheca were initially segmented automatically before final details were segmented manually. A resulting stereolithographic surface was then imported into Hypermesh 10.1 (Altair Engineering, Menlo Park, California, USA). Surfaces for each individual material (cortical bone, trabecular bone, suture and rhamphotheca) were meshed using the ‘shrink wrap’ function in Hypermesh. The result was a continuous mesh with the rhamphotheca directly attached to the cortical bone, with the sutures connected to cortical bone (and directly to the rhamphotheca in the places where the rhamphotheca overlies the suture), and many separate enclosed surfaces of trabecular bone located throughout the skull. A surface-based meshing approach was then used to generate 3D volumetric meshes of the cranium. Five meshes of between 1,333,804 and 1,749,533 elements (between 298,515 and 385,811 nodes) were produced for convergence testing (Bright & Rayfield, 2011b) for both Tet4 (four nodes per tetrahedral element) and Tet10 (10 nodes per element) models. The ‘shrink wrap’ function in Hypermesh was unable to produce usable meshes for element numbers outside of these limits for the ostrich cranium. When results at equivalent gauges sites for successive meshes were less than 10% different in maximum and minimum strain, and the difference in total strain energy was less than 5%, the mesh was considered converged. On this basis, the Tet10 mesh with 1,749,533 elements was selected for further analysis (see Supplemental Information for further information). Each model surface was cleaned by removing t-connections, holes, and free edges. Virtual gauges matching the location of the ex-vivo gauges were created using 2D membrane elements placed as close to the actual gauge locations as possible, in order to account for any discrepancy between the 3D strains of the models being projected in to the 2D plane of the gauge (Bright & Rayfield, 2011a). Load and constraint points equivalent to those experienced in the ex-vivo cranium were applied to the model. Two constraints prevented dorsoventral translation at the anterior of the beak, and four on each quadrate (at the jaw joint) preventing movement in all directions. The model was loaded with a force of 25 N on the right side and 26 N on the left in the orientation specified by the experimentally applied load (Figs. 2 and 3). These loads were distributed across 14 nodes at each muscle attachment site (four nodes at each of the three screw locations to replicate higher point loads induced by the introduction of a screw and two further nodes equidistant between the screw locations).

Nanoindentation data

Nanoindentation values for ostrich cranial bone, rhamphotheca and suture are listed in Table 2.

Ex-vivo experimental strain

The maximum (E_max; Table 3), minimum (E_min; Table 4) principal strains, strain ratios (Table 5) and strain orientations (Table 6) show similar trends between the two ex-vivo tests, but there is variation in magnitude between tests. The data, however, are significantly correlated for all metrics (see Figs. S3–S9).

As expected from the loads applied and bending observed by eye during the experiments, principal strains are highest in the premaxilla (G5), nasals (G7) and lacrimals (G8 and G9), and lowest in the braincase regions (G10–13), despite the latter area being closer to the applied muscle loading. Strain ratio is a measure of the E_max divided by the absolute value of E_min and is useful as it removes the effects of strain magnitudes. Where the strain ratio is between 0.67 and 1.5, shear is the predominant regime. When the strain is greater than 2, or less than 0.5, the regime is either primarily tension or compression respectively (Dzialo et al., 2014). Most of the cranium is in compression, which reflects the bending experienced on the dorsal surface, but there is tension in gauge 7 and a combination of shear and compression in gauge 12 (Table 5).

FEA

Material properties from the literature (models L1–L9)

Overall there is a high level of similarity between strains in all of the FE models with material properties taken from the literature (Figs. 4A and 4B). There are some similarities between FE tests and the ex-vivo experiments but also some key differences. In terms of absolute values of strain magnitude, the areas of closest numerical match are G9–G13, whilst G1–G4 strains follow similar trends to the experimental data but are all offset in E_max with the models being less stiff (i.e., higher strain magnitudes) than the experiment. The FE models represent the peak in E_min at G5 but do not represent the peak in E_max at G7. The FE model strain ratios show the same mode of strain (tension or compression) for G3–G7 (Fig. 5 and Table 5). However, most of the rest of the strain ratios in the models show a combination of shear or either tension or compression, whilst compression is found for the gauges in the ex-vivo experiments. Only G4 (on the posterior rhamphotheca) and G7 (on the nasal) have strain orientations within 10°of the orientations recorded in the FE models (Fig. 6 and Table 6). Beyond this there is a range of discrepancy between the experimental strain orientation data and the models (Table 6). Euclidean distance (ED) metrics show that the model with the closest match to the ex-vivo data for E_max, E_min and strain orientations was L3 but for strain ratios and L2 provided the best overall match (Table 7). Regressions show there is no significant relationship between the experimental data and the averaged data for the FE models (i.e., values for L1–L9 averaged) for E_max (Figs. S10–S11), and strain orientations (Fig. S16). There is a significant correlation between the experimental data and the averaged data for E_min (Fig. S12) and strain ratios (Fig. S14), but when possible outliers are removed (see Supplemental Information) the correlations in E_min (Fig. S13) and strain ratios (Fig. S15) become non-significant. The Euclidean distances also change if the possible outliers are removed, with L9 performing best for strain magnitudes instead of L3 (Table S3).

Table 6:

Principal strain orientations from grid one direction on the gauge, in degrees.

Gauge site	1	2	3	4	5	7	8	9	10	11	12	13
Ex-vivo
Test 1	34	60	17	55	−10	55	−6.5	48	16	−59	−46	−29
Test 2	48	68	78	59	−4.5	57	−6.0	45	12	−59	−46	−29
Literature
L1	57	23	53	54	56	45	38	37	78	70	55	19
L2	64	20	55	55	52	47	19	17	43	50	54	17
L3	62	26	52	55	51	47	18	18	44	51	54	17
L4	61	15	57	55	51	49	28	1.0	71	70	56	25
L5	60	23	54	55	51	49	27	2.5	71	70	56	25
L6	60	21	54	55	51	49	26	3.7	72	69	56	25
L7	61	18	55	55	51	49	25	3.5	71	70	56	25
L8	61	16	56	55	51	49	23	5.1	72	69	56	25
L9	59	20	55	54	56	45	43	39	72	69	56	25
Posthoc
PH1	60	26	53	55	52	47	28	11	67	70	56	24
PH2	59	23	54	55	52	49	30	4.0	72	69	55	25
PH3	59	24	53	55	52	49	31	7.3	72	69	56	25
PH4	59	23	54	54	56	45	40	37	72	69	56	25
PH5	61	21	55	55	52	47	32	11	68	70	56	24
PH6	61	25	53	55	50	49	11	3.6	68	70	56	24
PH7	63	25	53	55	50	48	8.4	5.0	67	70	56	23
PH8	62	26	53	55	51	48	16	4.8	67	70	56	23
Nanoindentation
N1	63	2.6	64	55	51	49	22	3.7	75	70	56	24
N2	63	6.0	63	55	52	49	28	1.5	73	70	56	25
N3	63	1.9	65	55	51	49	21	3.7	75	70	56	24

DOI: 10.7717/peerj.1294/table-6

Table 7:

Euclidean distances from ex-vivo experimental measures to FE model data.

Data in bold are the closest to the experimental data for each specific metric.

	Emax	Emin	Strain ratio	Orientation
Literature
L1	922	591	4.39	203
L2	893	657	4.04	180
L3	873	556	4.10	179
L4	890	648	4.89	207
L5	877	561	4.93	204
L6	874	557	4.91	204
L7	886	627	4.90	205
L8	884	627	4.89	205
L9	933	805	4.09	205
Posthoc
PH1	864	559	4.96	201
PH2	860	560	5.28	204
PH3	864	690	4.90	204
PH4	934	660	4.82	204
PH5	931	1,003	4.76	203
PH6	917	1,108	4.92	200
PH7	920	1,106	4.68	200
PH8	1,002	1,587	4.75	200
Nanoindentation
N1	1,171	1,912	4.75	210
N2	932	873	4.74	210
N3	3,773	8,526	4.55	210

DOI: 10.7717/peerj.1294/table-7

Posthoc values (models PH1–8)

For E_max, the results are similar to those for the literature models. Gauges 1–4 show a similar pattern of E_max values, but again, FE model magnitudes are greater than the ex-vivo data (Fig. 4C, Table 3). Some of the models (PH6–8) replicate peaks at G7, but do not achieve the magnitude of strain recorded in the ex-vivo data (Fig. 4C). Other models do not replicate this peak. For E_min (Table 4), the strain magnitudes at sites G2 and G3 in the FE models closely match the ex-vivo experiment but strains in G1 and G4 are higher in the FE models than experimentally (Fig. 4). The peak at G5 in the experimental data is replicated in the FE models, as is the decline in E_min magnitude in the more posterior gauges G8–G13.

The ex-vivo experimental data are broadly similar to the FE models with respect to strain ratios from G3–7 (Table 5). However, most of the rest of the strain ratios in the models show a combination of shear or either tension or compression, whilst compression is found for the gauges in the ex-vivo experiments. Again the strain orientations in all of the models (PH1–PH8) vary in how much they match experimental strain orientations, with only G4 and G7 falling within 10° of the ex-vivo values (Fig. 6 and Table 6). The Euclidean distances show that PH2 is most similar to the ex-vivo data for strain magnitudes whilst PH7 best matches the strain orientations and strain ratios (Table 7). The regressions show there is no significant relationship between the experimental data and the averaged data for the FE models for E_max (Figs. S10 and S11), and the strain orientations (Fig. S16). There is a significant correlation between the experimental data and the averaged data for E_min (Fig. S12) and strain ratios (Fig. S14), but when possible outliers are removed the correlations in E_min (Fig. S13) and strain ratios (Fig. S15) become non-significant. The Euclidean distances also change if the possible outliers are removed, with PH4 becoming best for strain orientations (Table S3).

Rhamphotheca

One pair of models (L4 and L6) possessed the same material properties except for the rhamphotheca (Table 1). As with the sutural tests, an increase in Young’s modulus of keratin leads to a decrease in the overall maximum principal strains within in the model (Fig. 8A and Table 3). For minimum principal strains, the peak values are found in the model with the lowest Young’s modulus (Fig. 8B and Table 4). The strain ratios are approximately the same between the models with differing properties, and there are no differences with respect to compression or tension at the gauge sites (Fig. 8C and Table 5). Strain orientations are all very similar between both models (Table 6).

Material property considerations

Within the models, each material is modelled as isotropic. It is unlikely that many biological bony structures are isotropic throughout, and probably have variable material properties and anisotropy depending on location (e.g., Zapata et al., 2010). This is supported by our nanoindentation data documenting the high variability in material property standard deviations across just an individual test location, and between different sites on the rostrum (Table 2). The use of isotropic and internally homogeneous material properties for bone, keratin and sutures may explain the lack of statistically significant correlation between our ex-vivo results and FE models, despite the fact that the models are broadly deforming as expected. Furthermore, the methods used here may fail to accurately model the complex bone/beak/suture interfaces. The use of anisotropic material properties may have improved model data fit, but this was beyond the scope of this study and would require more extensive nanoindentation sampling across the entirety of the cranium. It is possible that such concerns could be mitigated by restricting future validation tests to smaller regions of presumably similar function, for instance either the facial region (e.g., Soons et al., 2012a) or the braincase, particularly when model performance seems to be so affected by anatomical location, as seen here and in previous validation studies on pigs (Bright, 2012) and ostrich mandibles (Rayfield, 2011). There is also a chance that there is intraspecific variability with regards to the material properties; a different ostrich was used for the nanoindentation than for the FE modelling. In this case the animals were from the same stock and culled at the same age, so these effects were minimised as best possible.

Effects of sutures

The ex-vivo experimental data follows the results of many previous experimental studies showing that the strain environment can change dramatically between bones separated by a suture (Jaslow, 1990; Herring & Mucci, 1991; Jaslow & Biewener, 1995; Herring & Teng, 2000; Bright & Gröning, 2011; Bright, 2012). The maximum (Table 3) and minimum principal strains (Table 4), strain ratios (Table 5) and strain orientations (Table 6) change markedly between G5, G7 and G9 which are located on the premaxilla, right nasal and right lacrimal and each separated by a suture (Fig. 3). There is a rapid drop in Von Mises stress across the sutures in the rostra, particularly between the high stress area of the premaxilla and the lower stressed nasals (Fig. 7).

The presence of sutures and variation in their elastic properties is known to effect overall FE model stress in certain taxa (Reed et al., 2011; Bright, 2012; Moazen, Costantini & Bruner, 2013), but not in all (e.g., macaques—Wang et al., 2010). The sutural effects may vary depending on the taxon, volume of sutures or even the method used to model the sutures. The sutures in this study were modelled as tetrahedral elements with lower Young’s moduli than the surrounding bone, as for previous analyses (Moazen et al., 2009; Wang et al., 2010; Reed et al., 2011; Bright & Gröning, 2011; Bright, 2012). The presence of sutures in our models affected all measured metrics such that the models with sutures were better at replicating the experimental minimum principal strains, strain ratios and orientations than those that did not possess sutures. However, models L1 and L9, which do not have sutures, have values on each extreme for the Euclidean distances: L1 is very different in maximum principal strain to the ex-vivo data, and L9 is the most similar to the experimental maximum principal strain. Sutures alone therefore do not rectify all of the differences seen between the ex-vivo experimental values and the FE models.

Ostriches possess cranial sutures between all of their cranial bones, and through ontogeny many of these fuse. This fusion is not uniform; those across the braincase fuse first and those in the high strain regions of the premaxilla/nasal/lacrimal region stay unfused into adulthood (AR Cuff, pers. obs., 2011). This, combined with the experimental data, suggests that the sutures are playing an active role in strain mitigation; regions of low strain such as the braincase do not need sutural flexibility so they fuse, whilst sutures in regions of higher strain such as the beak remain patent and maintain a functional influence.

Effects of the rhamphotheca

The rhamphotheca is a complex, composite, keratinous structure that covers portions of the bony rostra in birds. Its extent varies between palaeognaths and neognaths, and within each group, but in all taxa its primary function is feeding. Three previous studies have tested the ability to model the beak using FEA in toucans (Seki, Mackey & Meyers, 2012), and finches (Soons et al., 2012b; Soons et al., 2012c). These studies focus only on the rostral regions containing the beak, but show that FE models can accurately replicate strain patterns across the beak. Indeed, they suggest that the keratinous rhamphotheca surrounding the bone (which may well be heavily trabeculated in the beak region) protects the bone from buckling during impact loadings and removal of the rhamphotheca leads to failure of the bone under equivalent forces (Soons et al., 2012b).

The rhamphotheca in ostriches was modelled as a homogenous structure overlying the bone with a lower Young’s modulus than the cortical bone it surrounds. The experimental data shows that the gauges on the dorsal surface of the rhamphotheca are less strained than the nearest gauges directly on the bone, particularly in minimum principal strain. Increasing the Young’s modulus of the keratin decreases the strains experienced both on the rhamphotheca and throughout the model. The mechanical influence of a rhamphotheca would be especially important across the areas of high strain and thin bone of the rostra and buccal margins, as was found in the beaks of herbivorous theropod dinosaurs (Lautenschlager et al., 2013). In ostriches, the rhamphotheca extends over the premaxilla, and the ventrolateral edges of the maxillae and anteriormost jugals. There may also be a gradational transition of keratinisation from the rhamphotheca sensu stricto to the skin (Genbrugge et al., 2012). This phenomenon was not modelled here, but warrants further investigation. Comparison to material property data obtained from previous nanoindentation studies on various areas of the beak suggest that ostrich beak keratin may have a far lower Young’s modulus (by up to a factor of two) than that seen in other birds (Seki, Bodde & Meyers, 2010; Soons et al., 2012a). It should be noted that the areas tested in our nanoindentation study were the internal regions of the rhamphotheca, which may have a far lower Young’s moduli than the tougher outermost layer. Further modelling will be required to test the overall effect this may have on strain patterns and orientations in FE models.

Mechanical behaviour of the beak versus the braincase

In general, the experimental and modelled strain ratios match more closely in the beak than they do in the braincase. This may be because the braincase is closer to the constraints emplaced on the quadrate. Yet, strain magnitudes in the braincase region are low, meaning a small change in strain magnitude could switch the braincase from experiencing primarily compression to primarily tension, or vice versa. These very different strain signals point to a functional signal represented in this data, suggesting a modular biomechanical response between the functions of the beak and the cranium.