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.

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.

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 L0 that occurs at foot-on and foot-off, and has stiffness k. The mass Mb 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 L0sinθ0 and vertical position(y) of L0cosθ0. The mass also has initial velocities Vxi and Vyi in the horizontal and vertical directions, respectively. The velocities change over time (a positive value of Vx denotes a forward velocity and a positive value for Vy 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 Fy=Fmaxsin(at), where Fy is the vertical force, Fmax 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, Tc, is therefore:
\[\ {{\int}_{0}^{T_{\mathrm{c}}}}F_{\mathrm{max}}{\,}\mathrm{sin}\left(\frac{{\pi}}{T_{\mathrm{c}}}t\right)dt=\frac{2F_{\mathrm{max}}T_{\mathrm{c}}}{{\pi}}.\]
(1)
It follows then that:
\[\ \frac{2F_{\mathrm{max}}T_{\mathrm{c}}}{{\pi}}=m\mathbf{\mathit{g}}(T_{\mathrm{c}}+T_{\mathrm{f}}).\]
(2)
This can be simplified as Tc and Tf(flight time) are related by the duty factor, β:
\[\ {\beta}=\frac{T_{\mathrm{c}}}{T_{\mathrm{c}}+T_{\mathrm{f}}}.\]
(3)
So the vertical GRF in terms of time t is:
\[\ F_{\mathrm{y}}=\frac{{\pi}m\mathbf{\mathit{g}}}{2{\beta}}\mathrm{sin}\left(\frac{{\pi}}{T_{\mathrm{c}}}t\right).\]
(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:
\[\ A_{\mathrm{y}}(t)=\frac{{\pi}\mathbf{\mathit{g}}}{2{\beta}}\mathrm{sin}\left(\frac{{\pi}}{T_{\mathrm{c}}}t\right)-\mathbf{\mathit{g}}.\]
(5)
The integral of acceleration with respect to time is velocity. Hence the vertical velocity (Vy) of the mass is:
\[\ V_{\mathrm{y}}(t)=\frac{-T_{\mathrm{c}}\mathbf{\mathit{g}}}{2{\beta}}\mathrm{cos}\left(\frac{{\pi}}{T_{\mathrm{c}}}t\right)-\mathbf{\mathit{g}}t+C_{1},\]
(6)
where C1 is integration constant. As Vy at time 0 s = Vyi (where Vyi is the initial vertical velocity and has a negative value), then:
\[\ C_{1}=V_{\mathrm{yi}}+\frac{T_{\mathrm{c}}\mathbf{\mathit{g}}}{2{\beta}},\]
(7)
and
\[\ V_{\mathrm{y}}(t)=\frac{-T_{\mathrm{c}}\mathbf{\mathit{g}}}{2{\beta}}\mathrm{cos}\left(\frac{{\pi}}{T_{\mathrm{c}}}t\right)-\mathbf{\mathit{g}}t+\left(V_{\mathrm{yi}}+\frac{T_{\mathrm{c}}\mathbf{\mathit{g}}}{2{\beta}}\right).\]
(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 Vyi can be written in terms of Tf, Tc orβ (which are all inter-related; see Eqn 4) so Vyi 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:
\[\ y(t)=\frac{-\mathbf{\mathit{g}}T_{\mathrm{c}}^{2}}{2{\pi}{\beta}}\mathrm{sin}\left(\frac{{\pi}}{T_{\mathrm{c}}}t\right)-\frac{\mathbf{\mathit{g}}t^{2}}{2}+\left(V_{\mathrm{yi}}+\frac{T_{\mathrm{c}}\mathbf{\mathit{g}}}{2{\beta}}\right)t+C_{2},\]
(9)
where C2 is an integration constant and is equal to the vertical height at foot-on, which occurs at y(0):
\[\ C_{2}=L_{0}{\,}\mathrm{cos}{\,}{\theta}_{0},\]
(10)
and
\[\ y(t)=\frac{-\mathbf{\mathit{g}}T_{\mathrm{c}}^{2}}{2{\pi}{\beta}}\mathrm{sin}\left(\frac{{\pi}}{T_{\mathrm{c}}}t\right)-\frac{\mathbf{\mathit{g}}t^{2}}{2}+\left(V_{\mathrm{yi}}+\frac{T_{\mathrm{c}}\mathbf{\mathit{g}}}{2{\beta}}\right)t+L_{0}{\,}\mathrm{cos}{\,}{\theta}_{0}.\]
(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 Ax(t) follows the curve:
\[\ A_{\mathrm{x}}(t)=-\mathrm{D}{\,}\mathrm{sin}(\mathrm{G}t),\]
(12)
where D and G are constants. It follows then that the horizontal velocity in terms of time Vx(t) would be equal to the integral of ax(t), i.e.
\[\ V_{\mathrm{x}}(t)=\frac{\mathrm{D}}{\mathrm{G}}\mathrm{cos}(\mathrm{G}t)+C_{3},\]
(13)
where C3 is an integration constant. Therefore horizontal position in terms of time x(t) will follow:
\[\ x(t)=\frac{\mathrm{D}}{\mathrm{G}^{2}}\mathrm{sin}(\mathrm{G}t)+C_{3}t+C_{4},\]
(14)
where C4 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– L0sinθ0 and the gradient at this point is the initial horizontal velocity Vxi. At midstance, (t=Tc/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=Tc the x position is L0sinθ0 and the gradient is again equal to Vxi.
In addition to this, the area under the horizontal velocity–time trace during stance phase will be equal to the total horizontal distance travelled:
\[\ {{\int}_{0}^{T_{\mathrm{c}}}}V_{\mathrm{x}}(t)dt=2L_{0}{\,}\mathrm{sin}{\,}{\theta}_{0}.\]
(15)
As the horizontal velocity at t=0 is Vxi, by substituting this into Eqn 13, then:
\[\ V_{\mathrm{xi}}=\frac{\mathrm{D}}{\mathrm{G}}+C_{3}.\]
(16)
As C3 is the gradient of the line connecting the horizontal position at t=0 and t=Tc,then C3 is equal to:
\[\ C_{3}=\frac{2L_{0}{\,}\mathrm{sin}{\,}{\theta}_{0}}{T_{\mathrm{c}}}.\]
(17)
As 2π/G is the period of a full sine wave, then:
\[\ \mathrm{G}=\frac{2{\pi}}{T_{\mathrm{c}}},\]
(18)
so D can be found by substituting G and C3 into Eqn 16:
\[\ \mathrm{D}=\mathrm{G}(V_{\mathrm{xi}}-C_{3}),\]
(19)
\[\ \mathrm{D}=\frac{2{\pi}}{T_{\mathrm{c}}}\left[V_{\mathrm{xi}}-\left(\frac{2L_{0}{\,}\mathrm{sin}{\,}{\theta}_{0}}{T_{\mathrm{c}}}\right)\right].\]
(20)
The constant C4 can be calculated as, when t=0,then x is equal to– L0sinθ0, hence:
\[\ C_{4}=-L_{0}{\,}\mathrm{sin}{\,}{\theta}_{0}.\]
(21)
These constants can now be incorporated to find the relationship between horizontal (x) acceleration, velocity, position and time:
\[\ A_{\mathrm{x}}(t)=-\left\{\frac{2{\pi}}{T_{\mathrm{c}}}\left[V_{\mathrm{xi}}-\left(\frac{2L_{0}{\,}\mathrm{sin}{\,}{\theta}_{0}}{T_{\mathrm{c}}}\right)\right]\right\}\mathrm{sin}\left(\frac{2{\pi}}{T_{\mathrm{c}}}t\right),\]
(22)
\[\ V_{\mathrm{x}}(t)=\left[V_{\mathrm{xi}}-\left(\frac{2L_{0}{\,}\mathrm{sin}{\,}{\theta}_{0}}{T_{\mathrm{c}}}\right)\right]\mathrm{cos}\left(\frac{2{\pi}}{T_{\mathrm{c}}}t\right)+\left(\frac{2L_{0}{\,}\mathrm{sin}{\,}{\theta}_{0}}{T_{\mathrm{c}}}\right),\]
(23)
\[\ x(t)=\left(\frac{V_{\mathrm{xi}}T_{\mathrm{c}}-2L_{0}{\,}\mathrm{sin}{\,}{\theta}_{0}}{2{\pi}}\right)\mathrm{sin}\left(\frac{2{\pi}}{T_{\mathrm{c}}}t\right)+\left(\frac{2L_{0}{\,}\mathrm{sin}{\,}{\theta}_{0}}{T_{\mathrm{c}}}\right)t-L_{0}{\,}\mathrm{sin}{\,}{\theta}_{0}.\]
(24)
Horizontal forces Fx can be calculated from horizontal acceleration using Newton's second law of motion, F=ma:
\[\ F_{\mathrm{x}}(t)=mA_{\mathrm{x}}(t).\]
(25)
Therefore total resultant force (F), from the vector sum, is:
\[\ F(t)=\sqrt{F_{\mathrm{x}}(t)^{2}+F_{\mathrm{y}}(t)^{2}},\]
(26)
and total energy E is the sum of kinetic, potential and elastic potential energy:
\[\ E(t)=[\frac{1}{2}mV_{\mathrm{x}}(t)^{2}]+[\frac{1}{2}mV_{\mathrm{y}}(t)^{2}]+m\mathbf{\mathit{g}}y(t)+\frac{1}{2}kc^{2},\]
(27)
where k is limb stiffness and c is limb compression, which can be calculated by:
\[\ c(t)=L_{0}-\sqrt{x(t)^{2}+y(t)^{2}}.\]
(28)
Thus mechanical energy (ME) is the sum of kinetic and potential energies:
\[\ ME(t)=[\frac{1}{2}mV_{\mathrm{x}}(t)^{2}]+[\frac{1}{2}mV_{\mathrm{y}}(t)^{2}]+m\mathbf{\mathit{g}}y(t).\]
(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 viaSimulink, 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–(variableanalytical/variablenumerical)]×100,where the values for each variable were either the maximum values (cases Vy, Ax, Ay, F) or minimum values (cases y, Vx). 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, Kvert, was calculated by:
\[\ K_{\mathrm{vert}}=\frac{F_{\mathrm{y}},\mathrm{peak}}{h_{\mathrm{peak}}}.\]
(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:
\[\ {\theta}_{0}=\mathrm{sin}^{-1}\left(\frac{V_{\mathrm{x}},\mathrm{stride}{\times}T_{\mathrm{c},\mathrm{stride}}}{2{\times}\mathrm{leg}{\,}\mathrm{length}}\right),\]
(31)
where Vx,stride is the horizontal velocity of the stride and Tc,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 Vxi, Vyi, θ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(Fx/mg) and dimensionless vertical force(Fy/mg). Dimensionless time was calculated as ω0t where ω0=(vertical leg stiffness/m)0.5.

Comparison of the numerical and analytical solutions

The required initial horizontal velocity (Vxi) and resulting stance time (Tc) for each of the computer simulations are shown in Fig. 3. In addition Kvert for the numerical and analytical solutions are shown. Kvert 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–1for 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, Vx, Vy, Ax, Ay, 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%, Vx 0.40%, Vy 0.56%, Ax 30.75%, Ay 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 Tc for the forelimb was 158 ms and the average Tf 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 hvalues 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.

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–1in 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 klegMb0.67, where Mb 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 Kvert 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.5kc2, 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.

     
  • Ax

    horizontal acceleration (m s–2)

  •  
  • Ay

    vertical acceleration (m s–2)

  •  
  • c

    total limb compression (m)

  •  
  • COM

    centre of mass

  •  
  • E

    total energy (J)

  •  
  • F

    total resultant limb force (N)

  •  
  • Fmax

    peak force=peak vertical force

  •  
  • Fx

    horizontal force

  •  
  • Fy

    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)

  •  
  • Kvert

    vertical stiffness (kN m–1)

  •  
  • L0

    initial leg length (m)

  •  
  • L

    leg length (m)

  •  
  • Lm

    leg length at midstance (m)

  •  
  • m

    mass (kg)

  •  
  • Mb

    body mass

  •  
  • ME

    mechanical energy (J)

  •  
  • t

    time (s)

  •  
  • Tc

    contact time (s)

  •  
  • Tf

    flight time (s)

  •  
  • Vx

    horizontal velocity (m s–1)

  •  
  • Vy

    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

Alexander, R. M. (
1992
). A model of bipedal locomotion on compliant legs.
Phil. Trans. R. Soc. Lond. B
338
,
189
-198.
Alexander, R. M. and Jayes, A. S. (
1983
). A dynamic similarity hypothesis for the gaits of quadrupedal mammals.
J. Zool.
201
,
135
-152.
Alexander, R. M., Maloiy, G. M. O., Hunter, B., Jayes, A. S. and Nturibi, J. (
1979
). Mechanical stresses during fast locomotion of buffalo (Syncerus caffer) and elephant (Loxodonta africana).
J. Zool. Lond.
189
,
135
-144.
Blickhan, R. (
1989
). The spring-mass model for running and hopping.
J. Biomech.
22
,
1217
-1227.
Blickhan, R. and Full, R. J. (
1993
). Similarity in multilegged locomotion: Bouncing like a monopode.
J. Comp. Physiol. A
173
,
509
-517.
Bullimore, S. R. and Burn, J. F. (
2002
). An approximate analytical solution for the planar spring mass model of locomotion.
Proc. 4th World Cong. Biomech.
2002
,
216
.
Cavagna, G. A., Saibene, F. P. and Margaria, R.(
1964
). Mechanical work in running.
J. Appl. Physiol.
19
,
249
-256.
Cavagna, G. A., Heglund, N. C. and Taylor, C. R.(
1977
). Mechanical work in terrestrial locomotion: two basic mechanisms for minimizing energy expenditure.
Am. J. Physiol.
233
,
R243
-R261.
Cavagna, G. A., Franzetti, P., Heglund, N. C. and Willems,P. (
1988
). The determinants of the step frequency in running,trotting and hopping in man and other vertebrates.
J. Physiol.
399
,
81
-92.
Donelan, J. M., Kram, R. and Kuo, A. D. (
2002
). Simultaneous positive and negative external mechanical work in human walking.
J. Biomech.
35
,
117
-124.
Farley, C. T. and González, O. (
1996
). Leg stiffness and stride frequency in human running.
J. Biomech.
29
,
181
-186.
Farley, C. T., Glasheen, J. and McMahon, T. A.(
1993
). Running springs: speed and animal size.
J. Exp. Biol.
185
,
71
-86.
Ferris, D. P., Louie, M. and Farley, C. T.(
1998
). Running in the real world: adjusting leg stiffness for different surfaces.
Proc. R. Soc. Lond. B
265
,
989
-994.
Full, R. J., Blickhan, R. and Ting, L. H.(
1991
). Leg design in hexapedal runners.
J. Exp. Biol.
158
,
369
-390.
Gatesy, S. M. and Biewener, A. A. (
1991
). Bipedal locomotion: effects of speed, size and limb posture in birds and humans.
J. Zool.
224
,
127
-147.
Geyer, H., Seyfarth, A. and Blickhan, R.(
2005
). Spring-mass running: simple approximate solution and application to gait stability.
J. Theor. Biol.
232
,
315
-328.
Hutchinson, J. R., Famini, D., Lair, R. and Kram, R.(
2003
). Are fast-moving elephants really running?
Nature
422
,
493
-494.
Kram, R. and Dawson, T. J. (
1998
). Energetics and biomechanics of locomotion by red kangaroos (Macropus rufus).
Comp. Biochem. Physiol.
120B
,
41
-49.
Kram, R. and Taylor, C. R. (
1990
). Energetics of running: a new perspective.
Nature
346
,
265
-267.
Lee, C. R. and Farley, C. T. (
1998
). Determinants of the center of mass trajectory in human walking and running.
J. Exp. Biol.
201
,
2935
-2944.
McGowan, C. P., Baudinette, R. V. and Biewener, A. A.(
2005
). Joint work and power associated with acceleration and deceleration in tammar wallabies (Macropus eugenii)
J. Exp. Biol.
208
,
41
-53.
McGuigan, M. P. and Wilson, A. M. (
2003
). The effect of gait and digital flexor muscle activation on limb compliance in the forelimb of the horse Equus caballus.
J. Exp. Biol.
206
,
1325
-1336.
McMahon, T. A. and Cheng, G. C. (
1990
). The mechanics of running: how does stiffness couple with speed?
J. Biomech.
23
,
65
-78.
Ruina, A., Bertram, J. E. and Srinivasan, M.(
2005
). A collisional model of the energetic cost of support work qualitively explains leg sequencing in walking and galloping, pseudo-elastic leg behaviour in running and the walk-run transition.
J. Theor. Biol.
21
,
170
-192.
Schwind, W. J. and Koditschek, D. E. (
2000
). Approximating the stance map of a 2-DOF monopod runner.
J. Nonlinear Sci.
10
,
533
-568.
Seyfarth, A., Geyer, H., Günter, M. and Blickhan, R.(
2002
). A movement criterion for running.
J. Biomech.
35
,
649
-655.
Taylor, C. R. (
1985
). Force development during sustained locomotion: a determinant of gait, speed and metabolic power.
J. Exp. Biol.
115
,
253
-262.
Whittackler, E. T. (
1904
).
A Treatise on the Analytical Dynamics of Particles and Rigid Bodies
, 4th edn. New York; Cambridge University Press.
Witte, T. H. and Wilson, A. M. (
2005
). Accuracy of WAAS-enabled GPS for the determination of position and speed over ground.
J. Biomech.
38
,
1717
-1722.
Witte, T. H., Knill, K. and Wilson, A. M.(
2004
). Determination of peak vertical ground reaction force from duty factor in the horse (Equus caballus).
J. Exp. Biol.
207
,
3639
-3648.