## Abstract

Finite element modelling is well entrenched in comparative vertebrate biomechanics as a tool to assess the mechanical design of skeletal structures and to better comprehend the complex interaction of their form–function relationships. But what makes a reliable subject-specific finite element model? To approach this question, we here present a set of convergence and sensitivity analyses and a validation study as an example, for finite element analysis (FEA) in general, of ways to ensure a reliable model. We detail how choices of element size, type and material properties in FEA influence the results of simulations. We also present an empirical model for estimating heterogeneous material properties throughout an elephant femur (but of broad applicability to FEA). We then use an *ex vivo* experimental validation test of a cadaveric femur to check our FEA results and find that the heterogeneous model matches the experimental results extremely well, and far better than the homogeneous model. We emphasize how considering heterogeneous material properties in FEA may be critical, so this should become standard practice in comparative FEA studies along with convergence analyses, consideration of element size, type and experimental validation. These steps may be required to obtain accurate models and derive reliable conclusions from them.

## 1. Introduction

Although finite element analysis (FEA) is a computational tool to quantify deformation in complex structures such as bones, it has been widely used to infer function from the form of extant and even extinct taxa. But how complex does a finite element (FE) model need to be in order to address hypotheses with complex form (such as heterogeneous bones) and function (such as locomotion)? Sources and magnitudes of error are a continuing concern in applying the FEA method to skeletal elements because of their complex geometry, heterogeneous materials and assumptions concerning loading boundary conditions. Simplifications and assumptions with regards to these variables must necessarily be made in order to make a model computationally tractable. Nonetheless, the quantitative assessment of the degree of simplification and the discretization error of an FE model using convergence tests, sensitivity studies and validation analyses is an essential prerequisite in order to ensure model reliability. Convergence tests are conducted by decreasing the element size until the results are repeatedly obtained. Sensitivity studies indicate how much the results change in response to small perturbations in certain parameters, such as material properties or loading conditions. Validation analyses compare FEA results with experimental results to maximize model accuracy (e.g. [1]).

A convergence study involves varying the size and type of elements in the FE model's mesh. The resolution of the mesh (size and type of elements) affects the results of the analysis. As the mesh resolution increases (i.e. smaller elements), the solution of the model converges. After convergence is reached, more elements will not change the solution. Thus, it is critical that FE models have a sufficient number of elements and/or the appropriate type of elements, so that the FEA results do not change between models and the model converges. At this convergence point further increases to the mesh density will result in the same solution. Once an FE model is determined to be converged, it can be compared against experimental (validation) results to assess its overall accuracy.

Convergence studies are commonly conducted in mechanical engineering studies [2,3] and in clinical studies of human bone [4–9], yet are widely overlooked in zoological applications. Many zoological studies involve arbitrary selection of element sizes or types for their FE model, sometimes apparently just the default values in the mesh generation software. Such practice, coupled with simplifications in the assignment of material properties such as modelling bones as homogeneous, can lead to inaccurate conclusions in studies such as fracture mechanics, safety factors and remodelling fields, in which quantitative accuracy is critical for obtaining specific, rigorous results.

Ethical limitations may prohibit invasive experiments when planning *in vivo* validation, so an *ex vivo* validation is often the only viable alternative. *Ex vivo* experimentation involves the application of predefined loading boundary conditions to a cadaver that can be replicated in the FE model, and preferably should be conducted on the same specimen. Limitations in terms of accessibility of animal cadaveric material and computed tomography scans, particularly when the model animal is a rare and large specimen such as an elephant, leave two options. One may either not conduct validation, assuming that the degree of error between the FE model and the experiment is minimal (i.e. assume that the model is reliable enough to draw conclusions from its quantitative results), or to conduct a validation analysis using a different specimen of similar age and geometry. A critical caveat here is that unvalidated models may have critical flaws that could cause unpredictable errors in the accuracy and reliability of biomechanical models [1].

In this study, we seek to construct the first reliable FE model of the femur of an adult Asian elephant (*Elephas maximus* Linnaeus 1758) in order to move towards a better understanding of its form–function relationship and re-examine previous ideas about bone strength, scaling and athleticism in these giant land animals [10,11]. We first test our model's precision by conducting a series of convergence tests with regards to the size and type of elements used in the geometric mesh. We then perform a sensitivity study on the model by determining the effects of varying its material properties. Finally, we compare our FEA results against the experimental data obtained under prescribed loading conditions in the laboratory (using laser speckle interferometry). We hope that our study will raise the methodological standards in zoological applications of FEA by encouraging the inclusion of convergence tests, sensitivity studies and validation analyses to ensure model reliability.

## 2. Materials and methods

Our finite element analysis was based on the left femur of an adult skeletally mature female Asian elephant (*E. maximus*; 3432 kg; 24 years old; 95 cm femur length measured between proximal and distal articular surfaces) (subject no. 1). The specimen was euthanased, for reasons unrelated to this study, at Whipsnade Zoo (Bedfordshire, UK). The elephant femur was frozen (−20°C) prior to our experiments. Our *ex vivo* experiment using laser speckle interferometry was conducted on the left femur of a different adult skeletal mature female Asian elephant (3920 kg; 40 years old; 98 cm femur length) (subject no. 2) from Paignton Zoo (Devon, UK). Ideally, the validation study should have been conducted on the same femur with the one used in FEA. We performed the *ex vivo* experiment on a different elephant to the FEA specimen because subject no. 1 was previously destructively analysed for mechanical testing of bone material properties [12]. We used BoneJ [13] software to quantify structural strength variations between the two elephant femora by quantifying the polar section modulus (Zpol) (figure 1 and table 1). The polar section modulus is a geometric measure of how strong an element is in torsion, specifically the polar moment of area divided by the maximum radius from the centroid. We normalized the polar section modulus by the length of the femur and compared values at every 5 per cent of total bone length. The femur used in our FEA (subject no. 1) has lower Zpol values than the femur used in our *ex vivo* experiment (subject no. 2) in the proximal epiphysis and the mid-diaphysis. Decreased Zpol values in these areas are an indication of the decreased regional strength of subject no. 1 compared with subject no. 2. Nevertheless, there is only a minimum (approx. 4%) difference between the Zpol values around the proximal femur (30% length)—the area sampled for this study's validation test. Thus, it is reasonable to infer that the effect of this difference on the FEA results is minimal.

### 2.1. Finite-element analysis

FEA is a numerical technique for solving the function of independent variables and their partial derivatives using linear ordinary differential equations. With FEA, a continuum structural system is divided into discrete sub-regions of finite size, the elements. Elements are systems with linear algebraic equations that model their response to loads and are bonded with each other at their vertices, called nodes. Each node has a displacement vector which represents 3 d.f. corresponding to displacements in the *x*, *y*, *z* directions. The number of nodes in an element of a given number of sides is determined by the order of interpolation used to interpolate the field variable (e.g. displacement) between nodes. More information on FEA is available in earlier reviews [14–16]. The three main steps involved in FEA are: model creation (figure 2*a–g*), model solution (figure 2*h*) and model validation (figure 2*i*).

### 2.2. Model creation

The creation phase involves the assignment of the model's geometry, transformation of the geometric model into a set of discrete finite elements, and the assignment of its material properties, loading and constraints. The skeletal geometry of the femur of subject no. 1 was CT-scanned on a Picker International Inc./PQ 5000 scanner (Cleveland, OH) at the Royal Veterinary College (100 mA, 120 kV, slice thickness 3 mm, 1.347 pixels mm^{−1}, 512 × 512 pixel images). The femur was scanned along with calibration phantoms (Gammex, Inc.) of known radiodensity (*ρ*) and apparent density (*D*_{app}). The femur CT dataset consisted of 193 slices with an isometric slice thickness of 0.742 mm. This dataset was processed to extract three-dimensional volume datasets that were converted into meshed files of elements using Simpleware image processing software (ScanIP v. 3.2; Simpleware Ltd, Exeter, UK).

### 2.3. Convergence analyses: element size and type

After segmentation in ScanIP, masks were exported to ScanFE software (ScanFE v. 3.2; Simpleware Ltd, Exeter, UK) to create finite-element models using the FE Grid algorithm. We generated mesh files of different uniform element sizes ranging from 3 to 7 mm and four different combinations of element types (four-node tetrahedra, 10-node quadratic tetrahedra, a mixture of four-node tetrahedra and eight-node hexahedra, eight-node hexahedra and 20-node hexahedra). A linear tetrahedral element has four nodes (one at every corner) and field variables within the element are interpolated from the nodal values using linear functions. A quadratic tetrahedral element has six extra nodes, one located at each of the six edges, in addition to the four at the corners of the element. The quadratic shape functions allow the element to have both quadratic (i.e. curved) geometry as well as quadratic displacement fields. The quadratic displacement field enables the element to more accurately respond to bending loads. Eight-node hexahedral or brick elements have eight corners (one node at every corner) and six faces. Twenty-node hexahedra have 12 additional side nodes. The mixed mesh file was composed of four-node linear tetrahedra elements to get smooth faceted surfaces and eight-node linear hexahedra elements for thick sections and flat surfaces.

### 2.4. Sensitivity study: material properties

We constructed an empirical model of Young's modulus (*E*) of the femur of subject no. 1 as a function of the apparent density, *D*_{app}, which is itself inferred from the radiodensity in Hounsfield units (HU), *ρ*, obtained via the CT scan data. The relationship between the radiodensity and the apparent density described below is scanner-specific and cannot be used to infer apparent densities from radiodensities in other scanners. Nor can the relationship between radiodensity and Young's modulus be applied to all species or specimens. Researchers who wish to use the empirical model described below should determine their own relationship between radiodensity and apparent density and validate their FE models if the presented relationship between apparent density and Young's modulus, which depends on the species, type of bone and sample location, is used. We modelled the relationship between radiodensity in HU and apparent density using CT phantom data (calibration data of known radiodensity to apparent density). We quantified the relationship between apparent density and radiodensity using a linear model comprising two segments (electronic supplementary material, figure S1; also appendix I). Such a model has four parameters: *K*_{1}, *K*_{2}, *C*_{1} and *C*_{2} that were estimated via least-squares fitting (using Matlab's fminsearch, Nelder–Mead simplex method, tolerance (10^{−4} HU)).

Having obtained a relationship between radiodensity and apparent density, we then estimated the Young moduli (electronic supplementary material, figure S2) for the femur. In Zioupos *et al*. [12], the Young moduli values of the same elephant femur (subject no. 1) were measured along three orthogonal directions (*E*_{x}*, E*_{y} and *E*_{z}). In addition the apparent and material densities were measured. Using Zioupos *et al*.'s [12] Young moduli data, we constructed an effective Young modulus (see electronic supplementary material, appendix I). We then modelled the relationship between the apparent density and the effective Young modulus by fitting a power law with three parameters (*α*, *β* and *γ* defined in the electronic supplementary material, appendix I and figure S2).

These procedures established our empirical model to calculate heterogeneous and isotropic Young modulus values from radiodensity. We then applied our model to the FE model of the femur of subject no. 1. We applied defined loads to the FE model (see below) and then analysed the difference in surface strain magnitudes, patterns, ratios and orientations between the heterogeneous model and a homogeneous model of the same femur with *E* = 19 GPa and *v* = 0.3. The homogeneous and isotropic value of 19 GPa for Young's modulus is a volumetric average of the heterogeneous values assigned to the heterogeneous FE model. We compared both the heterogeneous and homogeneous FE models to strain results derived from the *ex vivo* experiment (below) under the same loading conditions.

### 2.5. Loads and constraints

In the convergence tests, we applied a static compressive force of 5000 N magnitude at the most medial node of the femoral head to all FE models with various element size and type. All external nodes on the femoral condyles of the distal epiphysis were constrained in all three directions.

In order to validate the model with experimental tests, all FE models with various material properties and element types were loaded in a non-physiological manner to simulate the loading conditions of the *ex vivo* experiment (see below). To prevent rigid motion, all external nodes on the femoral condyles of the distal epiphysis were constrained in all three directions.

### 2.6. Model solution

We solved all FE models using a direct linear equation solver in ABAQUS/CAE 6.9-EF1 software (Dassault Systèmes Simulia Corp., Providence, RI, USA). The solver uses a sparse, Gauss elimination method to solve a large system *KU* = *F*, where *K* is the stiffness matrix, *F* is the external force with constraints and *U* is the final displacement.

### 2.7. *Ex vivo* validation using laser speckle interferometry

We obtained strain measurements for the *ex vivo* experiment using full-field laser speckle interferometry (Dantec Q-100. Dantec Dynamics, Ulm, Germany). The laser speckle interferometer is a non-contact system for measuring component displacement from which principal strains can be measured [17–21]. Prior to measurements, we treated the experimental specimen's surface for strain measurement with isopropanol and abraded an area of 25 × 33 mm^{2} with pumice powder, covering the most caudodistal aspect of the proximal femur, halfway between the femoral head and the greater trochanter, to remove any remaining fascia. We attached the ring of the Dantec Q100 to the bone surface using a glue mixture. We then positioned the femur in a mechanical testing machine (Lloyds Instrument Materials Testing System EZ50; electronic supplementary material, figure S3) and loads were applied on the most proximomedial aspect of the greater trochanter (electronic supplementary material, figure S3). The distal femoral epiphysis was fixed in a frame and constrained along all three anatomical axes (electronic supplementary material, figure S3). We then attached the Dantec sensor to the Dantec ring, the area of interest was laser-illuminated and a reference image of the sampling field in all three directions was taken. We checked the measurement conditions by applying a small load to the greater trochanter to assess the quality and stability of fringes and to calculate phase maps of the area of interest. We then sequentially applied compressive loads of 5000 and 10 000 N to the femur. We computed the three-dimensional strain component measurements qualitatively and quantitatively using IstraQ 100 2.7 software (Dantec Dynamics, Ulm, Germany).

The FE model's validity was checked by comparing strain magnitudes, ratios and orientations of the FE models of subject no. 1 with various material properties and element types with those of the *ex vivo* femur (subject no. 2) under the same loading conditions. We did not expect exact matches between the strain magnitudes and ratios for the *ex vivo* experiment and the valid FE model, not least because of the assumptions of isotropy used to model the material properties of the femoral bone in FEA. Differences in strain magnitudes were also expected owing to intraspecific variation because the FE femur (subject no. 1) and the *ex vivo* femur (subject no. 2) used for the experiment were from different individuals and thus presumably have variations in material properties and bone geometry (see below). Also, although both femora come from skeletally mature adult individuals, their age and size differences (in addition to many other aspects of individual variation) could also influence their bone material properties and thus give different strain magnitudes [22]. However, both FE models were expected to give similar strain patterns and orientations when compared with the experimental data.

For comparisons with our *ex vivo* experiment, we computed the mean maximum principal strain (*ɛ*_{1}) and mean minimum principal strain (*ɛ*_{3}) magnitudes and orientations from the centre of 60 surface nodes at 13 homologous locations on the proximal femur, covering a surface area of approximately 33 mm^{2} at each location (electronic supplementary material, figure S4). The limited field of view of the laser interferometer meant it was not possible to experimentally estimate strain magnitudes from the entire surface of the elephant femur. The sampling region was as large as possible and included areas expected to have substantial strain variations (i.e. localized changes in cortical thickness and subcortical geometry).

Pair-wise comparisons between the *ex vivo* experiment and each of the FE models under sequential axial loads of 5000 and 10 000 N magnitudes were conducted by calculating the normalized Euclidean distance (the square root of the sum of the square differences), root mean square error and standard deviations of the residuals for the strain values. To estimate whether the FE models and the experiment deformed in a similar pattern; i.e. whether at a particular location, the FE models and the experiment experienced compression or tension; we also computed strain ratios.

The interpretation of our validation test is based on the assumption that if surface strains and strain ratios in the sampling region *ex vivo* match those of the FE model under approximately identical conditions, then the model validation test is successful, building confidence in the assumption that strain values elsewhere in the femur would have accurate matches between actual and FE model results.

## 3. Results

### 3.1. Convergence analyses: element size and type

For our convergence analyses, we varied element size within the range of mesh densities chosen for the femur of subject no. 1 (figure 3). The FE models with element sizes 3–5 mm had similar strain magnitudes, forming a horizontal straight line on figure 3 and thus converged to the same solution. The FE mesh files with element sizes of 6 and 7 mm display greater strain magnitudes than the rest of the FE models. The FE mesh file with an element size of 4 mm and 425 934 total elements was used for further analysis because it produced results similar to the 3 mm model FE solution, yet required less computational resources. The computational times needed for the solution of the 3 and 4 mm element size models using a Dell core 2 duo quad-core extreme edition processor were (in seconds): 3 mm—user time 5064; system time 146; total CPU time 5210.9; 4 mm—user time 891; system time 6.7; total CPU time 897.5.

We also ensured that models with different types of elements converged. Figure 4 shows that the strain results for the tetrahedral and the mixed tetrahedral and hexahedral element types were very similar (approx. 3% strain difference) and form a line. Hexahedral elements increase strain magnitudes by 20 per cent. To determine which element type gives more accurate strain results for the scope of our study, we compared the FE models with the element types that converged (tetrahedral and the mixed tetrahedral and hexahedral element types) against the *ex vivo* experiment. Strain results from the *ex vivo* experiment and the FE models are presented in figure 5 and table 2, depicting alterations of element type and material properties under comparable loading of 10 000 N. Pair-wise comparisons between each FE model and the experiment are also expressed in terms of root mean square errors, Euclidean distance and the standard deviation of the residuals (table 3). Slopes of regressions and *R*^{2} values between the experiment and the FE models with various element types and material properties are, respectively, displayed in figures 6 and 7. Differences in strain magnitudes between the *ex vivo* experiment and the FE models under a load of 10 000 N magnitude were also expressed in terms of percentage differences (figure 8). Strain magnitudes of both the FE heterogeneous model with four-node tetrahedra and eight-node hexahedra element type and the *ex vivo* experiment decreased proportionately, as expected, when we reduced the applied force to 5000 N magnitude (electronic supplementary material, table S1).

## 4. Discussion

In this study, we assessed the effect of various element sizes, types and material properties on the reliability of biomechanical models of an elephant femur using a combination of subject-specific finite element models, convergence analyses, sensitivity tests and *ex vivo* experimentation. Our convergence analyses that varied elements' size and type showed that FE meshed files with 3–5 mm elements and tetrahedral and a mixture of tetrahedral and hexahedral interpolations converged to the same solution.

Our validation study showed that the FE meshed files with four-node linear tetrahedral interpolations and a mixture of four-node linear tetrahedral and eight-node linear hexahedral interpolations displayed the minimum error (approx. 7% difference in strain magnitudes) when compared against the *ex vivo* experiment in terms of root mean square error, Euclidean distance and the standard deviation of the residuals, as well as slopes and *R*^{2} values of regression. The 20-node and eight-node hexahedral interpolations displayed different strain magnitudes from the experimental strains by approximately 100 per cent. Hexahedral interpolations may produce FE solutions with minimal error in strain magnitudes compared with tetrahedral interpolations owing to their anisotropy and angle regularity [23,24], but are not suitable for surface meshing because they are not smooth enough, unless the available mesher can produce very high resolution meshes.

For the assignment of heterogeneous material properties in our FE model, we constructed an empirical model (i.e. mathematical equation) to estimate Young's modulus of the femur of subject no. 1 from radiodensity. We calculated the relationship of Young's modulus as a function of the apparent density, which itself was inferred from the radiodensity obtained via the calibrated CT scans. We developed our model from empirically determined material properties, which had been previously measured using destructive mechanical testing [12] for the same specimen. Relatively few studies have conducted validation of an empirical model against experimental data measured with mechanical testing, and none have done so for very large skeletal structures such as an elephant femur. In the event that material properties can be determined through mechanical testing of the same specimen, it is of course preferable to use those rather than any empirical model. Nevertheless, in the eventuality that destruction of samples is either undesirable or impossible owing to ethical/experimental constraints, a validated empirical model as presented in this study has inherent utility. If it is assumed that our femur specimen's relationships of material properties with density are reasonably representative of elephant bone, then our model would be applicable to at least other adult elephant femora and potentially any elephant limb bones. Care must be taken to ensure that any effective Young's modulus is representative of the material properties of the region of interest.

To test the sensitivity of heterogeneity and homogeneity in an FE simulation, we compared our FE model strain results with *ex vivo* experimental strain data under the same loading conditions. Our resulting validation test is novel, in addition to its first application to an elephant's long bone, in that the elephant femur enabled us to directly determine material properties in different regions (electronic supplementary material, figure S2), a procedure that cannot be applied to a mouse bone and is rarely done for human bones. Our validation test showed a high degree of congruence between the *ex vivo* experiment and the FE models for strain patterns. The heterogeneous model produced strain magnitudes of only 5 per cent difference from the experimental strain results for the proximal femur, which had relatively identical geometry between the two specimens in this study. On the contrary, the homogeneous FE model typically displayed a 60 per cent difference in strain magnitudes when compared with the experimental strain values. Thus, we judge the heterogeneous model to be surprisingly accurate because its validation test revealed approximately 5 per cent accuracy, whereas the homogeneous model was surprisingly inaccurate in having over an order of magnitude greater error in strain estimation.

The results of our and other validation studies [25] have shown that homogeneous models produce relatively similar strain patterns as heterogeneous models and can be indeed assumed for certain questions (e.g. simple comparisons of strain patterns and/or magnitudes) without producing unrealistic and biologically inaccurate FE results. Nevertheless, when compared with experimental data, heterogeneous models produce smaller errors in strain magnitudes (and ratios) than homogeneous models. Thus, to address questions on fracture mechanics, bone remodelling fields, safety factors and failure criteria, when absolute stress or strain magnitudes are needed, heterogeneous FE models are essential. The means to assign heterogeneous material properties along the orthogonal axes using either experimental techniques and/or empirical models are available, hence heterogeneity should not be taken lightly. In some cases, it may not be clear *a priori* whether heterogeneous material properties are essential or not, so future analyses should quantitatively test the assumption that homogeneity does not affect the reliability of their conclusions.

In this study, we combined experimental loading, computer simulation techniques, convergence studies and empirical models of material properties to test the reliability of an FE model of the femur of an adult Asian elephant. Validation is an essential procedure through which we can assess the inaccuracies of our computer models and should be ideally conducted on the same skeletal specimen with the one used in FEA. There will be cases when access to the same skeletal specimen is not possible (as per this current study). Such cases though should not be used as an excuse to not conduct a validation analysis but researchers must quantify how different the two specimens are and cautiously consider how such differences can affect their ultimate conclusions. In our case, the strain results between the experiment and the heterogeneous FE model displayed less than a 5 per cent difference. Bearing in mind that our experimental and the FE models are geometrically very similar in terms of polar section modulus in the region compared, similarities in strain magnitudes suggest that individual variations in morphology and size between the two elephant femora had little effect on measured and predicted strains. There will always be some mismatch between experimental data and FE models because the latter can only approximate reality (additionally, the former is seldom if ever free of its own errors). When the error in strain magnitudes is small, as in our case, then we can securely assume that our FE model does not display any critical flaws, at least within the sampling area, and therefore should be sufficiently accurate to draw reliable conclusions from it.

To this extent, we hope that our study can set an example for future biological/zoological studies to follow in order to develop FE models with sufficient accuracy and resolution to address the biological questions of interest. We urge that convergence and sensitivity analyses with regards to element types, use of accurate approaches to estimate material properties from radiodensity (where possible), comparisons of the assumptions of heterogeneity versus homogeneity, experimental validation and similar cautious approaches to model formulation and validation be critically and regularly considered in comparative FE studies. An adherence to these principles would not only create a foundation of stronger conclusions on which future studies can continue to build, but also bolster confidence in FE modelling studies in the broader biomechanical community.

## Acknowledgements

Thanks to Michael Fagan University of Hull for allowing access to the mechanical testing machine and the laser speckle interferometer. We thank Peter Njuguna, Michael Doube, Sue Taft and Andy Greenhalgh for technical support and Peter Zioupos and Callum Ross for useful discussions. Thanks also to Michael Doube and Stephanie Pierce for their help in the use of Bone J. We thank Whipsnade Zoo and Paignton Zoo for supplying cadaveric material. This study was funded by a grant from the BBSRC to J.R.H.

- Received May 24, 2011.
- Accepted June 22, 2011.

- This journal is © 2011 The Royal Society