March of the Titans: The Locomotor Capabilities of Sauropod Dinosaurs

Skeletal Capture

The initial stage in building the simulation is construction an appropriate musculoskeletal model. The first step is to acquire a digital model of the skeleton of the target species. In this case, our aim is to explore the locomotor capabilities of the largest of the sauropod dinosaurs and we chose the to use Argentinosaurus huinculensis, as reconstructed by the Museo Municipal Carmen Funes, Plaza Huincul, Argentina, which also houses the original fossil material. Permission was granted by Museo Municipal Carmen Funes, Plaza Huincul, Argentina to scan their reconstruction. The reconstruction was performed in-house at the museum. This reconstruction is shown in Figure 1. It is 39.7 m long and stands 7.3 m high at the shoulder. The reconstruction is based on rather fragmentary material [45] but includes well preserved fibula and vertebral elements that have allowed mass estimates to be obtained of between 60 and 88 tonnes depending on the regression equation used [46]. The reconstruction was scanned using a Z+F Imager 5006i LiDAR scanner from multiple locations in the gallery. The individual scans were aligned by Z+F Germany, using the multiple printed targets placed around the gallery as automatically detectable shared reference points. The tail, torso, neck and head and the individual limb bones and girdles were segmented out and decimated using of Geomagic Studio (www.geomagic.com) and the resultant 3D objects posed using 3DS Max (www.autodesk.com). The quality of the scan is variable due to limitations on where the scanner could be placed. Therefore limb bones on the side that had been better scanned were mirrored to produce a completely symmetrical model and the torso was moved slightly so that its centre of mass was exactly in the midline. This produced the reference pose illustrated in Figure 2. It was not possible to raise the scanner above floor level so the quality of the scan for dorsal elements such as neural spines is relatively poor. However the limb bones and girdles are well digitised and these are the most important in terms of subsequent modelling steps.

Segmental Mass Properties

Once the skeleton has been captured it is necessary to define the body segments that are used in the simulation. In common with nearly all locomotor analysis, the body is treated as a series of rigid, linked segments [47]. As in all modelling exercises it is necessary to decide on the level of complexity that is going to be used. It is perfectly possible to model every single bone as a separate segment but doing so greatly increases the calculation time for the simulation and having a large mass difference between body elements tends to cause numerical instability. For the sauropod model, 3 segments were defined for each limb representing the stylopodium, zeugopodium and autopodium. The head, neck, torso and tail were considered a single combined segment. Each segment is a six degree of freedom rigid element that has a position and orientation as well as a mass and inertial tensor. In the reference pose, the position is defined as the position of the centre of mass of the segment, and the orientation is set to a rotation of zero, with the inertial tensor calculated at this orientation. In the palaeontological literature there are two approaches for generating mass properties. Firstly these can be scaled from experimentally derived data of similarly shaped modern species and this is probably the commonest approach among hominoid workers (e.g. [32], [35]) with reference data from humans [48], [49] or chimpanzees [50]. Secondly these can be obtained from volumetric models of the target animal [51]–[53]. The modern locomotor analogues for dinosaurs have very different body shapes so the scaling approach is probably less useful than the volumetric approach. However whilst these are based on external body measurements when used with living animals, for fossil animals these soft-tissue measurements cannot be measured directly. This leads to an undesirable subjective element to these reconstructions and in an attempt to improve on this we have developed an objective technique based on convex hulling [54]. In its original form, this technique produced a mathematically unique minimum wrap around the individual skeletal components to estimate body mass. However since these are simply closed 3D shapes, all the other mass properties can also be calculated. The only difficulty is that our previous analysis found that approximately 20% of the mass was lost in the minimal wrap and this needs to be recovered. Figure 3AB shows the results of convex hulling the skeletal elements. The main place where the segments are clearly far too small is the thigh and upper arm and so the missing mass was added to these segments by using an appropriate scale factors orthogonal to the long axis of the bone. Figure 3CD shows the effects of this scaling. This choice of where to put the extra mass is somewhat arbitrary but it is believed that at low speeds, the choice of mass properties in the limbs is relatively unimportant [55]. The calculated mass properties for each segment in the reference pose are shown in Table 1. The total calculated body mass for the reconstruction using convex hulling approach [54] is 83,230 kg which is within the range previously predicted for this species [46] and certainly helps us have confidence in the reconstruction. However it must be remembered that these values are necessarily estimates. We do know how much soft tissue was associated with the skeletal segments and these estimates are means based on a limited dataset of modern animals. However we also know that the choice of mass parameters has relatively little effect on experimental [55] or simulation outcomes [33], [56].

Muscle and Joint Locations

From the reference skeleton it is now possible to define the joints and muscle paths, although there will always be ambiguities in specific cases. As with the choice of segments, it is necessary to simplify these to prevent undue model complexity. The joints were therefore all considered to be hinge joints operating in various parasagittal planes (i.e. with hinge axes directed laterally), with the joint centre measured from the skeleton. This is probably reasonably accurate for all the joints except the shoulder and hip joints, which should be ball-and-socket joints. However it is likely that there is very little abduction/adduction or axial rotation in normal walking so this is a reasonable approximation for a model of straight line walking and greatly simplifies the control processes. The joints chosen are listed in Table 2. It is also necessary to define contact points on the skeleton which are simply the parts of the feet that make contact with the ground. The foot contact points chosen are listed in Table 3. We also define contact points on the head and the tail but these are simply used to abort the model if the simulation falls over. Muscles are another area where simplification is necessary. It is actually very straightforward to simulate a large number of muscles and this causes very few problems, and relatively little simulation computational cost. However, each muscle needs to have its activation level controlled and therefore each additional muscle increases the dimensionality of the optimal control search space. This causes a huge additional cost in terms of search and it is therefore important to have as few functional muscles as possible. Since we also have the problem that we do not know the sizes of the individual muscles even if we can infer their probably identity using an extant phylogenetic bracket [57] it makes sense to reduce the model’s complexity by using a more idealised set of muscles that represent the functional actions that are likely to be available. These muscles can be defined with arbitrary paths and moment arms as long as they produce equivalent actions to anatomical muscles. The muscles chosen are listed in Table 4, including their origin and insertion points, and illustrated in Figure 4. Most muscles are not implemented as simple point-to-point muscles. This is because they need to wrap around bones to maintain their moment arms throughout the range of movement. This effect can be achieved using multiple via points but this approach often leads to unrealistic muscle paths at the extremes of joint action. It is also possible to define the muscle as a chain of linked segments and to calculate how these would slide over the bone morphology (and even other muscles). This is very computationally expensive and can cause numerical instability issues. Instead we define cylinders or pairs of parallel cylinders that allow a wrapping path to be calculated as needed with relatively minimal cost. The radius of the cylinder is chosen to match the effective moment arm of the muscle as it wraps around the condyles of the long bones.

Muscle Properties

As has been shown on several occasions [43], [56], [58], the most important property to estimate correctly in locomotor simulations is muscle mass. This is because the power available is proportional to muscle mass, and the force available, which is proportional to muscle area, is therefore proportional to the (muscle mass/muscle fibre length). Limb muscle mass as a fraction of total body mass is known for a number of animals and it is usually assumed that a value of 50% is an absolute maximum [58] and with values of 25 to 35% found more typically [59]. From the limited current data an approximate partitioning can be estimated with ∼60% of the muscle found around proximal joints, ∼30% around the intermediate, and ∼10% around the distal joints. Similarly muscle is split approximately ∼60% extensors to ∼40% flexors and ∼45% forelimb to ∼55% hindlimb [59]. Comparative data for greyhound, hare and reindeer are shown in Figure 5 and it can be seen that there is a relatively consistent pattern even for quadrupeds of different sizes and locomotor specialisations. Knowing these patterns it is therefore possible to calculate the masses of the individual muscles in the model based on their actions. This procedure works with any number of muscles as long as we assume that the mass is distributed evenly. Multiple joint muscles are simply divided among their multiple actions. To do this we need to use the model parameters listed in Table 5. Muscle density used is 1056 kg m⁻³ [60]. Force per unit area was chosen to be 300,000 Nm⁻² [61] but there are other values in the literature: Umberger et al [62] uses 250,000 Nm⁻², Alexander [63] reports an in vitro maximum value of 360,000 Nm⁻² for frog and 330,000 Nm⁻² for cat for parallel fibred leg muscles. Zheng et al. [64] recommend a value of 400,000 Nm⁻² for human quadriceps, and Pierrynowski [65] suggests 350,000 Nm⁻². There is a similarly large range for maximum contraction speed. Winter [47] suggest values from 6 to 10 times the muscle’s resting length per second for humans. This value is clearly highly dependent both on the fibre type composition of the muscle and on the temperature. Westneat [66] reports a range of values for fish from 3 to 10 s⁻¹ for different fibre types and Umberger et al [62] recommends values of 12 s⁻¹ for fast twitch and 4.8 s⁻¹ for slow twitch. A value of 8.4 s⁻¹ was chosen to represent a mixed fibred muscle. However it should be noted that there is data to suggest that this value reduces with body size [67] although there is very little data for large bodied animals and there is considerable scatter. The activation K value used is the recommended value for the muscle contraction and energetics model used [68].

Muscle maximum contractile force is determined by its physiological cross section area, which is calculated by dividing the muscle volume (obtained by dividing the mass by the muscle density) by the mean fibre length [47]. Unfortunately muscle fibre length is problematic to estimate. It is usually estimated by scaling from related species. This scaling can work well if there is a good modern analogue as is probably the case for early hominin musculoskeletal models [34], [35], but is considerably less reliable for morphologically more distinct species such as dinosaurs [43], [61]. This is particularly problematic if muscles with a similar action are being combined together to provide a more abstract joint driver since in that case there is no single muscle that can be used as a homologous reference. However there is a possible solution to this difficulty that can be derived from what we know about how vertebrate muscle contracts. Muscle can only generate force from approximately 60% of its resting length to about 160% [69]. Since the force follows an inverted U shaped curve we would expect most muscles to operate well within these limits in normal use, and since muscle physiology appears to be well conserved among the vertebrates, that this useful fraction of muscle length to be similar for different species. The length a muscle shortens depends on the change in angle at the joint multiplied by the moment arm [70]. So if we know the likely range of motion at a joint and the moment arm then we can predict the likely change in muscle length, and hence predict the muscle fibre length.

To test this prediction that vertebrate skeletal muscles exhibit a preferred length change, a literature survey was performed to identify suitable experimental data. What was required were studies that reported muscle fibre length and where length change could be calculated from moment arm and range of motion data. Since many muscle show changes in moment arm with joint angle this restricted studies to those where moment arm was measured over a range of joint angles. It was also decided that only studies that reported a reasonably large number of muscles should be included otherwise there would be bias associated with large numbers of studies on a relatively few specific muscles. There were relatively few suitable studies found, and of these several were of closely related primate species (hominoids including humans) and it was felt that including all these would produce a taxonomic bias. In the end the following species were chosen: chimpanzees [71], greyhound [72], [73], ostrich [74], [75] and horse [76]. For the chimpanzee, ostrich and horse the literature gave the best-fit polynomials for the tendon travel during joint rotation so that the length change of the muscle could be calculated directly. For the greyhound, the moment arm data was integrated over the range of angles presented to calculate length change. The chimpanzee and greyhound datasets included both fore- and hindlimbs whereas the ostrich and horse were hindlimb only. Ideally for this study the joint range of motion should match that seen in vivo for a range of movements. This is difficult to duplicate in cadaver studies since dead bodies tend to stiffen up which can restrict movement. Conversely as muscles are dissected away the joints become more mobile and this can lead to excessive movements at joints. In the case of the ostrich the joints were only moved through the range of movement associated with running and particularly for the hip and knee this was felt to be rather restricted. The analysis was repeated using a nominal, much larger range of movement for the ostrich data but this had no effect on the results and the conclusions remained unaltered so only the data as calculated directly from the paper is reported here.

Figure 6 shows the (extension/fibre length) ratios for the 121 muscles assessed subdivided by action and location. The modal value in the pooled case is 0.4–0.6, and only in two of the subdivided cases is the mode less clearly defined (0.2–0.6 in both cases). This suggests that assuming that muscle extends 50% of its resting fibre length (or conversely, that the resting fibre length is double the extension distance) is a reasonable assumption for most muscles. Very low values are probably due to one of two of factors. Firstly these are muscles whose prime action is neither flexion nor extension and therefore do not change length appreciably during this movement at the joint. Secondly these are muscles that cross more than one joint but whose action is mainly over a different joint. Very high values are more interesting because muscles cannot generate active force over these large extension ratios. Again there are two possibilities. Firstly these represent muscles that do not extend over the observed in vitro range in vivo. This includes two joint muscles where the full range of movement is not possible at both joints simultaneously. The human hamstrings are a good example of this where full hip flexion is not possible if the knee is extended. Secondly these represent muscles where part of the joint movement is accommodated by tendon stretch. The crural part of the camel m. plantaris is perhaps the most extreme example [77].

We can thus calculate the fibre length of the muscle by calculating the length change of the muscle which is equal to the joint range of motion multiplied by the moment arm. Moment arms are not necessarily easy to obtain for extinct species since exact points of attachment can be difficult to define. Furthermore, moment arms themselves depend on the presence of other soft tissue elements and exact instantaneous joint centres which are also unknown and need to be estimated (e.g. [42]). However if we use length change to define muscle fibre length, then the choice of moment arm does not actually matter in the simulation. If we choose a small moment arm, then we get a small length change, and hence a small fibre length. Since the volume of the muscle is defined by the mass which we have calculated a priori, a small fibre length leads to a large physiological cross section area which allows greater force production. Since all these relationships are directly proportional, the greater force production exactly compensates for the reduced moment arm in terms of the eventual torque around the joint. The contraction velocity is similarly exactly compensated: shorter muscle fibre, slower contraction velocity, but smaller moment arm leads to faster angular velocity around the joint. This is exactly as would be predicted from simple lever theory.

The key parameter then becomes joint range of motion. However there have been very few studies that have systematically looked at joint ranges of motion, and whilst some joint limits can be identified from skeletal features, others depend on soft tissue to limit the movement and thus can not. Ren et al. [17] compared elephant joint ranges of motion to cats, dogs, and humans and contrary to expectations did not find any body size related patterns. We thus created models with a range of different joint ranges of motion based on (1) estimation of joint range of motion from the skeleton; (2) range of motion matched to the functional range of motion for an elephant; (3) range of motion based on the previous two versions but with a restricted ankle range of motion. These ranges of motion are shown in Table 6. Using each of these ranges of motion allows us to calculate the length change of the individual muscle groups using the attachment points and wrapping cylinders previously specified. The tendon length is simply chosen so that the muscle tendon unit is slack when the joint is halfway between its maximum and minimum excursion. The calculated values for the muscles under the different range of motion conditions are shown in Table 7. Again there is no good comparative data on slack lengths and it is difficult to obtain since there is appreciable post mortem shrinkage and stiffening so that measurements taken from cadavers are probably not useful. Measuring passive elastic moments [78], as has been done for human models [79], might allow this to be calculated but the data would have to be taken from anaesthetised animals which would make it much more difficult to collect.

One useful side effect of calculating muscle fibre length from joint range of motion is that you can calculate the minimum muscle mass needed for joint extensors to be able to support a particular load. This is easiest to see for the ankle or wrist but is applicable for all the joints in each limb. If we consider Figure 7 which represents the ankle joint supporting the body weight of the animal (or some fraction thereof for multi-legged animals), we can see that the torque around the ankle (T) must be equal or greater to the ground reaction force (F) multiplied by the moment arm (M). This torque is generated by the ankle extensors, and using the methodology for specifying muscle fibre length outlines above we can show that:

(1)

Where K is the peak force generated per unit cross section area (N.m⁻²) as specified in Table 5; k is the (extension/fibre length) ratio chosen (0.5); m is the mass of the muscle (kg); Δθ is the joint range of motion (radians); and ρ is the muscle density (Kg.m⁻³).

Since T=FM we can rearrange this equation to calculate the minimum extensor mass:

(2)

Where B is the effective body mass (kg); and g is the acceleration due to gravity (m.s⁻²). Effective body mass is the body mass that would need to be supported by this leg alone. This would be equal to the body mass for a biped but would equal 1/3 of the body mass if we assume that 3 legs were on the ground at all times.

Of these values, only Δθ is unknown for a fossil animal and thus the muscle mass is directly proportional to the joint range of motion chosen. In fact the effect of joint range of motion may be greater than that because a larger range of motion may lead to a larger horizontal moment arm too. We performed this calculation for the Argentinosaurus model for all the joints using the maximum possible moment arm, as calculated by the maximum horizontal distance from the foot centre of pressure to the joint centre at either full extension or full flexion, as a way of checking that the model had adequate muscle to function.

Gait Simulation

Once all the muscle, joints, segments, and contacts have been defined it is necessary to find an appropriate activation pattern for the muscles that produces effective walking. To do this we use a feed-forward control system where a central pattern generator sends out muscle activation signals. This is a very simple approach but it is effective in a simulation environment which is entirely uniform. For these simulations we have adopted boxcar functions for the activation patterns [35]. A boxcar function is a rectangular function that has a zero value for a specified time and then a non-zero value for another specified time before falling back to zero. A boxcar function can thus be specified by 3 parameters: a delay, a width, and a height. This is a very concise way, in terms of control parameters, of specifying an activation pattern. If more precise control is required then two or more boxcar functions can be summed which rapidly allows very complex activation shapes to be generated, although single boxcar functions are the only ones that have been used in these simulations. The boxcar functions are duration normalised so that they work in a time interval from 0 to 1, and wrap around. The cycle time for all the functions is specified by a single master cycle time. The gait is assumed to be symmetrical so the left hand size drivers are identical to the right hand side drivers but are half a cycle out of phase. For these experiments the cycle phase was fixed externally. Since the model has 19 muscle groups per side, this equates to 57 unknown parameters to control the model.

We need to do two things: (1) find a good set of values for these parameters to allow high quality locomotion; (2) find a set of starting conditions that allow the simulation to work in a cyclic steady state. We do this using our now standard procedure of starting our simulant in its reference pose with all segments set at zero velocity, and using a genetic algorithm multiparameter optimisation procedure to find a pattern that maximises the forward distance moved by the model in fixed time. Once we have found a pattern that manages a good degree of forward movement, we use the segment poses and velocities from the middle of this simulation as a new set of starting conditions, and use the solution set as a best estimate solution set for a new optimisation run. This time the optimisation criteria is the maximum distance forward for a given amount of metabolic energy as calculated by the simulation. Once a good solution has been found, we repeat the process of selecting a mid-simulation set of velocities and poses, and reusing the solution set for a new optimisation run. In this way we bootstrap our start conditions, and eventually we end up converging on a largely steady state simulation that minimises the cost of locomotion since this is commonly considered the major goal of low speed locomotion [34], [35].

The simulation was performed using our in-house open source simulator, GaitSym. The software and the model specification files can be downloaded from www.animalsimulation.org. The simulation runs at about half real time on a modern processor, so a typical simulation run takes about 30 seconds of CPU core time. A single optimisation run requires 100,000 repeats of the simulation run, and typically 30 repeats of the bootstrap process are needed to get convergence. This equates to about 25,000 CPU core hours for each run condition tested. We had access to the HECToR, the UK National Supercomputer Service (www.hector.ac.uk) and were able to access up to 32,768 CPU cores at any one time. Our previous traditional genetic algorithm implementation [44] was very successful up to 512 cores but did not scale well for use with larger numbers of cores. Traditional genetic algorithms are highly synchronised [80], effectively because they use a seasonal breeding model. We re-implemented the algorithm using a continuous breeding and therefore asynchronous model and achieved excellent scaling up to 32,768 CPU cores (see Figure 8) which allowed us to explore considerably more options in terms of gait generation in a reasonable length of time.