Maximum running speed is an important locomotor parameter for many animals—predators as well as prey—and is thus of interest to palaeobiologists wishing to reconstruct the behavioural ecology of extinct species. A variety of approaches have been tried in the past including anatomical comparisons, bone scaling and strength, safety factors and ground reaction force analyses. However, these approaches are all indirect and an alternative approach is to create a musculoskeletal model of the animal and see how fast it can run. The major advantage of this approach is that all assumptions about the animal's morphology and physiology are directly addressed, whereas the exact same assumptions are hidden in the indirect approaches. In this paper, we present simple musculoskeletal models of three extant and five extinct bipedal species. The models predict top speed in the extant species with reasonably good agreement with accepted values, so we conclude that the values presented for the five extinct species are reasonable predictions given the modelling assumptions made. Improved musculoskeletal models and better estimates of soft tissue parameters will produce more accurate values. Limited sensitivity analysis is performed on key muscle parameters but there is considerable scope for extending this in the future.
Chasing down prey is a vital factor in the lives of extant predators, as is the avoidance of being captured for prey animals. It is therefore of little surprise that speed estimation is of such interest to palaeobiologists who study dinosaurs. The range of predicted speeds is as variable as the methods chosen with some authors favouring high speeds (Paul 1988, 1998) while others prefer moderate (Farlow et al. 1995) or low speeds (Alexander 1989; Hutchinson & Garcia 2002). Current analysis techniques are based on anatomical comparisons, bone scaling and strength, risk factors and ground reaction forces, and a recent review (Hutchinson & Gatesy 2006) summarizes the current state of the art and concludes, among other things, that a ‘rigorous dynamic simulation of a moving dinosaur, one encompassing all motions and forces, cannot yet plausibly be done’. However, such an approach is clearly the best option because it both explicitly requires a complete set of modelling assumptions and is conceptually simple. It requires a musculoskeletal model to be constructed necessitating assumptions about skeletal geometry, body mass and mass distribution, together with muscle and tendon properties. These properties all have a considerable effect on locomotor performance and all that happens in any quantitative technique that does not make explicit assumptions about their values is that they will be implicitly scaled from whatever reference species are employed. This is likely to be suboptimal given that some of these values can often be inferred from fossil evidence directly, and unknown values need to be clearly identified so that their importance to the final result can be assessed.
It is now possible to produce highly detailed musculoskeletal models and these can be used for fossil locomotor reconstruction (Nagano et al. 2005) and such could be constructed for dinosaur species. However, such complex models require a great deal of work to construct and good predictive results can be obtained from rather simpler models (Yamazaki et al. 1996; Ogihara & Yamazaki 2001; Sellers et al. 2003, 2004, 2005). These latter models have the great advantage of requiring no knowledge of locomotor kinematics and allow the generation of a range of gaits de novo, which maximize or minimize some global optimization parameter such as speed or energy consumption. These models use an evolutionary robotics approach where so-called evolutionary algorithms are used to find muscle activation patterns (Nolfi & Floreano 2000). However, such models are still computationally demanding. It is perfectly possible to solve the mechanical constraints of the system (internal forces generated by muscles and the spring recoil of tendons, segment movement constraints imposed by joints, external forces generated by gravity and ground reaction through contacts with the environment) in more or less real time on a modern computer using freely available software. However, finding the activation pattern of the muscles that produces high-quality gait is extremely challenging. Even a simple model such as ours with 12 muscles and 5 activation levels per step (half a gait cycle) leads to 61 dimensions in the search space, which is consequently far too large to search exhaustively. Fortunately, using a parallel implementation of a suitable non-exhaustive, evolutionary search procedure such as a genetic algorithm means that even this is a tractable problem that can be solved in days or weeks on a modern supercomputer.
It is important that any model attempting to predict the behaviour of extinct species be tested on extant animals. Thus, a model attempting to predict the top speeds of a range of fossil bipedal dinosaurs should also be tested using equivalent data from living bipeds. However, there is remarkably little high-quality top-speed data for any living species other than those actively involved in racing, such as humans, horses and greyhounds (Alexander 1989) and even elephants (Hutchinson et al. 2006). The values used for comparison in this paper (Alexander et al. 1979; Patak & Baldwin 1998) are based on anecdotal observations and need to be treated with some caution since there is no way of knowing whether these animals were actually running as fast as they could. The anatomy, posture and gait of bipedal dinosaurs is unique, with no equivalent modern analogue available for the comparison of locomotor abilities. The successful testing and validation of the computational approach using extant species indicates that it is possible to generate a robust model applicable to extinct species.
2. Material and methods
Our previous models (Sellers et al. 2003, 2004, 2005) used the Dynamechs simulation library (McMillan et al. 1995). However, this system has very limited support for contacts between the simulation and the environment so, for this model, we switched to using the open dynamics engine (ODE; http://www.ode.org) to provide the physics simulation. This is again a C++ library so we were able to modify our existing GaitSym code to use the new simulator. The models themselves were specified in a custom XML format that defined the necessary segments, joints, muscles, tendons and contacts. To perform the optimization, we modified our existing distributed genetic algorithm system to support the new format. While it is difficult to compare the results using the two systems, it is certainly our impression that ODE is both faster and more numerically stable than Dynamechs in this application. Communication between the optimization code and the simulation code was performed using sockets via the PTypes abstraction library (http://www.melikyan.com/ptypes) to allow easy portability over the variety of computer systems available to us, and to allow the code to take advantage of multiple, remote computer clusters, using over 300 processor cores when available. Graphical display and user interaction used the GLUI OpenGL widget library (http://glui.sourceforge.net/) to again allow easy portability. The muscle model was derived from Minetti & Alexander (1997) to generate both eccentric and concentric velocity-dependent force characteristics extended with the addition of Hill-style linear serial and parallel elastic elements (Hill 1938). The combined musculotendinous unit is solved analytically and exported as a standard C function using Mathematica (http://www.wolfram.com) that greatly aids numerical stability. A cylindrical wrapping operator was used where necessary to maintain the muscle and tendon path around joints.
A set of three extant and five extinct bipeds were modelled. These are two-dimensional models with a rigid trunk, and left and right thigh, shank and composite foot segments. These segments are linked using three hinge joints per limb. The segment properties are based on published data (Hutchinson 2004a,b) with the single composite foot segment combining the metatarsus and foot with all floor contact occurring at the distal end of the metatarsus. The species were chosen to cover a reasonable size range and are listed in table 1. The published dataset does not include moments of inertia so these were calculated by modelling the segments as geometric shapes chosen to match the published lengths, masses and centres of mass. Conic segments were used for the limb segments and back-to-back circular cones were used for the trunk. Muscle fibre lengths, physiological cross-section areas (PCSAs) and moment arms were available for limb extensors (Hutchinson 2004a,b), and these were used to create appropriate muscle paths on the model. Cylindrical wrapping operators were used for the knee and ankle extensors. Limb flexors were assumed to have the same properties, but only 59% of the mass was based on human proportions (Pierrynowski 1995) but not dissimilar from the values found in other species (Alexander et al. 1979; Maloiy et al. 1979; Bennett 1996).
The extensor muscle mass was limited to 5% of the body mass per joint (Hutchinson 2004b) and all muscles were considered to act over a single joint, with each joint having a single flexor and extensor. Muscle volume was calculated using the standard value of 1056 kg m−3 for muscle density (Winter 1990) and PCSA was calculated by dividing this volume by the fibre length. Force per unit area was chosen to be 300 000 N m−2 (Hutchinson 2004b), but there are other values in the literature: Umberger et al. (2003) use 250 000 N m−2, and Alexander (2003) reports an in vitro maximum value of 360 000 N m−2 for frog and 330 000 N m−2 for cat for parallel-fibred leg muscles. Zheng et al. (1998) recommend a value of 400 000 N m−2 for human quadriceps and Pierrynowski (1995) suggests 350 000 N m−2. There is a similarly large range for maximum contraction speed. Winter (1990) suggests values from 6 to 10 times the muscle's resting length per second for humans. This value is clearly highly dependent on both the fibre-type composition of the muscle and the temperature. Westneat (2003) reports a range of values for fishes from 3 to 10 s−1 for different fibre types and Umberger et al. (2003) recommend values of 12 s−1 for fast twitch and 4.8 s−1 for slow twitch. A value of 8 s−1 was chosen to represent a mixed-fibred muscle. Joint limits were set to extreme values for all species since these limits are not generally available and the range of motion permitted by the muscle and tendon lengths specified should be sufficient to limit joint excursion. Serial and parallel elastic constants were set so that the serial element strain was 6% at the maximum isometric contraction and parallel element strain was 60% for the same force as used in a previous model (Sellers et al. 2005), and the lengths of the tendons were scaled based on average human proportions for the hip, knee and ankle (Pierrynowski 1995). Muscle attachment points were arranged so that the muscle plus tendon was at its resting length in the published mid-stance positions (Hutchinson 2004a,b), so that the moment arm was approximately correct. This procedure was performed automatically using a custom Matlab program (http://www.mathworks.com/). The full specification for each of the models is included as human-readable XML files as electronic supplementary material.
Gait is generated using a distributed genetic algorithm optimization system using a genome that represents the gait cycle duration and the muscle activation levels at 10 time periods through the gait cycle. The genome contains 61 parameters representing the cycle time and the 5 activation levels for 12 muscles for half a gait cycle. The left–right activation levels are simply swapped for the other half of the gait cycle. This is implemented as a client–server architecture with the simulators running on multiple client machines and a central server that gathers the results from the multiple simulations. The starting conditions used published mid-stance positions with a trunk forward velocity giving a Froude number of approximately 1.5 (which equates to a medium running speed) using the formula velocity=1.5×(9.81×(thigh length+shank length+metatarsus length)) following Alexander (2003) but using leg length as a proxy for hip height. Forward velocities for the support leg segments decreased to zero depending on the height of the segment from the ground, and the velocities of the swing-leg segments increased to double the trunk forward velocity again depending on height. The fitness criterion used for the genetic algorithm optimization was the maximum forward distance achieved in a fixed time (3 s for most species but 5 s for the Allosaurus and Tyrannosaurus to allow a reasonable number of complete gait cycles). Thus, runs where the animal fell over scored very badly and the runs with the highest average speed would score the highest. The population was 1000 and up to 1000 generations were allowed in each run unless a steady maximum average forward velocity was achieved earlier. This procedure was repeated at least five times until a good quality run was obtained. Runs were judged to be good quality when the animal did not fall over within the time limit and managed at least 15 m forward movement. The best run was then used as the basis for a gait morphing procedure (Sellers et al. 2004), where the best results for a previous run were used to generate the starting conditions for subsequent runs. This procedure was repeated at least 20 times to obtain the highest speed estimate for each species. An individual simulation ran in approximately real time but at least 1 000 000 repeats were needed to generate the optimized running gait for each species.
All the models generated high-quality running gaits that were stable over the whole simulation period as shown in figure 1. This figure shows a series of snapshots of the generated gaits at 1/10 cycle time intervals. The lines connect the joint centres, except for the foot line that is drawn between the ankle joint and the contact point on the metatarsal head and the trunk line that is drawn from the hip joint to the position of the centre of mass of the trunk segment.
The top speeds achieved for each species are shown in table 1 and there is reasonably good correspondence for the extant species between the speeds generated by the model and those cited in the literature: Dromaius 14 m s−1 (Patak & Baldwin 1998); Struthio 17 m s−1 (Alexander et al. 1979); and Homo 10 m s−1 (Alexander 1989). To investigate whether mass alone is a good predictor of top speed, we plotted the top speed against body mass in figure 2, which shows a monotonic reduction in speed for the extinct bipeds but not for the extant species. The extinct species have a consistent 15% of body mass as muscle mass for each limb, whereas the value for Dromaius is 15.5%, for Struthio is 12.2% and for Homo is 10.6%, so this alone does not explain the speed differences seen in the extant species.
Similarly, there is a well-known empirical relationship between speed, leg length and stride length based on Froude numbers that is often used for trackway analysis: λ/h≈2.3(u2/gh)0.3, where λ is the stride length; h is a characteristic length (in this case, leg length); u is the forward velocity and g is the acceleration due to gravity (Alexander 1976). This information is shown in table 1 where the characteristic length used is leg length (the sum of thigh, shank and metatarsus lengths). Figure 3 shows a plot of λ/h against u2/gh for the models to compare the simulations against the empirical relationship demonstrating extremely good agreement.
A key area of uncertainty in our model is the choice of muscle parameters, so a simple sensitivity analysis was performed. The maximum force available per muscle is directionally proportional to the mass of muscle (Fmax=δM/ρl, where δ is the force per unit area; M is the muscle mass; ρ is the muscle density and l is the muscle fibre length). Of these the greatest uncertainty is likely to be the value of muscle mass (Hutchinson & Garcia 2002); hence, this was varied over a range of 2.5–7.5% per joint, which should encompass the probable variation. Similarly, Vmax is known to vary widely and it may scale negatively with body mass (Medler 2002), although there is little information about contraction velocities in large animals. Values between 4 and 12 s−1 were tested to both encompass the probable range and give the same degree of variation (±50%) as used for muscle mass. The results of the sensitivity analysis on the Tyrannosaurus rex model are shown in figure 4. Both parameters have an approximately linear effect over this range, and the effect of muscle mass is approximately double the effect of contraction velocity. However, the model was unable to sustain locomotion when the muscle mass was reduced to 2.5% per joint.
Multibody dynamic simulations using evolutionary robotic optimization approaches appear to provide reliable estimates for the maximum running speeds of extant animals. Multiple simulations with small changes in both starting conditions and muscle activation patterns produced highly consistent estimates. Maximum running speed is a highly variable character and difficult to estimate in any animal. The observed speed is always a lower bound estimate of maximum speed since the animal might not be running as fast as it can. In addition, many quoted values for running speed are based on observations made in less than ideal conditions, which may lead to considerable errors (Garland 1983; Alexander 1989). Even for humans, the situation is not straightforward: while 200 m sprint averages are in excess of 10 m s−1, the peak speed reached can be in excess of 12 m s−1 (Brown et al. 2004), and these are values for elite athletes who have considerably greater leg muscle mass than the average values used in our simulations. Tests on more general female athletes from other sports give typical speeds of approximately 6 m s−1 with short bursts of less than 8 m s−1 (Brown et al. 2004). Our estimates are broadly in line with other biomechanical estimation techniques, which predict 18 m s−1 for ostrich and 13 m s−1 for emu (Blanco & Jones 2005). It is self-evident (and has been demonstrated in various models; Hutchinson & Garcia 2002; Sellers & Paul 2005) that changes in muscle mass will affect maximum speed and this is a major source of uncertainty in these predictions. However, we would propose that limits can be set using functional bracketing to set minimum and maximum values. The 15% value used here is similar to that found in large extant bipeds and as our sensitivity analysis shows, an appreciably smaller value does not even allow the animal to walk. An upper limit would be harder to estimate but one approach that is possible (even if computationally expensive) is to calculate the bone loading during high-speed locomotion and compare the safety factors with those known for extant animals (Alexander 1997). For these models, the general decrease in top speed as body size increases is certainly in line with predictions made elsewhere (Hutchinson 2004b) and the predicted top speeds, with the possible exception of 17.8 m s−1 for Compsognathus, are not exceptional. However, we would expect that the actual top speeds are likely to be somewhat higher than those given here since it is probable that improvements could be made by altering the distribution of the leg muscle and optimizing the fibre and tendon lengths. In addition, all the models are relatively simplistic and lack multiple-joint muscles and accessory elastic storage structures, which when combined should increase the maximum speed. Our initial sensitivity analysis shows that changing our assumptions about the muscles has a considerable effect on our estimates. A 50% increase in muscle mass leads to a 60% increase in top speed and a 50% decrease stops the model working at all. Similarly, a 50% increase in maximum contraction velocity leads to a 30% increase in top speed with a 50% decrease leading to a 20% decrease. However, there is considerable scope for further sensitivity analyses such as variation in muscle attachment points (Sellers & Crompton 2004) and the effects of perturbation (Wilson et al. in press).
The kinematics generated by the models show a wide range of variation. It is very hard to judge the probable accuracy of these values without subject-specific kinematic matching, and in any case with such simple models it may be optimistic to expect high-quality kinematics. In particular, we would expect better estimates of muscle and tendon lengths (and to a lesser extent joint limits) to have a considerable effect on kinematics. However, as shown in figure 3, the mechanical accuracy of the simulations is high. The models represent a wide range of Froude numbers and follow the predicted stride length relationship very closely indeed. Overall, simulations such as ours illustrate how an animal could have moved given its physiological and morphological constraints, and perhaps also indicate probable movement patterns, but we are still some way off saying that this is how it must have moved.
There are a number of possible future developments in animal gait simulation. Our models already allow us to calculate the energetic costs of gaits in some detail. This includes estimates of metabolic and mechanical power, including active contraction and passive elasticity contributions on a muscle by muscle basis. However, these predictions need better validation and this can be achieved using a combination of species-specific modelling and in vivo physiological measurement. As computers get faster and cheaper, it would be useful to increase the biofidelity of the models by adding more muscles, including the forelimb in the simulation, and having considerably greater detail in the feet and the interactions with the substrate. In addition, reconstruction improvements such as new estimates of inertial parameters (Hutchinson et al. 2007) and better comparative anatomical and physiological data will certainly help. In particular, very little comparative data are available on the elastic properties of muscles. Similarly, the addition of a feedback-based control system will increase the range of locomotor activities that can be encompassed and allow this technology to look at non-continuous activities, such as starting, stopping, cornering and complex movements such as sitting and standing up.
We would like to thank Prof. Robin Crompton for allowing access to his Beowulf Cluster to perform the required simulations; Cliff Addison for arranging access to the NW-Grid; funding from BBSRC, NERC and the Leverhulme Trust; and two anonymous referees for their helpful comments on the original draft of the manuscript.