An intermittent control model of flexible human gait using a stable manifold of saddle-type unstable limit cycle dynamics

Stability of human gait is the ability to maintain upright posture during walking against external perturbations. It is a complex process determined by a number of cross-related factors, including gait trajectory, joint impedance and neural control strategies. Here, we consider a control strategy that can achieve stable steady-state periodic gait while maintaining joint flexibility with the lowest possible joint impedance. To this end, we carried out a simulation study of a heel-toe footed biped model with hip, knee and ankle joints and a heavy head-arms-trunk element, working in the sagittal plane. For simplicity, the model assumes a periodic desired joint angle trajectory and joint torques generated by a set of feed-forward and proportional-derivative feedback controllers, whereby the joint impedance is parametrized by the feedback gains. We could show that a desired steady-state gait accompanied by the desired joint angle trajectory can be established as a stable limit cycle (LC) for the feedback controller with an appropriate set of large feedback gains. Moreover, as the feedback gains are decreased for lowering the joint stiffness, stability of the LC is lost only in a few dimensions, while leaving the remaining large number of dimensions quite stable: this means that the LC becomes saddle-type, with a low-dimensional unstable manifold and a high-dimensional stable manifold. Remarkably, the unstable manifold remains of low dimensionality even when the feedback gains are decreased far below the instability point. We then developed an intermittent neural feedback controller that is activated only for short periods of time at an optimal phase of each gait stride. We characterized the robustness of this design by showing that it can better stabilize the unstable LC with small feedback gains, leading to a flexible gait, and in particular we demonstrated that such an intermittent controller performs better if it drives the state point to the stable manifold, rather than directly to the LC. The proposed intermittent control strategy might have a high affinity for the inverted pendulum analogy of biped gait, providing a dynamic view of how the step-to-step transition from one pendular stance to the next can be achieved stably in a robust manner by a well-timed neural intervention that exploits the stable modes embedded in the unstable dynamics.


A.2 The model of ground reaction forces
The ground reaction force is modeled using nonlinear dampers and springs. By defining the sagittal positions of left-toe, left-heel, right-toe and right-heel, respectively, as (X l,t ,Y l,t ), (X l,h ,Y l,h ), (X r,t ,Y r,t ), (X r,h ,Y r,h ), the vertical ground reaction forces acting at these four points in this order are defined as follows.
The horizontal ground reaction forces acting at the four points in the order of listed above are defined as follows.
The parameters used in the model of ground reaction forces are summarized in the following Table. B Numerical evaluation of Jacobian matrix and Floquet multipliers The state space representation of the biped model as a non-autonomous dynamical system with the vector field f (x,t) is defined by Eq. 10 in the main text, and its the Jacobian D ϕ 0 (t) for its linearized equation Ground reaction force parameter 0.3 Table 1: Parameter values used in the model of ground reaction forces is defined by Eq. 14. The numerical evaluation of D ϕ 0 (t) was performed by numerical partial derivative of f (x,t) using the double side finite difference method. That is, the i-j element of D ϕ 0 (t) was obtained as ···,19, and ∆x j is the difference of the i-th element of the state vector x. Throughout this study, we decided to use ∆x j =10 −3 . The size of ∆x j should be determined with a care, because it should balance with the time step ∆t, which was 10 −5 in this study. Our choice of ∆x j was based on the fact that changes in x r (t) for the short duration of time ∆t along the one gait cycle is about between 10 −3 and 10 −6 . Hence the value of ∆x j was the lower band of this variation. We validated the use of ∆x j =10 −3 by examining that the numerical evaluation of D ϕ 0 (t) using ∆x j = 10 −3 and loci of FMs as the function of the PD-gains, as in Fig.4 of the main text, for various values of ∆x j ranging from 10 −1 to 10 −6 . The result of this examination showed that was the loci of FMs were qualitatively and quantitatively the same quite robust for a wide range of ∆x j between 10 −2 and 5×10 −6 , and ∆x j =10 −3 is the middle of this valid range.

C Impedance and dynamic impedance
Here we summarize several definitions of dynamic impedance (stiffness and viscosity). Remind the biped motion equation as J(q)q+B(q,ω)+K(q)+G(q,ω)=U ff (q(t),q(t),q(t))+U fb (q,ω,q(t),q(t)) (1) where with P =diag{0,0,0,P a ,P k ,P h ,P a ,P k ,P h } and denote the total joint torque as U =U ff +U fb . The joint stiffness K d and viscosity B d are usually defined, respectively, by the derivatives of the total joint torque with respect to the position and the velocity, which are equal to P and D in our biped model. That is, Thus, we considered simply the PD-gains of the feedback controller as the joint impedance in this study. It is also worthwhile to consider a different type of joint impedance, we call it total impedance, which is more directly related to stability of the steady state solutionq,q as the limit cycle. Considering a perturbed solution as q=q+q, ω=q+ω,q=q+q. We have and We call K total and B total total stiffness and total dynamic viscosity, respectively. This can also be interpreted as the dynamic balance on each timing on the limit cycle.
In Eq.9, if the perturbation intends to cause the deviation away from the limit cycle, K total and B total will counteract the diverging torque and drive the trajectory back to the limit cycle. So it is also natural to define K total and B total as joint stiffness and viscosity.
Jacobian matrix could be obtained by differentiating the vector field of Eq.10 as follows: For the Jacobian matrix evaluated around the limit cycle, which is denoted by D ϕ 0 in the main text, from the comparison, we can easily see that In the calculation of Floquet matrix, the D ϕ 0 (t) has already been calculated numerically. So it is very prompt to obtain the dynamic stiffness and dynamic viscosity by multiply −J to ∂F 2 /∂q and ∂F 2 /∂ω block of D ϕ 0 . If we only care about the leg joints and are not interested in the correlated influence between the joints, we select only the diagonals of K total and B total , and consider them as dynamic stiffness and dynamic viscosity of leg joint. The non-diagonal element could be interpreted as inter-joint dynamic impedance. Dynamic stiffness and dynamic viscosity of leg joint during steady-state gait with the feedback PD gains of 1500 and 10 are illustrated in Fig.1 and Fig.2. It is noted that the peaks are due to the foot impact.