The Tarsometatarsus of the Ostrich Struthio camelus: Anatomy, Bone Densities, and Structural Mechanics

Background The ostrich Struthio camelus reaches the highest speeds of any extant biped, and has been an extraordinary subject for studies of soft-tissue anatomy and dynamics of locomotion. An elongate tarsometatarsus in adult ostriches contributes to their speed. The internal osteology of the tarsometatarsus, and its mechanical response to forces of running, are potentially revealing about ostrich foot function. Methods/Principal Findings Computed tomography (CT) reveals anatomy and bone densities in tarsometatarsi of an adult and a young juvenile ostrich. A finite element (FE) model for the adult was constructed with properties of compact and cancellous bone where these respective tissues predominate in the original specimen. The model was subjected to a quasi-static analysis under the midstance ground reaction and muscular forces of a fast run. Anatomy–Metatarsals are divided proximally and distally and unify around a single internal cavity in most adult tarsometatarsus shafts, but the juvenile retains an internal three-part division of metatarsals throughout the element. The juvenile has a sparsely ossified hypotarsus for insertion of the m. fibularis longus, as part of a proximally separate third metatarsal. Bone is denser in all regions of the adult tarsometatarsus, with cancellous bone concentrated at proximal and distal articulations, and highly dense compact bone throughout the shaft. Biomechanics–FE simulations show stress and strain are much greater at midshaft than at force applications, suggesting that shaft bending is the most important stressor of the tarsometatarsus. Contraction of digital flexors, inducing a posterior force at the TMT distal condyles, likely reduces buildup of tensile stresses in the bone by inducing compression at these locations, and counteracts bending loads. Safety factors are high for von Mises stress, consistent with faster running speeds known for ostriches. Conclusions/Significance High safety factors suggest that bone densities and anatomy of the ostrich tarsometatarsus confer strength for selectively critical activities, such as fleeing and kicking predators. Anatomical results and FE modeling of the ostrich tarsometatarsus are a useful baseline for testing the structure’s capabilities and constraints for locomotion, through ontogeny and the full step cycle. With this foundation, future analyses can incorporate behaviorally realistic strain rates and distal joint forces, experimental validation, and proximal elements of the ostrich hind limb.


Introduction
The ostrich (Struthio camelus) is the largest and fastest extant ratite, with great capacity for long distance locomotion [1,2,3,4]. The morphology of the ostrich hindlimb has been the subject of numerous studies [1,5,6,2,4,3,7]. High-powered muscle output is channeled through a multi-jointed system of interconnected bi-and tri-articular muscles [6,3] (Table 1). The distal tibiotarsus and tarsometatarsus are elongated and lightened to provide more efficient locomotor capability, with the muscle forces transmitted to the pes via long tendons [8,9,7]. Distal limb length is further extended by elevating the metatarsophalangeal joint above ground [7]. This articulation enhances elastic storage shock absorption during locomotion [2,4]. Ostriches have further reduced distal limb mass by eliminating the second digit and associated musculature, and are thus the only known didactyl birds [3].
The internal anatomy of the ostrich tarsometatarsus, and its distributions of compact and cancellous bone, have yet to be illustrated in detail or assessed in relation to their structural response during locomotion. We describe the external and internal anatomy of the ostrich tarsometatarsus based on specimen observations, reviews of the literature, and computed tomographic (CT) scans. Pertinent to ostrich biomechanics, we illustrate mineral densities of adult and juvenile tarsometatarsi, and present a finite element (FE) model of an adult's tarsometatarsus, which is useful for simulating its response to external loads. Our goals are 1) to present osseous anatomy and densities of the tarsometatarsus; 2) to relate anatomy and densities to stresses and strains under a sample loading regime for the adult, and 3) to identify refinements for locomotor finite element analysis (FEA) of individual limb bones. Our goal is not to claim definitive anatomical/functional correlations, when the accuracy of the most thorough models of ostrich locomotion is unknown [10], even when carefully compared to in vivo experimental data. Instead, the anatomical data and FE model will be useful for future analyses of structural response through the step cycle (during both running and turning: [10]). Forces for such studies will be derived from experimental data and motion simulation through methods of multibody kinematics and dynamics [5,11,2,4,12,13]. The current study is a first and necessarily simplified stage of correlating anatomy with FEA, which will ultimately combine simultaneous multibody dynamics (MBD) and FEA, as Snively, Kumbhar et al. (2013) [14] have applied to chewing pigs. The current ostrich loading regime is a "snapshot" (single time-increment) from manual MBD calculations, which are presented here fully for replication and as a guide to coding in programs such as MATLAB and Mathematica. This approach is a transparent complement to off-the-shelf MBD programs, including MSC Adams and the open-source GaitSym and OpenSim, which researchers can use as efficient black-box solutions for calculating muscle and reaction forces in many poses.

Interpretation of CT bone densities in juvenile and adult Struthio camelus
We compare densities of juvenile and adult ostrich tarsometatarsi to assess changes in ontogenetic and locomotor adaptation. Ostriches grow quickly, with males growing from 1 to 100 kg within a year as their tarsometatarsi lengthen approximately 8-fold [15]. Their limb bones experience shifts in histology and presumably in density, due to conditions necessary for fast growth in juveniles (including calcified cartilage interdigitated with long marrow tubes for bone deposition: [16,17]) to denser, stiffer elements essential for rapid locomotion in adults. CT scans enable us to assess densities of mineralized tissues. These densities are directly proportional to stiffness ( [18,19,20], and thus resistance to deformation from locomotor forces. Densities are often expressed as Hounsfield units (HU) [19], indexed from x-ray attenuation values, which enable comparison of tissue density to that of water. A Hounsfield value of 0 is that of water, corresponding to a density of 1 g/cm 3 . For air the HU value is -1,000, and bone ranges from 700 HU for dense cancellous to 3000 HU for highly dense compact bone.
We scanned juvenile and adult tarsometatarsi on medical scanners (see Materials and Methods), and further scanned proximal and distal ends of the juvenile specimen on a micro-CT scanner. Most medical-resolution scanners differentiate a 12-bit range of densities, 4,096 Hounsfield units wide, and the micro-CT has 16-bit resolution of 65,536 density levels (not indexed as traditional HU). Because the juvenile TMT had a narrow range of density, the micro-CT scan revealed gradations of its density not possible with the medical scanner. To visualize these gradations and highlight density variation, we use full-color palettes ( [21,22]) Table 1. Material properties assigned to FE model of the ostrich tarsometatarsus, and yield and ultimate values for comparison with FEA results. E = Young's modulus (stress/strain), G = shear modulus, ᴠ = Poisson's ratio, σ yield and σ ult = yield and ultimate stresses, ε ult and ε ult = yield and ultimate strains, ε ult * = strain up to which cancellous bone retains some load-carrying ability. Moduli are in GigaPascals (GPa), and σ yield and σ ult in MegaPascals (MPa). Yield and ultimate stresses are reported along the z (long) axis because experimental test samples are typically oriented along this axis, in uniaxial tension or compression tests, and in bending tests which cause compression and tension along different sides of the long axis. The rationale behind this materials testing practice is that most in-vivo loads of limb elements are assumed to be primarily oriented longitudinally. Sources: R = Reed  instead of solely using a greyscale range. To best visualize densities, each color scale is centered on an average HU level with a range around it. HU color scales differ for the adult and juvenile because the average density is lower and range of densities smaller in the juvenile.

Finite element analyses: properties, loadings and constraints
Bone density and stiffness are testably correlated with stresses and strains from locomotor force [23], using finite element analysis [23,24,25]. Table 2 lists stiffnesses of and other properties of bone used in this study. FEA calculates mechanical response of structures modeled as a mesh of small elements, connected at nodes and assigned stiffness and other properties of the original material. The model is constrained to prevent rigid body motion, subjected to external forces, and the resulting displacement of nodes enables calculation of stresses (internal force/ area), strains (change in dimension/initial dimension), and reaction forces. We constructed a finite element model of the adult Struthio tarsometatarsus in a quasistatic, mid-stance pose in which ground reaction force of running would be vertical and at its greatest magnitude. The FE model was constrained at the ankle joint by ligaments and contact with the tibiotarsus. Components of the ground reaction force F GRF on the distal end of the tarsometatarsus and muscle force magnitudes and directions necessary to counteract F GRF , were calculated in the "ground" coordinate system, with the vertical z-axis in line with the ground reaction force (see Materials and Methods for details). Because forces were applied to the FE model in its "anatomical" coordinate system (with z along the TMT proximodistal axis), we shifted force directions using coordinate frame rotation common in robotics [26].
Forces from other muscles and ligaments (including the reduced M. fibularis brevis) not necessary to extend the tarsometatarsus were accounted for (Table 2) through the applied ground force, derived from reaction forces at finite element constraints, or applied as stays that would counteract bending. Proximal muscles (not inserting on the TMT) held the femur and tibia stiff to hold the body over the ankle joint, contributing to the ankle's simulated instantaneous constraint (and resulting reaction force) in this simulated quasi-static pose. The phalangeal flexors and extensors were considered as momentary stabilizers of the TMT-phalangeal joints, and thus contributed to the applied loading at these joints both subsumed into the vertical joint reaction force and with additional applied horizontal force in a sensitivity analysis.
We undertook several such sensitivity analyses, to examine the effects of model resolution, material properties, and loading regimes including varying muscle forces. Table 2. Muscle forces applied to the ostrich tarsometatarsus (TMT), and subsumed into finite element reaction forces at the mesotarsal joint. Forces for ankle extensors are calculated as necessary to counteract the ground reaction moment, and in one analysis digital flexor and extensor forces are applied to stabilize the TMT-phalangeal joints. Effects of muscles acting on proximal limb elements emerge from FEA as joint reaction forces at the proximal surface of the TMT.

Bone:
Axes E GPa 1. We ran convergence analyses to test how closely our initial model approached the peak stress results of high-resolution models, and thus to test its suitability for testing the effects of many different loads and material properties.
2. In addition to force from the phalanges supporting the tarsometatarsus in the vertical direction, in a separate analysis we applied forces from phalangeal flexors and extensors to simulate the effects of stabilizing loads from these muscles.
3. Unusually, our dissections of a young ostrich showed the M. gastrocnemius tendon inserting partly on a structure called the hypotarsus, along with M. fibularis longus that normally inserts there. To test the effects of both possibilities on TMT stress and strain, we ran analyses with both insertion configurations.
4. As a further test of varying biological parameters, we tested the effects of uniform material properties versus our initial model with both compact and cancellous bone.

Interpreting FEA results
We apply finite element analysis to discover how the adult tarsometatarsus would respond to the applied loads. Stress (σ: internal force/area) and strain (ε: proportional deformation) results from FEA enable us to visualize distribution and magnitude of these quantities, and to estimate how close a structure is to breaking (strength, or ultimate stress or strain: σ ult or ε ult ) or permanent deformation (yield: σ yield or ε yield ). FEA differentiates between axial (compression or tension) and shear stress and strain (Appendix 1), and calculates principal stresses and strains as eigenvectors, similar to principal components in multivariate statistics. To assess damage in bone, we use von Mises stress and strain, σ vM and ε vM , which are functions of principal values that correlate well with failure in experimental tests (see Appendix 1 of the current paper; [25]). Under the von Mises criterion, ultimate stress σ ult of dense compact bone is about 180-200 MegaPascals (MPa; N/m 2 ) in compression, 150 MPa in tension, and 80-100 MPa in shear [23]. von Mises ultimate and yield values are consistent for bone across vertebrate taxa, and are therefore reasonable assumptions for ostrich bone. Dividing the ultimate or yield von Mises stress by an element's experienced σ vM , for example, gives a safety factor under the given loading regime [21]. For example, under a compressive load, bone experiencing σ vM of 20 MPa would have a safety factor of about 10 against breaking, if σ ult is 200 MPa. Tables of stress and strain at sampled points, and color-coded illustrations of these FE results, enable assessment of safety factors throughout the TMT. Full constraints in FEA give artificially high stresses and strains, and reliable interpretations of safety factor are possible at characteristic distances from the constraint. For example, stresses and strains within a cylinder constrained across the entire surface at one end can be safely interpreted only within the part of the cylinder that is separated from the constrained end by a distance greater than the cylinder's diameter. Constraints applied to smaller surface areas result in higher (artificial) peak stresses, but enable safer interpretation closer to the constraint.

Review of tarsometatarsus external osteology
Anatomical descriptions are from our dissections and observations, primarily following terminology of Gangl [6] and Smith [9]. Figs 1 and 2 present the external anatomy and bone densities of the adult Struthio tarsometatarsus, as rendered from CT scans; labels for Figs 1-3 also associate features with forces and constraints for FEA. As in other birds, the ostrich tarsometatarsus is comprised of fused metatarsal (MT) bones II, III, IV, and the distal tarsals at the mesotarsal joint. Unique among known avian species, MT II does not articulate with phalanges and is externally lost in adults.
Anteriorly, the ostrich tarsometatarsus is broad proximally and slender distally. Proximally, the concave and oval cotyla medialis and cotyla lateralis articulate with the tibiotarsus to form the intertarsal joint. Inferior to the cotyla, the fossa infracotylaris forms a central depression anteroproximally. The crista tibialis cranialis sits within the fossa infracotylaris. Projecting directly posterior to the mesotarsal articular surface is the hypotarsus, where the extensor M. fibularis longus inserts. The hypotarsus grades distally into the Crista medianoplantaris. Ridges of the Cristae plantares lateralis and medials flank the hypotarsus and Crista medianoplantaris, and continue distally. Narrowing and tapering distally, the entire tarsometatarsus shaft forms a cylinder-shaped barrel before expanding and terminating into trochlea metatasi III and IV. MT IV diverges laterally from the main shaft at the foramen vascular distale to form the wedge-shaped trochlea metatarsi IV, attachment for phalanges of digit IV. MT III continues to expand distally past MT IV, becoming trochlea metatarsi III for phalanges of digit III, forming the metatarsophalangeal joint. M. tibialis cranialis inserts onto the Crista tibialis cranialis, slipping beneath the Retinaculum extensorium tibiotarsi. Osteological correlates for attachment of the retinaculum, such as Impressiones retinaculi extensori seen in some birds [27], are ambiguous in Struthio.

Muscle and ligament attachments
M. fibularis brevis originates distolatetal of the fibula on the tibiotarsus. While fully developed in other birds, the M. fibularis brevis is reduced to a tendon functioning as a ligament in Struthio. When fully extended the M. fibularis brevis runs transversely and crosses the origin of the Ligamentum collaterale laterale. The muscle inserts on the proxinal plantar surface on the tarsometatarsus. When moving from extension to flexion, the M. fibularis brevis crosses the Ligamentum collaterale laterale at full flexion ( Fig 5).
M. fibularis longus originates primarily from the collective tendinofacial sheet associated with the distal femur and proximal tibia, and an aponeurosis from the lateral cnemial crest at the knee. This muscle inserts proximally on the tarsometatarsus, and functions in ankle extension.
M. extensor brevis digiti IV originates dorsally on the distal third of the tarsometatarsus. Passing through the Canalis interosseous tendineus at the distal end of the tarsometatarsus, this muscle inserts medially on the first phalanx of the fourth toe.
M. extensor brevis digiti III originates dorsally on the tarsometatarsus and distomedially of the origin of the M. extensor brevis digiti IV. Attaching dorsal to the joint capsule of the metatarsophalangeal joint of the third toe, this muscle inserts on the dorsal process of the articular cartilage of the proximal articular surface of the first phalanx of the third toe.
M. lumbricalis originates on the dorsal tendon of the M. flexor digitorium longus near the beginning of the distal third of the tarsometatarsus. This muscle divides into two Crura, finally inserting on the Ligamenta plantaria of the metatarsophalangeal joints of the third and fourth toe.
On lateral and medial surfaces, pits occur for collateral ligaments connecting the metatarsals with proximal phalanges. Proximal to the collateral ligament pits, there is a transverse structure of spiculated bone. This appears to be ossified ligament connecting the diverging free shafts of MT III and MT IV. CT results: bone densities and internal tarsometatarsus osteology Tables 3 and 4 list relative densities in Hounsfield units on the surface of the adult tarsometatarsus, depicted in Figs 1 and 2, for discrete osteological features and specified locations along shaft (Fig 6). The internal structure of the tarsometatarsus (Fig 7) is complicated by fusion of the distal tarsals and metatarsals II, III, and IV. According to position criterion of homology, MT II is located medially on the shaft, with MT IV lateral, and MT III in the median plane. A high concentration of low-density trabecular bone infuses the distal-and proximal-most articular ends, with little superficial or deep cortical bone present in these regions. Larger struts of bone occur just distal and proximal to the element's extremities. Trabecular bone gives way to compact bone in the tarsometatarsus shaft. Proximally, three distinct internal cavities are present, separated by partitions of compact bone. In the longitudinally central 50% of the tarsometatarsus shaft in the adult, these partitions are obliterated and there is only one central cavity as MT II, III, and IV fuse completely. MT II narrows distally, tapering to the medial aspect of the tarsometatarsus, where on all other bird species the trochlea metatarsi II would diverge from the main shaft to form an articular surface. Cortical bone thins at the distal end, with trabecular bone increasing in frequency and density. Splitting laterally from the main shaft, MT IV forms trochlea metatarsi IV, comprised nearly exclusively of trabecular bone. MT III continues to flare distally, forming trochlea metatarsi III. Like MT IV, MT III is composed primarily of trabecular bone.

Juvenile tarsometatarsus: morphological comparison and bone densities
In contrast with an adult, the juvenile ostrich tarsometatarsus (Fig 8) is relatively more flared and broad at the proximal end. The juvenile exhibits an underdeveloped hypotarsus, and the intercondylar fossa is virtually absent. The juvenile TMT's diaphysis is significantly shorter relative to the element's width, and ankylosis of MT II, III, and IV is more clearly observed on the external surface in the juvenile. CT imaging reveals internal tripartite division of metatarsals II, III, and IV throughout the entire length of the element in the juvenile (Fig 9). This is not seen in the adult, though present proximally, this division is lost before the midpoint of the shaft. In most adults, both internal and external evidence of MT II is lost after the midpoint of the shaft. A small spur of MT II is present in some adult specimens. This contrasts with the juvenile, as MT II splits from the main shaft and terminates as a vestigial but distinct element at the point where MT III and IV begin to flare into the articular condyles. Extensive cartilage is present in the juvenile tarsometarsus in both proximal and distal articular ends. Compared to the adult, the juvenile tarsometatarsus is considerably less dense, with extensive cartilage, less densely ossified compact bone, and relatively more cancellous bone (Figs 7 and 9). Table 3. Hounsfield attenuation units (HU), stresses in the tarsometatarsus's coordinate system (+x anterior, +y medial, +z proximal), and von Mises stresses σ vm , at anatomical features and positions along the tarsometatarsus shaft. Brick numbers are sampled tetrahedral elements at the specified locations (Fig 7). Abbreviations are as in Figs 1-4; CLP = collateral ligament pit, and ITM IV = the centerpoint of Trochlea metatarsi IV. Positive values indicate tension, and negative indicates compression. Except at constraints (LCML = ligamentum collaterale mediale longum, CL = Cotyla lateralis), the greatest compressive stresses occur along the bone's long axis (σ ZZ ) on the medial and anterior surfaces, and tensile stresses posteriorly and laterally at midshaft and just distal to this ("Mid-distal"). σ vm are also highest at these positions. Note artificially high stresses where compact and cancellous elements meet in the collateral ligament pit of MT IV, low shear stresses (σ XY , σ YZ , σ XZ ) except at constraints, and low stresses at the trochlear applications of the ground reaction force (TM III lateral and medial; ITM IV).

Finite element results
Proximal region. Tables 3 and 4 present stresses and strains from forces simulated as experienced during a fast run, at discrete anatomical features (Figs 1 and 2) and at locations along the tarsometatarsus shaft (Fig 6).  Table 4. Hounsfield attenuation units (HU), strains in the tarsometatarsus's coordinate system (+x anterior, +y medial, +z proximal), and von Mises strains ε vm , at anatomical features and positions along the tarsometatarsus shaft. Brick numbers are sampled tetrahedral elements at the specified locations (Fig 7). Abbreviations are as in Figs 1-4; CLP = collateral ligament pit, and ITM IV = the centerpoint of Trochlea metatarsi IV. All notable patterns are the same as for stresses, explained in the caption for Table 3.  stresses. Proximal values are near constraints, and must be interpreted cautiously. Artificial peaks of von Mises stress σ vm and strain ε vm , up to 291 MPa and 5.91%, respectively, where cancellous and compact bone come together in our model, which we consider to be an artifact of assigning material properties. Stress is moderate otherwise, even at the proximal constraints against long-axis displacement (Figs 10-12). Stresses diminish rapidly away from the constraints; for example, a node at the Cotylus medialis experiences σ vm of 5.03 MPa and ε vm of 1.01%. Close proximity to constrained nodes gives fairly high stress and strain at proximal entheses. Von Mises strains at the attachments of M. gastrocnemius to the hypotarsus (ε vm = 2.05%) and of M. fibularis brevis (ε vm = 1.85%) approach maximum safe values for cancellous bone (1.65-2.11%; Keaveny [28] reports these and lower values). However, stress and strain become much lower deep to the surfaces of these attachment sites (Fig 12), and tendon-tobone gradations within muscle attachment enable greater strains than bone can withstand alone. σ vm stresses on the hypotarsus, with a strong pull by M. gastrocnemius, are in the range of 10 MPa. Main tarsometatarsus shaft. On the tarsometatarsus shaft distal to "mid-proximal" (  [28,29], and much higher compared with values at which cancellous bone becomes non-functional (ε ult Ã in Table 4; [28]).
Stresses and strains with variants of muscle attachments. Fig 13 depicts a loading scenario (very likely with most dissection descriptions of Struthio) with separate muscle attachments for ankle extension, through the M. gastrocnemius tendon only inserting posteriorly on the tarsometatarsus, and M. fibularis longus inserting alone on the hypotarsus. This loadcase distributes the force between the muscles, reducing the force at either attachment. This array of muscle insertions trivially decreases peak von Mises stresses on the TMT at the proximal constraints, although stresses in the distal part of the TMT shaft are similar to those of the main analysis.
Another variant on muscle forces in the original analyses simulates a high posterior force on the TMT from the digital flexors pulling back on the phalanges. A high, 800 N force from the digital flexors increases tension anteriorly and compression posteriorly on the TMT (Fig  14), the opposite of the pattern with ankle extensors alone.
Convergence analysis. Away from the constraints, all models experienced peak von Mises stress of 50 MPa, plus or minus 2 MPa. The discrepancy of less than 5% between all models suggests that our initial model had sufficient resolution for further analyses. Fig 15 highlights the similarities, depicting the highest and lowest resolution models. There was no discernable pattern to the peak stresses at the constraints. As expected, the lowest resolution model (78,000 nodes; 385,000 tetrahedra) had the low von Mises stress at the constraints (84.6 MPa). However, the highest resolution model (180,000 nodes; 931,000 elements) had a lower constraint Effects of varying material properties on stress and strain. Fig 15 depicts models with compact bone properties applied to the entire element. Regardless of node and element numbers, these models had greater stress but lower strain proximally at the constraints, as expected from the much greater elastic modulus of compact bone compared to cancellous. However, the discrepancy was not as great distally at the force applications. The all-compact bone model also        Convergence and material property sensitivity in highest (A) and lowest (B) resolution models. The lowest resolution model (77,722 nodes) has lower peak von Mises stresses (84.6 MPa) than the highest resolution model (180,029 nodes; 100 MPa). Peak stress occurs at the proximal constraints in both models. Stresses elsewhere in the models are similar. The lower peak value in the stress color histogram (B) makes anteriomedial stresses appear greater (yellower) than in the high resolution model, but sampled stresses are nearly identical throughout the TMT shaft. Both models have material properties of compact bone. Note the smoother transition from proximal to distal stresses than in models with both cancellous and compact bone properties in appropriate regions (Fig 11). doi:10.1371/journal.pone.0149708.g015 Tarsometatarsus of the Ostrich Struthio camelus proximal to slice B in Fig 7). Saint-Venant's principle [30,31] gives greater confidence in values farther from the proximal constraints, such as in the main shaft of the tarsometatarsus and distally near the application of reaction forces. Materially heterogeneous models potentially give more realistic results than homogeneous ones [31,32,33] (but also see [34]). One of our models has cancellous bone at the proximal and distal ends and compact bone elsewhere. Subsequent to our initial model construction, Eichenseer et al. [35] showed how a thin shell of plate elements with compact bone properties can be placed to surround the cancellous bone, which would be more realistic improvement to our model despite fairly low densities evident at these surfaces (Figs 1 and 2). A multi-material model with semi-automated assignment of elastic modulus [21,32,36] will likely improve on our results, especially with a smoother gradient between proximal cancellous and compact bone (Fig 11A).

The initial mesh has adequate resolution to predict results of validation studies
Mesh characteristics also influence the likely accuracy of an FE model. Our convergence analyses revealed small differences in peak, non-constraint stresses (<5%) between our initial model and one with 1.6 times the number of nodes. Bright and Rayfield [32,37] investigate how well    FEA approximates measured strains with varying mesh density, element size, and anatomical detail. The Struthio tarsometatarsus mesh consists of small (mostly 1 mm), four-node tetrahedra that Bright and Rayfield [37] report as good for approximating measured strain. In analyses of a pig cranium, Bright and Rayfield [32] discovered that FE results closely approach strain gauge results in flat regions of bone, and near force application. Conversely, results were less precise near finely contoured areas such as blood vessel grooves, and near constraints [32]. Based on these findings [32], we predict that the Struthio model will prove most consistent with strain gauge measurements along the tarsometatarsus shaft, and reasonably consistent with experimentally measured strains near an applied F GRF [32].

Tarsometatarsus densities and FEA: implications for ostrich locomotion
The adult Struthio tarsometatarsus has generally high density and stiffness of compact bone in the shaft of the tarsometatarsus, where it experiences the highest stresses. At a finer level, however, we found no specific correlations of density (as CT Hounsfield units) and stress or strain in the specimen's compact bone. For example, regressing log10 HU against percentage ε vm of sampled points in Table 2 yields an R 2 of 0.07746. Assessing detailed, adaptive correlations between CT densities and stress will require a more complex material model, that indexes Hounsfield units of bone to Young's modulus [18,19,20,21,36]. Even with only two material properties, σ vm stresses indicate safety factors of 4-5 relative to yield stress, and 5-7 relative to ultimate stress (Table 3; [38]), in the tarsometatarsus shaft and distally at the tarsometatarso-phalangeal joints. This is consistent with ostriches running fast enough to exceed F GRF of 2.5 times body weight. By the σ vm stress yield criterion, bone strength and yield in the tarsometatarsus do not appear to limit speeds in adult ostriches (as with isometric muscle force: [39]). The greatest ε vm strain values are closer to the ε yield and ε ult of compact bone (0.2-1%; [28,40]). Distal cancellous strains, however, have high safety factors relative to ε yield or ε ult . Articular cartilage is highly viscoelastic, and excellent at resisting compression at high strain rates and protecting underlying bone. Cancellous bone in the distal condyles, deep to this cartilage, probably had greater safety factors than modeled here.  [50] for components in the ground reference frame/coordinate system (ground c.s.), and rotation matrices [51] in the metatarsus's reference frame/coordinate system in Strand7 (Strand7 MT c.s.). Angles are depicted in Fig 15,  Our application of digital flexor forces, pulling the phalanges towards the tarsometatarsus, probably had too great a magnitude, inducing too great a moment on the TMT for realistic results (with a moment arm nearly as large as the entire TMT length). However, the reversal of the bending pattern (compare Figs 10 and 14) is enlightening. In life a lower force magnitude would reduce compression anteromedially and reduce tension (increasing compression) posterolaterally compared to the pattern evident in Fig 10. Because bone is stronger in compression than tension, a posterior component from digital flexors (and extensors) would increase safety factors in the ostrich TMT. When muscle firing patterns are unknown, exploring the effects of several hypothetical muscle force regimes can point the way to realistic picture of bone stress and strain in life.

Biologically relevant refinements for future simulation
Anatomy, densities, and FEA results for the ostrich tarsometatarsus suggest specific improvements on our biomechanical modeling, beyond expansion of loading regimes [2,4,7,8,9,11] and material properties [21,36]. In particular, further consideration of force of the digital flexors will refine distally applied forces at the tarsometatarso-phalangeal joints. Studies of human and cow feet that incorporate soft tissues [41,42,43,44] will guide models that include cartilage, ligaments, tendons, and perhaps more informative constraints. Experimental validation, especially with fresh specimens and loadings at behaviorally realistic strain rates, will arbitrate the accuracy of the modeling and contribute inputs for refined analyses.

Ontogenetic, selective, and evolutionary context of Struthio TMT functional morphology
Regardless of the n-th refinements possible with the FE model and its loadings, the descriptions of adult and juvenile tarsometatarsi in this paper are a stable and fundamental baseline for interpreting biomechanical results related to ostrich locomotion. CT-imaged internal structure and material densities are described here for the Struthio TMT for the first time, including fine density gradations in the juvenile specimen. Because ostriches grow so rapidly, these endpoints can anchor comparisons of ostrich anatomy and locomotion through their ontogeny, and investigations of selection and the ostrich TMT.
Bird species dependent on ground speed for locomotion rely on long legs and powerful muscles. Biometrics have been used successfully to explain ecological variability among species [45,46,47]. This study provides a novel means of studying form, function, and mechanics by combining FEA and CT data, useful for investigating the tarsometatarsus of Struthio or any other bird species.
The tarsometatarsus has high predictive power to infer behavior and ecological adaptation from morphology in birds [45,46]. Correlations of avian limb length to primary method of locomotion across broad groups have found that an elongated tarsometatarsus is not necessarily important in terrestrial environments [47]. However, an elongated tarsometatarsus does appear to be evolutionarily advantageous in Struthio. FEA revealed high safety factors, suggesting the tarsometatarsus has momentarily excessive construction [48] for ensuring sufficient foot strength when maneuvering to escape predators, delivering fatal kicks, or other less-usual but selectively critical demands.
Elongation of the tibiotarsus and tarsometatarsus compensates for limited femoral movement in birds, especially advantageous for fast running ground birds such as ratites. Ground reaction forces strongly load the femur in bending, which reduces the length of the element in cursorial birds [49]. This is reflected in ostrich morphology, as the femur is short and the tibiotarsus and tarsometatarsus are both elongated. This morphology results in a longer step length, hence increasing travel speed. An elongated tarsometatarsus may facilitate locomotion in their natural environment, as high speed would be ecologically advantageous when avoiding predators in a grasslands ecosystem.
Loss of the medial toe in Struthio is believed to reflect an evolutionary modification to reduce weight and increase speed. Our results of higher stresses on the anteromedial aspect of the distal tarsometatarsus may be a result of loss of the second digit.
CT data presented herein demonstrates that vestiges of an internal three-part division of the tarsometatarsus persists even in adults. Fusion and elongation of the TMT is supported by the fossil record in basal birds and closely related non-avian theropod dinosaurs (see Brusatte et al. [50] for a recent review). With both morphological and molecular analysis supporting a nonavian dinosaur origin for birds, modern cursorial birds offer a modern analogue to extinct theropod dinosaurs. Therefore, a detailed study of the tarsometatarsus in Struthio acts as a starting point for similar studies of extinct theropods.

CT scanning, imaging, and descriptive anatomy
The ostriches used in this study were salvaged from commercial breeders after dying from natural causes. They were donated to the University of Calgary Museum of Zoology in Calgary, Alberta and the Ohio University Vertebrate Collections in Athens, Ohio. No animals were harmed or sacrificed for the purpose of this research.
The tarsometatarsus, tibiotarsus, and femur of an adult ostrich (UCMZ) were CT scanned on a GE Lightspeed scanner (Canada Diagnostic Centres, Calgary, Alberta), at 140 kV and 175 mA, with 1.25 mm spacing and 1.25 mm overlap. A small juvenile ostrich (UCMZ) was scanned on a NewTom 3G orthodontic scanner (Aperio Services, Verona, Italy), at 110 kV and 6.19 mAs; slice thickness was 0.5 mm. Its tarsometatarsus was dissected, removed, and scanned on a SkyScan 1174 compact micro-CT, at 50 kV at a resolution of 20 μm. The juvenile specimen was placed in a carved depression in a block of low-density foam, which was then secured to the scanner's rotating stage. Because the juvenile's tarsometatarsus is taller than the Sky-Scan's detector, it was scanned twice with respective proximal and distal ends placed in the foam.
To evaluate internal anatomy, CT scans were read as DICOM data into OsiriX, and evaluated as slice data and volumetric reconstruction. In the 3D reconstructions, the tarsometatarsus data was isolated from other elements with the scissor tool. We visualized densities with the NIH color palate, which facilitates density comparisons better than a grayscale spectrum. Internal structures, external osteology, and soft tissue attachments of the adult specimen were described with reference to the literature. The terminology concurs with the Nomina anatomica avium [51] with supplementary information provided by Gangl's studies [6,52].

Finite element geometry and meshing
CT data for the adult Struthio tarsometatarsus were exported as DICOM series into Avizo (VSG, Burlington MA, USA), and surface meshes constructed. The initial surface mesh consisted of 3 million triangles, many of which had high aspect ratios which introduce numerical errors and artificially high strains [52] into FE simulations. Ideally all triangles will be isometric, but aspect ratios of 10-1 or less are adequate for accurate results. Correcting all aberrant triangles in a 3-million triangle surface mesh is time-and computationally prohibitive (often crashing the program when it attempts to correct intersecting triangles). To streamline the process and ensure adequate resolution, surface meshes were reduced to 50,000 triangles without loss of anatomical detail, and triangle aspect ratios reduced to 10-to-1 or less. These high-quality triangles were then subdivided to produce many more, smaller elements across the surfaces ("Refine faces" in Avizo), which produces better-quality triangles than in the initial large mesh with equivalent triangle numbers. A few of the new, smaller triangles had high aspect ratios, which were again corrected to 10:1 or less. The small surface elements ensured that when solid meshed, the bones' cortical walls and internal struts were modeled at least three elements deep, for accurate simulations of bending. The surface model was meshed in Strand7 (Strand7 Pty Ltd, Sydney, NSW, Australia) to produce an initial FE model of 112,630 nodes and 564,088 tetrahedral elements, with element sizes of approximately 0.8 mm internode distance. Although element numbers are typically reported for FE meshes, the number of nodes is a more realistic indicator of a model's resolution and computational complexity. (High element numbers on their own are theoretically misleading; infinite elements can converge on a single node.) The initial model has fewer elements than in many FE studies of vertebrate skulls. However, the tarsometatarsus is a geometrically simpler structure (closer in complexity to a single mandible), and our model has a comparable number of elements to an analyzed elephant femur  102°) and ω (5.006°) used for calculating y and z components, respectively (center and right equations). In these views the bones appear slightly shorter than their true lengths, because they are angled from the vertical.
( [33], although see [34]). Simpler structures require fewer nodes and elements; for example, simulation of a rectangular beam only three elements deep approximates bending stresses quite accurately.
To determine if our initial model had sufficient resolution to apply other loadcases for parameter-effect (sensitivity) analyses, we assessed convergence (how closely peak stress and strain results approximate those of a model with infinite resolution) by constructing a model with 443,248 nodes and 2,086,152 elements. We constructed five further models with 77,722 nodes (385,953 tetrahedra), 125,246 nodes (636,995 tetrahedra), 151,947 nodes (779,940 tetrahedra), and 180,029 nodes (931,000 tetrahedra). If any lower-resolution model's stresses at sampled points were within 5% of the highest-resolution model, we were confident to proceed with further loading regimes with the lower-resolution model.

Material properties
We applied material properties appropriate to their respective locations: compact bone that occurs in the main shaft of the tarsometatarsus, and cancellous bone that occurs proximally and distally. Assigned material properties ( Table 1) are those of compact and "bulk" cancellous bone (not of individual trabeculae, which have a stiffness closer to that of compact bone: [53,54]). Appendicular compact bone is orthotropic, with the highest elastic modulus parallel to the long axis. To the long axis of the element, we applied the elastic modulus E zz of emu (Dromaius novaehollandiae) femoral cortical bone under high strain rates [55] (15.86 GPa). Both emu and ostrich appendicular bone has E zz of 13-14 GPa under lower strain rates [55]. Using the high-strain-rate E zz value for emu is appropriate because bone stiffness is load-rate dependent, and loads from fast running would introduce high strain rates. Using an elastic modulus collected at high strain rates introduces some realism to a linear, quasi-static analysis, which necessarily simulates the load as steady-state. We assume that E zz a running ostrich's appendicular bones would be similar to that of emu, as they are at lower strain rates. Longitudinal modulus was multiplied by 0.57 to obtain both transverse moduli [39]. In the absence of data for ratites, assigned shear moduli were those of human compact bone as determined through ultrasound experiments [56]. Poisson's ratio of the compact bone was set to 0.3 [57]. Cancellous bone was considered isotropic, with an elastic modulus of 0.64 GPa and Poisson's ratio of 0.29 [58,59].

Finite element simulations: boundary conditions
Constraints are necessary in FEA to prevent rigid body motion and allow deformation of a structure. However, fully fixing all nodes against motion in x, y, and z is unrealistic, and finite element simulations can give artificially high stress and strain results at constrained nodes. We assumed that the tibiotarsus would constrain the tarsometatarsus from moving proximally (Zaxis of the universal coordinate system; Fig 16A and 16C; Fig 17; Fig 18), and that the lateral colateral ligament and M. fibularis brevis constrained it from slipping laterally. Varying the size of the constraint area resulted in trivial variation in stress and strain results.

Finite element simulations: forces
Ground reaction force F GRF and phalangeal joint reaction forces. Simulations neglected inertial forces. Simulations were run under the assumption of linear behavior, in which stress and strain scale linearly with force magnitude, strains are small, and deflections do not change the stiffness of the structure. Under these assumptions, deflections, stress, and strain for any magnitude of load can be calculated by multiplying results of an initial simulation (even with unit forces) by the ratio of the desired load magnitude to the initially simulated force. We applied a baseline magnitude as the ground reaction force (F GRF ) with a maximum magnitude of 2.5 times body weight, equivalent to a fast run. A regression of shank length to mass [60] yielded 100 kg for the specimen, and a maximum F GRF of 2495 N.
Note that "F GRF " in this case is really the force of the proximal phalanges on the distal condyles of the tarsometatarsus. This requires several simplifying assumptions, clarifications about their effects, and two loadcases that account for forces across the tarsometarsus-phalangeal joints. Applying both of these loadcases (numbers 4 and 5 below) shows how unknown force transmission from the phalanges to the tarsometatarsus would affect stress in the TMT.
1. The proximal phalages support the tarsometatarsus vertically, so that their reaction forces on the tarsometatarsus are in line with the vertical F GRF direction.
2. The true ground reaction force, from the subdigital pads through the body's center of mass, requires contraction of the digital flexors to counteract extension of the phalanges.
3. In one set of analyses, forces from the digital flexors are simulated as contributing to the vertical reaction force on the TMT, and exert negligible non-vertical force on the TMT Moments about the center of rotation (proximal dot) must equal 0. The ground reaction moment R GRF x F GRF is therefore equal and opposite to a balancing extensor moment R mz x F mz . When the vertical muscle force F mz is calculated, the resultant F m and its x and y components can be determined by the law of cosines [50] (Figs 15 and 16; Table 5). doi:10.1371/journal.pone.0149708.g017 condyles. These forces would be incorporated into the 2495 N value for the reaction force, and not further calculated for the analysis.
4. An another set of analyses, we applied a posterior component to the phalangeal reaction force to simulate what would happen when digital flexors pulled the proximal phalanges towards the joint, in addition to their vertical contact. The forces are transmitted across the joint the anterior surfaces of these condyles.
For the analyses described in 4) above, an 800 N posterior load was applied, distributed proportional to the maximum isometric forces exerted by the flexors of digits 3 and 4 [10], seen in Table 7. These were 29% to digit 4 (231 N), and 71% to digit 3 (569 N). The 800 N is an arbitrary value, but useful to examine the sensitivity of results to such forces.  Table 4, and listed in Table 5. Components of F m in the global and anatomical coordinate systems are listed in Table 6. FEA can determine many quantities listed as "Known", including F lig at the ligamentous constraints, and tibiotarus and tarsometatarsus joint reaction forces F ttRF and F tmtRF at the constrained proximal surface of the metatarsus. Moment arms r will vary with angles through the motion, producing moments when coupled with their forces; the FE simulations incorporate these quantities at one position. Gravitational force F g and its moments on the tarsometatarsus cm (incorporating mass tmt and mass moment of inertia I tmt ) contribute trivially, compared to effects of muscular, ground, joint forces. doi:10.1371/journal.pone.0149708.g018 Estimating force components was more complex than the magnitude of the ground reaction force. F GRF at the midpoint of the step cycle in a straight run [11] is at its greatest magnitude and simplest direction (vertical), neglecting friction. We therefore graphically derived joint angles from kinematic data of Jindrich et al. [11], at a point in stance phase when the ground force is vertical. The local coordinate system for the tarsometatarsus was set with the x-y plane at the distal end of MT III, and z along the long axis. In this pose, the tarsometatarsus was inclined at 82.1°to the horizontal in the sagittal plane (lateral view), and 81.3°in the transverse plane (anterior view; Fig 19). The z-axis is at an angle of 5°from the vertical. The resultant GRF and ankle forces were therefore vertical in the universal coordinate system, but angled relative to the z axis of the local system. Their components in the universal x (anteroposterior), y (mediolateral), and z (proximodistal) directions were calculated with the law of cosines (Fig 10; [60]), and applied to the distal or proximal joint for respective simulations.
Force to counteract F GRF moment. To calculate the total extensor force F m , (distributed between M. gastrocnemius and M. fibularis longus), moments about the ankle joint were set to 0 so that the ground reaction and extensor moments about the mesotarsal joint were balanced (Fig 17). The center of rotation in lateral view was estimated to be at the midpoints between anatomical and experimental axes determined by Jindrich et al. [11]. Moment arms were measured graphically in OsiriX, with the CT render positioned so that the angles between moment arms and forces were equal to 90 degrees. The vertical component of the extensor force, necessary to counteract the ground reaction moment, was 4202 N. Future analyses will assess contribution of posterior attachments of M. gastrocnemius, which will exert additional moments as long as they are posterior to the transverse axis of joint rotation.
Extensor tension F m . The extensor muscle tension F m was greater than its vertical component, because the resultant of M. gastrocnemius acts along the tibia is angled in 3D relative to the tarsometatarsus. Unknown components and the resultant of F m were first calculated with the law of cosines [60,61] in the reference frame of the ground, using known vertical component and angles of the tarsometatarsus (Table 5; Figs 15 and 16). The force vector F m in the global coordinate system was reoriented into the anatomical coordinate system (with the z-axis coincident with the tarsometatarsus long axis) with rotation matrices [27]. For the first rotation  Fig 15. A and B are for a lateral view, and C and D in anterior view. A and C are in the global reference frame and coordinate system of the ground (Fig 15). B and D are in the anatomical reference frame and coordinate system of the tarsometatarsus (TMT). Table 5 includes rotation angles and matrices for forces in this frame. Variable definitions are in Table 4. doi:10.1371/journal.pone.0149708.g019 Tarsometatarsus of the Ostrich Struthio camelus about the y-axis, using angles from Table 5 and Fig 18, the matrix (Eq 1) is: The second rotation about the x-axis has the matrix (Eq 2):  Table 6 presents components of F m in the global and anatomical coordinate systems. Table 7 lists muscle forces components applied in all loadcases, plus the ground reaction components.
Applying F m to attachments of the M. gastrocnemius tendon. Fig 20 diagrams all forces and constraints applied to the FE model, including distribution of F m . The magnitude of F m was assumed to be constant in the M. gastrocnemius tendon, and intuitively would be divided equally among all nodes of its attachments to the hypotarsus and the posterior surface of the tarsometatarsus. Directional components of F m on the hypotarsus were calculated as described above (Figs 15 and 16) but components on the Crista plantares lateralis and Crista plantares medials (not in line with the direction of pull on the hypotarsus) were not obvious. To simulate displacements of the tendon along its attachments to the posterior surface (Fig 3), we applied masking tape to the attachments on a Struthio tarsometatarsus (Ohio University Vertebrate Collections OUVC 4023), with a free end replicating the tendon as it approaches the hypotarsus. We marked the bone and masking tape with pencil, positioned the metatarsus at the proper angles, and pulled on the free "tendon" in the resultant direction of F m . Markings at the proximal portion of the hypotarsus became displaced in the direction of pull. Markings at attachments more distally on the hypotarsus, and to the posterior cristae, were displaced almost exclusively along the bone's long axis, with negligible lateral displacement. At the posterior nodes of M. gastrocnemius attachment, we therefore restricted tension of M. gastrocnemius to the z, long axis of the bone. Appendix 1: Stress and strain theory and notation Stress, the internal force per unit area of a material, in a three dimensional structure is accounted for with a stress tensor, which records stresses in three directions (x, y, z) on three orthogonal surfaces (x, y, z). These stresses are conventionally and conveniently represented in matrix form, in Eq 3. Normal stresses (on the diagonal of the matrix), σ xx , σ yy , σ zz , occur where compressive (crushing) or tensile (pulling) stress is perpendicular to a surface in the x, y, or z directions, and thus parallel with these axes. Shear stresses σ xy , σ xz , σ yz , are parallel to a surface but perpendicular to each other. For example, σ xy signifies that in the xy plane, stress is in the x direction relative to the y axis. This shear is analogous to having children's toy blocks stacked vertically (on the y axis), and pushing on a block in the x direction (perpendicular to the y axis). Because a shear stress such as σ xy is equivalent to σ yx , the stress tensor is symmetric and designated by sym in the matrix notation. von Mises theory is a common choice for reducing the stress tensor to a single value that can be compared to experimental results of a uniaxial tensile test, and reflects the likelihood of material failure. This theory recognizes that many materials fracture not due to stresses that change the volume (hydrostatic or volumetric stresses), but rather those that change the shape (distortional stresses). The von Mises equivalent stress is a scalar value that reflects these distortional stresses (Eq 4).
s vM ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 0:5½ðs xx À s yy Þ 2 þ ðs yy À s zz Þ 2 þ ðs zz À s xx Þ 2 þ 6ðs 2 xy þ s 2 yz þ s 2 xz Þ q 4Þ Strain is a measure of relative deformation (change in length divided by original length), which also takes a form of a tensor of 6 unique values, (ε xx , ε yy , ε zz , ε xy , ε xz , ε yz ). Von Mises strain can be calculated with the same equation above, substituting strain components for stress. The levels of von Mises stress σ vm and strain ε vm relative to a material's yield values σ yield and ε yield indicate how close a structure is to nonlinear deformation and therefore damage, and their approach to ultimate values σ ult and ε ult indicate how close the structure is to immediate fracture. Dumont et al. [24] and Rayfield and Milner [25] cogently discuss von Mises criteria as they relate to bone. We assessed the mechanical integrity of the ostrich tarsometatarsus by comparing simulated stresses and strains to σ ult and ε ult of compact and cancellous bone. σ ult and ε ult are unknown for bone of the ostrich tarsometatarsus, and were assumed to be similar to those of emu, human, and bovine cortical bone [55,36,40], and human and dog appendicular cancellous bone [62,63,64,28].
Supporting Information S1 Dataset. Strand7 finite element model of Struthio tarsometatarsus, with compact and cancellous bone properties. This FE mesh has 77,722 nodes and 386,953 tetrahedra. It is the maximum size for PLOS ONE supplementary files. Higher-resolution models are available from the second author. This model gives realistic strain results within cancellous and compact bone, but unreliable results where the materials meet.