A simple model of a quadruped discovers single-foot walking and trotting as energy optimal strategies

It is widely held that quadrupeds choose steady gaits that minimize their energetic cost of transport, but it is difficult to explore the entire range of possible footfall sequences empirically. We present a simple model of a quadruped that can spontaneously produce any of the 2500+ planar footfall sequences available to quadrupeds. Through trajectory optimization of a work and force-rate cost, and a large sample of random initial guesses, we provide stronger evidence towards the global optimality of symmetrical four-beat walking at low speeds and two beat running (trotting) at intermediate speeds. Using input parameters based on measurements in dogs, the model predicts the correct phase offset in walking. It also spontaneously reproduces the double-hump ground reaction force profile observed in walking, and the half-sine profile observed in trotting, despite the model’s lack of springs. However, the model incorrectly predicts duty factor at the slowest speeds, and does not choose galloping as globally optimal at high speeds. These limitations might point to missing levels of complexity that could be added to future quadrupedal models, but the present results suggest that massive legs, elastic components, head-neck oscillations and even multi-linked legs are not critical determiners for the optimality of natural quadrupedal walking and trotting. Author summary

But why are particular gaits energy minimizing for particular animals at particular 23 speeds? Some have pointed to elastic elements as a key aspect [14][15][16]. Gan et al. [15] 24 found that natural gaits of horses could be replicated by a passive dynamic model with 25 springs as legs. Xi et al. [16] found that a "four-beat walk" is work minimizing at slow 26 speeds and a two-beat run (trotting or pacing) minimizes work at high speeds. They 27 pointed to their use of elastic elements as an important reason why their model 28 exhibited a double-hump ground reaction force profile in walking. This is a surprising 29 finding, given that the double-hump profile emerges in bipedal walking as a collision 30 mitigation strategy [17,18], even in the absence of springs [19][20][21]. 31 If elastic elements are prerequisites for gait optimality, it has important implications 32 for the evolutionary history of gait. In particular, it suggests that appropriately-tuned 33 tendons must have evolved before the gaits we see in modern organisms first appeared 34 in their ancestors [50]. 35 We may do well in taking guidance from bipedal research using work-based 36 optimization methods. By incrementally increasing the complexity of the models used, 37 we can determine what parameters are necessary prerequisites for generating the 38 phenomenon of interest. Srinivasan and Ruina [20] used a "minimal" biped; effectively, 39 a point mass with one massless leg 1 . By minimizing work at a given step length and Importantly, both models used no elastic elements, suggesting that elastic morphological 48 features are not necessary prerequisites to the optimality of the solutions observed (and 49 more likely, evolutionary post-adaptations). 50 Further, as more degrees of freedom are added, the parameter space of possible gaits 51 tends to increase in scope and dimensionality, and so is necessarily more difficult to 52 explore. Xi et al. [16] relied on preset footfall sequences to seed optimization, bringing 53 into question whether other footfall sequences not explored in the study might have 54 been energy-minimizing. 55 Finally, while limb work was the only cost in the above optimization models, there is 56 increasing evidence from the human physiology literature that muscle activation costs, 57 in addition to muscular work, are an important determiner of gait. That is, there is 58 increased metabolic cost for changing muscular forces quickly, even if work is held 59 constant [22][23][24]. While the mechanism is currently unknown, rate penalties explain 60 much about gait selection, particularly the stride length and duty factor vs. speed 61 relationship [19,25,26]. 62 To what extent can a simple, minimally constrained energy-optimization model, 63 without springs, forced pendular exchange or forced stability, predict gait choice in 64 quadrupeds? Do natural gaits emerge from such a simplified model and are they energy 65 1 As it was constrained to perform symmetrical gaits with no double-stance, the single leg could be reflected into the other leg for the next step, effectively representing a biped.
March 17, 2019 2/31 minimizing? Does a force-rate cost substantially enhance empirical agreement? We 66 explore these questions with a planar, two-point-mass optimization model that can 67 perform any of the 5040 footfall sequences available to quadrupeds. We compare the 68 simulation results to empirical data on dogs only, but we anticipate the results are more 69 broadly applicable across quadrupedal mammals. 70 Methods

71
Trajectory optimization of a minimal quadruped 72 Generalized model of a quadruped 73 We use a simplified model of a quadruped based on the bipedal model in [20], but 74 extended to four legs. A rigid trunk connects two point masses. Each point mass has 75 two massless limbs that can extend or contract. The quadruped is two-dimensional and 76 lives in the sagittal plane (Fig 1).
The simple quadrupedal model used in this paper. Two point masses sit on massless legs that can extend and contract. The point masses are connected with a rigid trunk. Five morphological parameters are necessary as input: the mass of the fore and hind quarters (m F and m H , respectively), the trunk length (l b ), and the maximum length of the fore and hind limbs (l Fmax and l Hmax , respectively). Two kinematic variables are also necessary: the stride length (D) and the average horizontal speed (U ). From these seven constant inputs, the optimal control chooses footfall positions (f ij ) and over twenty independent temporal variables, including ground reaction forces and their rate of change (F ijk andḞ ijk ), center of mass position and velocity (x andẋ) and body pitch angle and angular velocity (θ andθ).
Each limb comes as a "Trailing" and "Leading" virtual limb pair. A "trailing" limb 78 acts through a rear footfall position, and a "leading" limb acts through a forward temporal ones). Since we are interested in periodic motion, we can always shift the 83 temporal frame of reference so that, for one limb, touchdown occurs before liftoff. For 84 this limb, only one virtual limb and one footfall position is necessary. 85 The organism travels in a vertical gravitational field with downward acceleration g 86 at an average horizontal speed U in stride length D. The quadruped has forequarter 87 mass m F and hindquarter mass m H . Its maximum limb lengths are l Fmax and l Hmax for 88 the fore and hindlimbs, respectively. 89 We use trajectory optimization techniques (sometimes referred-to as optimal control) 90 to discover optimal gaits (for a beginner's overview see [27]; for a more advanced 91 treatment see [28]). The basic idea is that the aspects of gait the organism can directly 92 control (e.g. footfall positions and forces in time) are chosen so as to minimize some 93 cost function within particular constraints. These constraints can include physical or 94 geometric limitations (e.g. mechanical dynamics or maximum limb length) or even 95 biological limitations (e.g. an upper bound on force output or muscle shortening rate). 96 Our hypothesis is that quadrupedal mammals choose whichever gait minimizes 97 metabolic cost per unit distance. The primary source to energetic cost is limb work, the 98 combined positive and negative work performed by each limb: where T is the stride period, F i is the absolute force acting through limb i and l i is the 100 length of limb i. l i is defined as the distance between the fore-or hindquarters and 101 footfall position f i (see Appendix A). Note that we make no distinction between 102 positive and negative muscle work efficiencies. While these efficiencies can alter optimal 103 solutions in unsteady gaits [29], in steady gaits on level ground an equal amount of 104 positive and negative work must be performed. Therefore, in the absence of other forms 105 of dissipation (like foot-ground collisions or viscous effects) or asymmetric elastic energy 106 recovery, the relative efficiencies of positive and negative work do not change the 107 optimal activation pattern.

108
The secondary cost is that of a force-rate penalty. Since the exact choice of rate 109 penalty has little effect on the agreement with empirical data [30] (some fitting is 110 required in any case, since the mechanism is unknown), we used a well-behaved squared 111 cost in force rate, where c 1 is a fitting parameter.
The resulting objective 3 of the optimization, reflecting the metabolic cost of 114 transport in the absence of a basal metabolic rate, is where D is the stride length.  contact sequence a priori, an advantage over the multiple shooting method used by [16]. 148 One disadvantage of this method is that solutions exhibit a discontinuity at ground 149 contact, which requires special treatment and makes the problem more difficult to solve. 150 The relative simplicity of our optimal control problem makes the implementation 151 tractable, despite this disadvantage.

152
Constraint P5 ensures that legs are always ventral to the body, a sensible limitation 153 for cursorial quadrupedal mammals. Curiously, if this constraint was removed, the 154 simulation frequently discovered bipedal running. This finding deserves further analysis 155 as it may have implications for the evolution of bipedal locomotion, but is beyond the 156 scope of the present study.
The rationale is that breaking this constraint would cause the animal to turn or tip. 165 Including the constraint lowered solve time, while removing the constraint did not seem 166 to change the optimal solution. 167 Variables in optimal control problems may be bounded in some way. For the present 168 case, the ideal bounds on variables and time are:  Numerically well-conditioned optimal control problem 184 Three numerical issues arise in the ideal implementation. The first is the presence of a 185 non-smooth objective due to the absolute-value function of the integrand (equation 1). 186 The second is due to non-smooth optimal solutions. The final issue is due to poor  Non-smooth functions in optimal control problems can be reformulated using slack 190 variables [28]. Two slack variables are required for every instance of the absolute value 191 function; thus 14 new controls are introduced (p i ≥ 0 and q i ≥ 0). In addition, new path 192 constraints are required, one set for redefining the operand of the absolute value 193 function, and another for enforcing orthogonality of the slack variables: Unfortunately, Eq 4 is a complementarity condition, as are several other constraints 195 (P3, P4 and P5) which require special treatment [28,32]. Complementarity constraints 196 are of the form, where a i (t) and b i (t) are arbitrary variables [28]. One way to manage these conditions 198 is to introduce "relaxation parameters" that effectively smooth the constraint, as These relaxation parameters decrease the accuracy of the solution (compared to the  We can further ratchet down the relaxation parameters by incorporating them into 203 the optimization problem. We make s i decision variables (controls), and augment them 204 with the penalty where c 3i is a constant chosen to be sufficiently large [32]. Then, the optimizer tends to 206 reduce the value of s i , making the solution more accurate.

207
While it is not necessary to use both the strategy 5 and 6, we found that 208 incorporating both reduced solve time and eased implementation. Our strategy was to 209 start by not enforcing complementarity conditions (c 3i were 0 and s i was allowed to 210 take on any positive value) and then to increased the coefficients c 3i as we increased the 211 resolution on each iteration. An exception to this was in the slack variable 212 complementarity condition itself (Eq 4), which, on every iteration except the first, was 213 not enforced except through the augmentation where c 2 = 1 E -3.

215
The resulting cost function for the optimization (prior to normalization), was The final issue, poor scaling, was partially helped by normalization (next section), 217 but further ameliorated through the automatic scaling feature 218 ('automatic-hybridUpdate') in GPOPS-II [33]. Automatic scaling requires all 219 unbounded variables to be redefined with finite bounds. Selected bounds are shown in 220 Supplementary Table 1, and were chosen to be much larger than the expected 221 magnitude of the variables, but small enough to improve scaling of the optimal control 222 problem (see Appendix B for further details).
In this paper, normalized variables are noted with a prime exponent (e.g. x → x ′ ).

228
Since the force rate coefficient c 1 (Eq 2) may have biological significance, it was also constantT D as documented in [34]. For all other cases, c ′ 1 (T ) was scaled as The augmentations to the cost function (Eq 8) were not scaled.

235
While the characteristic length for the computational problem was l b , we will report 236 non-dimensional speed as the square of the Froude number, following a common convention in other literature [1,34].

244
Numerical solutions to optimal control problems require an initial guess. Values at 245 each of the 16 initial grid points for each variable were chosen from a uniform 246 distribution within the bounds of each variable. An exception is the horizontal position 247 state variable, which was set as [0, D] for the start and ending conditions, respectively, 248 and interpolated linearly between those values for intermediate grid points.

249
The algorithm for finding the pseudo-global optimum (the optimum among all valid 250 solutions) is as follows. An initial randomized guess was passed to GPOPS-II with one 251 mesh refinement step. For this initial step, the slack variable orthogonality constraint 252 (Eq 4) was explicitly enforced with c 2 and all c 3i set to 0 (Eq 8). This step allowed the 253 program to more quickly approach a feasible region. The output was then  For exploring the effect of c 1 on the solution in one test case, 100 initial guesses were 263 used. For the three detailed test cases, 500 initial guesses were used. For comparison to 264 the large empirical dataset from [37], the number of guesses depended on the speed,  Each solution was determined as "valid" on the basis of (1) satisfactory mesh error 270 (lower than the tolerance of 1e-4) and (2) satisfactory (low) constraint violation. The 271 second criteria was especially critical for complementarity conditions, which were not 272 explicitly enforced and were violated more often than the other constraints.

273
Among valid solutions, the lowest cost solution was selected as the (pseudo-)global 274 optimum as determined by J tot (Eq 3). That is, since constraints were checked through 275 the above validation step, the augmented parts of the objective (regarding constraint 276 violation) were ignored as part of the relative optimality of one solution as compared to 277 the others.

278
Empirical data 279 The model requires five inputs: stride length, speed, body length, center of mass 280 position, and limb lengths. As validation of the model, duty factor, phase offset, and 281 ground reaction force profiles are also required. These empirical data were acquired for 282 single speeds from four studies [34,38,39], representing dogs moving at a slow walk, 283 medium walk, and trot, respectively. In addition, [37] provides curves for duty factor, 284 stride length, and phase offset as a function of speed from a large dataset. While all 285 these studies supply a length calibration (e.g. shoulder height, withers height), they do 286 not supply all necessary body proportions (except for [34], which supplies shoulder and 287 hip locations in Figure 12 of that study). These additional proportions were taken from 288 Table 1. Empirical data used for validation and as input (bold) for the model. PL is pair lag, the fraction of stride between left hind and left fore touchdown (sometimes referred to as pair offset or phase offset). DFh and DFf are mean duty factor of the hind and fore limbs, respectively. *Not reported at the original source [37]; derived from [38] Main Source: [ website [40][41][42]. The acquired parameters are summarized in Table 1.

290
Forelimb length was taken as the length from manus to shoulder joint in standing. Hindlimb length was taken as length from pes to hip joint in standing. Body length was taken as hip to shoulder joint in standing. Masses were assumed concentrated at shoulder and hip, following [38]. To determine center of mass location (or, equivalently, the anterio-posterior body mass distribution [43]), a balance of torque approach was used following [44]. Briefly, since the net moment about the center of mass must be zero during steady, periodic locomotion, the net torque produced by the hindlimbs in the sagittal plane must be equal and opposite to the forelimbs. Assuming the trunk is approximately horizontal during the stride, then where DFf and DFh are the duty factors of the fore and hind limbs, respectively,ŷ is 291 the vertical unit vector, F F and F H are the fore and hindlimb ground reaction forces, 292 and t s is "stance time" (time during stance for a given foot, normalized to stance  , the solution becomes more impulsive, as expected from a work-minimizing "bang-coast-bang" solution [31]. As c ′ 1D increases, the force peaks become more shallow, but in all cases the optimal solution maintains the double-hump profile characteristic of walking bipeds and quadrupeds under most situations. The morphological data is based on a Dalmation from [34] (see also Fig 4B) The simulation was run for the morphology, speed and stride length derived 302 from [34], as outlined in Table 1 for a number of rate penalties (c 1 in Eq 2). Empirically, 303 the prescribed speed corresponded to a moderate walk [9]. The simulation spontaneously 304 discovered a four-beat walk 5 as cost-minimizing regardless of rate penalty.

305
For all rate penalties, the simulation predicts a "double-hump" shape in ground 306 reaction force to be optimal (Fig 2), a common characteristic of both bipedal and 307 quadrupedal walking. In bipedal locomotion, this shape arises from the pre-heelstrike 308 pushoff [46], a dynamic "trick" for minimizing energetic losses at foot contact. The As force-rate penalties increase, duty factor increases to slow the generation and reduction of force. Optimal duty factors remain close, but exhibit some divergence at larger c ′ 1D . Phase offset decreases slowly with force rate penalty. At no point do optimal duty factors and phase offset simultaneously match the empirical values (dotted lines) from [34]. The morphological data is based on a Dalmation from [34] (see also Fig 4B).
do. The same "trick" seems to be discovered by the model; a second peak in the trailing 312 limb ground reaction force occurs prior to the first peak in the leading limb.

313
The amplitudes of the peaks are affected by the magnitude of the force-rate cost. A 314 purely work-based objective is expected to have impulsive foot-down and 315 push-off [20,21,30,31], but this would incur a large force rate penalty. The model 316 compromises on a solution with smooth GRF, with larger peaks if the force rate penalty 317 is small, and smaller peaks if the penalty is large.

318
Five multipliers on rate penalties are compared in Fig. 3. As rate penalty increases, 319 both fore and hind limb duty factors increase and remain relatively similar, but exhibit 320 increasing divergence. Duty factor is strongly dependent on rate penalty because mean 321 vertical force must remain constant (supporting weight), but longer actuation time 322 allows the required force to be supplied over a longer time. At the lowest rate penalties, 323 duty factor and phase offset approach the work minimizing solutions of 0.5 and ≈ 0.25, 324 respectively, for symmetrical walking [47]. Phase offset is fairly insensitive to rate 325 penalty, but exhibits a slight decrease with increased force-rate cost.

326
While either fore or hind duty factor can be tuned to be arbitrarily close to 327 empirical values, matching both is not possible through adjusting rate penalty alone in 328 this case. Furthermore, since phase offset is relatively insensitive to rate penalties, it is 329 likewise difficult to tune. The results of this analysis imply that, while adding a rate 330 penalty enhances agreement with empirical data, it alone cannot explain the 331 discrepancy between work-minimizing solutions and empirical results, at least for this 332 highly simplified model. This finding is in contrast to work on human locomotion, 333 where the simple addition of a force-rate penalty enhances empirical agreement 334 substantially [19,30]. 335 c ′ 1D = 3 E -3 matches the empirical GRF profiles (Fig 2) and duty factors (Fig 3) for 336 this test case fairly well. It is used as the coefficient of choice for further exploration in 337 the remaining test cases.

338
Detailed comparison to four test cases 339 Three cases are examined in detail in Fig 4. In the first two cases, the optimization 340 correctly predicts a symmetrical, four-beat walking solution with double-hump ground 341 reaction forces. In the third and fourth case, the optimization correctly predicts 342 two-beat run with near-simultaneous fore-hind contact and single-hump GRFs. The 343 optimization also correctly predicts the gait to be a "grounded run", in the sense that it 344 exhibits in phase potential and kinetic energy changes (Fig 5H), but no flight phase [48]. 345 For trotting, two empirical results from related studies [39,45] are compared to the same 346 optimization solution, demonstrating the range of "normal behaviour" in trotting for 347 this breed. Overall, the simulation matches either empirical result fairly well, 348 demonstrated by the substantial overlap of stance and similar shape of the GRFs, but 349 matches neither perfectly. For example, the simulation overestimates hind duty factor in 350 one and underestimates fore duty factor in the other.

351
The substantial overlap between optimization and empirical gait diagrams 352 demonstrates the success of the optimization approach. It must be emphasized that any 353 of the 2520 planar quadrupedal event sequences could have been selected by the 354 optimization, and even gaits excluding certain limbs. Indeed the optimizer discovered 355 many gaits as locally optimal (Fig 5 shows some examples), including asymmetrical 356 (Fig 5A,B,E,F) and symmetrical (Fig 5C,D,G,H), some that can be related to natural 357 gaits (Fig 5C,D,F-H), some that are clearly unnatural (Fig 5A,B), and others that used 358 fewer than four limbs (Fig 5E) Yet the optimal gaits are readily identified as those 359 chosen by the organism at those speeds (Fig 5D,H) and match the empirical solutions 360 fairly closely (Fig 4B-D).

361
Only the second case ( Fig 4B). was "tuned" with the force rate penalty; the same 362 penalty was used in all other cases. Yet even within this tuning, the optimization 363 always selected a gait that was qualitatively similar to the empirical data Fig 2 and 3 . 364 These results come despite fairly rough empirical estimates, varying morphological 365 parameters, and many simplifying assumptions in the model (avoiding springs, most 366 joints, and using the simplest possible mass distribution).

367
Forces are predicted fairly well, at least qualitatively in terms of shape and 368 magnitude. While shapes of the ground reaction forces are an optimization decision, the 369 relative magnitude of these forces are highly constrained. The mass distributions used 370 were taken from empirical vertical forces in the first place (Eq 11-13), and kinematic 371 periodicity ensures that net torque around the body is zero. Furthermore, mean net 372 vertical force must equal body weight. Assuming the body remains relatively horizontal, 373 then these two conditions completely determine the mean vertical force produced by 374 either pair of limbs. While these arguments emphasize that little merit should be seen 375 in matching relative magnitudes (apart from duty factor predictions), they do show that 376 the quadruped's center of mass must be biased towards one pair of limbs to get 377 fore-hind differences in (mean) GRF magnitude. Without bias, there can be no 378 difference in mean vertical forces produced by the fore and hind limbs.

379
While the mean magnitude of GRF is not very informative, the shape is. The 380 optimization predicts a double-hump GRF profile for walking, and single hump for 381 trotting. These predictions are born out in the empirical cases, although the magnitude 382 of the force peaks in the slowest case is substantially lower than predicted. The 383 empirical forces for the hind feet in the intermediate walking case (Fig 4B) are derived 384 using a questionable method from simultaneous contact of fore and hind limbs on a 385 force plate [34], and are not considered reliable.

386
The simulation forces appear realistic for two reasons. First, the force-rate penalty 387 ensures they are smooth, (in the same way averaged or filtered ground reaction forces 388 tend to be). Second, the work cost promotes double-hump walking and single hump 389 trotting GRFs. The double-hump profile is representative of a "pre-heelstrike" 390 trailing-limb pushoff, which redirects velocity prior to contact of the leading limb, 391 avoiding some negative work. The roughly symmetrical, single hump profile is 392 representative of a short, largely vertical contact of the limbs, which are work 393 minimizing in point-mass bipedal models [49]. Our interpretation is in contrast to [16], 394 who say all shape characteristics (double hump, single humps) are from "oscillation 395 modes" due to springs and resonant properties of the system. Despite not having any 396 springs, similar double-hump walking and single hump running GRFs emerge from our 397 energy-minimizing model.

398
While spring-mass systems are often used to approximate the emergent gait 399 behaviour, they are not the only system that describes the behaviour well and should 400 not be taken as a prerequisite for "spring-like" behaviour. Rather, spring-like behaviour 401 emerges from energy minimization [20,21,49,50], even in the absence of springs.

402
It is difficult to tell if differences between solutions in Fig 4 are due to differences in 403 morphological input or prescribed speed. To control for these differences, the model was 404 compared to a large dataset using Belgian Malinois dogs of similar sizes over a large 405 speed range [37]. 406 Comparison to large kinematic dataset 407 The results of the simulations are compared to empirical results in Fig 6. Optimal mean 408 duty factor (averaged across all limbs) matches empirical values very well for U ′ h ≥ 0.5, 409 (Fig 6a) following its slow decay as speed increases. For U ′ h < 0.5, the predicted duty 410 factor is almost constant at approximately 0.6. As with the other test cases (Fig 4), the 411 model predicts singlefoot walking as the optimal gait at low speeds, and trotting as 412 optimal at high speeds (Fig 6b). The phase offset is predicted well for all walking and 413 trotting speeds, with all phase offsets lying within 0. considering that gait transition speeds can be quite variable [37]. While a true gallop 417 never emerges from the optimization, even well past the natural trot-gallop transition, 418 slight asymmetry begins to emerge from the solution for U ′ h > 2 (Fig 6c).

419
Response of optimal Duty Factor to speed 420 Above U ′ h = 0.5, optimal values lie within one standard deviation of the empirical 421 trendline. The decaying behaviour of duty factor is at odds with an impulsive trot (DF 422 = 0), which would be work-optimal. Below about U ′ h = 0.5, the optimal duty factor 423 levels off at about 0.6, slightly above the work-minimizing prediction of 0.5 for a 424 walk [47]. A similar pattern of constant duty factor shifting suddenly to a decay with 425 increasing speed was also observed by Hubel and Usherwood [26] in their biped model 426 with an analogous work + rate penalty objective (though the rate penalty was 427 peak-power in their case, rather than force rate). increase in duty factor above work-minimizing predictions in our model is due to the 432 force rate penalty, which allows both these additional costs in order to avoid impulses. 433 However, as speed increases, enhancing stance duration involves ever-increasing 434 excursion angles, and so larger and larger fore-aft forces, increasing work costs.

435
Therefore, the optimal solution is to decrease duty factor as speed increases to manage 436 this tradeoff.

437
It is interesting that the simple addition of a rate penalty yields such high agreement 438 with duty factor at high speeds, especially since the model does not include massive 439 limbs, springs, or a compliant back, nor does it use the correct (galloping) gait at high 440 speeds. If rate penalties (whether force/time, power or activation penalties) are real 441 phenomena, we would expect their influence to be most pronounced at higher speeds, 442 where, since stride time is shorter, by necessity forces must be produced more quickly, 443 power is higher, and muscles must be activated more frequently.

444
The discrepancy at lower speeds can be partially explained by leg swing dynamics. 445 There is a cost to swinging legs during locomotion [52,53] that is minimized when the 446 leg is swung at its natural frequency [54,55]. As a result, quadrupeds (including dogs) 447 tend to decrease stance duration, rather than swing duration, as speed increases [37,56]. 448 If swing costs dominated, at all speeds the animal would want to adjust its duty 449 factor such that swing time was exactly half the natural swing period, i.e.
where f n is the leg natural frequency and T is stride time. Taking T as given (based on 451 empirical data), and using the allometric equations from [57] for dog leg natural and equal fore and hindlimb lengths [47,58] Previous studies have shown that galloping becomes optimal at high speeds when 483 elasticity is added to the spine [60][61][62]. Future work could add active compliance to our 484 model's trunk, to see whether true elastic return is necessary for the optimality of 485 galloping, or whether any form of compliance makes galloping more efficient.

486
Despite initializing the optimization from over 8000 randomized initial guesses, only 487 singlefoot walking and trotting / pacing (or slight deviations thereof) were discovered as 488 pseudo-global optimal solutions. This is strong evidence (but not definitive proof) that 489 these are globally optimal solutions for the quadrupedal configuration (Fig 1) with a 490 work-based cost and force-rate penalty.

491
Again, these results emerge despite the simplicity of the modelling approach, and the 492 absence of springs. This is further evidence that natural walking, trotting and pacing 493 are energetically optimal even without elastic elements. We are not suggesting that 494 tendons (or springs in robotic quadrupeds) do not potentially reduce energetic cost, 495 only that they are not necessary prerequisites for the optimality of natural gait. A gait 496 could emerge as a recurrent strategy within a species because of its low energetic cost, 497 and well-tuned tendons (such as the digital flexor and interosseus tendons in horses [63]) 498 can then arise through natural selection to complement those strategies.

499
The discovery of trotting (as opposed to tolting) as energetically optimal is 500 consistent with an analysis of Usherwood [64]. A dimensionless pitch moment of inertia 501 can be defined as [65] 502 or the moment of inertia relative to the radius of gyration of the body. To a first 503 approximation, forÎ > 1, tolting is energetically optimal; otherwise trotting is [64]. In 504 our case,Î = 0.93, not far from a more realistic estimate of 0.84 for dogs [34,66]. Since 505 these dimensionless moments of inertia are less than 1, we expect the model (and dogs) 506 to find trotting to be less expensive than a running walk.

507
However, some deviation from exact trotting (with perfect symmetry and fore-hind 508 offset of 0.5) was observed at higher speeds. Moreover, while a clear singlefoot contact 509 pattern was observed for U ′ h < 0.65, and a clear trotting pattern emerged for 510 0.65 < U ′ h < 1.2, at higher speeds the footfall pattern was less consistent. Pair 511 dissociation became more common, but was not always optimal (for example at 512 U ′ h = 1.8 and 2.2) and the optimal gait became asymmetric (but regularly increasing 513 asymmetry was not observed in the pseudoglobally optimal solutions).

514
One reason for the lack of a consistent pattern may be that there are more locally 515 optimal gaits at higher speeds. To explore this possibility, we plot phase offsets for all 516 valid solutions for four speeds from the walking, trotting, slow gallop and fast gallop between HL and FL). These solutions emerged at all speeds, but were never optimal, 546 reflecting the results of [16].

547
Why are solutions more variable at higher speeds? One reason for this may be that 548 phase and event sequence may have less of an impact on cost at higher speeds than at 549 lower speeds. Fig 8A shows how cost changes with increasing speed. Reflecting both 550 human [69] and animal studies [10], the cost of transport in walking is highly sensitive 551 to speed, while the cost of transport of running is less so; however, minimum and 552 median cost always increases with speed. Despite the increased costs in running, the 553 range in CoT is fairly small at higher speeds compared to walking. Relative to the 554 increases in cost, the CoT range of feasible solutions decreases considerably from lowest 555 speed to highest speed (Fig 8B) At high speeds, virtually all combination of phase are 556 explored (Fig 7), with little change in cost (Fig 8). Optimizing locomotion at higher 557 speeds appears to involve searching a "flat and fuzzy" cost landscape for the best gait 558 (Tad McGeer, pers. comm.).

559
This result from the simulation may somewhat reflect reality; dogs exhibit 560 substantial natural variation in gait at higher speeds. For example, empirical data 561 from [39,45] indicate that duty factor and pair dissociation can vary at the same 562 trotting speed in labradors (Fig 4C,D). Moreover, Maes et al. [37] found that Belgian 563 Malinois tend to use trots and rotary gallops at running speeds, but will also pace and 564 use transverse gallops. If the relative cost between these gaits are not substantial, then 565 factors apart from energetics might begin to dictate locomotory decisions.

566
The plots in Fig 8 represent a combination of symmetrical and asymmetrical   567 solutions. Hildebrand [51] studied the symmetrical gaits of domestic dogs in great detail, 568 and discovered that dogs (across a large range of speeds and breeds) choose a limited 569 subset of symmetrical gaits from the total "gait space" available to them (every 570 combination of duty factor and Phase from 0 to 1). Does the limited scope of C. lupus 571 domesticus gait represent a fluke of evolution, a physical constraint, or is it the result of 572 an energy optimization process?

573
In Fig 9 we plot all the symmetrical gaits preformed by dogs and recorded by 574 Hildebrand [51], as well as the "C-shaped" contour that he thought captured the scope 575 of locomotory behaviour most meaningfully. Overlaid on these data are all the feasible 576 convergent results 6 from every condition studied in Fig 6. Color represents relative cost, 577 that is the deviation of the cost for a given solution from the minimal cost at the same 578 speed.

579
The vast majority of symmetrical solutions, as well as the lowest cost, lie within a 580 similar "C-shape" to that observed by Hildebrand. Infrequently, sub-optimal but feasible 581 solutions are discovered outside the range observed by Hildebrand. Their presence 582 indicates that such solutions are possible-so there seems to be no plausible physical 583 constraint to prevent dogs from entering that region of gait space-but their relative 584 high cost indicates that they are unlikely to be chosen from an optimization process.

585
In contrast, the solutions along the C contour are frequently discovered by the 586 optimization process. Not all are global optima for the given condition, but they are 587 local optima; solutions that the optimizer converged on and halted when the underlying 588 cost landscape was sufficiently flat. The remarkable similarity between the simulations 589 and Hildebrand's empirical observations suggests that dogs are also optimizing cost of 590 transport when selecting gait.

591
Notably, Hildebrand's results combined observations from dog breeds of varying sizes 592 and shapes (from Basset hounds to Great Danes [51]), yet the simulations used only one 593 morphological set, based on Belgian Malinois. The consistency between simulations and 594 empirical data, despite the differences in size, is due to dynamic similarity. While a 595 small and large dog may move differently at the same absolute speed, they are moving 596 at different dimensionless speeds (the square root of Froude numbers, Eq 10) and so 597 have different dynamic constraints and opportunities. Yet when their dimensionless 598 speed is similar, their behaviour is similar [1]. As our simulation probed a large range of 599 Froude numbers, we explored the scope of dynamic regimes that dogs (of all sizes) 600 might be expected to experience in level, steady movement.

601
Harder to explain is the similarity between our results and Hildebrand's despite 602 differences in morphological shape. These results seem to suggest that general shape (at 603 least as far as limb length to body length ratios and mass distribution) has little effect 604 on the optimality of trotting, pacing, and single foot walking over other symmetrical 605 gaits. The effect of morphology on the optimality of gait is something that could be 606 explored with the present model by systematically varying morphological inputs.

607
One major discrepancy between the empirical and simulation results is the presence 608 in the latter of a lower "C" that is a reflection of the first, upper "C". In this planar 609 model there is no energetic distinction between a phase of p and p ± 0.5, so it should 610 come as no surprise that for every optimal gait in the upper region, there is an equally 611 optimal gait in the lower region with the same duty factor but with a phase increase of 612 0.5. However, the discrepancy opens the question of why dogs (and most other 613 quadrupedal mammals [7,70]) use symmetrical gaits with pair lag less than 0.6 as 614 opposed to higher pair lags. This is a question that a planar model cannot answer, but 615 may be explained by considering full three-dimensional dynamics.

617
We developed a minimally constrained quadrupedal model capable of any of the 2520 618 contact sequences available to quadrupeds. Of those possibilities, it found two basic 619 gaits-four-beat walking in singlefoot, and two-beat trotting or pacing-as energetically 620 optimal. It also transitioned spontaneously from walking to trotting at a realistic speed. 621 While limb work was the chief contributor to cost, an additional force-rate penalty led to 622 better prediction of duty factor. Despite no enforcement of gait symmetry, the optimal 623 gaits were highly symmetrical except at high speeds, as observed in natural gait.

624
The GRF profiles were double-hump in walking and single-hump in trotting despite 625 indicated by the NLP solver. features result from "oscillation modes of the legged system?s compliance" [16]. Instead, 627 they are the result of work-minimizing strategies; the "pre-heelstrike" pushoff in walking 628 and near-vertical "pseudo-elastic" contact in running [49].

629
Galloping did not emerge as optimal at high speeds in our model, but elastic-limb 630 models exhibit the same result if elastic elements are not included in the trunk. Future 631 work will examine whether galloping emerges as optimal if the torso is allowed to extend 632 or contract actively. This would test whether elastic elements of the spine are truly 633 necessary for the optimality of galloping.

634
Limb work with a small force-rate penalty has substantial explanatory power. affecting duty factor choice at slow walking remain a mystery, we believe that duty 639 factor at higher speed is determined by the competing costs of work (lowering duty 640 factor) and force rate (increasing it). However, little is known about the physiological 641 mechanisms for this rate cost, and it could be an exciting avenue for further research.

642
Singlefoot walking and trotting emerged as optimal gaits despite over 8000 random 643 guesses being used to seed optimization. This strongly suggests (but does not prove) 644 that these gaits are globally optimal strategies at their respective speeds, at least for the 645 modelling configuration used. It is unlikely that alternate gaits exist which would 646 further minimize quadrupedal cost of transport at slow and intermediate speeds.

647
During locomotion, mammals are likely optimizing a work-based cost function with 648 some form of force-rate related energetic penalty. Although elastic elements can lower 649 metabolic cost of transport, they were likely not prerequisites to the evolution of 650 modern gait in mammals, but rather complemented preexisting, metabolically optimal 651 strategies.
And the limb forces are As the body consists of only one rigid link, the dynamics are simplÿ The normalized cost of transport is  with the ground, its shortest possible duty factor is 1/10. If all four limbs are 671 operational, the shortest possible duty factor is 1/40. Dogs are never observed to have 672 DF < 0.2, so this force limitation is well above the empirical limit.

673
If |ẋ|, |ẏ| ≤ 4D ′ , then the COM speed is allowed to be four times the average 674 horizontal translation speed. If |θ| ≤ 4, then the body could conceivably more than fully 675 rotate in one gait cycle. The maximum possible power (represented by p ijk and q ijk ) is 676 for the force to be maximal when speed at either the hips or shoulders is maximal. The 677 speed at either joint could be as high as max(θ) + max(ẋ) 2 + max(ẏ) 2 = 4(1 + 2D ′ ). 678 D ′ was as high as 5, so in principle the instantaneous power could have been as high as 679 323 (much larger than the limit of 20 set).

680
However, increasing the bounds on p and q led to larger solve times, and in practice, 681 the peak power was much lower. Peak power was on occasion saturated, but only for  Similarly, the limit for |Ḟ ijk | was 100, which was much larger than any observed 688 valid solution. Among valid solutions, the largest |Ḟ ijk | was about 52. Other 689 complementarity slack variables (s ijk ) never exceeded 1 E -2 for valid solutions.

690
Footfall positions were bounded as |f ij | < max(l imax ) + D ′ . Geometrically, footfalls 691 at their set bounds would not allow actuation without exceeding limb length constraints, 692 so these limits cannot be exceeded. Left column shows an outline of the dog breed in question, and the points used as morphological measurements. Middle column shows that gait diagrams for the empirical (grey) and optimal solutions (light blue) have substantial overlap (dark blue). Right two columns demonstrate that empirical ground reaction forces (GRF, black lines) qualitatively match optimization results. The optimization GRF from left (solid blue line) and right limbs (dashed red line) are very similar, reflecting the symmetrical solutions that are discovered. At slow speeds (A-B), the simulation discovers a singlefoot walking gait, while at an intermediate speeds (C-D) it discovers trotting, matching natural gait choice in dogs. All optimizations use c ′ 1D = 3 E -3 . Empirical data and morphological measurements from [38,42] (A), [34] (B), [39,41] (C) and [41,45]   These locally optimal solutions are a sample of the diverse gaits possible with the model. Each stick figure is a freeze-frame of a different part of the gait cycle. The first still (leftmost, purple) is at touchdown of the left hind limb. Each still is separated from the next in time by one tenth of the stride time, culminating at the next left-hind touchdown. x marker represents center of mass location. At both a walking speed (left column) and running speed (right column) asymmetric (upper two rows) and symmetric gaits (bottom two rows) are possible. y-axis scale is equal to the x-axis. Pseudoglobal optimal solutions (markers) are compared to a large dataset for Belgian Malinois dogs [37]. (A) A slowly decaying duty factor with increasing speed is predicted well by the model for medium to fast speeds, using c 1D = 3 E -3. At slow speeds, the model remains constant, slightly above impulsive predictions, while the empirical data approaches 1. A swing period at half the natural pendular period (NPP) does not explain the discrepancy. 95% CI for empirical duty factor is shown as a thin dotted line. (B) The optimal phase offset matches empirical values well for walking and trotting speeds, and the optimal walk-trot transition is very close to the natural transition speed. However, the model does not transition to a gallop after the natural trot-gallop transition speed (C) Pair offset (time between contact of the hind and fore limbs) is shown against speed. A perfectly symmetrical gait would have pair offsets of 0.5 in both left and right limbs; walking and trotting in dogs is highly symmetrical both naturally [51] and in the model. Galloping is an asymmetrical gait, and past the trot-gallop transition, the optimal gait becomes asymmetrical. All local optima found for the model for each speed based on a Belgian Malinois morphology. Each row corresponds to a given speed condition, and each data point is a feasible local optimum arising from energy optimization from one random guess, with relative cost (CoT -min CoT for a given U Cost of Transport increases linearly for walking speeds (solid line), but exhibits a sharp change to a slower rate of increase after the walk-trot transition (dashed line), mirroring the response to speed of walking and running observed in human data [68]. The range of the costs of feasible solutions (whiskers) increases up to the walk-trot transition, and then settles into a near-constant range. The distribution of costs is heavily skewed, as indicated by the median value approaching the minimum at all speeds. This indicates that the solvers tended to discover solutions with costs close to the minimal value, but occasionally were "trapped" in local optima with costs far from the minimum. (B) As speed increases, the standard deviation of the distribution of costs of transport gets smaller, relative to the minimal cost, indicating that the variance in costs of local optima is relatively smaller at higher speeds than at lower speeds  , representing mean duty factor and Pair Lag for symmetrical gaits in dogs, are overlaid on all locally optimal symmetrical solutions (coloured circles) discovered by the model for the Belgian Malinois dataset (Fig 6). The empirical contour takes on a "C" shape in the upper (lateral sequence) region of the plot. The optimal solutions take on two "C" shapes, one in the upper half and one in the lower, as lateral and diagonal sequence gaits have equal cost in a planar model. While there is some discrepancy in lower duty factor, there is substantial overlap between Hildebrand's contour and the clustering of locally optimal solutions in the model. Each coloured data point represents one solution from a different random guess at any of the speeds represented in Fig 6. Colour represents cost relative to the minimal cost solution at that given speed.