Dynamics and stability of running on rough terrains

Stability of running on rough terrain depends on the propagation of perturbations due to the ground. We consider stability within the sagittal plane and model the dynamics of running as a two-dimensional body with alternating aerial and stance phases. Stance is modelled as a passive, impulsive collision followed by an active, impulsive push-off that compensates for collisional losses. Such a runner has infinitely many strategies to maintain periodic gaits on flat ground. However, these strategies differ in how perturbations due to terrain unevenness are propagated. Instabilities manifest as tumbling (orientational instability) or failing to maintain a steady speed (translational instability). We find that open-loop strategies that avoid sensory feedback are sufficient to maintain stability on step-like terrains with piecewise flat surfaces that randomly vary in height. However, these open-loop runners lose orientational stability on rough terrains whose slope also varies randomly. The orientational instability is significantly mitigated by minimizing the tangential collision, which typically requires sensory information and anticipatory strategies such as leg retraction. By analysing the propagation of perturbations, we derive a single dimensionless parameter that governs stability. This parameter provides guidelines for the design and control of both biological and robotic runners.


S1.3 Terrain model
The terrain is modeled as piecewise linear. This is achieved by first defining a one-dimensional grid with fixed grid spacing λ. Interpolating heights between the grid points k, located at x k to intermediate points (x t , y t ) in the kth terrain patch, yields a piecewise linear, continuous terrain profile, where terrain slope m k is discontinuous at the grid points, patch k : y t = m k x t + c k , (S1a) continuity condition: where m k and c k are constants within a patch. Terrain heights at all grid points are distributed according to ∼ U(−0.03, 0.03) (main text table 1). Our choice of the uniform distribution U is to improve simulation speed, even though beta distributions described in main text figure 2e most closely matched artificial terrain used in experiments [5,6]. The range of heights h ∈ [−0.03, 0.03] and grid spacing λ = 0.1 was chosen to match the artificially constructed rough terrains [5,6].
Step-like terrains with no slope distributions were simulated by picking a height from the probability distribution prior to landing. If the chosen landing height was above the apex height of any portion of the runner, we chose another landing height from the distribution. The probability of this re-sampling occurring is ∼ 10 −4 .

S1.4 Calculating ground contact point
The aerial phase ends when the runner collides with the ground. The landing position is determined by solving for the unknown intersection point x t of the runner's aerial phase trajectory (x G (t), y G (t)) with the condition for tangential contact between runner and ground, parabolic flight: touchdown: y G = y t + 1 where b 0 , b 1 , b 2 are constants that define the aerial phase trajectory. Equations (S1)-(S2) solved simultaneously yield a quadratic equation in x t , where The larger of the two roots of this quadratic is the true landing point x P , if the roots are real and the larger of the two roots is greater than x k . The other real root is always less than x k and corresponds to the intersection of the aerial phase trajectory with the terrain patch closer to the take-off point. On flat terrain, the smaller root is the location of the take-off point. Having solved for x P , the position of the center of mass at landing is determined using equation (S2b). However, if the runner lands on a grid point, the position of the center of mass would be indeterminate as the grid point x k is associated with two slopes, m k and m k+1 . In fact, we detect a corner collision if the larger root of equation (S3a) is less than x k , or if the roots are complex. Thus, we now know the position of contact point P, x P = x k , but cannot determine (x G , y G ) at contact using equation (S2b), since x k is associated with slopes m k and m k+1 . We regularize the slope at the point x k by accounting for the aerial phase trajectory. Substituting x t = x k in equation (S3a), we write equations (S3) as a quartic polynomial in unknown m k . We numerically find all the roots using the Jenkins-Traub algorithm [7] and pick the real root that corresponds to first contact between the ground and runner, i.e. when the parabolic trajectory describing the aerial phase is above the ground.

S2 Convergence of the Monte Carlo simulations
We ran Monte Carlo simulations with different ensemble sizes to test for convergence of the steps to failure distributions. We define the criterion for convergence as the first three moments of the distribution remaining unchanged to one significant decimal place as a function of ensemble size. This criteria is met with an ensemble size of 10 4 (figure S1). Hence we perform Monte Carlo simulations in the main text with an ensemble size of at least 10 5 .

S3.1 Passive collision
A runner with center of mass velocity v − G and angular velocity ω − collides with a terrain patch angled at θ with respect to the horizontal. By simultaneously solving equations describing the collision at the contact point P (main text equation (1a)), and using kinematic constraints associated  with a rigid body, the linear velocity v c G and angular velocity ω c of the runner after the collision are computed as, and and v G,t = v G,x cos θ + v G,y sin θ. (S4h)

S3.2 Relationship between t and tc
Based on a detailed consideration of how an animal may use its foot speed to control the tangential collision, we arrive at equation (2) in the main text. The velocity of the foot and center of mass are not constrained to be the same for a running animal. For example, the animal may retract the leg to S4 vary the foot's speed independent of the center of mass speed. We make the simplifying assumption that the foot comes to rest, or nearly so, upon colliding with the ground. In an animal, there is an additional degree of control corresponding to modulating the transmission of the collision between the foot and center of mass through changes in the leg posture and muscle co-contraction. We do not include this additional mode of controlling collisions in our consideration. The use of the foot enables the runner to reduce the effective collision independent of the body's overall linear and angular momentum. Thus the tangential coefficient of restitution t parametrizes the loss of speed of the foot relative to a point on a rigid runner that does not have the ability to modulate its foot speed. The foot speed of such a rigid runner is that of the contact point P in the rigid running model (main text figure 1b), one that translates with the center of mass and is also affected by any spin that the body may have. The tangential coefficient of restitution t at the moment of collision is therefore given by, These velocities are in a reference frame that is parallel to the tangential and normal directions of the terrain at the point of contact (t −n frame in figure S2). Point P refers to the contact point in the rigid runner that is coincident with the foot at the time of contact, but whose speed is solely a function of the body's linear and angular momentum. When the foot speed exactly matches ground speed as seen from the body-fixed frame, t = 1 and the tangential collision impulse would be zero. To an external observer, the foot's velocity vector would be exactly normal to the terrain at the point of contact and v − foot,t = 0. If the foot's speed equals that of point P of the rigid runner, then t = 0 and the runner experiences a significant tangential collision impulse.
The open-loop and anticipatory strategies differ in how they specify the foot speed. The intended coefficient of restitution tc for the open-loop strategy is specified under the assumption that the tangential speed of the contact point P equals the initial forward speed, v − P,t = v x0 . Furthermore, the open-loop strategy is only aware of the foot speed in a body-fixed frame, and its components relative to gravity, i.e. the components of v − P − v − foot in the x-y frame. These capture the three assumptions outlined in the main text (section 2.2). The intended coefficient tc for the anticipatory strategy makes no such assumption and simply follows equation S5. The respective equations for the anticipatory and open-loop tc are, By combining equation (S5) and equation (S6) we find : open-loop, tc : anticipatory. (S7) For small terrain angles θ, the tangential and horizontal components of the foot speed relative to that of point P are nearly equal, i.e.  Figure S2: The terrain-fixed push-off impulse is the laboratory-fixed push-off impulse rotated by the terrain slope angle θ about the contact point.
After the passive collision, the push-off impulse J imp is applied at the contact point P to propel the runner into the flight phase. Recall that J φ = 0 for the open-loop push policies that we consider (main text equation (1c)). The push-off impulse J imp leads to discrete changes in the linear velocity of the center of mass v imp , and angular velocity of the runner ω imp . As discussed in main text section 2.3, a terrain-fixed push-off policy or a laboratory-fixed push-off policy both satisfy the periodicity criteria on flat ground, namely v + G,steady = (v x0 , v y0 ) T and ω + = 0, but differ in their behavior on rough terrain. The linear velocity v + G and angular velocity ω + at take-off under these push-off policies are, where : terrain-fixed , and S4 Open-loop runners on rough terrain

S4.1 The open-loop passive collision
The tangential collision parameter t for open-loop runners varies from step-to-step on rough terrains (main text equation (2)) and is thus characterized by distributions (figure S3a) that evolve according to running dynamics described in main text equation (1). These distributions achieve stationarity as seen in the representative case of the t distribution with tc = 0.4, that appears to converge by 5 steps (  The stability benefits of higher mean t values (figure S3b) are mitigated by higher fluctuations in t at higher tc values (figure S3d) leading to larger step-to-step fluctuations in body angular momentum. Hence for open-loop runners, mean steps to failure shows a weak dependence on tc value. In main text section 5, we arrive at this same conclusion via an asymptotic analysis of the series expansion of the orientation change over a single step caused by an unexpected terrain slope perturbation.

S4.2 Additive noise in open-loop strategies
Additive noise in the tangential collision (main text section 3.5) of open-loop runners has little effect on mean steps to failure (figure S4) compared to the same noise intensity applied to anticipatory runners (main text figure 3b). For open-loop runners with human-like parameters, there is no optimal value of tc (figure S4) unlike the case with anticipatory runners where slight tangential collisions are optimal (main text figure 3b). The relative insensitivity of failure statistics to added noise for the open-loop runner is possibly because t varies significantly on rough terrains even in the absence of added noise (figure S3).

S5 Anticipatory runners on rough terrain
In main text section 3.4 we describe how eliminating the tangential collision impulse stabilizes body orientation, but deviations from this strategy quickly lead to a loss in stability. The dominant failure The mean steps to failure for an anticipatory runner with tc = 1 is captured by a single parameter v x0 /∆v x (figure S5b). The logic is similar to that in main text section 5 where the scaling analysis is for orientation failures, while here the failure mode is due to losing translational stability (main text figure 1c). For a runner landing on sloped ground with forward speed v − G,x = v x0 , the loss in forward speed over a single step due to the terrain perturbation is ∆v

S8
Following the logic outlined in main text section 5, the mean steps to failure N should scale as, where ∆v x = sin θ(v y0 (cos θ − 1 + n (cos θ + 1)) + v x0 sin θ(1 + n )). (S9b) Equation (S9c)   The linear stability analysis for the open-loop runner proceeds analogously to the analysis for the anticipatory runner discussed in main text section 4. We define a Poincaré section transverse to the runner's trajectory in phase space at the apex of the aerial phase (v y = 0), and a return map f ol that maps the state of the runner ψ n at the apex of the aerial phase at step n, to the state at the apex of the aerial phase on the following step ψ n+1 (main text figure 5). The state ψ, return map f ol , and its linearization T ol are defined by, The linearized return maps, T ol and T an are non-diagonalizable as discussed in main text section 4.
Thus an eigen-factorization of these matrices is not possible, so we perform a Jordan decomposition of T ol and T an . The Jordan decomposition of T ol is, and The linearized return map T ol has rank 4, with two stable eigenvalues; λ = 0 and λ = n (2 n − 1). The second stable eigenvalue corresponds to perturbations to the height of the apex of the aerial phase. The remaining eigenvalues are λ = 1 which are associated with eigenvectors ν 1 , ν 2 , and generalized eigenvector ν 3 , the same as the eigenvectors and generalized eigenvector associated with eigenvalues λ = 1 for T an (main text equation (7a)). Thus, a perturbation in the subspace spanned by these vectors ∆ψ 0 displays the same scaling with step number n as described in main text equation (8) for the anticipatory runner. These relationships for the open-loop runner are summarized as follows, and substituting a 1 = √ 2 n v y0 , and a 2 = which are the off-diagonal terms of J ol (equation S11b). While the growth of the instability shows the same dependence on step number n as the anticipatory runner (main text equation (8)), the value of a 2 differs, i.e. the projection onto ν 2 of a perturbation along ν 3 upon action of the return map T ol . The stability of the anticipatory runner is discussed in detail in main text section 4. The Jordan decomposition of T an when tc < 1 is, where and The anticipatory runner has a full rank return map T an , and stable eigenvalues λ = n (2 n − 1), which is identical to the open-loop runner, and λ = t which differs from the open-loop runner.
The off-diagonal elements of J an are a 1 = √ 2 n v y0 and a 2 = − √ 2 n v y0 (equation S13b), which correspond to the projection onto ν 1 and ν 2 respectively, of a perturbation along ν 3 upon action of the return map T an .
When tc = t = 1 for the anticipatory runner, the linearization of f an about ψ * (main text equations 6) yields a different form for T an . The Jordan decomposition of T an when tc = 1 is, where There is a single stable eigenvalue of λ = n (2 n − 1) which is identical to that of the openloop runner and the anticipatory runner with tc < 1. The remaining 4 eigenvalues are λ = 1 which are associated with eigenvectors ν 1 , ν 2 and generalized eigenvectors ν 3 and ν 4 (main text equations (9a), columns of V an from equation (S14c)). Perturbations along ν 3 project onto ν 1 , and perturbations along ν 4 project onto ν 2 , upon action of the return map T an (equation (S14b)). The off-diagonal elements of J an are a 1 = a 2 = 2 n v y0 .

S7 Scaling analysis of orientation failures
In main text section 5 we show that the orientation failure bound φ tol and the orientation change over a single step from encountering an unexpected terrain slope φ • predict the mean steps to failure S11 N , as N ∼ φ tol /φ • . The many parameters that define the runner and its gait do not separately predict mean steps to failure as accurately as the parameter φ tol /φ • ( figure S7). The change in orientation over one step can be calculated using the vertical velocity v + y,• and angular velocity ω + • at take-off as φ • = 2v + y,• ω + • (• = 'ol' for open-loop and 'an' for anticipatory). For a runner landing on a terrain patch angled at θ with respect to the horizontal, with center of mass velocity v − G = v x0 −v y0 T and angular velocity ω − = 0, the take-off velocities v + y,• , ω + • , and orientation change φ • , are given by The parametric dependence of φ tol /φ • (φ • defined in equation (S15e)) on tc and n is captured by a power series approximation of φ • to second order in θ (main text equations (11a)-(11b)). The contour plots of φ tol /φ • shown here using the truncated power series form of φ • (figure S8) display the same qualitative trends as the ones in main text figure 7 which use the complete expression for φ • (equation (S15e)). Therefore, we analyze failure statistics and runner morphology in main text section 5 and section 6 using the truncated power series approximation for φ • .

S8 Steps to failure statistics
All the numerically estimated distributions of the steps to failure show similar qualitative features. They are unimodal with a long tail, and are relatively insensitive to changes in the terrain height distribution (figure S1, main text figure 2e). The similarity of the distributions along with the results of the linear stability analysis (main text section 4) suggests that the stochastic dynamics underlying these distributions may be characterized by a small number of parameters. The linear stability analysis shows that perturbations from the ground do not directly affect linear or angular velocities, but affect position variables (main text section 4). For orientation variables, this implies that the rate of change of the angular velocity is zero. However, the terrain introduces perturbations at every step. When the correlation length of the random terrain is much smaller than a step length, like in the Monte Carlo simulations, the step-to-step perturbations may be treated as uncorrelated random perturbations. Irrespective of the distribution from which the random terrain is constructed, we make a heuristic appeal to the central limit theorem and posit that the net effect of the terrain perturbations is simply a Gaussian random forcing. Thus the stochastic step-to-step model for the orientation is that the discrete second derivative of the angle φ is forced by a Gaussian random variable: where w ∼ N (0, σ 2 ), (S16b) The random forcing w is drawn from a Gaussian with mean < w >= 0 and variance < w 2 >= σ 2 .
In the absence of noise, the solution to the second-order difference equation yields φ(n) = n(φ(1) − φ(0)) + φ(0). Thus the orientation of the runner φ grows linearly with step number n in response to an initial perturbation (φ(1) − φ(0)) as found in main text section 4.
Orientation failures occur when a runner's orientation exceeds φ tol . Thus, the probability of having failed by n steps, which we denote with p fall (n), is given by where erf( is the error function. The probability of failing at step number n, denoted by p steps (n), is given by the probability of taking n − 1 steps without failing, and then failing at step number n, i.e. p steps (n) = p(|φ| > φ tol , n |φ| < φ tol , n − 1). This is equivalent to the probability flux at φ tol at step n, given by =⇒ p steps (n) = erf φ tol The mean steps to failure N are thus given by We use the approximation f (n) = n 3 /3 and ζ(p) is the Riemann-zeta function defined by ζ(p) = ∞ i=0 1/i p for p > 1. The mean steps to failure depend on a single parameter, φ tol /σ, the ratio of the tolerance-angle φ tol to the intensity of the external noisy forcing w. We numerically solve for the infinite sum in equation (S20b) for figure S9c,d.
The probability of failing at step n, i.e. p steps (n) (figure S9a), and its cumulative distribution p fall (n) (figure S9b) are qualitatively similar to the numerically computed steps to failure distributions shown in figure S1 and main text figure 2e, and the cumulative distributions shown in main text figure 2a,b, in that they are unimodal with long tails.
We compare the mean of the numerically computed steps to failure distributions and the analytically derived mean from the stochastic model ( figure S9c,d). The standard deviation of the noise S14 S9 Effect of changing runner size relative to terrain grid spacing The size of the runner r relative to the grid spacing of the terrain λ determines the range of terrain slopes accessed by the runner (main text section 6, figure S10). For example, if λ/r = 1, the range of terrain slopes encountered by the runner is nearly the same as the slope distribution of the terrain itself (figure S10a). By decreasing this ratio to λ/r = 0.01 while keeping the terrain slope distribution constant, we find that the range of slopes encountered by the runner is significantly lower (figure S10b). We quantify this trend using the parameter σ runner /σ terrain where σ runner is the standard deviation of the slope distribution encountered by the runner, and σ terrain is the standard deviation of the terrain slope distribution. The parameter σ runner /σ terrain approaches 1 as λ/r increases (figure S10c). The mean slope encountered by the runner is always greater than zero (figure S10d) even though the slope distribution of the terrain has a zero mean. As the range of slopes encountered by the runner increases, the mean slope encountered by the runner deviates further from zero (figure S10d). This is consistent with the observation from the Monte Carlo simulations that runners slow down on rough terrain since the effect of the mean slope encountered is to redirect forward momentum Figure S10: Effect of varying the ratio of leg length r to terrain grid spacing λ. The probability density function of the terrain slope (tan θ) encountered by the open-loop runner (orange) with human-like parameters is plotted for grid spacing (a) λ = r and (b) λ = 0.01r . The height distribution of the terrain grid points was adjusted so that the terrain's slope distribution (blue) was the same in all simulations. The number of steps used to generate these probability density functions was ∼ 10 6 . (c) Shows the standard deviation of the slope encountered by the runner relative to the standard deviation of the terrain's slope distribution. (d) The mean slope encountered by the runner as a function of grid spacing λ.
into vertical momentum.

S10 Monte Carlo simulations with a laboratory-fixed push-off
The change in angular velocity due to the push-off impulse J imp applied as per the laboratory-fixed push-off policy destabilizes the runner as it exacerbates angular momentum fluctuations that arise due to variations in terrain slope angle θ. When θ > 0, the runner collides with lower tangential velocity compared to flat ground, resulting in a smaller clockwise pitching moment than if the runner were colliding with a flat terrain patch. Yet the clockwise pitching moment induced by the the push-off impulse J imp is larger than what it would be on flat ground (equation S8d), causing the runner to spin excessively in the clockwise direction. Conversely, when θ < 0, the runner lands at a shallower angle, increasing the tangential collision impulse and inducing a larger clockwise angular impulse on the runner compared to what it would experience when θ = 0. However, the push-off impulse induces a smaller counter-clockwise change in angular velocity at push-off ω imp relative to flat ground (equation S8d), causing the runner to spin in the anti-clockwise direction. In contrast, the angular velocity change due to the push-off impulse under the terrain-fixed push-off policy is independent of the terrain slope θ (equation S8d). As we have seen in the main text this push-off cannot stabilize the runner. However, it does not amplify angular momentum fluctuations caused by slope variations like the laboratory-fixed push-off policy. Thus, the laboratory-fixed push-off policy leads to orientation failures in fewer steps compared to the terrain-fixed push-off policy.

S10.1 Open-loop runners
Open-loop runners with a laboratory-fixed push-off and human-like parameters only fail due to losing orientation stability (figure S11a). Increasing energy dissipation in the direction normal to the terrain surface reduces mean steps taken by these open-loop runners. This trend contrasts that observed with the terrain-fixed push-off where increasing normal energy dissipation increases the steps taken (main text figure 2d). The increase in mean steps taken with the laboratory-fixed push-off is at most 1 for a 100% reduction in energy dissipation, from n = 0 to n = 1 at a fixed value of tc (figure S11b).
Reducing tc increases steps taken (figure S11b), unlike the trend observed with the terrainfixed push-off policy where changing tc had negligible effect on the mean steps to failure (main text figure 2d). Open-loop runners increase mean steps taken by 2 when tc is reduced from 1 to 0 while holding n fixed at any value (figure S11b).
Open-loop runners with a laboratory-fixed push-off policy perform best when n = 1, tc = 0 and worst when n = 0, tc = 1. The difference in mean steps taken at these extremes is just 3 steps. The maximum mean steps taken by the open-loop runner with a laboratory-fixed push-off is 9 steps (figure S11b), the same as the minimum number of steps taken by the open-loop runner with a terrain-fixed push-off policy (main text figure 2d). Thus, compared to the terrain-fixed push-off policy, open-loop runners with the laboratory-fixed push-off policy fare worse, and display differing trends in how stability is determined by the collision parameters, n and tc . However, similar to the terrain-fixed push-off policy, varying parameters governing the passive collision has a small effect on the mean steps to failure compared to similar changes in collision parameters in anticipatory runners.

S10.2 Anticipatory runners
Anticipatory runners with a laboratory-fixed push-off policy increase mean steps taken as energy dissipation in the direction normal to the terrain surface is reduced. This is in contrast to the trend observed for anticipatory runners with the terrain-fixed push-off policy, where increasing normal energy dissipation increases mean steps to failure (main text figure 3a). Anticipatory runners fail by losing orientation stability except when tc = 1 and n = 1 (figure S11a), in contrast to the terrainfixed push-off policy where anticipatory runners maintain orientation when tc = 1, and regardless of the value of n (main text figure 3a). At tc = 1, anticipatory runners with a laboratory-fixed push-off, increase mean steps taken by over 4-fold, from 7 to over 30 for a 100% reduction in normal energy dissipation from n = 0 to n = 1 (figure S11c).
Reducing the tangential collisional impulse at landing tc increases mean steps to failure for the anticipatory runner with a laboratory-fixed push-off if n is sufficiently high, while showing the opposite trend below that n value (figure S11d). This is in contrast to the trend observed in anticipatory runners with the terrain-fixed push-off policy where mean steps taken always increases when tc is reduced (main text figure 3a).  Figure S12: Effect of sensorimotor noise on stability, and scaling analysis of mean steps to failure. (a) Contour map of mean steps taken as a function of tc and ∆ t shows that mean steps to failure reduces with increasing sensorimotor noise ∆ t . Scuffing the ground upon landing is optimal as noise intensity ∆ t is increased. Red circles denote the optimal tc value for a given ∆ t . In these simulations n = 1. (b) For the anticipatory runner with tc = 1, n = 1, the mean steps to translational failure are predicted by the parameter v x0 /∆v x . The values of the remaining parameters used in this simulation are shown in figure S7. Mean steps to orientation failure are captured by the parameter (c) φ tol /φ ol for the open-loop runner, and by (d) φ tol /φ an , for the anticipatory runner. The parameter values used for these simulations are also shown in figure S7.