## SUMMARY

The spring-mass model is often used to describe bouncing gaits. Although at first inspection the mechanical system appears simple, the solution to the motion cannot be derived easily. An analytical solution would provide a fast and intuitive method to determine the kinetic and kinematics of the centre of mass of terrestrial animals during over-ground steady state locomotion. Here, an analytical approximation using sine wave simplifications for the motion is presented. The analytical solution was almost indistinguishable from the numerical solution across initial leg angles of 17.5–30°; percentage differences between the analytical solution and the numerical solution were less than 1% for total mechanical energy, centre of mass position, total limb compression and centre of mass velocity and less than 2% different for resultant limb force and vertical acceleration of the centre of mass. The solution matched the relationship between stance time and speed collected from a trotting racehorse and accurately characterised previously published biological data. This study has shown that a simple analytical solution can predict the kinetics and kinematics of a spring-mass system over the range of biologically relevant sweep angles and horizontal velocities, and could be used to further understanding of limb deployment and gait selection. Using this analytical solution not only the force profile but also the changes in mechanical energy can be calculated from easily observed morphological and kinematic data.

## Introduction

There are a variety of gaits used by terrestrial animals to locomote. Different gaits involve fundamentally different mechanics: forces, energies, metabolic and mechanical power requirements are all gait-dependent. Terrestrial gaits have been broadly classified as walking gaits or running gaits. Walking gaits are often modelled as an inverted pendulum where the centre of mass (COM) vaults over a relatively rigid leg and thus the kinetic energy and potential energy of the COM are out of phase (Cavagna et al., 1977; Lee and Farley, 1998; Donelan et al., 2002). Running is profoundly different: usually the trajectory of the COM decreases in both vertical height and horizontal speed from foot-on to midstance, before being accelerated forwards and upwards so that at foot-off an aerial phase is produced. Therefore kinetic and potential energies are in phase. Here, we compare a simple, reductionist technique for describing the mechanics of running gaits with measured data and numerical solutions. In the future, this will provide theoretical and experimental approximations for limb forces, and consequences in terms of mechanical energies. Throughout, we consider the limb of a running animal to act like a linear compression spring. During stance, the decrease in vertical height of the COM is achieved through flexion of the joints, which stretch the spring-like musculo-tendinous structures that span the joints (Taylor, 1985; Farley and González, 1996; Seyfarth et al., 2002), storing elastic strain energy (Taylor, 1985; Cavagna et al., 1988; Farley and González, 1996; Kram and Dawson, 1998; Seyfarth et al., 2002).

The spring-mass model has been shown to be an accurate tool for modelling trotting, running and hopping over a range of species differing dramatically in terms of body size and shape, including cockroaches, quail, dogs, humans, horses and kangaroos (McMahon and Cheng, 1990; Alexander, 1992; Blickhan and Full, 1993; Bullimore and Burn, 2002). Although the simplicity of the spring-mass system suggests that the mathematics of its mechanics would be straightforward, this is not the case: the motion is complex and cannot be solved using a simple analytical solution. Horizontal and vertical COM deflections, velocities and accelerations interact to make analytical solutions very difficult to describe. Indeed, the situation is formally classed as a `non-integrative Hamilton equation', which indicates that an explicit analytical solution does not exist (Whittackler, 1904; Schwind and Koditschek, 2000).

A number of approximations to the spring-mass system have been made. Schwind and Koditschek conducted a thorough mathematical investigation of this system and produced complex estimates of COM position in a two degrees of freedom monopod. Simpler analyses included consideration of specific points in the stance phase (Farley et al., 1993; Ferris et al., 1998), small angle sweep assumptions (Geyer et al., 2005) and numerical iterations that search for a symmetrical solution in terms of motion of the stance phase (Blickhan, 1989) or mapping of successive maximum heights (Seyfarth et al., 2002). Numerical step-by-step solutions, although computer-intensive, do generate accurate answers. However, they do not provide an intuitive relationship between the five input variables. The advantage of an analytical solution is that it is fast and intuitive; it also allows the inter-relationships between variables to be immediately identified.

A relatively simple solution that can elegantly predict the COM trajectory and limb kinetics produced by the spring-mass model is desirable. Arguably an exact answer is not required, as in the real world limbs are obviously not mass-less, perfectly elastic springs of fixed stiffness – even in steady, level locomotion, muscles are required to provide power to replace energy dissipated because of hysteresis in tendons or other loaded structures. The leg is not necessarily the same effective length (i.e. the length between the hip and the toe) at foot-on to foot-off due to plantar flexion at the end of stance. Therefore, we are interested in determining whether a simple spring-mass model and related mathematical approximations are adequate in modelling running in real animals.

In steady state bouncing gaits such as running and trotting, the motion of the COM can be described in two components. In the vertical component, the downward vertical velocity peaks at foot-on and gradually reaches zero by mid-stance. In the second half of stance the velocity is reversed and peaks again at foot-off. The differential of this motion, acceleration, would approximate to a half sine wave. As acceleration is directly proportional to force, then vertical ground reaction force (GRF) can be approximated to a half sine wave (Cavagna et al., 1964, 1977; Alexander et al., 1979; Full et al., 1991; Farley et al., 1993; Witte et al., 2004). In the horizontal component, the COM first decelerates from foot-on to a minimum velocity at midstance before being accelerated forward in the second half of stance. It is therefore reasonable to consider the horizontal acceleration curve of the mass during the stance phase as a full negative sine wave, since this approximates the shape of the horizontal forces measured using force plates.

We hypothesize that an analytical approximation to the spring-mass model based on sine wave assumptions matches both numerical solutions of the model across the biological range and data collected from real animals.

In this paper, the accuracy of a number of analytical solutions based on sine wave approximations is assessed in three ways. Firstly, they are compared to a numerical solution produced by a computer simulation. Secondly, they are compared to data collected in this study from a Standardbred racehorse trotting over a range of speeds. Thirdly, a comparison is made with previously published data of a man running, dog trotting and a kangaroo hopping from McMahon and Cheng (1990) based on the data of Cavagna (1988). The validity of the sine wave assumptions across the biologically relevant parameter space are investigated.

## Materials and methods

### An analytical solution to the spring-mass model

A spring-mass model can simulate the movements of a multi-legged animal if
limbs are used simultaneously, as in trotting, pronking, hopping and pacing,
or separately as in bipedal running. Only steady state locomotion on a level
surface is considered; this means that the horizontal impulse over the stance
period is zero. The motion of the mass is assumed to be symmetrical about the
vertical: the horizontal kinetic energy and the potential energy of the COM
(i.e. the height) are equal at foot-on and foot-off, the vertical velocity is
equal but opposite at foot-on and foot-off, and the stance phase is
symmetrical in terms of leg orientation. Finally, the leg is treated as a
mass-less, perfectly elastic compression spring of unloaded length
*L*_{0} that occurs at foot-on and foot-off, and has stiffness
*k*. The mass *M*_{b} of the animal acting upon the limb
is modelled as a point mass at the top of the spring and is represented in
Fig. 1 as a filled circle. As
the motion of the COM during the aerial phase is ballistic then the duration
of the aerial phase is solely dependent on the vertical velocity at
foot-off.

At foot-on the limb makes contact with the ground at coordinate position
(0,0) at an initial angle of –θ_{0} to the vertical. At
this time the mass has a horizontal position (*x*) of
*L*_{0}sinθ_{0} and vertical position
(*y*) of *L*_{0}cosθ_{0}. The mass also
has initial velocities *V*_{xi} and *V*_{yi} in
the horizontal and vertical directions, respectively. The velocities change
over time (a positive value of *V*_{x} denotes a forward
velocity and a positive value for *V*_{y} denotes an upward
velocity) and the spring compresses and then extends as the mass pivots over
the `foot' until the spring loses contact with the ground at `foot-off'. As
the motion of the mass during flight is ballistic, and so determined by the
final horizontal and vertical velocities of the stance phase, then only the
stance phase needs to be considered to determine the motion of the mass. For
clarity the motion will be first considered in the vertical direction
(*y*) and then the horizontal direction (*x*).

#### The vertical orientation

During the gait cycle in steady state locomotion, the impulse due to
gravity on the centre of mass (i.e. the product of body mass *m*,
gravitational constant * g* and stride time) is balanced over the
stride by the vertical impulse generated by the leg during stance. If the
vertical GRF of animals is modelled as a half sine wave
(Alexander et al., 1979), the
vertical force experienced by the COM can be calculated. The sine wave that
would represent the vertical GRF is of the form

*F*

_{y}=

*F*

_{max}sin(

*at*), where

*F*

_{y}is the vertical force,

*F*

_{max}is the peak force,

*t*is time and 2π/

*a*is the period of the wave. The area under this curve, or the vertical impulse produced by the leg during contact time,

*T*

_{c}, is therefore: (1) It follows then that: (2) This can be simplified as

*T*

_{c}and

*T*

_{f}(flight time) are related by the duty factor, β: (3) So the vertical GRF in terms of time

*t*is: (4) As force is the product of mass and acceleration (from Newton's second law of motion), the vertical acceleration of our mass would be: (5) The integral of acceleration with respect to time is velocity. Hence the vertical velocity (

*V*

_{y}) of the mass is: (6) where

*C*

_{1}is integration constant. As

*V*

_{y}at time 0 s =

*V*

_{yi}(where

*V*

_{yi}is the initial vertical velocity and has a negative value), then: (7) and (8) The aerial phase can be derived from the vertical velocity at foot-off (which is equal in magnitude but opposite in sign to that at foot-on during symmetrical, steady state locomotion). Therefore

*V*

_{yi}can be written in terms of

*T*

_{f},

*T*

_{c}orβ (which are all inter-related; see Eqn 4) so

*V*

_{yi}is not a new term but has been written in the above as a separate variable for clarity.

The integral of velocity is position, so the vertical position is:
(9)
where *C*_{2} is an integration constant and is equal to the
vertical height at foot-on, which occurs at *y*(0):
(10)
and
(11)

#### The horizontal orientation

The slowing of the mass during the first half of stance in the horizontal direction means that an analysis based on a constant horizontal velocity will underestimate stance time and the vertical impulse generated. The horizontal acceleration of the mass can be shown to approximate to a full negative sine wave. The horizontal velocity of the animal decreases from foot-on until a minimum velocity is reached at midstance before increasing throughout the second half of the stance phase. As velocity is the integral of acceleration, this is consistent with a full negative sine wave for acceleration. Since force is directly proportional to acceleration, and the horizontal force trace measured by a force plate during steady state locomotion resembles a full negative sine wave with a total impulse of zero, this also supports the approximation of the horizontal acceleration of the COM to a full negative sine wave.

So, in our analytical approximation, the horizontal acceleration trace in
terms of time *A*_{x}(*t*) follows the curve:
(12)
where D and G are constants. It follows then that the horizontal velocity in
terms of time *V*_{x}(*t*) would be equal to the
integral of *a*_{x}(*t*), i.e.
(13)
where *C*_{3} is an integration constant. Therefore horizontal
position in terms of time *x*(*t*) will follow:
(14)
where *C*_{4} is an integration constant. If we consider a
horizontal position-time trace, as shown in
Fig. 2, then at time,
*t*=0, the *x* position is–
*L*_{0}sinθ_{0} and the gradient at this
point is the initial horizontal velocity *V*_{xi}. At
midstance, (*t*=*T*_{c}/2) the horizontal position is 0
and the gradient at this point is at its shallowest as here the minimum
horizontal velocity is obtained. At foot-off, *t=T*_{c} the
*x* position is *L*_{0}sinθ_{0} and the
gradient is again equal to *V*_{xi}.

In addition to this, the area under the horizontal velocity–time
trace during stance phase will be equal to the total horizontal distance
travelled:
(15)
As the horizontal velocity at *t*=0 is *V*_{xi}, by
substituting this into Eqn 13, then:
(16)
As *C*_{3} is the gradient of the line connecting the
horizontal position at *t*=0 and *t*=*T*_{c},
then *C*_{3} is equal to:
(17)
As 2π/G is the period of a full sine wave, then:
(18)
so D can be found by substituting G and *C*_{3} into Eqn 16:
(19)
(20)
The constant *C*_{4} can be calculated as, when *t=*0,
then *x* is equal to–
*L*_{0}sinθ_{0}, hence:
(21)
These constants can now be incorporated to find the relationship between
horizontal (*x*) acceleration, velocity, position and time:
(22)
(23)
(24)

Horizontal forces *F*_{x} can be calculated from horizontal
acceleration using Newton's second law of motion, *F=m***a**:
(25)
Therefore total resultant force (*F*), from the vector sum, is:
(26)
and total energy *E* is the sum of kinetic, potential and elastic
potential energy:
(27)
where *k* is limb stiffness and *c* is limb compression, which
can be calculated by:
(28)
Thus mechanical energy (*ME*) is the sum of kinetic and potential
energies:
(29)

### Production of numerical solutions to the planar spring-mass model

We created simulations using a computer software model (MSC
visualNastran4D, MSC.Software, California, USA run through Matlab *via*
Simulink, The MathWorks Inc, MA, USA). The simulation consisted of a mass
above a perfect linear spring tethered to the ground at coordinates (0,0). The
simulation starts with the same initial conditions as for the above
mathematics and stops when the tension in the spring returns to zero at the
end of stance phase. Adjustments were made to spring stiffness or horizontal
velocity until the simulation produced a symmetrical trajectory for the mass.
The resulting stance time was used as the input for the mathematical
calculations. All simulations were run at an integration rate of step size
0.0002 s. Halving the step size to 0.0001 s did not change the required spring
stiffness or horizontal velocity, expressed accurate to 4 decimal places.

### Comparison of the analytical and numerical solutions

Absolute and relative differences for all kinematic variables (horizontal and vertical position, velocity and acceleration) as well as resultant leg force and mechanical energy were compared between 17 numerical simulation solutions and the corresponding analytical solutions. For the computer simulations, and those obtained from the biological data, the same initial leg length, mass, leg stiffness and initial vertical velocity were used. Initial leg angle was varied from 2.5° to 60° and the required initial horizontal velocity was found so that a symmetric solution was produced. The inputs for the analytical solutions were the constant input parameters stated above and initial horizontal velocity and contact times that were obtained from the output of the computer simulations.

To compare the results of the analytical and numerical solutions across the biologically relevant limb angles (15–30°), each variable was plotted with respect to stance time, where zero represents midstance.

Percentage differences for each variable were calculated using
[1–(variable_{analytical}/variable_{numerical})]×100,
where the values for each variable were either the maximum values (cases
*V*_{y}, *A*_{x}, *A*_{y},
*F*) or minimum values (cases *y, V*_{x}). In the case
of *x* position and total mechanical energy, this value was calculated
for every time point and the maximum value was taken.

Vertical stiffness, *K*_{vert}, was calculated by:
(30)

### Comparison to biological data

To determine how accurately the analytical and numerical solutions reflect
biology, data were collected from one trotting Standardbred racehorse. Data
collection took place during an exercise session on a training track, with the
horse pulling an exercise sulky driven by an experienced driver. Data were
collected across a range of trotting speeds. Footfall data (foot-on and
foot-off events) were measured using ±50 g solid-state capacitive
accelerometers (ADXL150, Analog Devices, Norwood, MA, USA) as described in
Witte et al. (2004). Each
accelerometer was encased with epoxy-impregnated Kevlar fibres (total mass 2
g) and attached to the dorsal midline of the hoof of each forelimb using hot
glue (Bostik Findley Inc., Stafford, UK) so that the sensitive axis was in the
proximo-distal direction. Each accelerometer was attached to a telemetry
transmitter and battery *via* a short fatigue-resistant cable. The
transmitter and battery were securely attached over a protective pad over the
lateral aspect of the third metacarpal bone using a custom modified exercise
bandage. Output signals were telemetered using custom programmed FM radio
telemetry devices (ST/SR500, Wood and Douglas Ltd., Tadley, Hampshire, UK) and
logged at 1000 Hz *via* a 12-bit A/D converter and PCMCIA card (DAQcard
700, National Instruments, Austin, TX, USA) into a laptop using custom-made
software (Matlab, Natick, MA, USA). Speed was measured using a WAAS-enabled
GPS receiver (G30-L, Laipac Technology, Richmond Hill, Ontario, Canada)
attached to the sulky and sampled at a frequency of 1 Hz. This system is
accurate to within 0.2 m s^{–1} for 57% of all samples, as
described in Witte and Wilson
(2005).

Forelimb leg length was measured using a tape measure and recorded accurate to the nearest cm. Limb length measurements were taken when the horse was standing square. Forelimb length was taken as the vertical distance between the ground and the approximate insertion of serratus ventralis. When standing the limb is loaded to approximately 30% body weight. This loading introduces a small error that is not significant.

Foot-on angle for each stride was calculated by assuming symmetry of the
stance phase about the vertical in terms of leg length (leg length) and angle,
such that foot-on angle, θ_{0} is:
(31)
where *V*_{x,stride} is the horizontal velocity of the stride
and *T*_{c,stride} is the contact time for the stride.

Trot is a symmetrical `2 beat' gait where the diagonal limbs can be
considered to operate as a single virtual leg. A computer simulation was made
by using average input variables for the horse trotting at 7 m
s^{–1}. The input variables were *V*_{xi},
*V*_{yi}, θ_{0}, mass and leg length. The spring
stiffness was adjusted using a custom made software (MSCvisualNastran4D run
through Matlab *via* Simulink) until a symmetrical stance phase was
produced in which the final force in the leg spring was 0 N and the vertical
height of the COM was equal at foot-on and foot-off. The resulting calculated
spring stiffness represents the combined leg stiffness of the fore and
hindlimb.

#### Comparison to previously published data

The mathematical results were compared to previously published data of
Cavagna et al. (1988) published
by McMahon and Cheng (1990).
These consisted of a step cycle from a running man, a hopping kangaroo and a
trotting dog. Absolute and percentage differences between the numerical and
analytical solutions are reported for peak vertical force, peak *h*
(decrease in vertical height from the foot-on position) and *k*. The
analytical solutions were also superimposed over the original data of Cavagna
et al. Data from the analytical solutions were normalised using the same
method as McMahon and Cheng: vertical accelerations were normalised by
dividing by accelerations due to gravity (9.81 m s^{–2}) and
forces divided by body weight to obtain dimensionless horizontal force
(*F*_{x}/*m g*) and dimensionless vertical force
(

*F*

_{y}/

*m*). Dimensionless time was calculated as ω

**g**_{0}

*t*where ω

_{0}=(vertical leg stiffness/

*m*)

^{0.5}.

## Results

### Comparison of the numerical and analytical solutions

The required initial horizontal velocity (*V*_{xi}) and
resulting stance time (*T*_{c}) for each of the computer
simulations are shown in Fig.
3. In addition *K*_{vert} for the numerical and
analytical solutions are shown. *K*_{vert} increased markedly
at higher contact angles. From 15° to 30°, values increased from 131
kN m^{–1} to 812 kN m^{–1} for the numerical
solution and from 130 kN m^{–1} to 808 kN m^{–1}
for the analytical solution.

Fig. 4 shows the results for
the numerical (blue line) and analytical solutions (red line) (the same colour
coding will be used throughout) for θ_{0} of 15°, 20°,
25°and 30° (i.e. the angles typically used by animals) to show how
*x, y, V*_{x}, *V*_{y}, *A*_{x},
*A*_{y}, *F* and *ME* change throughout the
stance phase for each condition. Here the stance phase has been adjusted so
that midstance occurs when time is zero. The analytical and numerical
solutions are indistinguishable for position data, vertical velocity, vertical
acceleration and resultant force as the red line superimposes on the blue
line. Slight differences are seen for total mechanical energy, horizontal
velocity and horizontal acceleration between the numerical and analytical
solutions, especially at the higher contact angles. As the analytical solution
considers the vertical and horizontal components separately, total energy is
not required to be constant, unlike the numerical solution that is based on
Newtonian mechanics.

To compare the analytical and numerical solutions, histograms were
constructed for all simulations to show the percentage difference in
positions, velocities, accelerations and resultant leg force
(Fig. 5). The analytical
solutions are most similar to the numerical solutions for initial leg angles
between 15° and 30° (shown by grey bars), which correspond to limb
angles commonly used by animals (the solutions outside that range are shown by
black bars). Within this range the maximum difference between the analytical
and numerical solutions in *x* position was 0.17%, *y* position
0.01%, total leg compression 0.34%, *V*_{x} 0.40%,
*V*_{y} 0.56%, *A*_{x} 30.75%,
*A*_{y} 1.95%, *F* 1.18% and *ME* 0.03%.

### Comparison to experimental data from a trotting horse over a range of speeds

The mass of the horse was 427 kg. The forelimb leg length was 1.43 m. At 7
m s^{–1} the average *T*_{c} for the forelimb was
158 ms and the average *T*_{f} was 139 ms (duty factor=0.36).
This flight time would require a final vertical velocity of 0.68 m
s^{–1}. A spring stiffness of 80.9 kN m^{–1} was
required for the numerical simulation to produce a symmetrical stance
phase.

Fig. 6 show that the simulations very closely matched the stance time and initial contact angle data obtained from the racehorse over the entire speed range.

### Comparison to previously published data

The analytical and numerical approximations produced results that were very
similar to previous published results for man.
Fig. 7 show the solutions of
McMahon and Cheng's method and the analytical solution for normalised vertical
acceleration, dimensionless horizontal force and dimensionless vertical force.
The two solutions are almost identical for normalised vertical acceleration
against vertical displacement and very similar for dimensionless force against
dimensionless time. Peak forces differed by less than 2% between the
analytical and numerical solutions to the published data (analytical
solution=1986 N, numerical solution=2022 N, published result=2020 N). The
vertical drop in height *h* was nearly identical for all methods, 6.25
cm for the analytical and numerical solutions *vs* 6.24 cm for the
published data. Leg stiffness was also very similar (less than 2% difference);
*k* for the analytical solution was 10740 N m^{–1}, 10930
N m^{–1} for the numerical solution and 10950 N
m^{–1} for McMahon and Cheng's simulation. Peak horizontal force
only differed by 26 N (6.7% difference) and 0.5 N (0.5% difference) for the
analytical and numerical solutions respectively, and was greater for the
analytical solution (420 N) compared to the published data (390 N).

Fig. 8 compared the normalised vertical acceleration against displacement traces for the trotting dog and the hopping kangaroo. There are small differences between the methods but the results are very similar.

Unfortunately force data for the trotting dog, and force and leg stiffness
data for the hopping kangaroo, could not be extracted due to the data
presented from McMahon and Cheng
(1990); however *h*
values could be measured from graphs. Only an *h* and a leg stiffness
value were extracted for the trotting dog. There was a difference of less than
3 mm in *h* between the two solutions and the published result, as
shown in Table 1. The
analytical solution produced a decrease in height of 1.2 cm whereas the other
study and the numerical solution had *h*=1.3 cm (percentage difference
of 7.7%). Smaller percentage differences were obtained for total leg
stiffness; 1420 N m^{–1}, 1490 N m^{–1} and 1490 N
m^{–1} for the analytical, numerical and McMahon and Cheng's
solution, respectively.

Only an *h* value was obtained for comparison with the analytical
and numerical solutions for the kangaroo. Values of 7.7 cm, 6.4 cm and 6.9 cm
(difference from the published data of 1.3 cm and 0.8 cm; 13% and 10.5%
difference, respectively) were obtained for the published data, the analytical
solution and the numerical solution, respectively.

## Discussion

The aims of this discussion are to explore the accuracy and implications of using sine wave approximations for a spring leg in terms of forces, energies and how these could contribute to understanding gait selection criteria.

### Sine wave approximations produce results similar to numerical simulations and biological data

This study has shown that a simple analytical solution, requiring a small
number of easily measurable input values, can predict the kinetics and
kinematics of a numerical solution to the planar spring-mass system over the
range of biologically relevant sweep angles and horizontal velocities. For
initial leg angles of between 17.5° and 30° the percentage differences
between the analytical solution and the numerical solution for the majority of
variables were less than 1% (total *ME*, COM position, total limb
compression and COM velocity) and a difference of less than 2% was obtained
for total limb force and vertical acceleration of the COM. Less accurate
results were obtained for horizontal acceleration, with up to 31% difference
between the analytical and numerical solutions.

The analytical solution was most accurate for leg contact angles of between 17.5° and 30°, which represent the range of angles used by running animals including the racehorse in this study (Fig. 6B). Running animals typically have an initial leg angle of approximately 17.5–20° when they make the transition from walking to running (Blickhan, 1989; Lee and Farley, 1998), and a typical contact angle of 30–35° at their maximum running speed (Farley et al., 1993; Alexander and Jayes, 1983; Gatesy and Biewener, 1991).

The analytical and numerical solutions matched the data obtained from our
trotting racehorse across the whole horizontal speed range. The racehorse in
this study used contact angles ranging from 15° at 3 m
s^{–1} to 32° at 11 m s^{–1}. Stance time and
duty factor data matched the predicted values. The value of combined limb
stiffness (i.e. the combined stiffness of the forelimb and hindlimb) required
for the numerical calculations was 80.9 kN m^{–1}. McGuigan and
Wilson found an overall forelimb leg stiffness of 55 kN m^{–1}
in horses of similar size to the one in this study
(McGuigan and Wilson, 2003).
If 57% of the mass of the horse were carried by the forelimb in a trot stance
phase (Witte et al., 2004)
then the total forelimb and hindlimb stiffness would be 96 kN
m^{–1}, a value approximately 20% higher than the value used in
this study. However using the isometric scaling relationship of
*k*_{leg}∝*M*_{b}^{0.67}, where
*M*_{b} is body mass
(Farley et al., 1993), which
was based on number of species ranging from kangaroo rat (0.112 kg) to horse
(135 kg), the combined leg stiffness for the racehorse in his study would be
57.9 kN m^{–1}. This is approximately 30% less than the value
than estimated in this study; however, the horse used in this inter-species
scaling study was considerably smaller (135 kg *vs* 426 kg), had a
shorter leg length (0.75 m *vs* 1.43 m) and trotted at a slower speed
(approximately 3 m s^{–1} against 7 m s^{–1}). Our
estimate of limb stiffness is therefore within the range of previously
published values.

Both analytical and numerical solutions predicted the changes in biological
vertical stiffness across the trotting speed range. A change in leg angle from
25° to 30° results in an increase in *K*_{vert} of
201%. This is of a similar magnitude to the 170% increase for a smaller
trotting horse reported by Farley et al.
(1993). Therefore the
analytical approximation not only matches the numerical solutions, but can
also represent the mechanics of a complex biological system – the
trotting horse, and can be expected to be appropriate for other symmetrical
bouncing gaits.

The solutions were also similar to other previously published models and
kinetic and kinematic measurements. Peak limb forces predicted by Alexander et
al. (1979) were very similar to
the peak forces produced by the analytical and numerical solutions. Our
solutions have similar accuracy to the analytical solution reported by Geyer
et al. (2005), who report a
less than 1% error for total leg compression and 0.6° difference for
angular motion compared to a numerical solution (based on a simulation with
input consistent with a running human). Our analytical solution produced
errors of less than 1% for total leg compression and a maximum difference in
leg angle less than 0.1° across the entire biological range; however, our
solution uses different assumptions and cannot be considered an equivalent
approach. Comparison of our analytical and numerical solutions with the
previously published results of McMahon and Cheng based on the biological data
of Cavagna et al. showed that the solutions were nearly indistinguishable for
the running man: maximum force and limb stiffness differed by less than 2% and
vertical height decrease (*h*) differed by 0.2%. The results for the
dog were similar: a difference of less than 5% in leg stiffness and a
difference of 0.1 cm for vertical compression (as measured from their figure)
were detected between the two methods and the published result. The kangaroo
data provided the poorest fit both for McMahon and Cheng's model and for the
analytical and numerical solution, and possible explanations are discussed
below.

While previous methods remain valuable, the advantages of our analytical solution are twofold. Like, and indeed following, Alexander's sine wave approximation, which just considers vertical GRF, it is very simple (Alexander et al., 1979), but our addition of the horizontal component allows for calculation of fluctuations in velocity and mechanical energies.

### Limitations

Both the numerical and sine wave analyses are limited to gaits that are similar to bouncing springs. In bipeds these are hopping and running (not walking and probably not skipping). In quadrupeds these are trot, pace and pronk (not walk and possibly not gallop). Both analyses are more accurate if the limbs behave as simple compression springs of constant stiffness. In the numerical solutions the leg stiffness (the gradient of the peak force against peak compression curve shown in Fig. 9) remains constant for all simulations. For the analytical solutions stiffness remains constant over the most of the solutions and the range of angles used by animals and only increases between the two most extreme initial conditions of 45° and 60° due to inaccuracies in peak force estimation at these extreme contact angles (Fig. 9).

Hopping macropodids, whilst clearly bouncing, deviate from our assumptions as the stance phase is decidedly asymmetric: leg angle and leg length are much greater at foot-on than at foot-off (McGowan et al., 2005), and the peak vertical GRF is significantly greater (36%) than that predicted using a sine wave (Kram and Dawson, 1998). In addition, kangaroos use large contact angles. A leg angle of 53° (Farley et al., 1993) has been measured. The angle of the virtual leg (i.e. foot to COM) would be significantly less, since the COM is cranial to the hip (McGowan et al., 2005). The required coefficient of friction at the foot is the tangent of the contact angle of the virtual leg. For instance a 45° contact angle would require a coefficient of friction of 1; improbably high for biological foot–ground interaction. Macropod locomotion therefore requires further consideration and the approach used here is not ideal for these species. Most animals bounce with less extreme leg angles (McMahon and Cheng, 1990; Alexander, 1992; Blickhan and Full, 1993; Farley et al., 1993); those studied are well represented by compression springs and would be modelled well by our analytical approach.

The analytical solution simplifies the motion of the spring mass model by considering the horizontal and vertical components of motion separately, thus the force vector does not have to remain in line with the virtual leg. However, the differences between the leg angles calculated using the numerical solution and those calculated using the analytical solutions were very similar across the biological range (Fig. 10A), with a maximum difference of less than 0.10°. A maximum difference in leg angle of less than 0.24° occurred between the two solutions for all simulations across the entire range of angles (Fig. 10B).

### Applications of the technique

Our findings extend Alexander's simplification to obtain vertical and horizontal forces and changes in mechanical energy from body mass, leg length, foot timings and horizontal velocity. This allows calculation of leg force and hence muscle and tendon forces and limb elastic energy storage, and mechanical work on the COM without the need for force plates and high speed kinematics. These measurements will allow novel insights into the compromises involved in locomotor strategy.

For instance, although the velocity amplitude through stance is greater in
the vertical than the horizontal components
(Fig. 4), the kinetic energy
fluctuations are much greater in the horizontal direction. This is because a
small change in velocity around a large mean velocity involves a large change
in kinetic energy as it is the difference in velocity squared. Therefore
horizontal components dominate changes in COM mechanical energy during
running. For an animal of a specific leg length and mass using a simple
bouncing gait the leg stiffness determines the relationship between speed and
initial leg angle. High stiffness legs use small contact angles. Conversely,
more compliant limbs require larger contact angles. High stiffness limbs are
beneficial since they result in small changes in mechanical energy and hence
low hysteresis losses, since biological springs do not return all energy
stored in them (the potential energy storage in the spring leg can be
calculated as 0.5*kc*^{2}, and can also be calculated using our
method). Stiff legs, however, result in brief contact times and high forces.
These high forces require a stronger leg (which may be difficult to protract
quickly) or a reduced safety factor and hence an increased risk of injury.
There may also be a detrimental effect on locomotor and muscle efficiency due
to shorter stance times (Kram and Taylor,
1990). If maximum contact angle determines maximum running speed
then stiff legs will also enable a higher maximum speed. This trade-off
explains why, in systems where force is a minor issue, high-stiffness,
vertically orientated limbs are preferred, whereas in biological systems,
where there is a structural and energetic cost to `force generation', the
compromise solution tends towards a more compliant limb. Therefore our simple
analytical approximation allows us to explore intuitively the consequences of
observed and postulated strategies of locomotion. Compliant limbs may, of
course, also carry benefits in control of locomotion under variable
conditions.

### Future directions

Apart from using these methods to determine forces and energies from the measurements as outlined above, we plan to develop these approximations further to provide insight into asymmetrical gaits, specifically skipping and galloping. We foresee that these approaches will provide alternative but convergent conclusions to the developing collision-based models (Ruina et al., 2005). Although understanding of galloping is improving quickly, the consequences and desirability of this gait are certainly not yet fully understood. Simple optimisations providing conclusions about gait parameters including footfall sequences will have to take account of both force and energetic costs. As discussed above, energy loss minimisation and peak force minimisation have quite different requirements. While collision-based models may be effective in determining energetic consequences of different gaits, spring-based models will be required if force consequences are to be understood.

### Conclusion

This study has shown that a simple analytical solution matches the kinematics and kinetics of the motion of a spring-mass system, accurately modelling the kinetics of running and trotting. Limb force, leg stiffness and changes in mechanical energy can be determined for these gaits from a few easily obtainable, morphological and kinematic observations. This will allow key energetic consequences of observed locomotion to be determined in the field and without the need for force plates. In addition this allows consideration of the consequences of postulated morphologies and gait strategies.

The authors wish to thank Henry Wilson and Thilo Pfau for their technical help, Jim Usherwood and Anna Wilson for help with manuscript preparation. J.R. is a BBSRC funded student and A.W. a BBSRC Research Fellow and holder of a Royal Society Wolfson Research Merit award.

## List of symbols

*A*_{x}- horizontal acceleration (m s
^{–2}) *A*_{y}- vertical acceleration (m s
^{–2}) - c
- total limb compression (m)
- COM
- centre of mass
- E
- total energy (J)
- F
- total resultant limb force (N)
*F*_{max}- peak force=peak vertical force
*F*_{x}- horizontal force
*F*_{y}- vertical force
- g
- gravitational constant
- GRF
- ground reaction force
- h
- decrease in vertical height of the COM, relative to foot-on (m)
- k
- leg stiffness (N m
^{–1}) *K*_{vert}- vertical stiffness (kN m
^{–1}) *L*_{0}- initial leg length (m)
- L
- leg length (m)
*L*_{m}- leg length at midstance (m)
- m
- mass (kg)
*M*_{b}- body mass
- ME
- mechanical energy (J)
- t
- time (s)
*T*_{c}- contact time (s)
*T*_{f}- flight time (s)
*V*_{x}- horizontal velocity (m s
^{–1}) *V*_{y}- vertical velocity (m s
^{–1}) - x
- horizontal position (m)
- y
- vertical position (m)
- θ0
- initial contact angle, relative to the vertical
- ω0
- (vertical leg stiffness/
*m*)^{0.5}

- © The Company of Biologists Limited 2005