## Abstract

Classic beam theory is frequently used in biomechanics to model the stress behaviour of vertebrate long bones, particularly when creating intraspecific scaling models. Although methodologically straightforward, classic beam theory requires complex irregular bones to be approximated as slender beams, and the errors associated with simplifying complex organic structures to such an extent are unknown. Alternative approaches, such as finite element analysis (FEA), while much more time-consuming to perform, require no such assumptions. This study compares the results obtained using classic beam theory with those from FEA to quantify the beam theory errors and to provide recommendations about when a full FEA is essential for reasonable biomechanical predictions. High-resolution computed tomographic scans of eight vertebrate long bones were used to calculate diaphyseal stress owing to various loading regimes. Under compression, FEA values of minimum principal stress (*σ*_{min}) were on average 142 per cent (±28% s.e.) larger than those predicted by beam theory, with deviation between the two models correlated to shaft curvature (two-tailed *p* = 0.03, *r*^{2} = 0.56). Under bending, FEA values of maximum principal stress (*σ*_{max}) and beam theory values differed on average by 12 per cent (±4% s.e.), with deviation between the models significantly correlated to cross-sectional asymmetry at midshaft (two-tailed *p* = 0.02, *r*^{2} = 0.62). In torsion, assuming maximum stress values occurred at the location of minimum cortical thickness brought beam theory and FEA values closest in line, and in this case FEA values of *τ*_{torsion} were on average 14 per cent (±5% s.e.) higher than beam theory. Therefore, FEA is the preferred modelling solution when estimates of absolute diaphyseal stress are required, although values calculated by beam theory for bending may be acceptable in some situations.

## 1. Introduction

In comparative biology and palaeontology, the relationship between diaphyseal cross-sectional properties of long bones and body mass (*M*_{b}) is frequently used as an indicator of skeletal strength and rigidity. In many instances, *relative* values of such properties as cortical cross-sectional area (*A*_{cort}), second moment of area (*I*) and polar moment of area (*J*) are used in comparative studies, from which skeletal posture and *in vivo* function may be inferred [1–4]. However, in other instances, diaphyseal cross-sectional properties are incorporated into equations to predict *actual* skeletal stress. Reliable estimates of diaphyseal stress are essential in determining maximum upper limits to terrestrial vertebrate body mass [5], estimating safety factors (the ratio of yield stress to maximum functional stress), and in reconstructing gaits via locomotion modelling [6]. While the preferred method for obtaining absolute values of stress must always be the *in vivo* application of strain gauges, this is often impractical, owing to ethical constraints associated with the study species, the sample size required or the fossilized nature of the sample.

When a particular biomechanical model is used to estimate stress in a range of specimens, we assume that the mechanical consequences of differing skeletal morphology will be illuminated [7]. Yet, the error magnitude involved in the calculation of stress and strain is also a function of the underlying geometry of the skeletal element (both external morphology and internal architecture). A model's suitability for estimating stress is dependent upon the extent to which each biological specimen in turn meets the conditions of the model. In a sample containing high morphological variability, application of classic beam theory to estimate stress may result in inconsistency in both the direction and the magnitude of errors, and may mask the functional morphological signal of interest.

Euler–Bernoulli beam theory [8] (hereafter referred to as ‘classic beam theory’) provides a means of calculating deflection of a beam and has been extensively applied to the estimation of stresses in vertebrate long bones, owing in large part to its simplicity [9–12]. Compressive loads acting through the centroid of the cross section generate normal stresses defined as
1.1where *σ*_{comp} is compressive stress, *F* is the applied force and *A*_{cort} the cross-sectional cortical area [13]. However, in instances when the beam possesses a degree of curvature, axial components of the applied force act longitudinally around the curvature and induce bending moments. In the case of long bones, the extent of induced bending is a function of the radius of curvature of the element, and resultant bending stresses may come to dwarf those of axial compression during dynamic loading [14]. Previous authors have sought to explain the scaling behaviour of long bone dimensions in terms of maintaining resistance to Euler buckling (failure of a thin-walled straight column under axial compression owing to elastic instability) [15,16], despite concern elsewhere that elastic deformation is unlikely to be an important factor in failure of mammal bones [17]. In addition, such scaling models have so far ignored the potential role of curvature in their calculations. Figure 1 illustrates the varying morphology characteristic of vertebrate hindlimbs included in this study, and highlights the divergence of long bones from the idealized straight beams considered in engineering.

A limited number of biomechanical studies have incorporated curvature into estimates of bending stress [18–20] or have investigated scaling of curvature to body mass [18,21,22]. This is particularly the case in the modelling of stress in primate mandibular symphyses [23,24], as they exhibit a particularly high degree of curvature relative to other skeletal elements. However, there has been no systematic application of curvature-corrected equations to long bones across the comparative anatomy and palaeontological literature, and the effect of ignoring this geometry on beam theory estimates has not been adequately explored.

The bending stress (*σ*_{bending}(*y*)) varies with position across a symmetric beam and is estimated as
1.2where *M _{x}* is the bending moment about the

*x*-axis,

*y*is the perpendicular distance to the neutral section and

*I*is the second moment of area about the

_{x}*x*-axis [13]. Application of this equation assumes that the cross section is symmetrical about the axis on which loading is occurring, that cross-sectional shape is maintained downshaft, and that plane sections remain undistorted and normal to the long axis following loading, i.e. ignoring shear deformation [13]. In beams possessing a low aspect ratio (length/diameter;

*l*/

*d*), warping owing to transverse shear contributes to the total stress experienced in a cross section. Standard engineering practice suggests that beam theory equations are reasonably accurate for objects only with an

*l*/

*d*ratio of 16 or greater [25], and many vertebrate long bones fall below this aspect ratio at which shear deformation could justifiably be ignored (table 1). Therefore, values of bending stress calculated using Euler–Bernoulli theory are likely to be underestimates in stout, long bones.

In the biological literature, the relationship between axial compression, bending and load vector has been described by combining equations (1.1) and (1.2)
1.3where *σ*_{combined}(*y*) is the sum of compressive and bending stresses, and *θ* is the angle between the loading direction and the longest principal axis [29]. Therefore, when *θ* = 0°, *σ*_{combined} is equal to *σ*_{comp}, while when *θ* = 90°, *σ*_{combined} is equal to *σ*_{bending} (figure 2). By combining equations (1.1) and (1.2), total stress can be calculated for long bones loaded neither in pure compression nor pure bending, but at some intermediate angle. However, this equation potentially suffers from the compounded problems of equations (1.1) and (1.2) when applied to irregular geometries.

The maximum shear stress owing to torsion (*τ*_{torsion}) in a hollow elliptical beam is calculated as
1.4where *T* is the applied torque, *r*_{ap} and *r*_{ml} are the radii in the anteroposterior and mediolateral directions, respectively, and *q* is the ratio of inner radius to outer radius [30]. Equation (1.4) makes the assumption that the endosteal and periosteal contours are similar concentric ellipses. Alternatively, when the cross section is characterized by possessing thin walls, the average *τ*_{torsion} in a hollow elliptical section can be approximated using an alternative ‘thin-walled ellipse’ model
1.5where *t* is the thickness of the cortical wall (assuming a uniform thickness across the section) [30]. If the highest torsional stresses are considered to occur where the wall thickness is at a minimum (*t*_{min}), a modified Bredt's formula may also be used to approximate *σ*_{torsion} in sections of varying cortical thickness [31]
1.6where *A* is the area enclosed by the median boundary (figure 3). This ‘minimum wall thickness’ model has been shown to be more suitable in estimating torsional stresses in asymmetric human tibial bones than the hollow ellipse model of equation (1.4) [32].

It is clear that all the beam formulae above are idealized approximations to the loading conditions actually experienced by long bones during dynamic loading, but the degree of error that they introduce is unclear. Classic beam theory remains the most practical and highly favoured modelling solution for very large comparative datasets, and for instances in which information regarding myology and material properties is lacking. However, finite element analysis (FEA) is increasingly becoming the preferred solution for estimating mechanical behaviour when a sample has an irregular and highly variable geometry [32]. Until recently, it was difficult to obtain an accurate model geometry for FEA of complex morphologies, hampering its widespread application to large biological datasets. However, access to computed tomography (CT) facilities capable of providing detailed and structurally faithful three-dimensional models is becoming cheaper and easier. As a result, the scope is broadening for generating a sample size, previously achievable only via simpler modelling techniques. Given the additional complexity of such an approach, it is reasonable to ask whether stress values predicted by CT-based FEA differ significantly from those of classic beam theory when applied to vertebrate long bones. Here, we set out to answer this question by comparing stress predictions of theoretical simple beam equations against those of FEA in a diverse sample of long bones.

We seek to test (i) *geometrical effects* of diaphyseal cross-sectional shape failing to meet the assumptions of beam theory formulae; (ii) *loading effects* of shaft curvature preventing solely compressional and torsional loading; and (iii) *shear stress effects* of incorporating shear stress components into stress estimates. While the incorporation of heterogeneous material properties into FEA is commendable and results in closer agreement between FEA models and *ex vivo* results [33], here we have sought solely to explore the consequences of incorporating the inherent complexity of long bone geometry into such models for the purpose of comparative anatomical studies.

In reality, a synergistic combination of FEA modelling and *ex vivo* experimental validation may provide the best means of reliably testing the mechanics of vertebrate long bones [34]. However, in this study, we provide recommendations for the application of FEA and improvements to existing beam theory equations for the majority of instances when destructive mechanical testing is not feasible.

## 2. Material and methods

### 2.1. Specimen identification

Specimens were taken from a pre-existing large dataset of CT scans and were selected in order to avoid any bias towards a particular group, skeletal element or body size, making the results of general application. Hindlimb bones from eight species of bird and mammal were acquired from various museum collections (National Museums Scotland, Edinburgh; Manchester Museum; and the World Museum, Liverpool). All specimens were skeletally mature (as determined by fusion of the epiphyses) and free of pathologies. When individual samples did not possess an associated *M*_{b}, typical values were assigned from the literature (table 1). External length measurements were taken using digital callipers (accurate to 0.1 mm), with the exception of the giraffe (*Giraffa camelopardalis*), which was measured with an anthropometer (accurate to 1 mm).

### 2.2. Finite element analysis

All but the largest specimen (*Giraffa*) were scanned in the Henry Moseley X-ray Imaging Facility, University of Manchester (X-Tek HMX 225 Custom Bay, Nikon Metrology Ltd, Tring, UK). Voxel size ranged between 64 and 119 μm, depending upon maximum bone length. Data were exported in unsigned 16-bit DICOM format (VG Studio Max v. 2.0, Volume Graphics, Heidelberg, Germany). *Giraffa* was scanned in a helical CT scanner at the University of Liverpool Small Animal Teaching Hospital (Siemens SOMATOM Volume, Erlangen, Germany) at a resolution of 391 μm and slice thickness of 3 mm, and reconstructed with Syngo (Siemens).

DICOM files were imported into OsiriX v. 3.8 [35], individually thresholded according to their greyscale values to accurately define the periosteal and endosteal surfaces, and surfaces exported as OBJ files, using the three-dimensional surface rendering function. Files were then imported into Geomagic Studio v. 12 (Geomagic, USA), the periosteal and endosteal contours isolated from one another, and exported as closed manifold OBJ files (available for download from the Dryad data repository: doi:10.5061/dryad.9ct2f). OBJ files were converted into SAT files, using Form·Z v. 6.1 (AutoDesSys, USA), and imported as solid parts into ABAQUS v. 6.10 (Simula, USA). In order to create a hollow bone part, the endosteal part was subtracted from the periosteal part, using a Boolean operation within ABAQUS. A homologous value for Young's modulus (stress/elastic strain) of 19 GPa and Poisson's ratio of 0.3 were assigned to all finite element models [36].

Hollow bone parts were meshed using a built-in Delaunay meshing algorithm, and for each bone, a one-way sensitivity analysis on the effect of changing element size was conducted. Total element number was progressively increased from approximately 200 000 elements to greater than 1 million elements and stress values were recorded at three locations at midshaft under a simple compressive loading regime. An optimal mesh size was considered to have been reached once stress values converged (i.e. formed a straight line on a plot of stress versus total element number) at all three locations. This occurred when element number exceeded roughly 800 000. However, convergence was reached at different mesh densities for each specimen, and the mesh sizes used for further analyses are detailed in the electronic supplementary material, S2.

A validation analysis carried out by Panagiotopoulou *et al*. [33] previously found four-node and eight-node tetrahedral, and mixed four-node tetrahedral and eight-node hexahedral FEA meshes to perform well, compared with *ex vivo* experimental data. By contrast, eight-node and 20-node hexahedral interpolations deviated significantly from recorded strain magnitudes. We, therefore, chose to carry out a comparison of C3D4 four-node linear tetrahedral and C3D10 10-node quadratic tetrahedral meshes. In mesh sizes beyond 200 000 elements, the difference in stress magnitudes between finite element models with 4- and 10-node tetrahedra was minimal (less than 5%; electronic supplementary material; figure S2). Meshes consisting of 10-node tetrahedra are computationally more expensive than those of 4-node tetrahedra, and C3D4 tetrahedra were, therefore, used throughout.

Each hollow bone model was loaded under combined compression and bending (0–90°), and axial torsion. For each loading regime, total applied load was calculated as 1 per cent of body mass (*M*_{b}; kg) multiplied by gravitational acceleration (*G*; 9.81 m s^{−2}; table 1). A force equivalent to 1 per cent of *M*_{b} was chosen in order to ensure that absolute strain magnitudes were small and deformation remained within the linear elastic region. The actual magnitude of the force is, therefore, largely unimportant in this study because there is a direct linear relationship between force and strain (stress).

For the combined compression–bending models, the condyles of the distal epiphyses were constrained in all three directions at 20 nodes on the distal surface using the ABAQUS ‘encastre’ boundary condition (figure 4). The applied force was spread across 10 adjacent nodes on the upper surface of the proximal epiphyses (figure 4). In order to minimize the extent of induced bending associated with off-axis application of compressive force, the orientation of the bones' principal axes was calculated using the ‘moments of inertia’ function of the BoneJ [37] plugin for ImageJ (US National Institutes for Health, MD, USA). To simulate compressive loading, the force was applied parallel to the longest principal axis, at nodes corresponding to the location at which the axis emerged onto the proximal epiphyseal surface. To simulate combined compressive and bending loading [29], FEA models for compressive loading were rerun while incrementally modifying the load vector from 10° to 90° from the principal axis.

To load the hollow bones under torsion, the condyles of the distal epiphyses were constrained in all three directions. A constraint control point (CP) was created on the proximal epiphyses, corresponding to the location at which the principal axis emerged onto the surface. The CP was constrained in all three directions, and a kinematic couple created between the load surface (defined as 10 nodes surrounding the CP) and the CP itself (figure 4). A torsional moment about the long axis was applied to the CP and transmitted to the load surfaces via the kinematic coupling. A linear elastic analysis was conducted on all models, and equations were solved using Gaussian elimination. Stress values used in this study were taken at a considerable distance from the constrained nodes. Stresses near constraints are known to be inaccurate in finite element modelling.

For bending and torsional loading regimes, the greatest value of principal stress (i.e. maximum principal stress, *σ*_{max}) was extracted from midshaft and used for comparison with beam theory. For models under compression, the most negative value of principal stress (minimum principal stress, *σ*_{min}) was recorded. For all loading conditions, the distribution of Von Mises stress (*σ*_{vm}) at midshaft was also noted. Von Mises stress combines the three principal stresses into one equivalent stress. *σ*_{vm} is signless and visually intuitive, and is, therefore, used in this study for illustrative purposes.

### 2.3. Simple beam theory

Basic morphometric properties were collected from the models to use in classic beam analysis. Cross-sectional geometrical properties of the hollow bone models were calculated at midshaft, using BoneJ. The coordinate system of both ABAQUS and ImageJ was in agreement, ensuring calculated values of radii and second moments of area corresponded to the load axis of the finite element models. Cross-sectional properties were entered into equation (1.3) to estimate maximum compressive and bending stresses at midshaft, under the same applied load as FEA models. Values of *t* (taken as mean cortical thickness), *t*_{min} and *A* required for equations (1.5) and (1.6) were calculated using a custom-written script in MATLAB v. 7.10 (The MathWorks Inc., Natick, MA, USA; see the electronic supplementary material, S1) before torsional stresses were estimated in the same manner.

Cross-sectional asymmetry was estimated as the ratio of two orthogonal measures of second moment of area (*I*_{max}/*I*_{min}). Normalized curvature lever arm (ζ) was calculated as the distance between the chord drawn between the proximal- and distal-most points of the epiphyses, and the location of the centroid at midshaft, divided by the radius (figure 3). Curvature was calculated in order to investigate the effect of irregular morphology on deviation between stress values predicted by FEA and beam theory. At this point, no attempt was made to correct the simple beam equations for curvature (but see §4 for further details). To test for relationships between bone morphology and deviation of FEA values from beam theory, type I ordinary least-squares regressions were carried out using the smatr package [38] of statistical software R (www.cran.r-project.org).

## 3. Results

### 3.1. Compression

When loaded in axial compression, all FEA models experienced significant levels of induced bending at midshaft. In figure 5*a*, a beam theory/FEA ratio of 1 suggests complete agreement between the models, whereas values greater than 1 indicate beam theory predictions fall consistently below those of FEA. At mid-cortex, beam theory and FEA are in agreement (figure 6*a*), whereas minimum principal stresses (*σ*_{min}) furthest from the neutral axis are on average 142 per cent (±28% s.e.) above predicted beam values. The model with the greatest normalized curvature lever arm (*Mustela putorius*, ζ = 1.63) also displays the greatest deviation of *σ*_{min} from predicted values (280%). Similarly, the model with the smallest degree of curvature (*Uria aalge,* ζ = 0.32) experienced the least deviation from predicted values (30%). There is a significant relationship between shaft curvature and variation between beam theory and FEA values (two-tailed *p* = 0.03, *r*^{2} = 0.56; figure 5*b*).

### 3.2. Bending

The location of maximum *σ*_{vm} (signless) varies between models, alternating between regions of the periosteal surface under maximum tension or compression (figure 7). Values of *σ*_{max} predicted by equation (1.3) are on average 12 per cent (±4% s.e.) different from those calculated in finite element models (figure 5*c*). For half the species modelled here, equation (1.3) underestimates *σ*_{max} compared with FEA models, whereas *σ*_{max} is overestimated in the remaining four species (figure 5*c*). No significant relationship is found between aspect ratio of the bone and deviation of FEA values from those predicted by beam theory. However, a significant relationship does exist between the asymmetry of the cross section and the ratio of FEA to beam theory values (two-tailed *p* = 0.02, *r*^{2} = 0.62; figure 5*d*).

Under Alexander's [29] model of combined compressive and bending loading, stress increases rapidly as the load vector shifts from parallel to perpendicular to the shaft long axis (figure 8). With the exception of *Giraffa*, the maximum deviation between FEA values and simple beam predictions occurs at *θ* = 90°, i.e. under compression (figure 8). As load vectors shift perpendicular to the bone long axis and bending dominates the loading regime, percentage deviation between FEA and beam theory is at minimum between 50° and 80°, and increases again towards 0° from vertical. A plot of *σ*_{combined} against bone orientation for FEA and simple beam models for all species is provided in figure 9.

### 3.3. Torsion

Comparing the values of *τ*_{torsion} predicted by equations (1.4)–(1.6) with those calculated using FEA, the hollow ellipse model (equation (1.4)) underestimates *τ*_{torsion} most (FEA values on average are 48% (±18% s.e.) above those of classic beam theory; figure 5*e*). The thin-walled ellipse model (equation (1.5)) and minimum wall thickness model (equation (1.6)) both provide reasonable estimates of *τ*_{torsion} (FEA values on average 21% (±16% s.e.) and 14% (±5% s.e.) above beam theory, respectively). It must be emphasized that the thin-walled ellipse model provides an estimate of *average τ*

_{torsion}in a cross section (using mean

*t*in equation (1.5)), while we are comparing these values with

*maximum*stress values extracted from FEA models. Applying either the hollow ellipse model or the thin-walled ellipse model, there is a significant relationship between cross-sectional asymmetry and variation between beam theory and FEA values (hollow ellipse model, two-tailed

*p*= 0.016,

*r*

^{2}= 0.65; thin-walled ellipse model, two-tailed

*p*= 0.004,

*r*

^{2}= 0.77; figure 5

*f*). However, when applying the minimum wall thickness model, no such relationship exists.

## 4. Discussion

When finite element models were loaded under axial compression, a significant bending stress resulted with values of *σ*_{min} averaging over twice as high as those predicted by beam theory. Previous *in vivo* studies have found that *σ*_{bending} contributes to total stress to a much greater degree than *σ*_{comp} during locomotion [39], in part owing to eccentric (off-axis) dynamic loading. The earlier-mentioned results highlight the additional importance of bending moments induced by on-axis forces acting about the longitudinal curvature of a static bone. Hence, pure compression as calculated by equation (1.1) is clearly an unrealistic loading modality for any long bone owing to shaft curvature. Here we show that a long bone, subjected to compressive axial loads, will probably fail, owing to the curvature-induced bending stress long before Euler buckling could be of concern. Therefore, it has been shown again here that in scaling studies there is no basis for the simplification of complex curved geometries down to idealized straight columns necessary to infer buckling as a viable failure mode. The extent to which finite element models deviate from simple beam predictions of compression is not constant, but is instead correlated with normalized curvature lever arm (*ζ*). This precludes the application of a single correction factor to account for shaft curvature and implies a measure of curvature must be incorporated into simple beam equations and calculations of ‘relative strength’ values, where previous studies have not done so [1,5,16,40].

This may be achieved simply by reusing equation (1.3) and assuming that a bone loaded in compression with a shaft curvature of 5° is equivalent to loading a straight bone at an angle of 5° from its longest principal axis (*θ*). Compressive loads are split into axial and transverse components as in equation (1.3)
4.1where is the angle between the chord drawn between the proximal- and distal-most points of the epiphyses, and the chord joining the proximal-most point with the centroid at midshaft (figure 3). Having subsequently corrected for curvature in our sample by applying equation (4.1), beam theory results overestimate *σ*_{comp} relative to FEA values in all but one instance, compared with the consistent underestimation of *σ*_{comp} when applying equation (1.1). In terms of percentage deviation between models, application of equation (4.1) brought the majority of beam theory predictions closer in-line with FEA results than when applying equation (1.1). As such, equation (4.1) provides a worst-case scenario estimate of *σ*_{comp} that is of particular interest when calculating safety factors of curved bones under compression.

The Euler–Bernoulli beam equation for estimating bending stress in cylinders (equations (1.2) and (1.3)) ignores potential transverse shear stresses, which may lead to an underestimation of the maximum stress in uniform circular cross sections of up to 12 per cent [41]. Indeed, the bending stress may even exceed normal stresses in some locations within long bones [42]. Here, we also find that FEA values and simple beam values differ on average by 12 per cent, with a maximum deviation of up to 35 per cent (*Giraffa*). Yet in half the models, beam formulae overestimated *σ*_{bending} relative to FEA values.

Interestingly, no relationship is found between aspect ratio of the bone and percentage deviation between the two estimates of *σ*_{bending}. Instead, this deviation appears to relate to the cross-sectional shape and the change in shape down-shaft. Those species with the greatest cross-sectional asymmetry (*Giraffa, Mustela, Erinaceus europaeus*) are found to have the greatest deviation between the two stress estimates (figure 5*c* and table 1). Therefore, in interspecific samples characterized by high morphological diversity, a measure of asymmetry should be incorporated into estimates of *σ*_{bending}. In engineering, it is standard practice to account for a lack of mirror symmetry by including an additional measure of product moments of area (*I _{xy}*)
4.2where

*M*is the bending moment about the

_{y}*y*-axis,

*x*is the perpendicular distance to the centroidal

*y*-axis and

*I*is the second moment of area about the

_{y}*y*-axis [43]. Inspecting the distribution of

*σ*

_{bending}values at midshaft across models (figure 7), it is clear that some bones experience a degree of twisting when loaded under bending. In

*Erinaceus*and

*Phoenicopterus ruber*, it is particularly noticeable that the neutral section (plane of zero stress where the bone is neither in tension nor in compression) has rotated away from the mediolateral axis. Such torquing of the bone under bending is not accounted for in equations (1.2), (1.3) or (4.2), and most likely results from irregularities in cross-sectional geometry away from the midshaft (i.e. trochanters, flanges, crests; figure 10). Unfortunately, there is no simple way of quantifying or accounting for down-shaft variation in geometry in classic beam theory, and FEA remains the preferred solution for bones of a complex shape.

The combined compression–bending model of equation (1.3) [29] is broadly supported by the results of our FEA models. The predicted increase in total stress as a result of shifting load vectors is present (figure 8). However, the combined model suffers from the compounded problems associated with equations (1.1) and (1.2). The decrease in effective stress as bones are orientated parallel to their load axes has been incorporated into the ‘effective mechanical advantage’ theory, by which safety factors in large mammals are maintained via adoption of increasingly erect postures during the stance phase [44]. The above results highlight the importance of incorporating curvature into such models in the future. As an example, assuming that an arbitrary value of 1000 kPa must be maintained in order to achieve a given safety factor in the tibia of *Giraffa*, equation (1.3) would predict the bone may be held at a maximum of 53° from vertical (figure 9*a*). By contrast, the stress curve produced under FEA predicts the tibia must not exceed 36° from vertical in order to maintain the same safety factor (figure 9*b*).

In agreement with results found elsewhere [32], the minimum wall thickness model is found to be most suitable in estimating *τ*_{torsion}. Equation (1.6) recognizes that minimum torsional strength occurs in regions where cortical thickness is least [30] and is supported by the distribution of cortical stresses noted in finite element models (figure 6*b,c*). While estimates of *τ*_{torsion} based on the hollow ellipse or thin-walled ellipse are still affected by the degree of asymmetry present in the cross section, no such relationship exists when applying the minimum wall thickness model, and equation (1.6) provides a reasonable estimate of *τ*_{torsion} as calculated by FEA models.

Figure 6*b*,*c* clearly illustrates the divergence of complex bone cross sections from the idealized toroidal cross sections familiar to engineers. Classic beam theory is a technique that rests upon certain geometrical assumptions that are clearly invalidated in the case of irregular long bone morphology, and the consequences are evident in the discrepancy between FEA and beam theory presented here. Both classic beam theory and FEA allow the biomechanical behaviour of long bones to be modelled. As such, the results outlined earlier afford a methodological comparison of two models of reality; however, they are unlikely to reflect *in vivo* loaded behaviour, particularly as both FEA and simple beam models assume uniform elasticity and material density. The results of a previous study, in which an FEA model of an elephant femur was validated by means of laser speckle interferometry [33], found strain magnitudes predicted by FEA models consisting of homogeneous isotropic material properties to differ 60 per cent from *ex vivo* experimental values. By contrast, FEA models incorporating heterogeneous material properties deviated only 5 per cent from experimental values. Incorporating the anisotropic behaviour of bone into FEA models has also been found to improve their accuracy in predicting regions of fracture [45], although further work is needed to clarify the errors associated with assumptions of bone isotropy [46].

Furthermore, the loading conditions considered in the present study are highly idealized, in that forces are applied only at articular surfaces. In reality, peak locomotory forces are generated by contraction of muscles attached at various locations along the length of the shaft, and may be a considerable distance from the centre of joint rotation. Although beyond the scope of this study, future work may focus upon comparing numerical models to *in vivo* strain gauge derived values of stress under average and extreme locomotor activity. Such a study would be advantageous in incorporating realistic muscle forces, and could illuminate the relative importance of compression, bending and torsion across taxa and locomotor type.

## 5. Conclusions

For interspecific samples of diverse morphology, reliance upon Euler–Bernoulli classic beam theory to address questions of comparative functional morphology in vertebrate bones will result in estimation errors of varying magnitude and direction. The utility of commonly used beam formulae is a function of the extent to which skeletal elements conform to the assumptions of classic beam theory. While beam-theory-derived values still probably correlate with actual stress values, care must be taken when interpreting mechanical function based on these values. In the small sample of vertebrate bones considered here, application of FEA leads to a change in the rank order of absolute stress values across species compared with beam theory predictions, most noticeably so in the case of compression. It is therefore of concern that real biomechanical signals within the sample may be lost, owing to the use of inappropriate beam formulae. Suggestions are made for alternative methods by which beam theory may be more reliably applied, such as the incorporation of curvature in estimates of *σ*_{comp}, including product moments of area into estimates of *σ*_{bending} and calculating *τ*_{torsion} at the location of minimum cortical wall thickness.

In spite of these extensions of beam theory, when absolute values of diaphyseal stress are required for biomechanical analyses, FEA remains the preferred solution. With improvements in computational power, user accessibility and the availability of CT facilities, it is now feasible to generate large datasets of finite element models for the purpose of comparative functional morphology. FEA benefits from the incorporation of whole-bone geometry into models and overcomes problems associated with curvature, asymmetry and shear deformations. With interest in longitudinal variations in bone strength indices increasing [47,48], FEA represents a technique by which muscle attachment sites and trochanters may be studied, when their irregular geometry might otherwise preclude the application of classic beam theory.

Future studies applying FEA to long bone stress estimation should proceed with caution, however, particularly when variables such as applied forces and material properties remain uncertain. This is necessarily the case in palaeontological studies, and, therefore, sensitivity analyses should be carried out in order to quantify the effect of the error introduced by these unknowns. It must be emphasized that applying an overly simplified FEA model to a complex biomechanical problem may result in incorporating just as many assumptions into the analysis as the application of classic beam theory.

## Acknowledgements

The authors are grateful to the staff of the Henry Moseley X-ray Imaging Facility (supported by EPSRC under nos. EP/F007906 and EP/I02249X), University of Manchester; Prof. Paul Mummery, Dr Tristan Lowe; the staff in Research Computing Services, University of Manchester; the University of Liverpool Small Animal Hospital: Martin Baker; Manchester Museum: Henry McGhie, Judith White; Liverpool World Museum: Tony Parker. The comments of four anonymous reviewers were extremely useful in refining the manuscript. This research was supported by the UK Natural Environment Research Council Doctoral Training Grant no. NE/1528134/1.

- Received October 8, 2012.
- Accepted October 31, 2012.

- © 2012 The Author(s) Published by the Royal Society. All rights reserved.