Phase resetting and intermittent control at critical edge of stability as major mechanisms of fractality in human gait cycle variability
Chunjiang Fu, Yasuyuki Suzuki, Pietro Morasso, Taishin Nomura
aa r X i v : . [ n li n . AO ] O c t Phase resetting and intermittent control at critical edge of stabilityas major mechanisms of fractality in human gait cycle variability
Chunjiang Fu a , Yasuyuki Suzuki a , Pietro Morasso b and Taishin Nomura ∗ a a Graduate School of Engineering Science, Osaka University, Osaka 5608531, Japan b Center for Human Technologies, Istituto Italiano di Tecnologia, 16152 Genoa, Italy
Abstract
The fractality of human gait, namely, the long-range correlation that characterizes scale-free fluctuations of gait descriptors, such as the stride intervals during steady-state walking,depends on the well-tuned organization of the sensorimotor controller. Gait fractality is ap-parent in healthy young adults but tends to disappear in elderly individuals and neurologicalpatients. Therefore, its partial loss may be indicative of pathological conditions. Despite itspotential to be used as a dynamical biomarker for fall risk assessment, the mechanistic ori-gin of gait fractality has been investigated by only a few studies, and even less attention hasbeen devoted to the link between gait fractality and gait stability. Here, we propose a novelcomputational model of gait by addressing the flexibility-stability trade-off first, and then byshowing that gait fractality is a natural consequence of the developed control mechanisms,including phase resetting and intermittent control, which supplement instability in a linearfeedback controller operated at stability’s edge. The results suggest that pathological gait,characterized by joint-rigidity and/or loss of fractality, may be caused by dysfunction in someof these mechanisms.
Gait cycle fluctuation during steady-state human walking, examined typically by a variation ina series of stride intervals between time instants of consecutive ground-contact events for either foot,exhibits a long-range correlation that decays in a power-law manner [1]. This variability is knownas gait fractality , described with the astonishing fact that a given interval is influenced statisticallyby the long-term history over thousands of past intervals, as in a class of scale-free physiologicalphenomena including beat-to-beat heart rate variability [2, 3] and behavioral activity [4]. Gaitfractality is characterized by the rate of power-law decay, referred to as the scaling exponent , whichinvolves critical information for quantifying pathological and age-related alterations in the gaitcontrol system [5, 6, 7, 8]: the exponent is close to unity, corresponding to 1 /f fluctuation inhealthy young adults, and it tends to be 1/2, corresponding to uncorrelated white noise in elderlyand patients with neurological diseases including Huntington’s and Parkinson’s diseases [9, 10].Moreover, the exponent could discriminate fallers and non-fallers among patients with high-levelgait disorder [11].Despite the high expectations of using the scaling exponent of gait fractality as a dynamicalbiomarker for the risk of falls [1], only a few studies have investigated the underlying sensorimotormechanisms that may be responsible for fractal gait [6, 12, 13, 14]. Moreover, there are even lessstudies that relate gait fractality directly to gait stability and the risk of falls [15]. Considering thetwo decades of growing researches suggesting the usefulness of gait fractality, a theoretical accountis keenly needed for gait fractality to become an established principle, in addition to being anacknowledged index in human gait research.A key concept here for understanding the deep rationale of the mechanisms leading to fractalityof behavior is to take into account the general principle that too little variability is not necessarilyfavorable for physiological systems, as they could imply rigidity of the system, whereas too muchvariability may indicate instability of the system [16, 17]. That is, the human motor controlsystem might be designed to avoid rigidity and instability, while aiming to achieve flexibility ∗ Corresponding author: [email protected] the intermittent control that exploits stable modes of the saddle-type unstableequilibrium of the passive dynamics of the body without active feedback control: the unstableupright equilibrium for standing and the unstable limit cycle oscillation for walking.In this study, we focuses on the gait model considering a multi-link rigid body system that is ac-tuated by a feedforward (FF) controller and a conventional time-continuous small-gain proportional-derivative (PD) feedback controller for tracking a prescribed periodic desired joint angle profile,with the addition of well-timed brief activations of time-discontinuous feedback controllers. Inparticular, two types of discontinuous feedback controllers have been proposed previously. One is an intermittent controller (INT), as mentioned above, that brings the state vector of the walkingbody close to the stable manifold of the unstable limit cycle to be stabilized during brief “on-periods” of the controller, and then leaves the state point to flow along the stable manifold towardthe unstable limit cycle during “off-periods” of the controller [25]. The other is a phase resettingcontroller (PR) that compensates timing of fluctuations of foot-contact-events due to endogenousand/or exogenous perturbations, in which the phase of the desired periodic joint-angle trajectoryis advanced or delayed impulsively in response to an occurrence of every early or late event, re-spectively [26, 23, 24]. As we show in the results section, by combining the three types of feedbackcontrollers (PD, INT, and PR), we can achieve at the same time flexibility and stability. Jointflexibility in the model is obtained by operating with small gains of the continuous PD controllerthat maintain joint impedance at appropriate low levels. By tuning the PD controller in such amanner, we would induce instability of the gait, but this danger is avoided in a robust way by thesubtle action of the two time-discontinuous controllers (INT and PR) that are activated brieflyonly at specific timings of the gait cycle: at double-stance phases for the intermittent controller[25] and immediately after heel-contact-events for the phase resetting controller [23, 24].A deeper motivation of the study was the expection that the mechanisms responsible for re-solving the flexibility-stability trade-off might lead naturally to gait fractality, since joint flexibilityis a basis of variability, provided that gait stability can be achieved even with such small jointimpedance. Based on this working hypothesis, we performed dynamic simulations of the gaitmodel with endogenous torque noise. In order to elucidate possible mechanisms necessary for thegait fractality, we systematically analyzed the relationship between gait stability and fractality inthe model with different combinations of the four types of neural gait controllers defined above,and determined the conditions for the genesis of fractal gait.
Results
We considered four combinations of four types of gait controllers to stabilize a limit cycle solutionrepresenting a periodic gait of the mechanical model with seven rigid-links (
Methods ), along witha set of periodic desired joint angle profiles prescribed for hip, knee, and ankle joints using motion-captured gait data [25, 23]. A FF controller (
Methods ) generating joint torques necessary for steady-state walking and a conventional PD feedback controller (
Methods ) at knee, ankle and hip joints totrack the desired profiles with the feedback gains of P = ( P k , P a , P h ) and D = ( D k , D a , D h ) werealways used for every combination. FF+PD controllers provide a basis of the gait, with whichthe steady-state gait can be attained stably if P -gains are large enough with a certain amount of2 -gains. FF+PD controllers.
Fig. 1A exemplifies a simulated time series of the stride intervalsfor the model with FF+PD controllers, when endogenous noisy torques of white Gaussian wereadditively superposed on the joint torques, exhibiting frequent stride-to-stride alternans betweenincrease and decrease in the intervals. As shown in the magnified time series with noise and asample of transient dynamics with a slow decay for the same model without noise (correspondingto the impulse response of the model), the alternans were an inherent property of the model . Asthe short period alternans with the slow decay in the stride intervals results in the anti-persistentcorrelation, the scaling exponent α estimated as the slope of linear regression for the plot of thedetrended fluctuation analysis (DFA for quantifying scale-free long-range correlations in fractionalBrownian motions; Methods ) was almost zero (Fig. 1A-right).The anti-persistency in gait cycle variability shown in Fig. 1A was common across all valuesof P -gain within the stability regions (the light/dark-gray-colored regions in Fig. 2A-1st-row), asindicated by the corresponding color-coded values of the scaling exponent α in Fig. 2B-1st-rowfor the model with FF+PD. That is, the gait mechanics stabilized by the conventional linear PDfeedback controller could never exhibit the fractal gait regardless of the P -gain values. Note thatthe P -gain used for Fig. 1A is indicated by the arrow-heads in the 1st rows of Figs. 2A and 2B. FF+PD+INT controllers.
The use of the INT controller (
Methods ) in addition to FF+FDbroadened the stability region substantially (Fig. 2A-2nd-row), as shown in [25]. However, it didnot contribute to the emergence of fractal gait (Fig. 1B and Fig. 2B-2nd-row). Indeed, the short-period alternans in the stride intervals was still present as in the model with FF+PD, leaving thecorrelation in the stride intervals anti-persistent with the scaling exponent α close to zero. FF+PD+PR controllers.
The stride interval fluctuation became white-noise-like when thePR controller (
Methods ) was appended to the FF+PD controllers (Fig. 1C), where the determin-istic impulse response was non-oscillatory and the degrees of alternans in the stochastic dynamicswere substantially reduced, leading to an exponent α close to 1/2 (Fig. 1C-right). Moreover, theuse of PR controller broadened the stability region in the P -gain parameter space more effectivelythan INT (Fig. 2A-3rd-row) as expected by our previous studies [23, 24], although the geometryof the limit cycle was altered as the P -gain approaches the marginal area of the extended stabilityregions, which is indicated by the light and dark gray colors of the stability regions in Fig. 2A.That is, the limit cycle in the light gray regions oscillates according to the prescribed referencedynamics (defined as γ in Methods ), whereas the one in the dark gray regions is qualitatively dif-ferent from the prescribed γ , referred to as γ ′ for convenience in this sequel. The emergence of γ ′ is due to destabilization of γ in some dimensions, representing a bifurcation phenomenon, but theoverall gait dynamics remained stable as long as the P -gain is located in the dark gray regions.In this way, the border between the light and gray regions represents the bifurcation curve or thebifurcation set in the P -gain parameter space.In the extended stability regions of the model with FF+PD+PR, the scaling exponent α valuesof about 1/2 were widely distributed as shown by the light-blue and green regions of Fig. 2B-3rd-row. Interestingly, the exponent exhibited large values between 0.8 and 1.0 for specific areas D = (10 , ,
10) Nm · s/rad was used throughout the study. The alternans is associated with the fact that the dominant mode characterized by Floquet multipliers of thelimit cycle of the model with FF+PD is a pair of complex conjugates with an amplitude close to unity and anargument greater than π/ π ) generates a perfect alternans, which was not the casefor this model. The mechanism of short-period alternans in the model with FF+PD+INT is not the same as that in themodel with FF+PD, which can be seen from the non-oscillatory deterministic impulse response. For the P -gainssmaller than those used for the model with FF+PD, the limit cycle of the model with FF+PD (as the model withFF+PD+INT during "off-period" of INT) is unstable, for which the pair of complex conjugate multipliers existedfor the model with large gain FF+PD moves outside the unit circle and located typically near or on the real axisin the right-half plane (see Fig. 4 of [25]), associated with an unstable manifold of the limit cycle, along whichthe model during "off-period" of the INT controller falls away from the limit cycle. The alternans here is due tothe random positioning of the state point on either side of the unstable manifold that is separated by the stablemanifold at every timing of the onset of "on-period" of the INT controller. igure 1: Stride interval time series of the gait model actuated by different combinations of feedback controllers withand without endogenous white Gaussian torque noise ( σ = 0 .
003 Nm) and graphs of DFA. (A) FF+PD controllersfor P = (1300 , , P = (1000 , , P = (1300 , , P = (700 , , P = (300 , , F ( n ) against the scales from 10 . to 10 stride number n for the tenlong-run trials, wherein the thickness of the red curve represents the standard deviations of the ten estimates of the F ( n ). The scaling exponent value of α was estimated from the linear regression (dashed line). F + PDFF + PD + INTFF + PD + PRFF + PD + INT + PR P k P a P h = P h = P h = P h = P h = P h = P h = P h = P h = P h = P h = FF + PDFF + PD + INTFF + PD + PRFF + PD + INT + PR P k P a AB P h = P h = P h = P h = P h = P h = P h = P h = P h = P h = P h = Figure 2: Stability and scaling exponents as the functions of the proportional gains ( P -gains) of the time-continuousPD feedback controller for four combinations of four neural controllers (FF, PD, INT, PR). All D -gains were fixedat 10 Nm · s/rad. (A) Stability regions are indicated in light and dark gray in the P -gain parameter space in eachpanel, which means that the gait model could not walk stably in black regions. The panels aligned in each rowrepresent the P k - P a planes for 11 different P h values of 0, 100, · · · , 1000 Nm/rad. Any stable limit cycle in the lightgray region was close enough to the limit cycle γ defined by Eq. (7) representing the prescribed desired periodicgait, whereas the limit cycles in the dark gray regions were differ substantially from the desired trajectory. Thus,the border between light and gray regions in each panel, if exists, represents a bifurcation set at which the stabilityof the limit cycle altered and a new limit cycle appeared. Stability of limit cycle (either γ or the one bifurcatedfrom γ ) was lost completely at another border between (light/dark) gray and black regions. Both types of bordersrepresent the critical edges of stability. (B) Scaling exponent α as the functions of P -gain values, where values ofthe exponents between 0 and 1 are color-coded defined at the bottom of the panels. All simulations for (B) wereperformed with the endogenous white Gaussian torque noise of σ = 0 .
003 Nm. For both (A) and (B), the four rowsof panels were obtained for the model with different combinations of the controllers. 1st-row: FF+PD, 2nd-row:FF+PD+INT. 3rd-row: FF+PD+PR. 4th-row: FF+PD+INT+PR. The arrow-heads in the 1st, 2nd, 3rd, and4th rows indicate the P -gain values used for the simulations shown in Figs. 1A, 1B, 1C, 1D and 1E(and also F),respectively. 1 /f -like fluctuation was observed at the grids with orange-yellow colors, which were distributed nearthe critical edges of stability. P -gainplane in Fig. 2A-3rd-row (Fig. 1D and the orange-yellow arcs in Fig. 2B-3rd-row). That is, thelong-range power-law correlation was observed when the P -gain for the model with FF+PD+PRwas located at the critical edge of stability of the prescribed limit cycle γ . FF+PD+INT+PR controllers.
The combined use of INT and PR controllers dramat-ically broadened the stability region in the P -gain parameter space (Fig. 2A-4th-row), by whichthe model could establish stable walking for very small P - D gains, implying stable gait with veryflexible (compliant) joint dynamics. In particular, the stability regions with the prescribed γ (thelight gray regions) occupied larger areas than those for the model with FF+PD+PR, as comparedwith the light gray regions between the 3rd and 4th rows of Fig. 2A. As a result, the border betweenthe light and dark gray regions marched into the small P -gain regime, wherein the exponent α values between 0.8 and 1.0 were widely distributed (the orange-yellow regions in Fig. 2B-4th-row)along the edge of stability (the bifurcation curve). Fig. 1E is such an example of stride intervaldynamics showing the long-range correlated 1 /f -like fluctuation with noise and the deterministicimpulse response with a slow decay approximated by a power-law function n − .Here, by reference to a clinical report showing a larger stochasticity in the amount of PR fornon-Parkinsonian patients with a symptom of freezing of gait (FoG) than for Parkinsonian patientswithout FoG [27], we examined how the 1 /f -like dynamics in the model with FF+PD+INT+PRcould be influenced by a stochasticity introduced on the amount of PR ( Methods ). It was shownthat the stochasticity of the PR could lead to a loss of the long-range correlation (Fig. 1F). Thisresult suggests that the randomness of the PR is one of the possible mechanisms that cause a lossof gait fractality.
Discussion
Summary.
Based on the numerical simulations of the gait model with different types of con-trollers, we showed that the model with endogenous torque noise can exhibit the fractal gait, ifthe PR mechanism is incorporated as a feedback controller, and the gains of time-continuous PDfeedback controller for tracking a reference gait motion are tuned at the critical edge of stability.Moreover, combined use of the INT feedback controller with PD and PR controllers remarkablylowers the critical gains of PD for stability and broadens the gain-parameter regions that exhibitgait fractality, by which the fractal gait becomes a commonplace phenomenon in terms of the areain the P -gain parameter space. Thus, these three conditions clarified in this study could be themajor mechanisms responsible for the genesis of fractal gait. It is logically apparent that the gaitdynamics that satisfied these conditions simultaneously naturally resolve the flexibility-stabilitytrade-off. The results suggest that pathological gait, characterized by joint-rigidity and/or loss offractality [9, 10], may be caused by dysfunction in some of these mechanisms. Stability versus flexibility.
A novel outcome of this study, besides the fractal issue, isthat the combined use of PR and INT with FF+PD controllers could establish gait stability withextremely small overall joint impedance (the small P -gain values). For example, the gain valuesof P a = P k = P h = 100 Nm/rad and D a = D k = D h = 10 Nm · s/rad corresponding to theparameter point that is located in the stability region of Fig. 2A-4th-row are quite comparablewith experimentally estimated joint stiffness of 200-300 Nm/rad for the lower extremities duringhuman walking [28, 29]. To the best of our knowledge, there might be no other human bipedal gaitmodels that can establish a stable walking with such a small joint impedance that allows flexiblejoint dynamics. Limitations.
The gait model used here is simplified from many points of view. In particular,it describes only two-dimensional sagittal motion during walking, unlike others [30], and is limitedto steady-state walking. The prescribed periodic desired joint-angle trajectory may have a limitedphysiological substance, though it can be conceptually considered as a set of motor output signalsproduced by central pattern generators [31] that exhibit PR in response to perturbations [26].6ecause the desired joint-angle trajectory used for the model is fixed as prescribed, and it cannotbe altered actively against perturbations, except its oscillating phase by the PR, the overall gaitstability of the model is probably more limited than that in a healthy human subject. For thisreason, the noise intensity used in this study could not be large enough, which made the varianceof stride intervals smaller than experimental reality. In perspective, however, the limitations ofthe model could be overcome by another computational layer, in addition to the control layerinvestigated in this study, that addresses the human capability to adapt to changing environmentalconditions by updating an original plan of action.
Related studies.
The mechanisms of stride-to-stride control that can generate the fractal gaitis not trivial; thus, it is not easy to hypothesize possible mechanisms. However, such mechanismsmight not be unique. Dingwell et al. have proposed a reconciling mechanism among optimality ofcontrol-signal (minimum intervention) and task-related error (errors in the walking velocity froma desired constant velocity), redundancy in stride length (distance) and stride interval (time) toachieve the desired constant walking speed (distance/time), and stochasticity by multiplicativemotor noise, with a concept of goal equivalent manifold (GEM) associated with the redundancy,which can generate the long-range power-law correlation of stride interval fluctuation along theGEM [14]. Unlike our neuromechanical modeling, their optimal feedback control model has nounderlying mechanics, and thus cannot be directly related to stability (falls) of the mechanicalbody dynamics. That is, the GEM is defined in the abstract space concerning the control strategy.Nevertheless, the gait model studied here and the GEM-based optimal control model could bedeeply inter-related with each other, as discussed in detail by [22] for a model of postural fluctuationduring standing, as the GEM, a kind of uncontrolled manifold defined for kinematic redundancyof the human body [32], might be associated with the stable manifold used for the INT controlas shown in our previous study [33]. That is, the redundancy for strategy and that for bodykinematics should be closely inter-related.Emergence of fractal gait near the edge of stability in this study is reminiscent of a recent reportshowing that energetic costs are minimized by the control at stability’s edge during expert stickbalancing [34]. Moreover, the best-fit model parameters for the stick balancing experiment usingthe INT control model in our previous study were also located at the critical edge of stability [35].Although the mechanical energy consumptions in our model did not show significant differencesamong four examined combinations of controllers (results not shown), which was probably becausethe fluctuation of the gait dynamics of the model was too small to make any difference, criticality instability and the long-range correlated fluctuation might generally be closely related. In this regard,clarification of the interrelationship among criticality, optimality, and fractality is the significantfuture challenge in motor control research.
Risk of falls and adaptability.
In view of importance of the fall risk assessment forprediction and prevention of falls in the aging and aged societies, the reliability of the gait indexavailable so far needs to be improved substantially [36]. Although some of the traditional measuresare not necessarily capable of characterizing motor functions properly, newly developed measures,such as those quantifying complexity and/or nonlinear dynamics of fluctuation are expected tocontribute to improving the quality of assessments [15, 37, 38]. A difficult, but yet scientificallyinteresting aspect of gait stability has been demonstrated well by Bruijn et al [39] who showedthe absence of a consistent pattern of correlations between orbital stability of the gait limit cycleevaluated by the local Lyapunov exponent and the size (the standard deviation) of variability ofgait dynamics. More importantly, alterations in a metric value can represent different meaningsdepending on the underlying mechanisms that alter the gait dynamics. For example, a recentstudy shows that increases in gait variability, which has been reported to be associated with fallhistory, may not imply impaired stride-to-stride control of walking in healthy older adults, butmay be due to increased neuromotor noise [40]. Gait fractality quantified by the scaling exponentmay also face a similar difficulty for practical use as a dynamical biomarker. As suggested bythe current study, the appropriate combination of feedback controllers can broaden the stabilityregions into the small-gain regime, which makes the fractal gait with flexible joint dynamics robustagainst changes in the gain-parameters. However, this fact per se does not imply high stability ofgait, but rather adaptability of the gait control system. Although local dynamics analyzed by any7ethodologies, including gait fractality, can never characterize the global stability of gait, such asstability determined by a largeness of the basin of attraction of the limit cycle for example [23, 24],we can speculate that a highly adaptable gait control system is more capable of exhibiting highglobal stability. Thus, the fractal gait measure should contribute to the assessment fall risks, withother measures characterizing different aspects of gait dynamics.
Methods
The equations of motion of the gait model, the controllers of FF, PD, INT, PR, and stochasticPR, and DFA are described below.
Equations of motion of the model.
The model with seven rigid-links in the sagittalplane from our previous work [25] was used in this study. The general coordinates of the modelcan be described as q = ( q , q , · · · , q ) T ≡ ( θ, x, y, θ la , θ lk , θ lh , θ ra , θ rk , θ rh ) T , (1)where θ , x and y are the tilt angle, horizontal position and vertical position of the head-arm-trunk(HAT) center of mass (CoM). θ ia , θ ik , and θ ih are the ankle, knee, and hip joint angles of the left( i = l ) and right ( i = r ) limbs. The equation of motion is represented as J ( q )¨ q + B ( q, ˙ q ) + K ( q ) + G ( q, ˙ q ) = U, (2)where J ( q ), B ( q, ˙ q ), K ( q ) and G ( q, ˙ q ) are the inertia matrix, centrifugal and Coriolis torque, grav-itational torque and the ground reaction force, respectively. U = (0 , , , U la , U lk , U lh , U ra , U rk , U rh ) T on the right-hand side represents the joint torques acting on the six joints. G ( q, ˙ q ) was modeledusing a nonlinear spring and damper acting forces on the heel and the toe of either foot [23, 25]. The desired joint angle profiles.
Joint angles of the model are represented by one partof the coordinate q , and we denote it byˆ q ≡ ( q , · · · , q ) T . (3)Time profiles of the desired joint angles or the desired joint angle trajectories, denoted by ˆ q d ( φ ) asa function of the gait phase φ or ˆ q d ( t ) as a function of time t , are defined for ˆ q asˆ q d ( φ ) ≡ ( q d ( φ ) , · · · , q d ( φ )) T (4)where the gait phase φ takes a value in [0 , T c ] with T c = 1 .
135 (sec) being the gait period of thesteady-state walking, or equivalentlyˆ q d ( t ) ≡ ( q d (mod T c ( t + φ )) , · · · , q d ((mod T c ( t + φ )) T , (5)for an appropriate initial phase φ of the desired joint angle trajectory, where mod T c takes a modulo T c of the argument. Note that φ and t are related as φ = φ ( t ) ≡ mod T c ( t + φ ). ˆ q d ( φ ) or ˆ q d ( t ) wasobtained as the Fourier expanded series of a human walking motion-captured data as detailed in[23, 25]. FF controller.
A benefit of the model used here with the desired joint angle trajectory ˆ q d ( φ ) isthat the model can walk stably when the time courses of all of the six joints in Eq. (2) are kinemat-ically constrained by ˆ q d ( φ ), by which only three kinematic degrees of freedom ( q , q , q ) and theirderivatives ( ˙ q , ˙ q , ˙ q ) are remained as dynamic variables [23]. A forward dynamic simulation ofEq. (2) with such kinematic constraints, after its transient dynamics, generates the correspondingtime profiles of ( q ss1 ( φ ) , q ss2 ( φ ) , q ss3 ( φ )). In this way, we obtain a complete sequence of kinematicsof the model during steady-state periodic gait as ¯ q ( φ ) = ( q ss1 ( φ ) , q ss2 ( φ ) , q ss3 ( φ ) , q d ( φ ) , · · · , q d ( φ )) T ,8nd the corresponding ground reaction force G (¯ q ( φ ) , ˙¯ q ( φ )). The inverse dynamics solution of Eq.(2) for the six joints, which are used as the outputs of the FF controller [25], was obtained as U ff ( φ ) = J (¯ q ( φ ))¨¯ q ( φ ) + B (¯ q ( φ ) , ˙¯ q ( φ )) + K (¯ q ( φ )) + G (¯ q ( φ ) , ˙¯ q ( φ )) . (6)The equation of motion with the FF controller is formulated as J ( q ( t ))¨ q ( t ) + B ( q ( t ) , ˙ q ( t )) + K ( q ( t )) + G ( q ( t ) , ˙ q ( t )) = U ff ( φ ( t )) . (7)In this study, we considered the periodic trajectory of the state point defined as ¯ x ( φ ) = (¯ q ( φ ) , ˙¯ q ( φ ))in the 18-dimensional state space of Eq. (7), referred to as γ in this section. By definition, γ is always a periodic solution of Eq. (7), regardless of its stability, although the FF controller alonecannot stabilize γ . PD feedback controller.
As in [25], the PD controller generates the feedback torquesdefined as U fb ( t ) = P (ˆ q d ( t + φ ) − ˆ q ( t )) + D ( ˙ˆ q d ( t + φ ) − ˙ˆ q ( t )) , (8)or equivalently U fb ( φ ) = P (ˆ q d ( φ ) − ˆ q ( φ )) + D ( ˙ˆ q d ( φ ) − ˙ˆ q ( φ )) , (9)where P = diag( P a , P k , P h , P a , P k , P h ) ,D = diag( D a , D k , D h , D a , D k , D h ) . The equation of motion with the FF+PD controller is formulated as J ( q )¨ q + B ( q, ˙ q ) + K ( q ) + G ( q, ˙ q ) = U ff ( φ ) + U fb ( φ ) . (10)Because γ is the solution of Eq. (7), it is also the solution of Eq. (10). FF+PD controllercan stabilize the limit cycle γ , if P and D are large enough. In this study, D was fixed at D a = D k = D h = 10 Nm · s/rad. P -gain values were altered systematically in the range of P h ∈ [0 , P a ∈ [0 , P k ∈ [0 , INT feedback controller.
The INT feedback controller provides feedback torques onlyduring an optimized brief time span, lasting for w (= 0 .
1) second, after the heel contact withinevery double support phase of the gait [25]. It is used together with FF+PD controllers, typicallywhen P - D gains of PD controller are small and located outside the stability regions. Because theINT controller is activated only in a sequence of brief time spans, it contributes to reducing theoverall impedance necessary for maintaining the walk.The INT controller exploits the stable manifold W s ( γ ) of unstable γ as the solution of Eq. (10)for small gains of PD controller, where W s ( γ ) was obtained locally around γ using the monodor-omy matrix (a linearized Poincaré map) with its stable eigenvalues (Floquet multipliers) and thecorresponding eigenvectors. Let γ ( φ ) be a state point on γ at phase φ , which means γ = ∪ φ γ ( φ ).Then, we can define W s ( γ ( φ )) as the intersection between W s ( γ ) and a Poincaré section of γ passing through γ ( φ ) for the linearized Poincaré map. The INT controller, when it is activated,aims to force the state point of the model at the phase φ toward W s ( γ ( φ )), not to γ ( φ ), using thefeedback torque defined as U int ( φ ) = P + (ˆ q s ( φ ) − ˆ q ) + D + ( ˙ˆ q s ( φ ) − ˙ˆ q ) , (11)where (ˆ q s ( φ ) , ˙ˆ q s ( φ )) ∈ W s ( γ ( φ )) with the gains of P + and D + defined as P + = diag( P + a , P + k , P + h , P + a , P + k , P + h ) ,D + = diag( D + a , D + k , D + h , D + a , D + k , D + h ) .P + a - D + a , P + k - D + k and P + h - D + h are the gains of the INT feedback controller acting on the ankle, kneeand hip joints, respectively, of the left and right limbs. U int is supplemented to the right-hand-side9f Eq. (10) intermittently only during on-period defined as [ φ on , φ on + w ] for w seconds. That is,the equation of motion with FF+PD+INT controllers can be described as J ( q )¨ q + H ( q, ˙ q ) = U ff ( φ ) + U fb ( φ ) + U int ( φ ) , (12)if φ ∈ [ φ on-L , φ on-L + w ] or φ ∈ [ φ on-R , φ on-R + w ], which is during the on-period of INT controller,and J ( q )¨ q + H ( q, ˙ q ) = U ff ( φ ) + U fb ( φ ) , (13)otherwise during the off-period of INT controller, where H ( q, ˙ q ) ≡ B ( q, ˙ q ) + K ( q ) + G ( q, ˙ q ). In thisstudy, φ on-L and φ on-R were set to the timings of 1.5 ms after every heel-contact event of the leftand right feet, respectively. See [25] for details. PR controller with FF+PD.
The PR controller used in this study resets the phase φ ofthe desired joint angle trajectory impulsively at every heel-contact of either foot. During steady-state walking of the model along the prescribed γ , heel-contact events of the left (right) foot occurperiodically with the gait period T c , for which the phases of left (right) heel-contact events arefixed at φ L ( φ R ). During transient gait or in the presence of noise that perturbs the state pointfrom γ , the phase of every heel-contact may deviate from the reference phase, i.e., φ L or φ R . ThePR controller shifts the phase φ of the desired joint angle trajectory to compensate such deviations.That is, when a heel-contact of the left (right) foot is detected, the phase φ of the desired trajectoryis reset to φ L ( φ R ). Thus, the amount of phase resetting is ∆( φ ) = φ − φ L for a left heel-contact,and it is ∆( φ ) = φ − φ R for a right heel-contact. If a heel-contact occurs earlier (later) than thereference timing, ∆( φ ) < φ ) > φ to φ − ∆( φ ), which is equal to either φ L or φ R for both of advanced and delayed resettings. Notethat the FF torques are also reset by the same amount when the desired joint angle trajectory isreset. It has been shown that the PR mechanism can increase gait stability and largely reduce theimpedance necessary for stability [24, 41].Stochasticity in the amount of PR was introduced using a truncated normal distribution. Thatis, for every heel-contact event at phase φ , the deterministic amount of PR ∆( φ ) was randomlymodified as ∆( φ ) + ξ where ξ ∼ N (0 , σ r ) with σ r was set to 10 ms and any outcomes beyond 2 σ r were discarded so that n took values within [ − σ r , +2 σ r ] only. FF+PD+INT+PR.
When the PR and INT controllers were used together with FF+PD,the sequence of executing two time-discontinuous control mechanisms was detecting a heel-contactfirst, then resetting the phase of the desired joint trajectory and the FF torques. Then, the INTcontroller was turned on 1.5 ms after the PR was performed.
Torque noise.
Stochastic fluctuation of the gait was induced by adding a set of six independentseries of white Gaussian torque noise (with an identical standard deviation σ = 0 .
003 Nm) actingon the 6 joints. Because the gait model in this study focused only on the recovery against smallperturbations when the desired joint-angle trajectory and thus the geometry of the limit cycle werefixed, the stride interval fluctuation could not be large enough in comparison with the real humandata. However, we intended to use large noise amplitude as much as possible.
DFA.
The standard procedure of DFA [42, 43] was applied for the stride interval time seriesdata { x ( i ) } N − i =0 of length N . We obtained an integrated time series y ( k ) = P k − i =0 ( x ( i ) − ¯ x ),where ¯ x is the sample mean of the series. Then { y ( k ) } Nk =1 was divided into equal segments oflength n without overlapping, so the segment number was M = N/n . In each segment of length n with an index l ( l = 1 , · · · , M ), a least square polynomial p ( l ) n ( k ) as the trend of the segmentwas fitted to the data. Square of deviations in the segment l were summed up to obtain f ( l ) n = P lnk =( l − n +1 ( y ( k ) − p ( l ) n ( k )) . The root of average square deviation of all segments was thencalculated as F ( n ) = [ M − P Ml =1 f ( l ) n ] / . For each n , the corresponding F ( n ) was calculated,finally the scaling exponent α was obtained by the linear fitting of log F ( n ) and log n .Longer data length is preferred to obtain reliable results. Typically in the literature, 1 hour walkincluding about 3,000 strides are tested experimentally [44]. Here we adopt the higher standard10f 3,000 strides for every P -gain parameter to quantify a scaling exponent. The orders of thepolynomial ranging from 1 to 3 were examined for many preliminary sets of P -gain, for each ofwhich we obtained 10 simulated trials of 3,000 strides to calculate a mean of 10 exponent values.We found that the exponent was not altered quantitatively for the different orders. Thus, we usedthe linear fitting (the order 1) for calculating the exponent for the wide-range parameter study,in which the exponent was calculated from only one trial of 3,000 strides for each P -gain due tosignificant computational time necessary for the simulations. Numerical simulation.
Numerical simulations of the model were performed by integratingequations of motion simply using the Euler-Maruyama method [45] with a time step of ∆ t = 10 − seconds. The robustness of this numerical simulation method was tested by time step ten timeslarger and smaller. For each simulation, an initial state ( q (0) , ˙ q (0)) and an initial phase φ werespecified, in which the initial state was set close to the limit cycle such that | q (0) − ¯ q ( φ ) | < ǫ wassatisfied for a given φ . ǫ was set as 10 − rad on the HAT tilt angle.In the model with PR and INT controllers in the presence of noise, we considered that a heelcontacted on the ground if the ground reaction force remained non-zero continuously more than 20time steps. PR was performed immediately after a detection of heel-contact. Then, with a latency(during which only FF+PD controllers were employed) of w = 1 . w = 0 . P a , P k , and P h werevaried between 0 and 1500 Nm/s with fixed D -gains at 10 Nms/rad) for the model with INT andPR controllers took about 10 days. Acknowledgements.
This work was supported in part by JSPS grants-in-aid 26242041 and26750147, MEXT KAKENHI Grant Number 16H01614 (Non-linear Neuro-oscillology), and MEXTgrant for Post-K project (primary issue 2).
Author contributions.
T.N. conceived and designed the study. C.F. and Y.S. performedthe modeling, simulations and data analysis. C.F, Y.S., P.M., and T.N. wrote the paper.
Competing interests.
The authors have declared that no competing interests exist.
References [1] Hausdorff, J. M.
Hum. Mov. Sci. (4), 555–589 (2007).[2] Peng, C.-K., Mietus, J., Hausdorff, J. M., Havlin, S., Stanley, H. E., and Goldberger, A. L. Phys. Rev. Lett. , 1343–1346 Mar (1993).[3] Goldberger, A. L., Amaral, L. A. N., Hausdorff, J. M., Ivanov, P. C., Peng, C.-K., and Stanley,H. E. Proc. Natl. Acad. Sci. USA (suppl 1), 2466–2472 (2002).[4] Gu, C., Coomans, C. P., Hu, K., Scheer, F. A. J. L., Stanley, H. E., and Meijer, J. H. Proc.Natl. Acad. Sci. USA (8), 2320–2324 (2015).[5] Pailhous, J. and Bonnard, M.
Behav. Brain Res. (2), 181–189 (1992).[6] Hausdorff, J. M., Peng, C. K., Ladin, Z., Wei, J. Y., and Goldberger, A. L. J. Appl. Physiol. (1), 349–358 (1995).[7] Frenkel-Toledo, S., Giladi, N., Peretz, C., Herman, T., Gruendlinger, L., and Hausdorff, J. M. J. Neuroeng. Rehabil. (1), 23 (2005).[8] Dingwell, J. B. and Cusumano, J. P. Gait Posture (3), 348 – 353 (2010).119] Hausdorff, J. M., Mitchell, S. L., Firtion, R., Peng, C. K., Cudkowicz, M. E., Wei, J. Y., andGoldberger, A. L. J. Appl. Physiol. (1), 262–269 (1997). PMID: 9029225.[10] Hausdorff, J. M. Chaos (2), 026113 (2009).[11] Herman, T., Giladi, N., Gurevich, T., and Hausdorff, J. Gait Posture (2), 178 – 185 (2005).[12] Ashkenazy, Y., Hausdorff, J. M., Ivanov, P. C., and Stanley, H. E. Physica A (1–4),662–670 (2002).[13] Ahn, J. and Hogan, N.
PLoS One (9), e73239 (2013).[14] Dingwell, J., John, J., and Cusumano, J. PLoS Comput. Biol. (7), 14 (2010).[15] Bruijn, S. M., Meijer, O. G., Beek, P. J., and van Dieën, J. H. J. R. Soc. Interface (83), 0(2013).[16] Stergiou, N., Harbourne, R., and Cavanaugh, J. J. Neurol. Phys. Ther. (3), 120–129 (2006).[17] Todorov, E. Proc. Natl. Acad. Sci. USA (28), 11478–11483 (2009).[18] Hogan, N.
IEEE Trans. Automat. Contr. (8), 681–690 (1984).[19] Burdet, E., Osu, R., Franklin, D., TE, M., and Kawato, M. Nature (6862), 446–449 (2001).[20] Bottaro, A., Yasutake, Y., Nomura, T., Casadio, M., and Morasso, P.
Hum. Mov. Sci. (3),473 – 495 (2008).[21] Asai, Y., Tasaka, Y., Nomura, K., Nomura, T., Casadio, M., and Morasso, P. PLoS One (7),e6169 (2009).[22] Suzuki, Y., Nomura, T., Casadio, M., and Morasso, P. J. Theor. Biol. , 55 – 79 (2012).[23] Yamasaki, T., Nomura, T., and Sato, S.
Biol. Cybern. (6), 468–496 (2003).[24] Yamasaki, T., Nomura, T., and Sato, S. Biosystems (1–2), 221 – 232 (2003).[25] Fu, C., Suzuki, Y., Kiyono, K., Morasso, P., and Nomura, T. J. R. Soc. Interface (101)(2014).[26] Schillings, A. M., van Wezel, B., Mulder, T., and Duysens, J. J. Neurophysiol. (4), 2093–2102 (2000).[27] Tanahashi, T., Yamamoto, T., Endo, T., Fujimura, H., Yokoe, M., Mochizuki, H., Nomura,T., and Sakoda, S. PLoS One (12) (2013).[28] Shamaei, K., Sawicki, G., and Dollar, A. PLoS One (3), e59935 (2013).[29] Shamaei, K., Sawicki, G., and Dollar, A. PLoS One (12), e81841 (2013).[30] Roos, P. E. and Dingwell, J. B. J. Biomech. (15), 2929 – 2935 (2010).[31] Minassian, K., Hofstoetter, U. S., Dzeladini, F., Guertin, P. A., and Ijspeert, A. Neuroscientist (6), 649–663 (2017). PMID: 28351197.[32] Latash, M., Scholz, J., and Schöner, G. Exerc. Sport Sci. Rev. (1), 26–31 (2002).[33] Suzuki, Y., Morimoto, H., Kiyono, K., Morasso, P. G., and Nomura, T. Front. Hum. Neurosci. , 618 (2016).[34] Milton, J., Meyer, R., Zhvanetsky, M., Ridge, S., and Insperger, T. J. R. Soc. Interface (119) (2016).[35] Yoshikawa, N., Suzuki, Y., Kiyono, K., and Nomura, T. Front. Comput. Neurosci. , 34(2016). 1236] Verghese, J., Holtzer, R., Lipton, R., and Wang, C. J. Gerontol. A Biol. Sci. Med. Sci. (8),896–901 (2009).[37] Toebes, M. J., Hoozemans, M. J., Furrer, R., Dekker, J., and van Dieën, J. H. Gait Posture (3), 527 – 531 (2012).[38] Zhou, J., Habtemariam, D., Iloputaife, I., Lipsitz, L., and Manor, B. Sci. Rep. (1) (2017).[39] Bruijn, S. M., van Dieën, J. H., Meijer, O. G., and Beek, P. J. J. Biomech. (10), 1506 –1512 (2009).[40] Dingwell, J. B., Salinas, M. M., and Cusumano, J. P. Gait Posture , 131 – 137 (2017).[41] Aoi, S., Ogihara, N., Funato, T., Sugimoto, Y., and Tsuchiya, K. Biol. Cybern. (5),373–387 (2010).[42] Peng, C.-K., Buldyrev, S. V., Havlin, S., Simons, M., Stanley, H. E., and Goldberger, A. L.
Phys. Rev. E , 1685–1689 Feb (1994).[43] Peng, C.-K., Havlin, S., Stanley, H., and Goldberger, A. Chaos (1), 82–87 (1995).[44] Hausdorff, J., Purdon, P., Peng, C.-K., Ladin, Z., Wei, J., and Goldberger, A. J. Appl. Physiol. (5), 1448–1457 (1996).[45] Higham., D. SIAM Rev.43