The validation of new phase-dependent gait stability measures: a modelling approach

Identification of individuals at risk of falling is important when designing fall prevention methods. Current measures that estimate gait stability and robustness appear limited in predicting falls in older adults. Inspired by recent findings on changes in phase-dependent local stability within a gait cycle, we devised several phase-dependent stability measures and tested for their usefulness to predict gait robustness in compass walker models. These measures are closely related to the often-employed maximum finite-time Lyapunov exponent and maximum Floquet multiplier that both assess a system's response to infinitesimal perturbations. As such, they entail linearizing the system, but this is realized in a rotating hypersurface orthogonal to the period-one solution followed by estimating the trajectory-normal divergence rate of the swing phases and the foot strikes. We correlated the measures with gait robustness, i.e. the largest perturbation a walker can handle, in two compass walker models with either point or circular feet to estimate their prediction accuracy. To also test for the dependence of the measures under state space transform, we represented the point feet walker in both Euler–Lagrange and Hamiltonian canonical form. Our simulations revealed that for most of the measures their correlation with gait robustness differs between models and between different state space forms. In particular, the latter may jeopardize many stability measures' predictive capacity for gait robustness. The only exception that consistently displayed strong correlations is the divergence of foot strike. Our results admit challenges of using phase-dependent stability measures as objective means to estimate the risk of falling.

JJ, 0000-0003-2525-0471; JHvD, 0000-0002-7719-5585; AD, 0000-0001-9107-3552; SMB, 0000-0003-0290-2131 Identification of individuals at risk of falling is important when designing fall prevention methods. Current measures that estimate gait stability and robustness appear limited in predicting falls in older adults. Inspired by recent findings on changes in phase-dependent local stability within a gait cycle, we devised several phase-dependent stability measures and tested for their usefulness to predict gait robustness in compass walker models. These measures are closely related to the often-employed maximum finite-time Lyapunov exponent and maximum Floquet multiplier that both assess a system's response to infinitesimal perturbations. As such, they entail linearizing the system, but this is realized in a rotating hypersurface orthogonal to the period-one solution followed by estimating the trajectory-normal divergence rate of the swing phases and the foot strikes. We correlated the measures with gait robustness, i.e. the largest perturbation a walker can handle, in two compass walker models with either point or circular feet to estimate their prediction accuracy. To also test for the dependence of the measures under state space transform, we represented the point feet walker in both Euler-Lagrange and Hamiltonian canonical form. Our simulations revealed that for most of the measures their correlation with gait robustness differs between models and between different state space forms. In particular, the latter may jeopardize many stability measures' predictive capacity for gait robustness. The only exception that consistently displayed strong

Introduction
Falling is a major threat, especially for older adults. About one-third of all adults older than 65 years fall at least once per year, often with serious injuries and fractures as consequences [1]. Without a doubt, there is an urgent need to identify individuals at risk, in particular, when it comes to customizing fall prevention [2]. Various measures have been proposed to quantify the risk of falling in humans. Given the apparent relation to dynamic stability, the local divergence exponent (maximum finite-time Lyapunov exponent) [3,4] and the maximum Floquet multiplier [5] are the first to mention. Lockhart & Liu [6] and Toebes et al. [7] suggested that the local divergence exponent may allow for discriminating fallers from non-fallers, with accuracy of 80% at best [8]. Although 80% accuracy is a fair achievement in view of the complexity of the dynamics involved in human walking, it may still lead to a large number of false positives ( predicted fallers that are not prone to fall) and-arguably worse-false negatives ( predicted non-fallers that are prone to fall) of the elderly population. Therefore, a higher accuracy in fall prediction is of great importance.
Most stability measures are obtained by looking at either stride-to-stride deviations of an isolated point (e.g. heel strike state) or at average deviations over all phases of gait cycles. Humans, however, have phase-dependent gait stability, as we display distinct responses to perturbations when applied at different phases of a gait cycle [9,10]. Estimating gait stability at an isolated point of or from an average over the gait cycle may thus fail. Norris et al. [11] identified the phase-dependency in local stability using a 'simple' walking model, and Ihlen et al. [12] reported intra-stride local stability variations in human gait. These observations and studies ascertained the potential for using phase-dependent local stability information to quantify gait stability.
While the term gait stability may be intuitively sound, it is often confused with gait robustness. Robustness can be operationalized as the maximum magnitude of perturbations that a walker can handle, while stability reveals whether or not a walker will return to its periodic motion after an infinitesimal perturbation. Infinitesimally implies locality (here close to a limit cycle). Stability, or better, local stability, does, in general, not relate to the maximum size of the perturbation that a nonlinear system can handle, i.e. robustness. As such, the relation between a given stability measure and gait robustness is far from straightforward. Take the maximum Floquet multiplier as an example: it is a perfectly valid stability measure for periodic systems [13] but modelling [14][15][16] as well as experimental studies [17] showed its limitations as proxy for gait robustness.
In a limit cycle system like a passive dynamic walker, locally stable phases and unstable phases may coexist. In that case, averaging local stability measures may lead to a loss of information regarding local stability tendencies. Perturbations applied during unstable phases may lead to larger deviations from the limit cycle than applied during stable phases. To address this, phase-dependent stability measures have been introduced by Ali & Menzinger [18]. By and large, they are based on a limit cycle's local stability quantified by the aforementioned responses to infinitesimal perturbations, but applied at different phases of the cycle. When it comes to human gait, phase-dependent stability measures have, in fact, been shown to potentially improve the accuracy of fall prediction over cycleaveraged parameters [19]. Still, only a limited number of studies [12,[19][20][21][22] on human gait used such phase-dependent stability measures, and further validation is needed.
The goal of the current study was to further validate phase-dependent stability measures for walking. Directly testing these stability measures on empirical data is hard, given the enormous complexity and variability of human walking. Therefore, rather than directly testing on empirical data, we asked: Can phase-dependent stability measures provide proper predictions of gait robustness in 'simple' walking models [23]? A positive answer will encourage application to experimental data, while a negative answer should be considered a call for alternative approaches. To answer our question, we devised several phase-dependent stability measures based on the stability analysis of Garcia [24], Goswami et al. [25] and Norris et al. [11] applied to a two-dimensional compass walker model with point feet [11,24,25]. To validate our stability measures, we determined their relations with the walker's gait robustness. We also tested whether these relations were stronger than those between the commonly used local divergence exponent and gait robustness. Since we aimed to test the theoretical merits of phase-dependent stability measures, we here circumvent the inclusion of time-series analysis method.
Hence, we did not apply the method proposed by Ihlen et al. [12,21]. Next, we asked whether the corresponding results transferred to a slightly more complicated model including circular rather than point feet. Finally, we investigated the extent to which our phase-dependent stability measures are invariant to different state space formulations, in particular, when formulating the walker dynamics in either Hamiltonian canonical or Euler-Lagrange form. A dependency of the choice of state space form would severely limit the usefulness of any measure for fall prediction in humans.

Compass walker with point feet
The first model we used is illustrated in figure 1. It is a passive dynamic compass walker with point feet [11,24,25]. In brief, the model consists of two massless legs connecting the hip point mass M and foot point masses m; we denote the mass ratio as β = m/M. The state variables are defined as s ¼ (u,w, _ u, _ w) T , where θ represents the angle of the stance leg with respect to the normal of the inclined plane, w is the angle between the legs, and _ u as well as _ w are the respective angular velocities. A gait cycle (one step) comprised a swing phase and an instantaneous double stance phase. At the double stance moment (= foot strike), we assumed a fully inelastic collision of the leading leg and exploited the conservation of angular momentum that gives rise to a linear mapping from pre-to post-collision state. The model walked down a surface with slope γ. Throughout our paper, we rescaled time by ffiffiffiffiffiffi ffi l=g p to ease legibility. A more detailed description including the Euler-Lagrange equations of motion can be found in appendix A, where we adopted the notation of Garcia [24] and Norris et al. [11].
We also analysed this compass walker model using Hamiltonian form. In that formulation, the state variables were defined ass = (θ, w, p θ , p w ) T , where p θ and p w denote the canonical conjugate momenta with respect to θ and w. The corresponding equations of motion can be found in appendix B. In order to find the general relations between stability measures and gait robustness of the model, we included a variety of configuration parameters by systematically varying the slope γ and mass ratio β (g ¼ 10 À3 . . . 1:2 Á 10 À2 in steps of Dg ¼ 2 Á 10 À4 and b ¼ 2 Á 10 À3 . . . 1 Á 10 À1 in steps of Db ¼ 2 Á 10 À3 ), yielding 56 × 50 = 2800 parameter combinations.

Compass walker with circular feet
We also simulated a two-dimensional compass walker model with circular feet. In this model, the state variables are s ¼ À w 1 ,w 2 , _ w 1 _ w 2 Á T , where the definitions of these angles are given in figure 2; more details of the configurations and equations of motions can be found in Wisse & Schwab [26]. To obtain the equations of motion and collision equations, we followed the approach by Casius et al. [27]. The circular feet walker also has instantaneous foot strike. The foot radius r was varied as r = 0.01… 0.49 in steps of Δr = 0.04 and the slope γ were varied like for the point feet walker model yielding 13 × 56 = 728 parameter combinations.

Numerical simulations
For the point feet walker, we used a Runge-Kutta (4,5) integrator, while for the circular feet walker, we used an integrator for non-stiff differential equations; all simulations were realized in Matlab (The Mathworks, Natick, MA, USA). The absolute tolerance and relative tolerance were both 10 −8 for the point feet walker and 10 −10 for the circular feet walker. To avoid foot scuffing, foot strike detection condition was defined by constraining θ to be less than −0.05 radians for the point feet walker, and the swing leg to be in front of the stance leg and the inter-leg angle to be larger than 5°for the circular feet walker. For every parameter combination, we searched for a period-one solution (i.e. a solution where the states at the end of the step are the same as those at the beginning of the step) using the Newton-Raphson method [26]. Whenever a period-one solution was detected, we determined the Floquet multipliers from the Poincaré map reflecting error multiplication factors from the step-to-step map following [26]. We ignored all the unstable solutions (and thus all the short-period gaits; [28]) in further calculations by removing all the solutions with a maximum Floquet multiplier with a modulus greater than 1.

. Gait robustness
In order to assess the robustness of the models, we determined the maximal step-up and step-down perturbation the two models could handle for every parameter combination. To this aim, the walker started from its stable period-one solution after which, in the first step, we applied a one-time step-up/ step-down floor height difference, see figure 3 for an illustration of these perturbations. To generalize gait robustness in the light of another type of perturbation, we also applied a one-time constant push or pull at the centre of mass of each leg. We performed this type of perturbation only for the circular feet walker. The push or pull was applied during the first step for 0.1 s, starting at the moment the hip angle was zero. Gait robustness was then quantified as the sum of the maximum allowable push and pull perturbation or the sum of the maximum step-up and step-down perturbation without falling. We gradually increased this perturbation size (with precision up to 10 −5 for the point feet walker, 10 −4 for the circular feet walker with step variation perturbation and 10 −2 with the push and pull perturbation) until, after the perturbation, the model was no longer able to complete 30 steps (which were empirically enough for the disturbance to be attenuated if the walker would not fall earlier).

The local divergence exponent
To determine the widely used local divergence exponent or short-term maximum finite-time Lyapunov exponent, we employed Rosenstein's algorithm [29]. In brief, per parameter combination, we performed five simulations of 200 steps where at each step, a floor height variation was added (zerocentred Gaussian white noise with a standard deviation of 2 · 10 −5 ). The time series of full states of the model were the direct input of the Rosenstein's algorithm for constructing the corresponding state   . Circular feet compass walker and its configurations. A snapshot is captured at the end of a step. For simulation purposes, gravity was tilted to simulate an inclined slope γ. w 1 and w 2 are the angles of the swing leg and stance leg with respect to the inclined slope; r is the radius of the foot.
royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 201122 space. Then, we applied Rosenstein's algorithm [29,30] to estimate the local divergence exponent for each of these simulations, and averaged over five simulations per configuration.

Swing phase
To quantify the local stability of the swing phase of our models, we assessed how their continuous dynamics, _ s ¼ f(s), behave in the presence of an infinitesimal perturbation δ(t). This behaviour can be approximated by linearization using the (time-dependent 1 ) Jacobian at every point along the periodone solution ð2:1Þ Note that J(t) is not explicitly dependent on time (both walkers are time-invariant systems), but implicitly through δ(t). The eigenvalues of J(t) quantify the rate at which the system returns to (negative eigenvalue) or moves away from ( positive eigenvalue) the period-one solution in the corresponding eigen-directions after infinitesimal perturbations. Here, we would like to add that these eigen-directions may be ( partly) along the period-one solution. In that case, the eigenvalues belonging to these eigenvectors quantify how perturbations yield a phase shift. Since such phase shifts in our model do not result in moving away from the period-one solution, they do not affect the stability of the period-one solution as a whole, and hence do not provide useful information about gait stability or gait robustness. One can circumvent this by performing a coordinate transformation to obtain eigenvalues that belong to eigenvectors that are orthogonal to the period-one solution.
Following Ali & Menzinger [18] and Norris et al. [11], we used a moving coordinate frame U(t), with one axis always remaining tangent to the period-one solution f(s(t)) T Á k f(s(t))k À1 , while the other axes span a hyperplane orthogonal to the period-one solution; cf. Figure 4. One can use this moving coordinate frame to transform perturbations from the global frame δ(t) to the moving coordinate Although local stability is defined in terms of time, for a period-one solution, this directly translates to phase.
royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 201122 By taking the derivative of equation (2.2) with respect to time and using equation (2.1), one can obtain in which J(t) describes how perturbations evolve. However, due to our coordinate transform, the perturbations δ 0 (t) in the tangent and orthogonal directions are uncoupled. That is, the J(t) matrix obeys the form where vanishing (zero-valued) l ? i k (t) imply that a tangent initial perturbation does not evolve onto the orthogonal hyperplane [11]. Given that phase shifts do not alter stability in our model, one creates a reduced Jacobian matrix J 0 (t) by removing the first column and top row of matrix J(t). The eigenvalues of J 0 (t) describe the rate at which infinitesimal perturbations return to (negative eigenvalue) or move away ( positive eigenvalue) from the period-one solution in the corresponding eigen-directions, which are all orthogonal to the direction of the period-one solution, i.e. they do not quantify a phase shift. Figure 5 illustrates the evolution of all eigenvalues of the reduced Jacobian J 0 (t) during the swing phase for the point feet walker with an arbitrary γ and β. These eigenvalues indicate phase-dependent local stability orthogonal to the period-one solution, and one may derive potentially useful measures from them. To do so, one can calculate the trajectory-normal divergence rate, Γ(t), as the sum of all eigenvalues of J 0 (t), indicating the mean contraction or expansion rate of all neighbouring perturbations rather than along a single eigendirection (figure 5, blue curve).
The advantages of using the trajectory-normal divergence rate instead of the maximum eigenvalue of the local Jacobian matrix to quantify local stability are threefold: (i) it quantifies the growth of all perturbations around the unperturbed trajectory rather than in a single eigendirection; (ii) it removes dimension that corresponds to the tangential; (iii) it is invariant against rotation of the normal coordinates since the growth of perturbations is evaluated only in the orthogonal space that is always normal to the time-varying tangential. From this trajectory-normal divergence rate, we calculated the following three measures.
Stability measure i. The maximum of trajectory-normal divergence rate during the swing phase, i.e. max t[Swing G(t), quantifies the point with the largest trajectory-normal divergent rate in the swing phase.
The higher the value, the more deviation an infinitesimal perturbation at this phase will yield from  . State space and moving coordinate frame. The period-one solution is displayed as the blue curve with a jump occurring at foot strike. We expressed our stability measures in the orthogonal hyperplane since we are only interested in the eigenvalues that belong to the eigenvectors that are orthogonal to the trajectory (i.e. the green arrows lying in the green plane), and can ignore eigenvalues that belong to eigenvectors that are along the period-one solution (i.e. red arrow perpendicular to the green plane). The blue straight line is due to the switch between swing leg and stance leg at foot strike.
royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 201122 the orbit. Since it is a nonlinear system, it remains to be shown if for the walker models analysed here, a higher value of this measure would result in a lower gait robustness. Stability measure ii. The minimum of the trajectory-normal divergence rate during the swing phase or min t[Swing G(t) quantifies the point with smallest trajectory-normal divergent rate in the swing phase. The higher this value, the more deviation an infinitesimal perturbation at this phase will yield from the orbit; cf. measure i above.
Stability measure iii. The integral over the trajectory-normal divergence rate over the swing phase, Ð t[Swing G(t)dt, provides the long-time multiplication factor of the volume of all neighbouring infinitesimal perturbations over the swing phase. Similar to the previous measures, the higher this value, the more the infinitesimal perturbations will cause on average a deviation from the orbit. Note that this measure was also studied in the physics literature where it is known as the trajectory-normal repulsion rate ρ T ; it indicates the normal growth of infinitesimal normal perturbations to the trajectory over the time interval [0, T ] [31,32].

Foot strike and gait cycle
A complete gait cycle (one step) of the compass walker contains a swing phase and an instantaneous double stance event (foot strike). Next to quantifying the local stability of the swing phase of the period-one solution of our models, we also quantified the local stability of the foot strike. We did this in a similar way as for the continuous dynamics. That is, we determined the Jacobian J FS of the foot strike event. Its eigenvalues quantify the degree of deviation a perturbation may cause during the discrete event (i.e. eigenvalue modulus greater than 1 indicates instability 2 ). Similar to the continuous case, where we eliminated the phase shift dimension before estimating the trajectory-normal divergence rate, we finally derived the divergence as product of all non-vanishing eigenvalues. Since the foot strike event can be considered an important phase in the gait cycle, we also considered its divergence as a phase-dependent stability measure.
Stability measure iv. The divergence of the reduced Jacobian of the foot strike, div J 0 FS , quantifies the multiplication factor of the volume of all infinitesimal perturbations surrounding the period-one solution before and after the foot strike event. The closer this measure is to zero, the more likely it is that the perturbations after foot strike is attenuated.

. Floquet multipliers
Stability measure v. The divergence of a full period-one solution can be obtained by taking the absolute value of the product of all non-vanishing Floquet multipliers. Although it is known that the maximum Floquet multiplier is not an optimal proxy for gait robustness [17], the divergence of a gait cycle could be an alternative for further validation. 3

Results
In order to evaluate how well these phase-dependent stability measures and the local divergence exponent predict gait robustness, we correlated them with gait robustness (= maximum perturbation the model could handle), given our sets of parameter combinations. Before selecting our correlation measure, the general rubric for a 'good correlation' maintains that (i) for a given stability measure value, there should be a small variation in the corresponding gait robustness; (ii) a monotonic increasing stability measure should be able to predict a monotonic increasing (positive correlation) or decreasing (negative correlation) gait robustness. Two simple examples of a 'good' and a 'bad' correlation are illustrated in figure 6b,e. Instead of using Pearson correlation coefficient which cannot quantify nonlinear associations, we used Kendall rank correlation as summarized in table 1. This nonparametric correlation coefficient (between −1 and 1) evaluates the similarities in the ordering of the data when ranked by two variables. A coefficient value of 1 (−1) satisfies the monotonic relation that when one variable increases, the other variable increases (decreases); a coefficient value of 0 implies the absence of any association between the two variables. We considered the Kendall rank correlation coefficients larger than 0.7 (or smaller than −0.7) to be indicative for strong correlations. For truly promising stability measures predicting the robustness of a human gait, however, we expect the correlation coefficients to be larger than 0.9 (or smaller than −0.9). For the point feet walker, two measures (measure iii in figure 6d and measure iv in figure 6e) showed good correlations with gait robustness (figure 6), better than those of the local divergence exponent, suggesting that these phase-dependent stability measures might be useful for predicting gait robustness. In particular, the divergence of swing phase (measure iii, the integral of the trajectorynormal divergence rate, figure 6e) and the divergence of foot strike (measure iv, figure 6f ) showed very good correlations.
Unlike for the point feet walker, for the circular feet walker model, the local divergence exponent and most of the phase-dependent stability measures were not correlated well with gait robustness (measure i, iii and v in figures 7 and 8). One exception was the minimal trajectory-normal divergence rate (measure ii Table 1. Kendall rank correlation coefficients of stability measures with gait robustness for the different models, different perturbation types and different state space forms. Green shading indicates negative correlations (higher stability measure value corresponding to lower robustness, as theoretically expected) and grey shading indicates positive correlations. We considered the Kendall rank correlation coefficients larger than 0.7 (or smaller than −0.7) to be strong correlations and displayed them in bold font in the table. All of these strong correlations are statistically significant ( p < 0.001).   688 (figure 8). Interestingly, applying different types of perturbations to quantify gait robustness (maximum step-up/step-down perturbation in figure 7 and maximum push/pull perturbation in figure 8) did not alter the general relations between stability measures and gait robustness, at least not qualitatively. The influence of state space choice (here the Euler-Lagrange versus Hamiltonian form of the point feet walker) on the correlations between our stability measures and gait robustness is shown in figure 9. Note that the effect of (maximal) perturbations is independent of the state space choice. Furthermore, note that the figure 9e,f is necessarily the same for these two different coordinate systems. As can be readily seen from the figure, for the majority of measures, there was a clear dependency on the choice of state space form-see also appendix C for more mathematical details. This means that the choice of the state space form (e.g. whether we choose angular velocity or angular momentum as our state variable) will affect the values of phase-dependent stability measures and how well they correlate to gait robustness.
All in all, our results showed that (i) comparing figures 6 and 7, the correlations between phasedependent stability measures and gait robustness highly depend on the walker model we considered; except for the divergence of foot strike (measure iv, figures 6e and 7e), all other measures show different correlations; (ii) comparing figures 7 and 8, robustness quantified by different types of perturbations (step-up/step-down versus push/pull) does not change the general relations between stability measures and gait robustness; (iii) comparing different state space forms in figure 9, the choice of state space form affects the values of phase-dependent stability measures and their correlations with gait robustness.

Discussion
We explored several phase-dependent stability measures, to assess if they may adequately and consistently predict robustness of simple walking models. To do so, we devised five phase-dependent stability measures of compass walker models and investigated their correlations with gait robustness. We also calculated the commonly used local divergence exponent to test whether it may be outperformed by phase-dependent measures. Our analysis revealed that some phase-dependent stability measures could predict gait robustness better than the local divergence exponent for the simple point feet walker model, but failed to predict gait robustness for the circular feet walker. It is worth mentioning that in Bruijn et al. [33], the local divergence exponent was found to be well correlated with gait robustness of an arced feet walker with hip spring, but there the manipulated parameters were different (hip spring stiffness and foot radius) and much fewer configurations were tested in that study. Therefore, the correlations between stability measures and gait robustness also depend on which configuration parameters are varied.
We considered two different types of perturbations that took place at different phases during the gait cycle. We used a step-up/step-down positional perturbation that occurred at foot strike while push/pull and a force ( pulse) perturbation that occurred around mid-stance phase. Even though the nature and timing of these two perturbations are completely different, there seems to be no qualitative difference in the correlations between stability measures and gait robustness quantified by these two different perturbations (figures 7 and 8). This may be explained by the fact that a step-down perturbation results in a larger stance leg velocity at the subsequent step [34], which is also the result of the push perturbations, and hence, the perturbation types may share similarities in their ultimate effects. In line with this, Chen et al. [35] applied acceleration disturbances at five different phases of a gait cycle for a planar-feet walker model. They found that the Pearson correlation coefficient of the maximum allowable perturbations added at different phases was higher than 0.65 for all pairs and greater than 0.8 for most pairs, suggesting that the measured maximum allowable perturbation (gait robustness) to be consistent along the gait cycle for simple walking models. Moreover, Bruijn et al. [33] reported similar correlations between local divergence exponent and gait robustness for a step-down perturbation and push/pull perturbation. Hence, when studying relationships between stability measures and robustness, using one robustness measure may suffice.
Our results showed that the phase-dependent stability measures depend on the state space form. This might be surprising, because the Lyapunov exponents are invariant to smooth coordinate transformation [4]. However, we note that all phase-dependent stability measures used here are subspace description of the original state space. They depend on the tangent vectors, which are not necessarily objective but may be coordinate-and norm-dependent. This may in fact be one of the reasons that different values royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 201122 of phase-dependent stability measures are obtained when using different state variables or different scales of them. Similar arguments have been put forward by Nave et al. [32] and Sternad et al. [36]. Interestingly, we also found the local divergence rate (divergence rate without eliminating the phase shift dimension), which is not a subspace description, to be dependent on the state space form. We analysed this in appendix C.1 and found that due to the fact that the mass matrix is state-dependent, the state space transformation from angular velocities to Hamiltonian momenta is thus nonlinear. While we share the appeal in finding 'objective' stability measures that are state space form invariant (more precisely, the measures that remained unchanged under any translation or rotation of reference frame), we must admit that formulating them is far beyond the scope of the current paper. Yet, we dare to suggest an alternative, namely to find the 'best' state space representation, for instance, using Hamiltonian canonical form or action-angle coordinates [37]. However, this requires further studies that we leave for future work.
We considered most of the phase-dependent stability measures that are based on conventional local stability to have limited predictive value about gait robustness. Local stability quantifies a system's linear response under the proviso of infinitesimal perturbations. It, hence, does not necessarily reveal information about responses to larger perturbations. Nevertheless, predictions of gait robustness may be improved by looking at local stability throughout a gait cycle and calculating (novel) phase-dependent stability measures more appropriately. Interestingly, Ross et al. [38] probed the 'stability frontier' (i.e. the border between unstable and stable regions in state space) in biomechanical time series using local stability measures. Unfortunately, the outcome of this approach would, for our model, consist of high-dimensional dynamical boundaries that we consider difficult to interpret and requires further research. Moreover, their use for estimating the risk of falling remains to be validated.
Only the divergence of foot strike appeared to be reasonably well correlated with gait robustness of our two walker models. Both of the models have an instantaneous double stance phase. While this double stance phase is hardly comparable to that in humans, it is interesting to note that Ihlen et al. [12] found that during heel strike, older adults had impaired gait stability compared to younger adults. Also, older fallers have been reported to have higher gait instability and a higher phasedependent entropy around 0 and 60% (heel strike and toe-off) of the gait cycle [19,21]. These findings do suggest that the stability of the foot strike could be of great importance in estimating/predicting gait robustness. A reason for that might be that (in our models) a dimension reduction occurs during foot strike, which constrains all post-impact perturbations in the state space into a two-dimensional plane. Apparently, a stability measure derived in that plane can contain valuable information about a model's robustness.
The divergence of foot strike does seem to be a promising measure, in particular when considering its consistent relation with gait robustness irrespective of state space form choices and numerical settings. Yet, note that we calculated this measure directly from the analytical expression of the walker model. When working with empirical data, a different approach is needed. For instance, one could use a method similar to Ihlen et al. [12] who calculated phase-dependent stability measures from experimental human data. Another interesting alternative, introduced by Nave et al. [32], is to calculate the trajectory divergence rate from vector fields obtained from experimental data. In this context, it is worth mentioning other potentially relevant variations of the trajectory-normal divergence rate, namely the normal infinitesimal Lyapunov exponent [39] that quantifies the instantaneous growth of perturbations transverse to the trajectory of the linearized dynamics, and finite-time Lyapunov exponents in the instantaneous limit [40], but both to our best knowledge have not yet been applied to time-series data. Finally, we would like to remark that there are of course also methods to calculate the entire Lyapunov spectrum in the presence of discontinuities, see, for instance, Müller [41], but supplemental numerical algorithms are far and few between and certainly not ready-for-use on empirical data. Moreover, since foot strike in human walking is not a discontinuity, it is still an open question whether such methods specifically designed for discontinuities should eventually be used.

Limitations of the current study
The phase-dependent stability measures used here are based on the idea that perturbations tangent to the period-one solution do not affect stability of the gait and will only result in a phase shift along the unperturbed trajectory. This idea has widely been used in controller design in engineering, for example, transverse linearization [42] and hybrid zero dynamics [43]. However, in human locomotion, royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 201122 phase shift may play an important role, indeed. Stumbling perturbations applied at different phases can result in a phase reset to improve dynamic stability of gait [44]. Therefore, we also investigated the local divergence rate of the Jacobian J(t) including both tangent and transverse local stability properties and found that the corresponding phase-dependent stability measures also did not strongly correlate with gait robustness (data not presented here but available via https://datadryad.org/stash/share/ sXS8kl3SiaW1yteinh7-Tjuq5UP5iU5-pIu0lTbNGl8). Furthermore, we looked at the time-varying eigenvalue belonging to the eigenvector that is tangent to the period-one solution and found it to correlate only poorly with gait robustness.
Our phase-dependent stability measures have been derived based on coordinate transformation of the conventional local stability of time-invariant systems. Here, we would like to note that for timevarying systems, in general, the use of the Jacobian's eigenvalues to quantify local stability can be limited, because they may not account for the explicit time-dependent changes of the dynamics. However, our coordinate transformation can compensate for such time dependencies. In appendix D, we included an example of a two-dimensional time-varying system [45] illustrating this. Yet, a more formal mathematical analysis showing that coordinate transformation can improve local stability analysis is pending.

Conclusion
Phase-dependent stability measures have been used to characterize stability changes in compass walker models and in human gait. However, our simulations revealed that they do not correlate with gait robustness in a consistent manner. In particular, they may deviate strongly for different walker models and-arguably more important-they change with the choice of state space form. In our analysis, we encountered only a single phase-dependent stability measure, the divergence of the (instantaneous) double stance phase that displayed relatively good and coherent correlations with gait robustness. This almost linear relation, however, appears model dependent. Overall, it seems that even in simple compass walker models, phase-dependent stability measures and gait robustness are not one-to-one related, at least not the ones evaluated in the current study, which poses challenges for applying these measures to empirical data.

Appendix D. A time-varying system example
Example: (adapted from Leonov & Kuznetsov [45]) Consider the linear system of non-autonomous differential equations with periodic coefficients A(t) ¼ 1 À 4cos 2 (2t) 2þ 2 sin (4t) À2 þ 2 sin (4t) 1 À 4sin 2 (2t) : ðD 2Þ A(t) defines a defective system with degenerated eigenvalues α 1 = α 2 = −1. That is, although A(t), its real-valued eigenvalues do not depend on time and are negative, which indicates local stability everywhere.  royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 201122 On the other hand, the function x(t) ¼ e t sin (2t) e t cos (2t) ðD 3Þ is an explicit solution of (D 1). One can readily see that the solution (D 3) is unbounded, i.e. lim t!1 x 2 (t) ¼ lim t!1 e 2t ! 1, which implies global instability of the dynamics (D 1). And, since (D 1) is linear in x, it also suggests local instability of (D 3). This example readily shows that a 'time-frozen' linearization of the Jacobian matrix may not suffice to determine a system's or its solution's stability, in particular, if systems, e.g. their coefficients, change as a function of time. By contrast, due to our coordinate transform, our approach does account for such time dependencies (through the time derivative of the rotation matrix; see equation (D 3) in our Methods section) rendering the transformed Jacobian an arguably better starting point for identifying local instability. To illustrate this, we simulated several trajectories of (D 1) over a time interval of t = 0…20, with five different initial conditions at t = 0 The corresponding trajectories are shown in figure 12 (only a small range of state values are shown). Subsequently, we applied the coordinate transform, which, interestingly, led to a new, time-invariant Jacobian matrix, J ¼ À2:2 À1:6 À1:6 0:2 ! , with two different, real-valued eigenvalues a 1 ¼À3 and a 2 ¼ 1.
The latter is positive, suggesting local instability of the corresponding solutions, i.e. this approach can identify local instability of the time-varying system. Of course, this statement remains a mere conjecture but we consider a rigorous proof beyond the scope of our current paper. royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 201122