Predictive simulations of running gait reveal a critical dynamic role for the tail in bipedal dinosaur locomotion.

Locomotion has influenced the ecology, evolution, and extinction of species throughout history, yet studying locomotion in the fossil record is challenging. Computational biomechanics can provide novel insight by mechanistically relating observed anatomy to whole-animal function and behavior. Here, we leverage optimal control methods to generate the first fully predictive, three-dimensional, muscle-driven simulations of locomotion in an extinct terrestrial vertebrate, the bipedal non-avian theropod dinosaur Coelophysis. Unexpectedly, our simulations involved pronounced lateroflexion movements of the tail. Rather than just being a static counterbalance, simulations indicate that the tail played a crucial dynamic role, with lateroflexion acting as a passive, physics-based mechanism for regulating angular momentum and improving locomotor economy, analogous to the swinging arms of humans. We infer this mechanism to have existed in many other bipedal non-avian dinosaurs as well, and our methodology provides new avenues for exploring the functional diversity of dinosaur tails in the future.

INTRODUCTION

As a key aspect of behavior in many animals, understanding locomotion is integral to deciphering the biology of extant and extinct species (), as well as broader-scale questions of ecological interactions, evolution, diversification, and extinction across major clades (, ). Inferring locomotion in extinct vertebrates typically draws upon the study of fossilized bones and footprints, which offer an incomplete view (), or locomotor behavior in extant relatives (), which may nevertheless deviate substantially in anatomy. Computational biomechanical methods offer an additional approach, by applying universal physical and biological principles, derived from the study of extant species, to quantitatively and mechanistically relate observed anatomy to behavior and performance (). Studies using these methods can provide important insight into posture, muscle function, bone loading, and limb control, and also provide estimates of maximum performance. Yet, locomotion is an inherently dynamic, complex behavior that can involve multiple parts of the body (), and static analyses of single limb poses are still a long way from deriving complete gait cycles for steady behaviors (), let alone more transient behaviors such as jumping or turning.

A dynamic, whole-organismal approach can help forge stronger links between fossilized locomotor anatomy, behavior, performance, paleoecology, and macroevolutionary patterns. Moreover, the fusion of biomechanics with methods in mathematical optimization can be used to generate simulations of behaviors de novo, by seeking to maximize some physiologically relevant performance criterion (). Optimization-based predictive simulations can be completely independent of any experimental data, where dynamic behaviors evolve naturally from the underlying physics of the system in question (). Being freed of the constraints of available empirical datasets, such as limited experimental data or scope for anatomical variation within extant species, they hence have formidable potential for investigating locomotion in the fossil record. Indeed, predictive simulations provide the only means to explicitly and rigorously test musculoskeletal function in anatomies that are not observed among extant species, and hence which cannot be studied in vivo ().

The application of predictive simulation approaches to questions of locomotion in extinct species has focused on hominins () or non-avian dinosaurs (, ). Yet, the highly nonlinear, hyperdimensional dynamic equations involved render these simulations computationally expensive, previously requiring supercomputing and/or drastic modeling simplifications (). Here, we leverage recent developments in rapid numerical methods of optimal control () to create three-dimensional (3D) simulations of gait in the extinct bipedal non-avian theropod Coelophysis bauri. This well-known Triassic taxon () represents the general bipedal and “cursorial” (i.e., with anatomies conducive to high-speed locomotion, such as long limbs) bauplan ancestral for dinosaurs, as well as most bipedal dinosaurs generally; it can therefore illuminate ancestral biomechanical traits for Dinosauria and Theropoda, and the evolution of fast-running bipeds (). Previous static or dynamic analyses of non-avian dinosaur locomotion were at best “2.5 D,” almost always ignoring the role of axial body segments (neck, thorax, trunk, and tail) by treating them as a single rigid entity [(, , ); but see 2D analysis of ()] and inferring the tail as just a static counterbalance to the cranial body segments. For the first time in an extinct species, we explicitly incorporated 3D axial body movements, seeking to determine how coordination of whole-animal movement interacted with morphology, neuromuscular control, and performance. In the process, we found a crucial and dynamic role of the tail that, after extensive testing, we infer to be a passive, physics-based mechanism for locomotor economy, providing insight into the potential role of the axial body in terrestrial dinosaur locomotion. While the tail of non-avian dinosaurs has frequently been assumed to have been important in various nonlocomotor behaviors, we here provide the first rigorous, mechanistic demonstration of its (previously unrecognized) dynamic role in steady-state terrestrial locomotion, broadening our perspective on the diversity of tail functions in the distant past.

RESULTS

Tinamou bird simulations

We developed a generalizable framework for producing Hill-type muscle-driven predictive simulations of steady-state, straight-line locomotion in a theropod (see Materials and Methods) and first evaluated its validity with a musculoskeletal model () of a tinamou bird (Eudromia elegans, 545 g), an extant homolog and analog for non-avian theropods. To facilitate comparison to prior empirical observations, we performed simulations constrained to move at the slowest (0.39 m/s) and fastest (2.62 m/s) reported trials, corresponding to slow walking and moderate aerial running, respectively (, ). Both optimal control problems (OCPs) were solved (converged) within 7 min using a single processing core. Without reference to any experimental data, the simulations spontaneously generated walking and running gaits that had a strong kinematic and kinetic match to empirical observations (Fig. 1 and movie S1). Key aspects of gait were qualitatively, if not quantitatively, replicated, including the patterns of ground reaction force (GRF) time histories [e.g., double-humped vertical GRF profile in walking, and early-skewed and single-humped vertical GRF profile in running ()] and whole-body center of mass (COM) mechanics, with vaulting in walking and bouncing in running [percent congruity () of 18 and 64% between kinetic and potential energies, respectively, for walking and running]. In terms of gross timing with respect to stance and swing phases, simulated muscle excitation patterns also generally showed good agreement with previously published electromyography data for avian walking and running (fig. S1). In particular, strong accord was observed for key hip extensors [e.g., iliotibialis lateralis pars postacetabularis (ILPO), flexor cruris lateralis pars pelvica (FCLP), and puboischiofemoralis medialis et lateralis (PIFML)], knee extensors [e.g., iliotibialis cranialis (IC) and femorotibialis intermedius (FMTI)], ankle extensors [e.g., gastrocnemius lateralis (GL) and fibularis longus (FL)], ankle flexors [e.g., tibialis cranialis caput tibiale (TCt)], and metatarsophalangeal (MTP) plantarflexors [e.g., lexor perforans et perforatus digiti 2 (FPP2) and flexor perforatus digit 4 (FP4)]. Most deep dorsal thigh muscles also showed strong agreement [iliofemoralis externus (IFE), iliotrochantericus medius (ITM), and iliotrochantericus caudalis (ITCa)], but the iliotrochantericus cranialis (ITCr) did not. As observed in previous inverse simulations (), excitations of the ambiens (AMB) and iliofibularis (ILFB) also did not show much similarity to experimental observations.

Fig. 1.

Tinamou simulations.

(A) Recovered gait cycle for walking at 0.39 m/s, showing limb kinematics, ground reaction forces (GRF) (aqua arrows), Center of mass (COM) position (green and black disc), and muscle recruitment. Muscles of only the right leg are shown for clarity. (B) Recovered temporal patterns of right-leg GRFs compared to experimental data () and scaling-based predictions derived from statistical models fitted to experimental locomotor data for 12 species of birds across a wide range of speeds and body sizes (). Note that scaling-based predictions are only available for F (fore-aft) and F (vertical). Some discrepancy in vertical force between simulation and empirical curves is explained by the empirical curves not necessarily being dynamically consistent, due to stride-to-stride variability for example (). (C) Comparison of key kinematic parameters describing gait for the simulations versus experimental data or scaling-based predictions: stance duration, duty factor, and stride length. (D) Recovered gait cycle for running at 2.62 m/s, showing limb kinematics, GRFs, COM position, and muscle recruitment. (E) Simulated temporal patterns of right-leg ground forces reasonably matching experimental data () and scaling-based predictions derived from locomotion data for 12 bird species (), for 2.62 and 3.4 m/s (maximal speed) running. (F) Correspondence of key kinematic descriptors of gait for the simulations versus experimental data or scaling-based predictions.

A few minor discords with prior empirical data were also evident for limb kinematics and kinetics. First, the model used a more upright (extended) hip posture than reality (, , ), and the digits barely cleared the ground during the swing phase, even during running. Similar effects have been observed in previous predictive simulation studies as well (, –, ) and potentially reflect the formulation of our objective function, which seeks to minimize a gross measure of effort (“energy expenditure”); in reality, tinamous may also prioritize other aspects, such as peak power demands () and/or stability (). The range of MTP angles used by the model was restricted compared to real motions, which may be a consequence of effort minimization in the objective function, or the use of a single contact sphere to model foot-ground interactions. It was deemed acceptable in this study to sacrifice some accuracy in the distal limb to facilitate examination of whole-limb and whole-animal mechanics, pending further advances in general foot-ground contact modeling methods for avian feet. One other notable discord is that the peak vertical GRF was higher in the running simulation compared to both experimental data and empirical predictions based on experimental locomotor data for tinamous and 11 other bird species (), which corresponds to the use of a slightly shorter stance duration and stride length.

In addition to verifying against prior empirical data, we also used the model and simulation framework to estimate maximal running speed in tinamous. The imposed target forward speed was progressively increased until the simulation could no longer converge, signifying the model’s maximum capable speed (3.4 m/s). Although it is difficult to empirically verify an organism’s true maximum performance capability in any behavior, our estimate compares well with the fastest reliably recorded running speeds previously reported for other bird species (Fig. 2).

Fig. 2.

Maximum speed estimate of tinamous.

The maximum speed obtained for the tinamou simulations compared with the fastest reliably recorded speeds previously reported for other species of ground-dwelling bird from (, ). Although it is difficult to ascertain an animal’s true maximum running speed, the estimate for the tinamou clearly compares favorably with data for similarly sized birds.

Simple muscle model

The use of Hill-type muscle models hinge on architectural parameters that are not preserved in the fossil record. To address this, we developed a simplified, architecture-free model of muscle contraction suitable for application to extinct species and identified the combination of muscle-tendon unit (MTU) parameters that replicated the maximal speed simulation (see Materials and Methods; fig. S2 and movie S2). When simulations were run using this simplified model, the minimally required muscle strength required to generate a 3.4 m/s run varied markedly in relation to MTU stiffness, kMTU (fig. S2B). We identified the “best-performing” combination of parameters as that which involved the lowest possible strength scaler [as a multiple of body weight (BW)]. This combination was found to be a model with kMTU = 2 and maximum force (Fmax) for every MTU set to 2.15 BW; the corresponding OCP converged in <4 min, less than half the time of the 3.4 m/s simulation generated using Hill-type muscles (9 min). (For comparison, in the original model with Hill-type muscles, the mean Fmax of the MTUs was 1.95 BW.) The resulting simulation kinematics and kinetics were comparable to those for the Hill-type muscle simulation (fig. S2C and movie S2), indicating that the simplified muscle model has merit for simulating at least some whole-organismal behaviors.

Coelophysis simulations

We applied our validated framework to a previously developed musculoskeletal model of Coelophysis () (13.1 kg; fig. S3). The OCP converged in <5 min, and we obtained a maximum speed of 6.65 m/s, comparable to previous static-based estimates and empirical scaling patterns generally. As with the tinamou, the qualitative temporal pattern of GRFs replicated empirical scaling-based predictions, although peak vertical GRF was higher than expected (fig. S4A) (). The models’ capability to sustain higher-than-expected vertical GRFs suggests that, if peak ground forces are a proximate limitation to maximum running speed () and GRFs in reality were more similar to empirical expectations, the current simulation framework may slightly underestimate maximum running speed.

The Coelophysis model used a strongly extended hindlimb posture during stance phase, consistent with previous assessments (), and simulated stride length was comparable to empirical scaling-based predictions (Fig. 3A, fig. S4B, and movie S3) (). However, stance duration and duty factor were lower than expected based on a scaling-based statistical model derived from experimental locomotor data for 12 extant bird species (); of course, predictions derived from extant birds may not be expected to be fully applicable to all extinct theropods, especially if craniocaudal COM locations differ (). As in the tinamou simulation, the swing limb minimally cleared the ground. Although the model kept its feet close to the midline, it used a wider step width than expected for its speed based on a previous study of Triassic theropod footprints ()—expected width: −0.021 m or −0.048 times hip height; simulation width: 0.058 m or 0.129 times hip height. Nonetheless, this amounts to a trivial difference in whole-limb posture: Assuming a midstance hip height of 0.45 m, the greater step width increases hip abduction by only 5°. Such a small change in global limb orientation, and positioning of center of pressure versus COM, would have minimal effect on mediolateral GRFs and whole-animal mechanics, as is seen in wide-gauge walking and narrow-gauge running in humans (–). Key hip extensor muscles such as the caudofemoralis longus (CFL; the only limb muscle that attaches to the tail) and adductors exhibited a single main burst of activation during stance, whereas deep dorsal thigh muscles such as the IFE and ITCa showed sustained bursts during both stance and late swing (fig. S4C). This supports inferences of the latter muscles as both important abductors and flexors of the non-avian theropod hip ().

Fig. 3.

Coelophysis simulation.

(A and B) Simulated maximum speed running gait, in lateral (A) and dorsal (B) views, showing limb and axial segment kinematics, GRFs, COM position, and muscle recruitment. Note that the black and white circles are only used to illustrate the relative cranio-caudal positioning of the left (L) and right (R) feet at a given instant in time; their true mediolateral positions are obscured, hidden underneath the body (see movie S2). (C) Time histories of limb and axial kinematics, showing temporal coupling between them. Prox., proximal; Dist., distal.

The recovered kinematic and kinetic patterns for Coelophysis were generally robust to variation in the location of the whole-body COM and the construction of the tail (fig. S4, D to G). With a caudally shifted COM, simulation kinematics and kinetics were extremely similar to the nominal simulation, although the neck + head segment was more ventrally depressed (extending its COM further cranially, slightly countering the caudal COM shift). A cranial COM position also produced a similar result, although the model was only able to attain a reduced maximal speed of 5.5 m/s, contrasting somewhat with previous findings (), where a simplified 2.5 D simulation of the much larger non-avian theropod Allosaurus was used and results were relatively insensitive to COM position. Altering the location of the proximal and distal tail joints in the model—but otherwise permitting the same mobility as in the nominal simulation—also reduced maximal attainable speed to 5.5 m/s, which is likely due to the effect of reduced moment arms of the CFL about the proximal joint (by 34.7 and 17.7% in dorsiflexion and lateroflexion, respectively, in the reference pose). Otherwise, the resulting kinematic and kinetic patterns were extremely similar to those of the nominal simulation.

Tail lateroflexion

Unexpectedly, the Coelophysis simulations involved pronounced movements of the neck and tail (Fig. 3, B and C, and movie S3). Dorsiflexion of the neck and proximal and distal tail joints was coupled with hindlimb protraction-retraction with two oscillations per stride cycle. There was also substantial lateroflexion of both tail joints (single oscillation per stride), where peak lateroflexion was closely coupled to peak protraction-retraction of the hindlimb; when the left hindlimb was retracted, the tail flexed to the left, and vice versa for the right hindlimb.

Lateral tail motions in straight-line terrestrial locomotion in dinosaurs had not been predicted before. We identified two explanations for why lateroflexion occurred in our simulations: enhancement of CFL force versus regulation of whole-body angular momentum; furthermore, the underlying mechanism may have been actively or passively driven. To test between these hypotheses, we ran 10 additional simulations wherein a certain aspect of the model or optimization problem was modified (Materials and Methods; Figs. 4 and 5 and figs. S5 and S6). Modeling the tail as a passive extension of the trunk, by removing all active control of tail joints, had a negligible effect on the simulation. Similarly, removing the potential for force enhancement via passive CFL stretch, or removing the CFL together, also had minimal effect. We therefore reject the possibility that the recovered tail motions were an actively controlled strategy related to increased force production in the CFL.

Fig. 4.

Further comparisons among Coelophysis tail lateroflexion tests.

(A) Key angles for axial body movement compared among each test and the nominal simulation. (B) Phase offsets between proximal tail lateroflexion and hip extension angles for each simulation, as well as between pelvic yaw kinematics for each variant and the nominal simulation; high values indicate out-of-phase, and low values indicate in-phase. n/a, not applicable. (C) Activation cost contributions from muscles (controlling limb joints) and torque motors (controlling axial joints, as well as supplementing the MTP joint) in each simulation. (D) Root mean square deviation for each coordinate (except forward translation) and each component of the GRF between each variant and the nominal simulation. l-flex, lateroflexion; d-flex, dorsiflexion. These results show that the synchronized tail simulation (and to a reduced degree, the rigid joint and heavier tail variants) by far had the greatest deviation in joint kinematics and their relative phasing.

Fig. 5.

Coelophysis tail lateroflexion tests.

(A) Mean activation of the right leg muscles throughout the gait cycle, illustrating the markedly lower economy of the synchronized (Sync.) tail simulation. (B) Illustration of yaw (Y) and pitch (Z) axes of axial body movement. (C) First principal component (PC1) describing segmental contributions to (and cancelations of) angular momentum about the yaw and pitch axes over the stride; bars of opposite sign signify cancelations between respective segments. met., metatarsus. (D) Time histories of yaw and pitch angular momentum (H and H, respectively) for the nominal and synchronized tail simulations showing the key role the tail has in regulating whole-body angular momentum.

Restricting the neck, back, or tail joints had a negligible effect on limb kinematics, a minor effect on the kinematics of unrestricted axial joints and marginally increased required effort (Fig. 4). When the tail joints were restricted (“rigid axis” and “rigid tail” simulations), the phase of pelvic yaw with respect to hindlimb protraction-retraction was strongly altered compared to the nominal simulation (Fig. 4B); this enabled the tail (as an effectively rigid extension of the pelvis) to lateroflex with respect to the hindlimbs in the same temporal pattern as the nominal simulation. When the tail was removed and its mass transferred to the trunk, the phase of pelvic yaw again was strongly modified compared to the nominal simulation (Fig. 4B), suggesting that the (now more massive) trunk segment was functionally compensating for the missing tail. These results indicate that the relative timing between tail, body, and limb movements—and in turn angular momentum generation—is critical. Encouraging the model to use tail lateroflexion–hip extension kinematics of opposite phase with respect to the nominal simulation (“synchronized tail”; Fig. 4B and movie S4) produced starkly differing results. The resulting simulation used markedly greater lateroflexion at the proximal tail joint, as well as increased pelvic yaw and roll (Fig. 4A). Furthermore, required muscular effort was massively increased, at least 2.6 times greater than in any other simulation (Fig. 4C), corresponding to both generally higher activations throughout the stride and a novel burst of heightened activity during the swing phase (Fig. 5A). The originally recovered kinematic pattern of tail lateroflexion is therefore a passive, economical behavior that evolves naturally from the underlying physics of the system.

Computation of segmental angular momentum about the whole-body COM (see the Supplementary Materials) revealed that the nominal simulation’s tail kinematics were strongly linked to fluctuations in whole-body angular momentum over the gait cycle. Principal components analysis (PCA) of temporal fluctuations () demonstrated that angular momentum about the vertical axis (yaw, H) was dominated by contributions from both tail segments (Fig. 5, B and C, and fig. S5). The tail’s contribution to whole-body H not only completely countered that of every other body segment combined but also completely reversed the temporal pattern of total H fluctuations (effectively flipping fluctuations about the abscissa; Fig. 5D). Over a whole stride, the hindlimbs were the principal contributors to fluctuations in angular momentum about the mediolateral axis (pitch, H), with H of the right leg countering that from the left leg; the two fluctuated out of phase with one another. Because angular momentum is disproportionately influenced by a segment’s distance from the point of reference in comparison to the mass of the segment (see equation S10), this explains the disproportionate influence of the distal hindlimbs on H compared to that of the more massive thighs (Fig. 5C). Over only half a stride, H was dominated by contributions from the neck + head and tail segments, fluctuating out of phase with one another. Extremely similar results were obtained in most of the other simulations tested. Changes to tail mass (within the ranges explored here), particularly a decrease, generally had minimal effect on simulation behavior; this again can be explained by the greater influence of distance than mass on angular momentum, an aspect highly relevant to a long tail.

Angular momentum patterns in the synchronized tail simulation starkly contrasted with those of the nominal simulation (Fig. 5, C and D, and fig. S5). Here, H of the tail did not counteract that of the other body segments but rather constructively added to it to create an oscillation of total H whose maxima exceeded that of the nominal simulation by 77%. Furthermore, maxima in total H occurred during swing (rather than stance), coinciding with the second burst of heightened limb muscle activity (Fig. 5A). Many of the muscles that experienced this second burst were hip flexors or long-axis rotators (e.g., iliotibialis 1, puboischiofemoralis internus 1, IFE, ITCa, and puboischiofemoralis externus 1 and 2); notably, during this phase, the CFL was recruited less compared to the nominal simulation. We infer that increased muscular effort in the synchronized tail simulation largely derived from needing to more actively control the disposition of the swing limb, probably for altered maneuvering to maintain dynamic stability (). Vertical angular momentum patterns were also strongly modified in the “no tail” simulation, where the more massive trunk segment dominated temporal fluctuations (fig. S6). Despite this, the temporal profile of total H remained qualitatively similar to that of the nominal simulation (albeit with peak magnitudes reduced by 48%); this further suggests that, by altering pelvic yaw kinematics, the trunk was functionally compensating for the role of the missing tail in this simulation. Fluctuations in angular momentum about the pitch and roll axes remained largely unaltered.

DISCUSSION

Predictive simulations are uniquely positioned in paleontological enquiry, because they permit explicit and mechanistic investigation of anatomies that have no extant equivalent. Here, we focused on a striding, parasagittal, obligate biped with a pronograde trunk and a long, heavy tail—a body plan that died out ~66 million years (Ma) ago. Unexpectedly, our Coelophysis running simulation involved marked tail lateroflexion that was coupled with hindlimb protraction-retraction. Dinosaur tails have been extensively discussed in terms of antipredator defense (), intraspecific display (), and the execution of unsteady behaviors such as jumping () or body repositioning (), and recently, their role in aquatic propulsion has been reconsidered (); most such considerations have been qualitative and phenomenological only, with limited testing. In analyses of steady, straight-line terrestrial locomotion, previous studies have almost always treated the tail as a static, caudal extension to the pelvis [but see ()], merely acting as a counterbalance to the body cranial to the hips. Here, we move beyond assumptions or vague speculation of the tail’s importance and mechanistically demonstrate a previously unrecognized, crucial, and 3D dynamic role of the non-avian dinosaur tail.

By accounting for movements of each of the axial body segments, our study provides the first insight into the dynamic, 3D role of the tail during terrestrial locomotion in extinct dinosaurs. Our results support the hypothesis that the tail was the principal regulator of whole-body angular momentum about the vertical (yaw) axis, acting as an inertial damper. Recovered lateroflexion kinematics did not serve to minimize total vertical angular momentum () but rather constructively modified it to produce a stereotyped (i.e., evidently “desirable”), nonzero temporal pattern, one that can be exploited for economical coordination of dynamic stability and limb control during gait. Although our current simulation framework may result in the use of an overly extended limb posture during stance, or limited ground clearance during swing, these would be expected to influence temporal fluctuations of whole-body angular momentum about the mediolateral axis, but not that about the vertical axis, the axis of primary relevance to the current study.

Our findings collectively demonstrate that lateroflexion of the tail of Coelophysis, and by inference many other bipedal non-avian dinosaurs, played a role analogous to the swinging arms of humans during walking and running (, , ). Both movements are largely passively driven, primarily modulate fluctuations in vertical angular momentum, and are coordinated with other body movements to minimize required muscular effort. When the tail was removed, our model still sought to achieve the same stereotyped temporal pattern of H, but adopted modified axial kinematic patterns to do so, and incurred a slightly (18%) higher muscle effort cost (although the caveat of a concomitantly heavier trunk segment in this model should be recognized). This implies that while having a tail may not be critical to the achievement of steady, straight-line locomotor performance per se (as fast-running birds readily demonstrate; Fig. 2), when the tail is present, it plays a fundamental role in governing whole-animal dynamics and functional integration: Improper use of the tail will substantially impede locomotor function and economy. If this is the case, why was the tail not reduced in other bipedal dinosaur clades apart from the tetanuran theropods? The answer is likely twofold. First, the tail housed a key hindlimb retractor muscle (the CFL) and provided a counterbalance that helped position the whole-body COM near the hips. Second, it likely played an important role in multiple other behaviors (locomotor and otherwise) as well, as noted above. Its reduction was therefore probably not a simple affair, and it took ~80 Ma for this to be accomplished in the avian stem lineage, in concert with marked changes in body proportions and limb posture (, ). Regardless of the proximate reason for its persistence in most bipedal dinosaurs, we contend that the presence of a long tail would have had a critical influence on the execution of straight-line locomotion such as fast running.

Locomotion studies of bipedal non-avian dinosaurs frequently consider just the hindlimbs, but we caution that this reductionist approach may miss important contributions from other body parts. Our results indicate that the coordination of movements throughout the whole body is important, and even a relatively “normal” tail as in Coelophysis can play a critical dynamic role in straight-line locomotion, in addition to more unsteady behaviors (, ). Much future work is therefore needed to understand how the locomotor role of the tail and other body segments (e.g., forelimbs) may have varied between different bipedal dinosaur clades. The diversity of tail construction observed among different clades, including short tails in oviraptorosaurs (), sustained caudal vertebral modification along the avian stem lineage (), and trellis-like ossified tendons in ornithopods (), suggests a diversity in dynamic tail functions that remains to be fully understood. Differences in other body proportions, such as very long arms in therizinosaurs () and deinocheirids (), also invite investigation of whole-animal functional integration across Dinosauria more broadly.

In the context of extant tailed saurians, the important role of the tail in locomotion has frequently been noted [e.g., (, )], but most such demonstrations are phenomenological only, stopping short of illuminating the biomechanical basis underpinning tail use. Few studies have mechanistically tested the underlying dynamics involved, and this has so far been principally concerned with movements in the sagittal plane (, ). Nevertheless, it is interesting to note that the timing of tail lateroflexion (with respect to hindlimb protraction-retraction) recovered for Coelophysis appears to be qualitatively opposite to that illustrated for locomotion in extant quadrupedal crocodylians () and quadrupedal and bipedal lizards (, , ), although this remains to be quantitatively investigated. This timing suggests that the patterns observed in extant saurians are not necessarily reflective of the dynamics of all extinct saurians, and further highlights the need for more mechanistically targeted analyses of this phenomenon. Furthermore, these differences hint that changes in mechanics associated with the adoption of parasagittal bipedality in dinosaur ancestors may have led to partial decoupling of the tail from the hindlimb, paving the way for their separation into distinct locomotor modules later in theropod evolution ().

MATERIALS AND METHODS

Tinamou musculoskeletal model

A previously published 3D musculoskeletal model of the tinamou was used (). This model is formulated in OpenSim (, ) and has a mass of 545 g, and each hindlimb is actuated by 36 MTU actuators, covering all the important muscles of the leg. Initially, each MTU was modeled with a Hill-type formulation of excitation-activation and activation-contraction dynamics (), incorporating previously derived architectural parameters, stiffness, and activation-deactivation time constants. Ground contact was modeled with a single contact sphere fixed to the digits segment with contact stiffness set to 150,000 N/m and dissipation set to 0.0387 N∙s/m, assuming that the latter scaled isometrically with body mass with respect to a previous human model (). Abduction-adduction and long-axis rotation degrees of freedom (DOFs) at the knee and ankle joints were omitted in the current study (i.e., fixed at the reference pose of 0°), because these DOFs play a secondary role in limb kinematics and are mostly controlled by passive forces (), making the simulations more tractable and comparable to those for Coelophysis. The head, neck, trunk, and (diminutive) tail of the tinamou were collectively modeled as a single rigid entity. This is justified, because birds have a very rigid trunk, as a result of partial vertebral fusion dorsally (notarium) and a massive rib-sternal complex ventrally, as well as reduction in the number of dorsal vertebrae compared to non-avian theropods. Moreover, during fast running, birds [including tinamous ()] extend their necks out from the body and do not appreciably move their head relative to the trunk. The final model as used in the current study had 18 DOFs.

We used recently developed numerical methods to represent the musculoskeletal model in a format conducive to rapid OCP solving. Variation in MTU length, velocity, and moment arms with respect to joint angles and velocities was represented with algorithmically differentiable polynomial functions (). Polynomials were fitted to the output from OpenSim’s MuscleAnalysis tool applied to 5000 randomly varying limb postures throughout the limb’s full range of motion (). System dynamic equations were represented using algorithmically differentiable implicit formulations (). The skeletal (rigid body) dynamics component of the musculoskeletal model was transcribed to an algorithmically differentiable OpenSim C++ source file (, ), where the states are generalized coordinates and velocities; we introduced additional “slack” controls that are the time derivatives of the generalized velocities (accelerations), and then imposed the nonlinear dynamic equations as algebraic constraints in their implicit form. As part of this representation, foot-ground contact was modeled with a smoothed implementation of OpenSim’s Hunt-Crossley formulation (, , ). Excitation-activation and activation-contraction dynamics of the MTUs were implemented with previously published models (, ), where the states are muscle activations and tendon forces. Here, a change of variables (, ) meant that slack controls are the time derivative of muscle activations and tendon forces. The real controls of muscle excitations are back-calculated from the solution (they are not involved in its determination) and are uniquely determined by the activations and their time derivatives.

Tinamou simulations

The OCP of generating a steady gait cycle was set up using direct collocation methods (, ). Full details of OCP formulation and solving are given in the Supplementary Materials. Briefly, the states and controls for a half-gait cycle (assuming left-right symmetry and stride cyclicity) that satisfied various dynamic and path constraints were derived, which minimized the value of an objective function incorporating muscle activations as an “effort” term. The OCP was constrained to satisfy some target forward speed, allowing us to explicitly control the speed of locomotion (, , ). Thus, we could not only simulate slow (walking) and fast (running) gaits in a controlled fashion, but also ascertain the model’s maximum capable speed by progressively increasing the target speed (in increments of 0.05 m/s) until the OCP was no longer able to converge. Our simulations used a “cold start” initial guess, where no prior knowledge of avian locomotor kinematics or kinetics was provided to the optimizer; the resulting simulation behavior therefore evolved entirely out of the underlying model dynamics and the objective function.

To assess how well our predictive simulations emulated reality, we ran simulations for walking and running at target speeds corresponding to those for previously published experimental data of overground locomotion in tinamou. The only difference between the OCPs for the two simulations was the imposed target speed. Our walking simulation moved at 0.39 m/s, corresponding to the walking trial of () for the same individual bird that the model is based on. Our running simulation moved at 2.62 m/s corresponding to the fastest running trial of (); note that, in that study, the fastest speed was actually 2.78 m/s but this was for a larger individual (790 g), and so, speed and other measurements were scaled to the model by virtue of mass and hip height, assuming isometry (). As very limited data (12 trials in total) exist for any speed of aerial running in tinamous, it is not currently possible to rigorously infer what a “representative” running trial is, nor to plausibly derive statistical confidence intervals; we hence chose to focus on the fastest recorded run, because this would better circumscribe the validity of the simulation framework insofar as speed was concerned. In addition to comparison of kinematics and kinetics between simulation results and empirical data, we also compared the time histories of simulated muscle excitations with previously published electromyographic data for birds, including tinamou ().

Typical Hill-type models of muscle contraction () rely on architectural parameters that are not preserved in the fossil record (e.g., fiber length). Various approaches to estimating such parameters have been proposed, but recent testing of their underlying assumptions has raised doubts as to how well these approaches work for all muscles across all species (, ). Here, we explored a novel approach to whole-muscle modeling, which is not contingent on architectural parameters and could be applied to extinct species. This simplified model represents MTU force as an active muscle force along a line of action (determined by activation level, a), supplemented with force due to passive stretching of parallel and in-series collagenous tissues (e.g., tendon and epimysium), determined as a function of MTU length (fig. S2A)whereis the component of total force due to passive stretch. Equation 2 expresses passive stretch as a linear function (with slope or “MTU stiffness” of kMTU) of the relative MTU length (lMTU*), modulated by a hyperbolic tangent function so as to produce zero force at relative MTU lengths less than one. Relative MTU length is the ratio of the current (instantaneous) MTU length to the “resting” MTU length, taken as the length when all limb DOFs are set midway between their upper and lower bounds on movement. This simple model therefore represents MTU behavior in a phenomenological fashion, and it does not explicitly attribute passive stretch to any specific part of the MTU; there is no literal (anatomical) interpretation of what the value of kMTU means. Note that if kMTU is set to zero, the model degenerates to just the active force along a line of action.

We explored how well this simple muscle model performed in the maximal speed simulation by setting different values of kMTU and iteratively testing whether the OCP could solve at different values of Fmax. Here, Fmax for every muscle was set to a single, constant multiple of BW, thus ascertaining the minimal strength required for the simplified model to replicate the fastest running gait. Arguably, strength scalers of lower multiples of BW are better: We conservatively did not give the model any extra strength (Fmax) than it needed. This assumed that maximum speed running is the most demanding behavior in terms of muscle force production, which is defensible (). In addition to these tests, the resulting simulation kinematics and kinetics were also compared to those for the maximal speed simulation using Hill-type muscles.

Coelophysis musculoskeletal model

Once benchmarked with an extant theropod, the above framework (including the simplified muscle model) was applied to the extinct, non-avian theropod C. bauri to generate a maximal speed running simulation. This relatively small [10 to 25 kg; ()], Late Triassic (~220 Ma) taxon is the best-known early theropod dinosaur, known from innumerable complete and partial skeletons. It is a good representative of the general bauplan ancestral for Dinosauria and had many anatomical features common to most bipedal dinosaurs in general, including long hindlimbs, short forelimbs, digitigrade stance, a long tail, a relatively long neck, and a small head (). Its relatively small mass (in dinosaurian terms) is closer to that of early dinosaurs including theropods () and is within the range observed in extant bipeds, thus avoiding the potentially confounding (allometric) effect of much larger body size as would occur with some other well-studied theropods (e.g., Tyrannosaurus). Coelophysis, therefore, is an ideal taxon for investigating locomotor biomechanics and its evolution in early bipedal dinosaurs, as well as the evolution of “cursoriality” in bipeds in general.

We used a previously developed model of Coelophysis (), with a few changes (Fig. 3, A and B, and the Supplementary Materials). The model has three DOFs at the hip (flexion-extension, abduction-adduction, long-axis rotation) and one each (flexion-extension) at the knee, ankle, and MTP joint, as well as 33 MTUs actuating each hindlimb. The forelimbs, which were free to move in the original model, were programmatically locked in a “tucked” position with respect to the body (shoulder extension angle fixed at 45°; see Fig. 2); their very small size (<0.5% of total body mass, shoulder-fingertip distance about a third of hip-toetip distance) indicates that they did not play an important role in locomotion and thus their movements could be safely ignored. Indeed, the forelimbs could have been used to hold food during locomotion. The neck and head formed a single segment, and its joint with the chest was assigned two DOFs (dorsiflexion and lateroflexion). The tail of the original model was split in two at midlength, and both proximal and distal tail joints were assigned two DOFs (dorsiflexion and lateroflexion). In addition, as Coelophysis and other basal theropods have dorsal vertebrae with pre- and postzygopophyses that are strongly inclined laterally (, –), this suggests capability for a modest amount of body lateroflexion, which may have been important. Lizards can exhibit marked levels of body lateroflexion during bipedal gait, as judged by published still images in dorsal view (). Thus, the originally singular body segment was split into a cranial “thorax” and caudal “trunk” segment, with the division parallel to the coronal (Y-Z) plane and positioned at the first intervertebral joint cranial to the pubes. These two segments were joined by a back joint that permitted lateroflexion only; dorsiflexion was probably limited given the vertebral morphology and the dorsoventrally deep extent of the body. Thus, the revised model used here comprised 25 free DOFs and had a mass of 13.1 kg. As with the tinamou, ground contact was computed using a single large contact sphere affixed to the middle of the digits segment. On account of the marked difference in size between Coelophysis and the tinamou, the values of several parameters of the system were scaled up from those of the tinamou, assuming isometry (): muscle activation and deactivation time constants, contact stiffness and dissipation, and joint damping.

As with the tinamou model, variation in MTU length, velocity, and moment arms with respect to joint angles and velocities was represented with polynomial functions, calculated from 5000 randomly varying postures for the hindlimb (six DOFs) and proximal tail (two DOFs). The bounds on the postures sampled were more restrictive than the “bones-only” ranges of motion reported previously () and were set based on what was considered plausible for in vivo running poses given the kinematics of extant bipeds [e.g., (, )], as well as MTU path behavior in more extreme poses. We therefore favored using a more restricted range of limb poses, rather than modifying the MTU paths to behave properly at (likely nonphysiological) extreme poses. For both the ambiens and fibularis longus muscles, the part of the MTU representing the inferred secondary tendon to the digital flexors was excluded for the sake of model simplicity.

Our approach to simulating Coelophysis gait used the simple muscle model developed above and closely mirrored the approach for the tinamou, with additional nuancing. The extra DOFs in the neck, back, and tail joints necessitated additional cyclicity and symmetry path constraints. Furthermore, as muscles actuated only the hindlimbs in our simulations, the active control of the axial DOFs necessitated the use of ideal torque motors, whose applied torque was proportional to activation, and whose excitation-activation dynamics was described by a simple first-order differential equation (). The use of torque motors was chosen primarily for expediency: They are simpler to implement, which reduces dimensionality and nonlinearity in the OCP (facilitating quicker convergence) and means that their effects on the resulting simulation behavior are more straightforward to interpret. Involving a more complex, MTU-driven back and tail would not just greatly increase OCP complexity, but might only be reasonable to do so if each intervertebral joint was modeled separately (as opposed to the coarse segmented approach used here) because axial musculature will cross multiple successive joints; it would also involve additional assumptions about postcervical axial musculature in non-avian theropods, which needs further study. The simulation framework developed here can be built upon by future studies that more finely target vertebral mechanics, either by incorporating MTU actuators explicitly or by using physiologically inspired torque actuators comprising active and passive elements [e.g., ()]. It is difficult to predict what refinements such nuanced studies may gain, and this may depend on whether each intervertebral joint of the back and tail was modeled separately.

The proximal tail joint was also actuated by the CFL muscle, which ran from near the tail base to the fourth trochanter of the femur. Whereas the hip, knee, and ankle DOFs were actuated solely by MTUs, the MTP joint also had a torque motor appended to it, acting as a “reserve’ on account of only two MTUs crossing the joint due to uncertainty in distal limb muscle reconstruction (, ). The objective function used sought to penalize both muscle and torque motor activations (see the Supplementary Materials) and additionally kept the head and foot placements close to the body midline (). Maximum speed locomotion was derived by progressively increasing the target speed (in increments of 0.05 m/s) until the OCP would no longer converge to a striding gait cycle, and the MTP torque motors started contributing more than 50% toward the external joint moment at the middle of the stance phase. As there is uncertainty in modeling actuation of the MTP joint—would the “missing” musculature not modeled here have contributed at most 50% of joint moment?—the maximum speed obtained here should be viewed as an approximation only.

Following the nominal simulation, we explored the sensitivity of the simulation behavior to changes in the craniocaudal position of the whole-body COM and location of the tail joints, because the former may at least influence limb posture and running performance (, , ). We altered the craniocaudal position of the COM of the trunk and thorax, as these segments are both relatively massive and likely have the greatest degree of error in flesh (and air space) volume estimation in extinct sauropsids (). For both segments, the COM was shifted craniocaudally by ±50% of their original location with respect to the segment origins (hips for trunk segment and back joint for thorax segment; fig. S4D). In addition to this, we investigated simulation behavior when a different tail configuration was used (fig. S3C). Rather than being located at the base of the tail (between the last sacral and first caudal vertebra), the proximal tail joint was moved distally, to between the fourth and fifth caudal vertebrae, emulating a tail whose base is more firmly restricted in mobility by its connection with the pelvis. The mass of the proximal-most part of the tail was therefore incorporated into the trunk segment. The distal tail joint was also moved distally so as to remain halfway between the proximal tail joint and the tail tip.

Coelophysis tail lateroflexion tests

In the predictive simulation of fast running in Coelophysis, the tail underwent appreciable lateroflexion throughout the gait cycle; given that previous discussions and analyses have treated the tail as a static counterweight caudal to the pelvis, this was an unexpected result that prompted investigation. We identified two hypotheses that could explain why sizeable tail lateroflexion may occur:

H1. If lateroflexion oscillations were appropriately timed with respect to hip flexion-extension oscillations, the CFL muscle might lengthen beyond its “rest length” at some key point in the gait cycle. According to the simple muscle model used (with kMTU > 0), passive stretch would increase CFL force output, enabling it to produce a greater hindlimb extensor moment, facilitating better support or propulsion of the animal.

H2. Lateroflexion of the tail helps regulate whole-body angular momentum, helping to counter that produced by the cumulative motions of other body segments such as the hindlimbs or neck and head.

In addition, we identified a third hypothesis:

H3. Regardless of the proximate reason for lateroflexion oscillations, these oscillations were actively controlled to maximize their effect (whatever that may be).

To test the above hypotheses, 10 variations on the nominal OCP were conducted:

1) “Passive tail”: All active and (non-inertial) passive forces were removed from the tail such that the tail moved passively as a consequence of inertial transfer via pelvic motions. The CFL remained in the OCP, but its relationship with the tail was removed, effectively modeling the muscle as though it originated from a rigid caudal extension of the pelvis.

2) “No stretch CFL”: The potential for force enhancement via passive stretch was omitted in the CFL by setting kMTU of just this muscle to zero.

3) “No CFL”: The CFL was removed entirely from the OCP so that the proximal tail joint was actuated solely by torque motors. However, the force of the CFL was “transferred” to the caudofemoralis brevis (CFB; a similar muscle that also runs to the femur but from the caudal pelvis) so that hip extensor capabilities were not markedly altered. Except for its strength, no other aspect of the CFB’s action was modified.

4) “Lighter tail”: The proximal tail segment was made 0.75× heavier with attendant changes also made to its inertia tensor. If the proximal tail was made any lighter, the MTP torque motor was recruited more than the allowable 50% threshold; thus, ×0.75 was the greatest viable decrease.

5) “Heavier tail”: The proximal tail segment was made 1.8× heavier with attendant changes also made to its inertia tensor. If the proximal tail was made any heavier, the MTP torque motor was recruited more than the allowable 50% threshold; thus, ×1.8 was the greatest viable increase.

6) “Rigid axis”: All neck, back, and tail joint angles were constrained to bounds of ±1°, effectively locking these joints to create a single axial body segment as in previous analyses of theropod locomotion.

7) “Rigid back”: The back joint angle was constrained to bounds of ±1°, effectively locking the joint to create a rigid body segment and more closely emulating the stiffened condition of birds.

8) “Rigid tail”: All tail joint angles were constrained to bounds of ±1°, effectively locking these joints to create a rigid tail that moved as a rigid caudal extension of the pelvis.

9) “Synchronized tail”: An additional term was introduced into the objective function that encouraged hip extension and proximal tail lateroflexion to oscillate with respect to one another in a fashion opposite to that observed in the nominal simulation (see the Supplementary Materials).

10) “No tail”: The tail segments were removed from the model. However, simulation behavior may be markedly affected just by the attendant reduction in model mass and cranial shift in whole-body COM. To compensate for this, the tail’s mass was added to the trunk segment (and moments of inertia were rescaled appropriately), and the trunk’s COM was modified such that the whole-body COM remained unaltered in the reference pose (all coordinates set to zero). In addition, the CFL’s origin was transferred to the trunk segment, but its location in the trunk segment coordinate system remained unaltered; thus, it effectively originated from a (massless) rigid caudal extension of the pelvis.

Variant 1 examined whether the tail’s observed behavior was a consequence of the objective function encouraging some preferred pattern of active torque production (muscle or torque motor), thus testing hypothesis H3. Variants 2 and 3 helped test hypothesis H1. Variants 4 to 9 helped test hypothesis H2 and, to a lesser degree, hypothesis H3. In particular, variants 4 and 5 determined the influence of the amount of available inertia that the tail could use to generate angular momentum; variants 6 to 8 identified the consequence of the axial segments being less able to generate whole-body angular momentum through active, independent movements (although torque motor actuation of the intervening joints was still required for moment balance); and variant 9 tested the importance of the timing of angular momentum generation by the tail with respect to that by other body segments. In emulating a tail-less Coelophysis, variant 10 provided a further test of hypotheses H2 and H3, as well as insight into the importance of the tail itself for locomotor mechanics in a bipedal dinosaur.

For each variant as well as the nominal simulation, the resulting whole-animal kinematic trajectories were examined, and peak pelvic, back, and tail oscillations, in particular, were recorded. In addition, four further aspects were computed:

1) The phase lag (offset) between proximal tail lateroflexion and left hip extension, the latter taken as a proxy for limb protraction-retraction motion. This was computed using cross-correlation analysis of the time histories of both angles and normalized to a full cycle of 360°; 0° indicates that the oscillations were perfectly in phase, while 180° indicates that the oscillations were perfectly out of phase with respect to each other. Right hip extension could equally be used, but then the meanings of “in phase” and “out of phase” would be offset by 180°.

2) The phase lag between pelvic yaw kinematics in a given variant and those of the nominal simulation, computed using cross-correlation analysis as above.

3) The mean activation of muscles (for one leg) at each instant in the gait cycle, as well as the individual contributions of the muscle and torque motor activation terms to the objective function. The latter represents the “effort” term in the objective function and is broadly comparable to the “integrated activation” metric of ().

4) Contributions of each segment to the total angular momentum about the whole-body COM, measured in the global coordinate system (see the Supplementary Materials). This was investigated both in terms of the time histories of angular momentum about each axis and through PCA to quantitatively describe angular momentum cancelation from different body segments ().

Sometimes, different configurations of the OCP meant that the model’s athletic performance was reduced: It was unable to achieve the top speed of the nominal simulation. As preliminary tests showed that faster running models tended to use more “vigorous” tail oscillations, all the variants (and the nominal OCP) were conducted at the same speed to remove the potentially confounding effect of different speeds producing different “intensities” of kinematics. This speed, 5.7 m/s, was the fastest speed that all variants were able to achieve, and although it is a modest (14%) decrease from the maximal speed of the nominal OCP (6.65 m/s), it still corresponded to a fast running gait (Froude number ≈ 6.6).