Ankle Push-Off Based Mathematical Model for Freezing of Gait in Parkinson's Disease

1. Introduction

Parkinson's disease results from the loss of neurons in the substantia nigra pars compacta of the basal ganglia (BG) (Davie, 2008), which has projections toward the motor, and cognitive areas (Albin et al., 1989; Alexander and Crutcher, 1990). Freezing of gait (FoG) is a motor disability in PD patients where subjects experience an “episodic absence or marked reduction of forwarding progression of the feet despite the intention to walk” (Nutt et al., 2011, p 734). This debilitating symptom occurs during the late-stage PD (Giladi et al., 2001) and is known to be very difficult to predict and control. The physiology of this symptom is not yet established conclusively, consisting of both neural and mechanical components. A set of correlations between the neural inputs (e.g., Dysfunction of Visuomotor and occipito-parietal pathways) and mechanical variables (e.g., gait pattern generation and automaticity) of PD-FoG have been studied (Heremans et al., 2013) but causality is not well-established. Apart from freezing, abnormal gait patterns in PD consists of high stride time variability with less reduction in stride length (Heremans et al., 2013). The relationship between these abnormalities and freezing is also not well-understood.

Gait has been studied fundamentally from two different perspectives. One that of robotics and control, and the second, biophysical. Mathematical models of passive gait have been studied extensively by several authors to understand their stability (e.g., Goswami et al., 1996; Manchester et al., 2011; Dai and Tedrake, 2013; Sadeghian and Barkhordari, 2020), and the effect of external conditions such as ramps (McGeer, 1990) and bifurcations has been investigated (e.g., Mahmoodi et al., 2013; Iqbal et al., 2014; Fathizadeh et al., 2018, 2019; Znegui et al., 2020). Impulsive dissipation at heel strike is studied for a multidimensional biped model in (Ros et al., 2015). There are other approaches to motor control using optimal control which demands an arbitrary or learned cost-functionals (e.g., Flash and Hogan, 1985; Pekarek et al., 2007; Parakkal Unni et al., 2017). These models are not sufficient to understand human locomotion in PD patients as these papers have focused on the stability behaviors and control of robots. Some of these cost functional/error minimization based models, even though they assume the existence of such a cost, have the advantage of being useful for extracting parameters easily from the data (Delp et al., 2007; Wu et al., 2019). However, they do not address explicitly how the external inputs result in high variability and freezing observed in PD gait (Heremans et al., 2013).

On the other hand, a biophysical model proposed in Taga (1995) considers the interaction with the Central Pattern Generators (CPG). The aim of the model is primarily to demonstrate walking as a stable limit cycle that emerges from the dynamic interaction between neural oscillation originating in CPG and the pendulum oscillation of body linkages, rather than involuntary stoppage of gait and variability. CPG-based complex model, which depends on several parameters such as the strength of neural connections, the magnitude of the coefficients in the rhythmic force controller, and strength of sensory inputs, has its significance. However, one drawback of such CPG-based complex model is often the dependence on an excessive number of parameters as described above to be determined for achieving a desired locomotor pattern over a large search space. To identify and tune such parameters for attaining involuntary freezing and variability of gait behavior for a wider population of patients is rather an arduous trial-and-error or learning and optimization based task. Involuntary stoppage of gait and variability is the key detail that is necessary to show the model's ability to display PD walking behavior. The effect of opposing forces generated by the plantar flexors observed in PD, as reported in Nieuwboer et al. (2004), is yet another detail for understanding the PD gait.

The model by Muralidharan et al. (2014) successfully captures the neural dynamics of basal ganglia (BG) but does not focus on the mechanics. A model which combines the chaotic region of the Lorentz system with the passive dynamic walker by Montazeri Moghadam et al. (2018), adds chaos externally, which makes it less relevant biophysically. However, these authors have established a need for explaining the variable nature of PD walking. As chaos is known to be absent in the basal ganglia (BG) of a PD patient (Mandali et al., 2015) the neuro-mechanical interactions need to be studied to find out its underpinnings. Another way to look at gait biophysically is through the equilibrium point hypothesis (Feldman, 1986; Duan et al., 1997), which suggests movements are the result of active changes in the equilibrium state of the motor system. Torque length characteristics of the muscles can be changed by a neural controller to achieve motion. It has been shown that the muscles act synergistically to reduce the variability in the targeted action. For example, the uncontrolled manifold hypothesis by Latash et al. (2002) explains the variability in muscle recruitment as “good” and “bad” regions of variability, depending on whether they achieve the targeted action or not. The muscle recruitments which achieve targeted action are considered “good” regions of variability, whereas the muscle recruitments which does not achieve targeted action are considered as “bad” regions of variability. The Equilibrium point approach (Feldman, 1986) makes the electromyogram (EMG) activity implicit. Another limitation of the equilibrium point hypothesis is in its difficulty in testing. The empirical determination of invariant characteristics (Sainburg, 2015) such as torque-length characteristics (Feldman, 1986) is necessary for validating the equilibrium point hypothesis. A way in which it is achieved is by asking the subjects “not to intervene” (Feldman, 1986; Sainburg, 2015) while doing a task such as unloading and assuming this results in stabilization of central commands to muscles (Sainburg, 2015). But this assumption is not necessarily true as there could be involuntary responses. The neuromuscular system is over-actuated with redundancies, as it contains more actuators than the degree of freedom. Use of muscle synergies (Latash, 2010) in models is one way to address redundancies. These models assume co-activation of a set of muscles as motor primitives to address the redundancy associated with muscle activation (Aoi et al., 2019; Tamura et al., 2020). The idea of muscle synergy is still debated and is considered difficult to refute by any empirical evidence or falsify (Popper, 2002; Olszewski and Sandroni, 2011) by some authors (Tresch and Jarc, 2009).

There is a need for a model to explain the empirical observations in PD gait, such as the high coefficient of variability and freezing near narrow passages (Snijders et al., 2008). Such a model will help in understanding the essential aspects of neural and mechanical systems contributing to PD gait, also shedding light into the future experimentations required. In this work, the relationships between the high variability and freezing will be studied by deriving a set of forces acting on the stance leg. They phenomenologically represent the EMG (Electromyogram) signals and therefore the activity of the CPGs. Kinetics of both swing leg and stance leg will be studied to better understand their roles under the action of these forces.

In summary, the model is built with two aims. The first aim is to explain the empirical observations that are seen in PD-Gait with a minimum number of variables. These include (1) a high coefficient of variation in PD subjects (Heremans et al., 2013), (2) a pattern of reduction of step lengths before freezing (Nutt et al., 2011), (3) the ability of sensory and visual cues to help reduce freezing (Rochester et al., 2005; Young et al., 2014; Amini et al., 2019), (4) the difficulty of freezing prediction, and (5) the occurrence of freezing near obstacles and narrow passages (Snijders et al., 2008). Secondly, the model aims to show the role of the swing leg as a supplier of ankle push of force as well as one that determines the time of heel strike. Hence, a bipedal model and a reduced low dimensional model resembling an inverted pendulum are studied upon the action of the ankle push-off force. The movement of the CoM under the action of the ankle push-off force is depicted in Figure 1. Hence, the hypothesis investigated in this study is that the variability and the motor symptoms associated with PD (Heremans et al., 2013) can be explained by the experimentally observed premature activation of plantar flexors observed in PD (Nieuwboer et al., 2004).

FIGURE 1

Figure 1. Anatomical representation of the stance phase from the ankle push-off to the heel strike is shown. The position of the center of mass (CoM) is shown as a red circle which is assumed to be rotating with respect to a pivot point S. The location of the plantar flexors, approximate region generating ankle push-off force and heel strike region are noted. The initial angle and angular velocity are represented by θ₀ and ω₀, respectively.

2. Materials and Methods—Modeling

2.1. Physiology

Walking is a complex process which involves the interaction of the brain, spinal cord and the musculoskeletal systems (Nutt et al., 2011). The typical gait cycle associated with walking involves “stance” and “swing” phases. The stance phase begins with a crucial heel strike phase, which is the initial contact that occurs instantaneously. As soon as the stance phase ends, the swing phase begins. The plantar flexor muscles of the trailing leg, supply energy to “push-off” the contra-lateral leading leg (Zelik and Adamczyk, 2016). Once the push-off occurs, the trailing leg enters the swing phase. The soleus and gastrocnemius muscles are the most notable plantar flexors, of which the significant role of the latter one in PD freezing/walking is established (Nieuwboer et al., 2004). Even though physiologically there is a non-linear relationship between the EMG signals and the torques generated (Genadry et al., 1988), a linear relationship can be assumed (Hof and Van Den Berg, 1977) between the envelope of the EMG (CPG firing) and the torques generated about the joint. Several other muscles are involved in walking, but the present study investigates only the effect of plantar flexors as these muscles supply most of the energy required for walking. In this work, the “freezing step” is defined as the step at which the legs do not have sufficient angular momentum to progress walking forward. When the physiology is modeled as an inverted pendulum-like system, the freezing results in backward motion of the stance leg. In a real-life scenario, this implies the patient either falls or stops movement. The remark 1 defines freezing and related terms used in this work.

REMARK 1. In this study, “freezing” or “freezing event” is defined as the condition where there is no forward motion of the stance leg. “Freezing episode” is defined as the events happening in the time interval between the heels strike phase the freezing event. Hence, the “start of the freeze” is defined to be at the heel strike phase after which a freezing event occurs.

Here, we carry out a systematic stability analysis, including unstable regimes, of the model in contrast to the stable limit cycle behavior studied in robotics (Grizzle et al., 2001) and passive walking dynamics literature (McGeer, 1990). Even though the complex freezing behavior can be explained through several possible routes (Nutt et al., 2011) (some of them purely based on neural control) an attempt is made here to explain it in the simplest possible way and to understand the effect of neuromuscular inputs in generating unstable and chaotic walking behavior as observed in PD (Heremans et al., 2013).

2.2. Dynamics of Walking

The dynamics of walking involves the coordinated action of neural input and muscles of the limbs. It consists in a continuous movement of the limbs as well as state reset at heel strike resulting in discrete dynamics (Sinnet et al., 2011). Plantar flexors of the swing leg supply the necessary torque to push the CoM forward. During walking, the CoM is supported by leg in stance phase for the majority of the time (80-90%) as the double support phase is approximately 10-20% of the overall gait cycle (Wahde and Pettersson, 2002; Kharb et al., 2011). The motion of the CoM with single support under the action of the plantar flexors is modeled in this work. The heel strike is modeled using discrete dynamics.

Traditionally the dynamics of walking is often modeled as biped model (Taga, 1995). In this section, a simplified biped model is presented. Further, a reduced, low dimensional, inverted pendulum system is considered as a special case relevant to PD. It is assumed that CoM is at the tip of the pendulum, and the links and the swing leg are massless. Therefore, the model generates the motion of the CoM of the human body. Running and jumping gaits are not considered in this model as the links are assumed to be rigid. It is also assumed that sufficient friction exists to avoid any slip. The angular displacements are assumed to be small enough (< 0.5 rad.) (Usherwood, 2005; Ranavolo et al., 2011; Polese et al., 2012) to allow for first/second-order approximations during the stance phase. Kane's method (Kane and Levinson, 1985) is used to derive the equations of motion (EoM). Kane's dynamical equation is of the form ${\overline{F}}_r+{\overline{F}}_r^{*}=0$, where ${\overline{F}}_r$ and ${\overline{F}}_r^{*}$ represents generalized active forces and generalized inertia forces, respectively, as described in Kane and Levinson (1985) (chapter 6, page 159). The equations in the form necessary for simulation is obtained using python libraries, the details of which are given in Gede et al. (2013).

Symbols m₁, m₂ represent the mass of the body and swing leg, respectively, as shown in Figure 2. The length of both the legs is represented by l. The variables associated with the system are the components of the vector $x={\left[{\theta }_1,{\theta }_2,{\omega }_1,{\omega }_2\right]}^T$ which are the angles and angular velocities (w.r.t inertial frame for stance leg and w.r.t. stance leg frame for swing leg) as indicated in Figure 2. There are two types of angular velocities. The one which corresponds to the rotation of the rigid body with respect to its center of rotation is called spin angular velocity and one which corresponds to the revolution of a point with respect to an origin is called orbital angular velocity. In this work, spin angular velocities with respect to the center of rotation of the rigid links are considered as they rotate about the center of rotation.

FIGURE 2

Figure 2. The toe-off, mid stance, and heel strike instances of the two connected links of the biped with its CoM while walking in the forward direction. CoM indicated as a circle represents a point mass at the tip of the pendulum. The swing leg is indicated in red. The point of collision during a heel strike instance is circled in red. The terms θ₁ and θ₂ are the angles that the stance and swing legs subtend w.r.t. the vertical (in the inertial reference frame), respectively. The quantities ω₁, ω₂ and m₁, m₂ are the corresponding angular velocities and masses of the body and swing leg, respectively. The torques that are acting on the stance leg are indicated as ζ₁(t) (Ankle push-off force) and ζ₂(t) (Torque due to the activation of the plantar flexors of the stance leg). States immediately before and after heel strike is indicated with “-” and “+” superscripts, respectively.

Hybrid systems (Lunze and Lamnabhi-Lagarrigue, 2009) are a class of dynamical systems, which exhibit both continuous states and discrete mode dynamics often associated with events such as resets, jumps, and switching. The continuous behavior is typically governed by a system of differential equations (similar to Equation 1) and the discrete part is governed by a vector-valued function (similar to Equation 2) (Lunze and Lamnabhi-Lagarrigue, 2009)(chapter 1). The transition between the discrete and continuous governing equations is determined by the state of the system in the overall phase space. The dynamics in this work is governed by the general hybrid dynamical system of the form,

$\begin{array}{ll}\stackrel{.}{x}=q(x),x^{-}\notin S& (1)\end{array}$

$\begin{array}{ll}{x}^{+}=\Delta (x^{-}),x^{-}\in S& (2)\end{array}$

where,

$\begin{array}{ll}S:=\{x\in \chi |g_{reset}(x)=0\}& (3)\end{array}$

where, q(x), Δ(x⁻), and g_reset(x) are continuous vector valued functions of x. In the absence of external torques acting on the leg, the term q(x) is,

$\begin{array}{ll}q(x)=\left[\begin{array}{c}{\omega }_1\\ {\omega }_2\\ \frac{1}{2l(m_1+m_2{\text{sin }}^2({\theta }_2))}(2g{m}_1\text{sin }({\theta }_1)-g{m}_2\text{sin }({\theta }_1+2{\theta }_2)\\ +g{m}_2\text{sin }({\theta }_1)+l{m}_2{\omega }_1^2\text{sin }(2{\theta }_2)-2l{m}_2{\omega }_2^2\text{sin }({\theta }_2))\\ -\frac{\text{sin }({\theta }_2)}{l(m_1+m_2{\text{sin }}^2({\theta }_2))}(g{m}_1\text{cos }({\theta }_1)+g{m}_2\text{cos }({\theta }_1)\\ -l{\omega }_1^2(m_1+m_2)+l{m}_2{\omega }_2^2\text{cos }({\theta }_2))\end{array}\right]& (4)\end{array}$

The function, Δ(.) is the reset map, χ ⊂ ℝ⁴ the state space, and g_reset(.) is the function that defines the heel strike. The set “S” defines a surface where the heel strikes the ground and the states change abruptly according to the reset map. The x⁻ and x⁺ indicate the states immediately before and after the heel strike, respectively. The functions Δ(.) and g_reset(.) are described in sequel.

As the body mass is considerably larger than the mass of the leg, the case where m₂ goes to 0 has only been considered. Further, small-angle approximation leads to the following equation for ${\stackrel{.}{\omega }}_1$ and ${\stackrel{.}{\omega }}_2$

$\begin{array}{ll}{\stackrel{.}{\omega }}_1=g\frac{{\theta }_1}{l}& (5)\end{array}$

$\begin{array}{ll}{\dot{\omega }}_2=-\frac{(g-l{\omega }_1^2){\theta }_2}{l}& (6)\end{array}$

The ankle push-off forces of the stance leg supply majority of the energy needed to propel the leg forward. When this neuromuscular forcing Γ(t) is added to the stance leg, the Equation (5) becomes,

$\begin{array}{ll}{\stackrel{.}{\omega }}_1:=g\frac{{\theta }_1}{l}+\Gamma (t)& (7)\end{array}$

The forcing term Γ(t) is derived in the following section.

2.3. Derivation of the Forcing Terms

The torques acting on the stance leg are derived in this section to generate the neuromuscular forcing term Γ(t) for the stance leg. Torque produced by the plantar flexors on the trailing leg is defined as G_r(t) and that on the leading leg as G_l(t). These torques are assumed to be linearly related to the envelop of the EMG signals which are positive functions of time (Nieuwboer et al., 2004). The torque G_r(t) generates the ankle push off force F(t) and is assumed to be in phase with the heel strike. In the proposed model, the force F(t) and the torque G_l(t) are assumed to be

$\begin{array}{ll}{G}_l(t):={\tau }_l(\text{sin}(2\pi f_{r_1}t+\varphi )+1)& (8)\end{array}$

$\begin{array}{ll}F(t):={\tau }_r(\text{sin}(2\pi f_{r_2}t)+1)& (9)\end{array}$

where τ_l and τ_r are constants. The variables, f_r1, f_r2 and ϕ represent frequencies and the phase difference between torques on the leading and trailing leg, respectively. Both frequencies are assumed to be unity. The ankle push off torque acting on the leading leg (in stance phase) can be calculated using the free body diagram shown in Figure 3. The pivot points of the trailing leg and leading leg are “O” and “S,” respectively. Trailing leg and leading leg subtends the same angle θ₁ w.r.t. the normal to the ground as the trailing leg and leading leg together with the ground is assumed to form an isosceles triangle. By balancing the moments about the point “S” in Figure 3 yields,

$\begin{array}{ll}I{\stackrel{.}{\omega }}_1=G_l(t)-lF(t)\text{sin}({\theta }_h)+mgl\text{ sin}({\theta }_1)& (10)\end{array}$

When angle θ₁ is small (sin(θ) ≈ θ) and since $I=m_1{l}^2$, Equation (10) is rewritten as,

$\begin{array}{ll}{m}_1{l}^2{\stackrel{.}{\omega }}_1=G_l(t)-lF(t)\text{sin}({\theta }_h)+m_{1}gl{\theta }_1& (11)\end{array}$

Substituting F(t) and G_l(t) from Equation (8) and (9) in Equation (11), one obtains,

$\begin{array}{l}{m}_1{l}^2{\stackrel{.}{\omega }}_1={\tau }_l(\text{sin}(2\pi f_{r_1}t+\varphi )+1)-l{\tau }_r(\text{sin}(2\pi f_{r_2}t)+1)\text{sin}({\theta }_h)\\ \begin{array}{l}+m_{1}gl{\theta }_1\end{array}& (12)\end{array}$

Rearranging Equation (12), angular acceleration of the leading/stance leg is,

$\begin{array}{ll}{\stackrel{.}{\omega }}_1=\stackrel{\text{Plantar flexors of the leading leg}(={\zeta }_2(t))}{\overbrace{\frac{{\tau }_l(\text{sin}(2\pi f_{r_1}t+\varphi )+1)}{l^2{m}_1}}}& (13)\end{array}$

$\begin{array}{l}-\stackrel{\text{Ankle push off from trailing leg}(={\zeta }_1(t))}{\overbrace{\frac{{\tau }_{r}l\text{ sin}({\theta }_h)(\text{sin}(2\pi f_{r_2}t)+1)}{l^2{m}_1}}}+\stackrel{Gravity}{\overbrace{\frac{g{\theta }_1}{l}}}\\ \begin{array}{l}:=\frac{{\zeta }_2(t)-{\zeta }_1(t)}{m_1{l}^2}+\frac{g{\theta }_1}{l}\end{array}& (14)\end{array}$

where $\Gamma (t):=\frac{(-{\zeta }_1(t)+{\zeta }_2(t))}{m_1{l}^2}$ (as shown in Figure 2) is the time varying neuromuscular forcing.

FIGURE 3

Figure 3. The diagram visualizes the forces and moments on leading and trailing leg that enable the movement of the center of mass (CoM) forward. Symbols S and O indicate the points on leading and trailing leg about which the torques G_l and G_r are applied. G_r is the torque generated by the plantar flexors of the trailing leg, which results in a force F (ankle push-off) acting on the leading leg, about the point “S.” The distance between the pivot point to the point of action of the force is l_f. The plantar flexors of the leading leg generate a torque G_l in the leading leg, in the opposite direction. The angle the leading leg subtends with the vertical axis at S is θ₁. θ_h represents the hip angle. The moments are balanced about “S” to get the equations of motion (EoM). An approximate position of the starting stance phase is shown in the background in gray color.

The initial velocity of the swing leg is assumed to be constant in every step. As the mass m₂ is assumed to be zero, the corresponding angular momentum is also equal to zero. Therefore, the angular velocity of the stance leg is reset to conserve its angular momentum and the initial angular velocity of the swing leg after reset is assumed to be a positive constant to account for the impulse during the ankle push-off. This is a valid assumption as the definition of freezing in this work is independent of the swing leg movement. A reset is carried out when θ₁ + θ₂ = 0 and θ₁ < 0. Using Equations (6) and (14), the equations of motion of a biped can be written as in Equations (1)–(2), where the functions q(x), Δ(x), g_reset(x) and the set S are written as follows,

$\begin{array}{ll}q(x):=\left[\begin{array}{c}{\omega }_1\\ {\omega }_2\\ \frac{{\zeta }_2(t)-{\zeta }_1(t)}{m_1{l}^2}+\frac{g{\theta }_1}{l}\\ \\ -\frac{(g-l{\omega }_1^2){\theta }_2}{l}\\ \end{array}\right]& (15)\end{array}$

$\begin{array}{l}S:=\{{\left[\begin{array}{c}{\theta }_1,{\theta }_2,{\omega }_1,{\omega }_2\end{array}\right]}^T\in \chi |g_{reset}\left({\left[\begin{array}{c}{\theta }_1,{\theta }_2,{\omega }_1,{\omega }_2\end{array}\right]}^T\right)\\ =0\wedge {\theta }_1<0\},\\ g_{reset}(.):={\theta }_1+{\theta }_2,\\ \Delta (.):=\left[-{\theta }_1,-{\theta }_2,{\omega }_1\text{ cos}({\theta }_h),{\omega }_2^0\right]& (16)\end{array}$

with initial conditions ${\omega }_1(0)={\omega }_1^0,{\theta }_1(0)={\theta }_1^0,{\omega }_2(0)={\omega }_2^0,{\theta }_2(0)={\theta }_2^0$.

REMARK 2. It may be noted that the torque generated by the plantar flexors is assumed to act about the point “O” as the pivot point. Balancing the moments due to the ankle push off force, F(t) and G_r(t) about the point O, the ankle push-off force can be determined in terms of G_r(t) as,

$\begin{array}{ll}F(t)=\frac{G_r(t)}{l_f\text{cos}({\theta }_1)}& (17)\end{array}$

Here, the distance, l_f is taken between the heel to the pivot point on the foot for calculating the moments, as shown in Figure 3. Therefore, implicitly, the following assumption has been made while prescribing F(t).

$\begin{array}{ll}\frac{G_r(t)}{l_f\text{ cos}({\theta }_1)}:={\tau }_r(\text{sin}(2\pi f_{r_2}t)+1)& (18)\end{array}$

2.4. Rationale for Using a Low Dimensional Model for Analysis

The angular velocity of the stance leg contributes directly to the angular acceleration of the swing leg. A higher absolute angular velocity of the stance leg leads to lower acceleration of the swing leg. However, the dynamics of the stance leg in Equation (14) is uncoupled from the dynamics of the swing leg, and hence resembles the dynamics of an “inverted pendulum system.” It may be noted that the term “$l{\omega }_1^2$” can be approximated to a constant as in the physiological range of low angular velocities (especially in PD patients) $g>>l{\omega }_1^2$ (typically quantity $l{\omega }_1^2=0.448\text{ m}.{\text{rad}}^2.s^{-2}=0.7\times 0.8^2$ is of order “0” while g = 9.8 m.s⁻² is of order 1). This results in a condition where the swing leg acts independently to the stance leg, effectively determining the step length. Therefore, an inverted pendulum walking model for PD subjects is valid when constant step length is assumed. Figure 1 depicts the physical rationale behind the use of an inverted pendulum model. The constant step length assumption is general enough to explain the variability in stepping as this leads to variability in stepping angular velocities rather than step lengths. In summary, in the following sections we present analysis of the stance phase walking model in light of the PD walking behavior at a constant step length. Physiologically, the hip applies torques on the swing leg and controls its initial angular velocity. The hip torques acting on the swing is not relevant in propelling the CoM forward, as most of the torque required for that is supplied by the ankle (Zelik and Adamczyk, 2016). Therefore, an assumption made on the swing leg angular velocity will not affect the applicability of the model to the freezing problem as freezing is related to the inability of the legs to propel the CoM forward in the case of walking. Hence, swing leg angular velocity is reset to ${\omega }_2^0$ in every step. Furthermore, the low dimensional model helps to avoid making any assumptions on the initial angular velocity of the swing leg ${\omega }_2^0$ as it doesn't involve a swing leg.

3. Analysis of the Reduced System

When considered independently of the swing leg, the dynamics has states corresponding only to the stance leg, i.e., $x={\left[{\theta }_1,{\omega }_1\right]}^T$. The terms defining the Equation (1)–(2) for the inverted pendulum case are given below.

$\begin{array}{ll}q(.):=\left[\begin{array}{c}{\omega }_1\\ \frac{{\zeta }_2(t)-{\zeta }_1(t)}{m_1{l}^2}+\frac{g{\theta }_1}{l}\end{array}\right]& (19)\end{array}$

$\begin{array}{ll}S:=\left\{{\left[\begin{array}{c}{\theta }_1,{\omega }_1\\ \end{array}\right]}^T\in \chi |g_{reset}\left({\left[\begin{array}{c}{\theta }_1,{\omega }_1\end{array}\right]}^T\right)=0\right\}& (20)\end{array}$

$\begin{array}{ll}{g}_{reset}(.):={\theta }_1-{\theta }_{reset}& (21)\end{array}$

$\begin{array}{ll}\Delta (.):={\left[\begin{array}{c}-{\theta }_1,{\omega }_1\text{cos}({\theta }_h)\end{array}\right]}^T& (22)\end{array}$

As the inverted pendulum model is analyzed independently, θ₁, m₁, ${\theta }_1^0$ and ${\omega }_1^0$ will be referred here as θ, m, θ₀ and ω₀, respectively. These equations are solved to produce the motion trajectory during the stance phase of the stepping cycle. The sequence of model evolution is depicted in Figure 4, with the beginning and end of the stance positions, initial angular position (θ₀), initial angular velocity (ω₀) and the angle at reset (θ_reset). Step length is defined to be equal to |θ_reset| where |.| denotes the absolute value.

FIGURE 4

Figure 4. Gait cycle for the low dimensional (inverted pendulum system) system is shown. Terms ${\left[{\theta }_0,{\omega }_0\right]}^T$ and ${\left[{\theta }_0,{\omega }_0\right]}^{T-}$ indicate the initial and final angular position and velocity of the leading/stance leg, respectively, at the beginning and end of the stance phase. θ_reset is the angular displacement at which the angle is reset. The CoM (center of mass) which is assumed to be acting as a point mass at the tip of the inverted pendulum is shown as a circle. States before and after heel strike is indicated with “-” and “+” superscripts, respectively.

3.1. Gait Cycle

The proposed model considers only the “stance” phase of the gait cycle. Therefore, “gait cycle” in this study has been defined as the process, where the model states evolve from an initial condition of a step (“double support phase”) until the reset condition (where the heel of the swing leg is assumed to collide with the ground or “heel strike condition”) is met and the initial condition of the next step is computed. Here, the state of the system moves through three different states (beginning of the “stance” (double support), end of the “stance” (before collision of the contra-lateral leg), heel strike (after collision of the contra-lateral leg) whose notations are given below (Equation 23) ^1,2.

$\begin{array}{ll}{\left[{\theta }_0,{\omega }_0\right]}^T\mapsto {\displaystyle \underset{\text{Before collision}}{{\displaystyle \underbrace{{\left[{\theta }_0^{-},{\omega }_0^{-}\right]}^T}}}}\mapsto {\displaystyle \underset{\text{After collision}}{{\displaystyle \underbrace{{\left[{\theta }_0^{+},{\omega }_0^{+}\right]}^T}}}}={\left[{\theta }_1,{\omega }_1\right]}^T& (23)\end{array}$

Here ${\left[{\theta }_0,{\omega }_0\right]}^T$, ${\left[{\theta }_0^{-},{\omega }_0^{-}\right]}^T$, ${\left[{\theta }_0^{+},{\omega }_0^{+}\right]}^T$ correspond to the states at the initiation of the step, states at the end of the flow “immediately” before collision, and states “immediately” after collision, respectively. The states immediately after collision form the initial condition for the next step ${\left[{\theta }_1,{\omega }_1\right]}^T$. The superscripts (“-,” “+”) need not indicate the relative sizes of the states but the chronological order in which they appear, that is “-” superscript represent before collision variables and “+” represents after collision. But it should be noted that the transformation from ${\left[{\theta }_0^{-},{\omega }_0^{-}\right]}^T$ to ${\left[{\theta }_0^{+},{\omega }_0^{+}\right]}^T$ happens instantaneously in the model. Counterclockwise angles are defined as positive. In a typical walking simulation this results in θ⁻ < 0 < θ⁺. That is, the stance phase ends at a negative value for the angle and resets to a positive value before beginning the next stance phase.

Subscripts (indicating the step number) will be dropped from before and after collision state symbols when the step number is not relevant for the derivation (Equation 24). The same superscript will be used while referring to other parameters which change during collision ³.

The states [θ, ω]^T evolve as a function of time except at the collision point, where the same time point maps to two different state values. ⁴

$\begin{array}{ll}{\left[{\theta }_n,{\omega }_n\right]}^T\mapsto {\displaystyle \underset{\text{Before collision}}{{\displaystyle \underbrace{{\left[{\theta }^{-},{\omega }^{-}\right]}^T}}}}\mapsto {\displaystyle \underset{\text{After collision}}{{\displaystyle \underbrace{{\left[{\theta }^{+},{\omega }^{+}\right]}^T}}}}={\left[{\theta }_{n+1},{\omega }_{n+1}\right]}^T& (24)\end{array}$

3.2. Heel Strike Condition

A heel strike is defined as the state at which the swing leg (trailing leg) collides with the ground. This is modeled using an appropriate reset condition. At heel strike, both the angle and angular velocity are reset from the “before collision” to “after collision” state as described in Equation (24). The collision of the swing leg (trailing leg) at heel strike is modeled to be inelastic with angular momentum conserved. Therefore, the magnitudes of the angular momentum about the point of collision after and before collision are equated in the following way to generate the transition rule for angular velocity

$\begin{array}{}\\ \\ lm(v^{+}\text{ sin}(90))=lm(v^{-}\text{ sin}(90+{\theta }_h))& (25)\end{array}$

$\begin{array}{ll}lm(l{\omega }^{+})1=lm(l{\omega }^{-})\text{ cos}({\theta }_h)& (26)\end{array}$

$\begin{array}{ll}{\omega }^{+}={\omega }^{-}\text{ cos}({\theta }_h)& (27)\end{array}$

at the n^th iteration(step)

$\begin{array}{ll}{\omega }_{n+1}={\omega }_n^{+}={\omega }_n^{-}\text{ cos}({\theta }_h)& (28)\end{array}$

where θ_h is the hip angle. The angle, on the other hand, will be reset from θ⁻ to −θ⁻. This results in the following transition rules at ${\theta }^{-}={\theta }_{reset}$

$\begin{array}{ll}{\omega }_n^{+}={\omega }_n^{-}\text{ cos}({\theta }_h)& (29)\end{array}$

$\begin{array}{ll}{\theta }_n^{+}=-{\theta }_n^{-}& (30)\end{array}$

Rearranging we obtain Δ(.) as

$\begin{array}{ll}\Delta \left({\left[\begin{array}{c}\theta ,\stackrel{.}{\theta }\\ \end{array}\right]}^{T-}\right)={\left[\begin{array}{c}-\theta ,\stackrel{.}{\theta }\text{ cos}({\theta }_h)\\ \end{array}\right]}^T.& (31)\end{array}$

3.3. Analytical and Numerical Solution of the Equations of Motion

The differential equation Equation (1), was solved using the definition of the vector field given in Equation (19) analytically to obtain the flows given below.

$\begin{array}{l}\theta (t)=f_{\theta }(t,{\omega }_0,{\theta }_0):=N_1/D_1\\ \omega (t)=f_{\omega }(t,{\omega }_0,{\theta }_0)=\frac{d}{dt}(f_{\theta }):=N_2/D_2(32)\end{array}$

where N₁, N₂, D₁, D₂ are given in Supplementary Section 1. The analytical solution is intended to be used for the bifurcation analysis as numerical solutions may not always be able to detect the chaotic behavior (Lozi, 2013). The analytical solution will therefore be used to generate the discrete map governing the motion in the following sections. A numerical solution of the Equations (1)–(2) using definitions given in Equations (19)–(22) with the appropriate reset conditions (in the physiological range) are solved to show the freezing behavior and dynamics in the phase plane. PD subjects freeze intermittently, and the amount of time the subject walks until the freeze is an important measure to quantify the transient walking behavior. A simulation for a 10 s- window is carried out for different values of the parameters τ_l and τ_r (for a constant initial condition). The total time for which transient walking behavior occurred is computed numerically as a function of the parameters τ_l and τ_r. Numerical methods are also used in solving boundary value problems to gain further insights into the system as given in the remark 3.

REMARK 3. The parameters τ_l and τ_r determines the amount of energy supplied to the system apart from gravity. To understand how they influence the kinetic energy of the system, the difference in speed between the initial and final states are compared for the boundary value problem with boundary conditions θ₀ = 0 rad. and θ_0.5 = −0.1 rad. with definitions given in Equations (19)–(22) unchanged. Here the boundary conditions are chosen from the physiological range.

The quantities τ_l, τ_r, ω, ϕ, θ_reset and step length has units N m, N, rad. s⁻¹, rad., rad., and rad., respectively, when not specified.

3.4. Derivation of a Map to Describe Successive Stance Phases

The evolution of the flow (given by Equation 32) is terminated when the swing leg meets the ground. In other words, when there is sufficient energy in the system for forward motion, there exists a “reset time” T(θ₀, ω₀) such that f_θ(T(θ₀, ω₀), ω₀, θ₀) = A_reset(θ₀). Here, A_reset(.) is a function of the joint angle and the ground that determines the angle of the stance leg while the foot strikes the ground. The arguments associated with the reset time T(θ₀, ω₀) will be dropped and will be referred to as T from here on. Accounting the transition rules in Equation (29) for reset and conservation of angular momentum,
one defines⁵

$\begin{array}{ll}{\theta }_1={\theta }^{+}(T):=-f_{\theta }(T,{\omega }_0,{\theta }_0)& (33)\end{array}$

$\begin{array}{ll}{\omega }_1={\omega }^{+}(T):=f_{\omega }(T,{\omega }_0,{\theta }_0)\text{ cos}({\theta }_h)& (34)\end{array}$

Following an induction hypothesis, for an arbitrary initial condition (θ_n, ω_n) the map is

$\begin{array}{ll}{\theta }_{n+1}=-f_{\theta }(T,{\omega }_n,{\theta }_n)& (35)\end{array}$

$\begin{array}{ll}{\omega }_{n+1}=f_{\omega }(T,{\omega }_n,{\theta }_n)\text{ cos}({\theta }_h)& (36)\end{array}$

The following definitions are made to make the notations compact for further analysis

$\begin{array}{ll}{\theta }_{n+1}={\tilde{f}}_{\theta }(T,{\omega }_n,{\theta }_n):=-f_{\theta }(T,{\omega }_n,{\theta }_n)& (37)\end{array}$

$\begin{array}{ll}{\omega }_{n+1}={\tilde{f}}_{\omega }(T,{\omega }_n,{\theta }_n):=f_{\omega }(T,{\omega }_n,{\theta }_n)\text{ cos}({\theta }_h)& (38)\end{array}$

where ${\tilde{f}}_{\omega }(T,{\omega }_n,{\theta }_n,\varphi )$ and ${\tilde{f}}_{\theta }(T,{\omega }_n,{\theta }_n,\varphi )$ are T parametrized family of maps for (ω_n, θ_n) ↦ (ω_n+1, θ_n+1).

To investigate the condition of same step lengths and to generate a 1D map for further evaluation, A_reset(θ, t) is set to be θ_reset. Here θ_reset is an arbitrary angle in the physiological range at which the swing leg meets the ground. Then for an intermediate step, (when there is sufficient energy to move forward) there exists a $\tilde{T}$ s.t. ${\tilde{f}}_{\theta }(\tilde{T}({\theta }_n,{\omega }_n),{\omega }_n,{\theta }_n)={\theta }_n=-{\theta }_{reset}$. When there is not enough energy and therefore momentum to move forward, the model behavior is defined as freezing.

To find the $\tilde{T}$ at which θ_n maps to itself the following minimization problem is solved ⁶ using Newton's method (Wolfram Research Inc., 2019). This detects implicitly the time at which the swing leg collides with the ground.

$\begin{array}{ll}\tilde{T}({\theta }_n,{\omega }_n):={\displaystyle \underset{T}{\text{arg min}}}({\tilde{f}}_{\theta }(T,{\omega }_n,{\theta }_n,\varphi )-{\theta }_n)& (39)\end{array}$

Substituting $\tilde{T}$ from Equation (39) in Equation (38) the following map is obtained.

$\begin{array}{l}{\omega }_{n+1}={\tilde{f}}_{\omega }(\tilde{T}({\theta }_n,{\omega }_n),{\omega }_n,{\theta }_n)\\ :=H({\omega }_n)\text{ when }{\theta }_n={\theta }_0\forall n\in \mathbb{N}\end{array}$

When used as a 1D map,

$\begin{array}{ll}{\omega }_{n+1}:={\tilde{f}}_{\omega }({\omega }_n)& (40)\end{array}$

The argument $\tilde{T}$ in the function will be dropped from here-on. The function $\tilde{T}$ acts on the same input ω_n and θ_n = θ₀. This map has been analyzed to show the freezing behavior and variabilities associated with PD walking. The map has been analyzed for a particular parameter value to show the presence of horseshoe in the Supplementary Section 2.

REMARK 4. Equations (1)–(2) represent a general hybrid system. When the hybrid bipedal system's solution is sought these equations are solved using the definitions given in Equations (15)–(16); and, Equations (19)–(22) are used for hybrid inverted pendulum system.

4. Results

Numerical simulation of the PD gait and associated freezing behavior is described in this section. The change in the angular velocity from negative to zero is a property of any solution containing freezing by definition. Typically in this model, the angular velocity changes to a positive value under the action of gravity during a freeze. The effect of variation of the parameters τ_l, τ_r, ϕ, θ_reset are also investigated. The work aims to show that, the two opposing torques modeled to be generated from the plantar flexors could elicit freezing and chaotic behavior. The ability of these torques to generate freezing behavior has been shown first in a simplified biped model described using Equations (15)–(16), and then in the inverted pendulum model generated by Equations (19)–(22). As argued previously the inverted pendulum dynamics sufficiently captures the PD walking scenario. The results are presented in the sequel to support this hypothesis. Also, walking is the process of moving the CoM by pushing the stance leg forward, and the inverted pendulum model helps to study the effect of the stance leg independent of other variables. A range of values for the constants τ_l, τ_r, and ϕ have been analyzed, such that the trend in behavior is clear to understand. The range in which the behavior of the map ${\tilde{f}}_{\omega }$ changes the number of periodic orbits from “0” to “more than one” in lower absolute value of angular velocity conditions, is given Table 1. Simulations are carried out to understand the behavior of the system over and above this range. But it may be noted that the maximum value of the torque for l = 0.6 m., θ_reset = −0.1 rad. is approximately 0.23|τ_r| N m and 2|τ_l| N m in forward and backward directions, respectively. Hence, in this case, forward pushing plantar flexors has to generate 8.7 times the “premature activation of plantar flexors” to nullify the effect if the phase is matched exactly. Physiologically the minimum value of these torques is zero and the maximum is subject-specific.

TABLE 1

Table 1. Summary of qualitative behavior of the map.

4.1. Freezing in a Biped Model

The hybrid system (Equations 1–2) defined by the Equations (15)–(16) are simulated numerically and the results are shown in Figures 5A,B. The figure shows normal walking for the first few steps and then freezing afterwards (highlighted). The gradual reduction in step length observed experimentally prior to freezing (Nutt et al., 2011) is also observed in the model.

FIGURE 5

Figure 5. Results of numerical simulation of the biped during freezing. Parameters are chosen to be τ_l = 2.3 N m, τ_r = 15.74 N, ϕ = −π/2 rad. and initial conditions ${\theta }_2(0)=-0.1\text{rad}.,{\theta }_1(0)=0.1\text{rad}.,{\omega }_2(0)=2\text{rad}.s^{-1},{\omega }_1(0)=-0.6\text{rad}.s^{-1}$. The change in the θ_reset in every step and gradual reduction in θ₁ nearer to a freezing event is evident. (A) Simulated time series of the states θ₁ and ω₁. Region of slowing down and freezing is highlighted (B) Shows numerical simulation of biped in the phase plane.

There are two dissipative forces in this model; they are, the opposing torques due to the plantar flexors and the dissipation at the heel strike. Long-range walking will be achieved when the speed gain in every step compensate these two effects. The dissipative effect of the heel strike can't be controlled by the neuromuscular system, but the effect of plantar flexors can. Figures 6A,B illustrates the effect of the plantar flexors in this regard.

FIGURE 6

Figure 6. (A) A simulation was carried out for a time window of 10 s, and the total time walked before freezing was computed. This time is plotted as a function of τ_l and τ_r. The initial conditions are set to be ${\theta }_2(0)=-0.1\text{rad}.,{\theta }_1(0)=0.1\text{rad}.,{\omega }_2(0)=1\text{rad}.s^{-1},{\omega }_1(0)=-0.6\text{rad}.s^{-1}$ and ϕ = −π/2 rad. The yellow region forms the optimal region of the parameters τ_l, τ_r where walking is achieved. (B) Contour plot of the difference in speed(|ω(0)| − |ω(0.5)|) as a function of τ_l and τ_r determined by solving a BVP numerically with ${\theta }_0^0=0$ and ${\theta }_1^{0.5}=-0.1\text{rad}.$ for the biped model with the other initial conditions being ${\theta }_2^0=-0.1\text{rad}.,{\omega }_2^0=1$ (and ϕ = −π/2 rad.). This illustrates a decrease in the speed and therefore kinetic energy when there is an increase in τ_l and a decrease in τ_r.

A simulation was carried out for a 10 s window and the time difference between, the start of the simulation and the time of the last heel strike before the freezing event (as defined in remark 1), has been computed. This is shown as a function of the parameters τ_l and τ_r in Figure 6A. The blue shades indicate eventual freezing and shorter walking time and the yellow region is the safer non-freezing region. There is an intermediate region of parameter τ_l and τ_r in which the walking happens without freezing in the 10 s window. A higher value of τ_l necessitates a higher τ_r for walking. But a very high τ_r doesn't necessarily produce balanced walking as it can result in a lack of coordination between the swing leg and the stance leg. Although the initial value problem (IVP) in Figure 6A and boundary value problem (BVP) in Figure 6B can't be directly compared, they show analogous qualitative results. That is, to achieve the same speed gain [or kinetic energy (KE) gain] a higher τ_l demands a higher τ_r. A key aspect of PD freezing is, therefore, the inability of the two plantar flexors to coordinate to produce the required energy. From an energy point of view, the role of the swing leg is mainly in the generation of ankle push of force. In the following sections, the dynamics of the stance leg is studied independently using an inverted pendulum model, reducing the role of the swing only as a supplier of the ankle push-off force.

4.2. Freezing in an Inverted Pendulum Model

Freezing is defined as the condition where there is no more forward motion of the leg. Numerical simulation of such a scenario in the inverted pendulum model is shown in Figure 7A where there is a freezing episode after 18 s. A gradual reduction in step length observed in the biped model translate to the increased time taken in making the final few steps before freezing. The simulation in the phase plane for the last three steps is shown in Figure 7B. The dissipative torques due to the opposing plantar flexors act in the same way in the case of the inverted pendulum model. Figures 8A,B illustrate this similarity, where, an increased τ_r generates higher speed gains and, an elevated τ_l results in lower speed gains and lower total walk times. This is because increasing parameter τ_r heightens the forward ankle push-off while larger τ_l amplifies the dissipative torque. A key difference between the inverted pendulum model and the biped model is that a higher τ_r will not result in an imbalance in the former as there is no swing leg in that model, while there is a lack of balance in the latter. The contour plot of the speed differences as a function of τ_l and τ_r is shown in Figure 8B. The figure shows that a higher value of τ_l and a lower value of τ_r result in negative speed gain (reduction in KE). Numerical simulation of the total time of the walk, defined as the difference between the time in which the first step is taken and the last step before freezing in 10 s is shown in Figure 8A. More than 9 s of walking is indicative of the fact that there is no freezing in that parameter range in that time frame. A higher τ_r and lower τ_l results in better walking performance as in the case of a biped. Therefore, energetically, PD related behavior that is of interest is analogous in the case of inverted pendulum and biped model. Therefore, the analytical solution of the inverted pendulum model is investigated further to understand the consequence of the change in parameters τ_l, τ_r and ϕ. These parameters are controlled by the neural system while others such as mass of the body and length of the legs are not. The quantity θ_reset differentiates the inverted pendulum model from the biped model. Hence, the effect of this parameter is also studied.

FIGURE 7

Figure 7. Simulation of the inverted pendulum dynamics for the parameter values $\varphi =-\pi /2\text{rad}.,{\omega }_0=-1\text{rad}.s^{-1},{\tau }_l=2\text{N m},{\tau }_r=11\text{ N}$. Freezing occurs after 18 s. (A) States θ and ω as a function of time (10–20 s). Region of slowing down and freezing is highlighted. (B) Numerical simulation of the inverted pendulum states in the phase plane for 16–20 s.

FIGURE 8

Figure 8. (A) A simulation was carried out for a time window of 10 s, and the total time of walking before freezing was computed. This time is plotted as a function of τ_l and τ_r. The angle is reset when θ(t) = −0.1 rad. and ϕ = −π/2 rad. The colors indicate the duration of the walk (see the legend). (B) Contour plot of difference in speed (|ω(0)| − |ω(0.5)|) as a function of τ_l and τ_r determined by solving a BVP with θ₀ = 0 rad. and θ_0.5 = −0.1 rad. (and ϕ = −π/2 rad.) in Equation (1)–(2) using definitions in Equation (19)–(22). This illustrates a decrease in the speed and therefore KE when there is an increase in τ_l and a decrease in τ_r.

4.3. Parameter Exploration of the Inverted Pendulum: Study of the Map ${\tilde{f}}_{\omega }$

The neural control on the muscles alters the magnitude and the phase of the control signals. Exploration of the parameters τ_l, τ_r, and ϕ, therefore, reflects the effect of the neural control on walking dynamics. One of the hypotheses that are investigated through the model is that of the generation of variability through the premature activation of the plantar flexors. We have quantified the phase difference of the “premature” activation using the parameter ϕ in the model. Figure 9 shows the bifurcation diagram of the parameter ϕ in the range 0 to −2π for constant values of τ_l and τ_r. A period-doubling route to chaos can be observed when ϕ is varied between −5π/8 and −π/4. The map ${\tilde{f}}_{\omega }$ is iterated for 500 walking cycles and the last 50 walking cycles are used to compute the equilibrium points. The Feigenbaum bound is found to be at ϕ = −1.37 rad. at which walking becomes fully chaotic. This indicates that the premature activation (or lack of coordination between the muscles) can generate highly variable behavior in the system despite deterministic neural signals. The region of chaotic ϕ is sandwiched between the periodic orbits and freezing region. This suggests a higher variability in walking likely arising from a shift in ϕ (early activation of plantar flexors) must be treated with caution. Figure 9 shows the presence of chaos in the system for carefully selected parameter values. Its presence and stability are illustrated for other parameters values and initial conditions using a set of maps in Figures 10–12B and bifurcation diagrams in Figures 13A–D. A summary of the insights obtained from the maps are given in the Table 1. The presence of a period 3 orbit in a one-dimensional map is indicative of other periodic orbits and chaos. The presence of horseshoe in any of the period 1, 2,..., n maps also indicates chaos. An illustration of the presence of horseshoe for a set of parameter values is given in Supplementary Section 2. The intersection of the ${\tilde{f}}_{\omega }^1$, ${\tilde{f}}_{\omega }^2$, ${\tilde{f}}_{\omega }^3$ maps with ω_n = ω_n+1 indicate period 1, 2, 3 orbits, respectively. Figures 10–12B illustrate how the maps change with respect to the change of parameters.

FIGURE 9

Figure 9. Stable ω is shown as a function of the parameter ϕ for ${\omega }_0=-0.433\text{rad}.s^{-1},{\tau }_l=5\text{N m},{\tau }_r=35\text{ N}$. The Feigenbaum bound was found to be at ϕ = −1.37 rad. where walking becomes fully chaotic. The period-doubling cascade has been highlighted and enlarged. This chaotic region forms only a small part of the overall parameter space of ϕ. This region is sandwiched between the walking and the freezing regions indicated in red.

FIGURE 10

Figure 10. Maps obtained by varying the parameter τ_l and fixing τ_r = 35 N, ϕ = −1.57 rad., θ_reset = −0.1 rad. The black, green, and red curves represent ${\tilde{f}}_{\omega }^1$, ${\tilde{f}}_{\omega }^2$, ${\tilde{f}}_{\omega }^3$, respectively, and the ω_n = ω_n+1 is shown in blue. The curves intersect the blue line at a higher absolute value of angular velocity forming an attractor, this is not shown in the figure. There are no periodic orbits for the low velocity regimes for τ_l = 0 − 2 N m but they appear afterwards. Units: τ_l, τ_r, ω, ϕ, θ_reset and step length has units N m, N, rad. s⁻¹, rad., rad., and rad., respectively, when not specified.

FIGURE 11

Figure 11. Varying the parameter ϕ and fixing τ_r = 35 N, τ_l = 5 N m, θ_reset = −0.1 rad. The black, green, and red curves represent ${\tilde{f}}_{\omega }^1$, ${\tilde{f}}_{\omega }^2$, ${\tilde{f}}_{\omega }^3$, respectively, and the ω_n = ω_n+1 is shown in blue. The curves intersect the blue line at a higher absolute value of angular velocity forming an attractor, this is not shown in the figure. Creation of the periodic orbits and its coexistence is observed. Units: τ_l, τ_r, ω, ϕ, θ_reset and step length has units N m, N, rad. s⁻¹, rad., rad., and rad., respectively, when not specified.

FIGURE 12

Figure 12. Variations of the parameters τ_r and θ_reset are depicted in Figures 12A,B, respectively. Units: τ_l, τ_r, ω, ϕ, θ_reset and step length has units N m, N, rad. s⁻¹, rad., rad., and rad., respectively, when not specified. (A) Varying the parameter τ_r keeping τ_l = 5 N m, ϕ = −1.57 rad., θ_reset = −0.1 rad. fixed. The black, green, and red curves represent ${\tilde{f}}_{\omega }^1$, ${\tilde{f}}_{\omega }^2$, ${\tilde{f}}_{\omega }^3$, respectively, and the ω_n = ω_n+1 is shown in blue. The curves intersect the blue line at a higher absolute value of angular velocity forming an attractor, this is not shown in the figure. Increasing τ_r has an analogous behavior as decreasing τ_l. (B) Varying the parameter θ_reset keeping τ_r = 35 N, ϕ = −1.57 rad., τ_l = 0.5 N m. The black, green and red curves (overlapped) represent ${\tilde{f}}_{\omega }^1$, ${\tilde{f}}_{\omega }^2$, ${\tilde{f}}_{\omega }^3$, respectively, and the ω_n = ω_n+1 is shown in blue. The curves intersect at a high absolute value of angular velocity. The unstable region moving to the higher absolute value of angular velocities (moving to the left) can be observed while |θ_reset| is increased.

FIGURE 13

Figure 13. Period one orbits are shown by varying τ_l and τ_r for two different values of initial angular velocities in (C,D). Period one orbits found by varying τ_l and ω₀ is shown in (A). Period one orbits found by varying τ_r and ω₀ is shown in (B). The parameter values used are given in the respective figures. The green region shows the stable region where ${\stackrel{.}{\stackrel{~}{f}}}_{\omega }<1$. The stable periodic regions are the ones where the green region overlap the curves of the periodic orbits. Units: τ_l, τ_r, ω, ϕ, θ_reset and step length has units N m, N, rad. s⁻¹, rad., rad., and rad., respectively, when not specified. (A) ϕ = −π/2 rad., τ_r = 40 N, θ_reset = −0.1 rad. (B) ϕ = −π/2 rad., τ_l = 5 N m, θ_reset = −0.1 rad. (C) ϕ = −π/2 rad., ω₀ = −0.4 rad. s⁻¹, θ_reset = −0.1 rad. (D) ϕ = −π/2 rad., ω₀ = −1 rad. s⁻¹, θ_reset = −0.1 rad.

Variation of the parameter τ_l or the magnitude of premature activation (as ϕ is set to -1.57 rad.) results in a set of rich dynamic behaviors as shown in Figure 10. The presence of the periodic orbits starts appearing approximately around τ_l ≈ 3 N m, where, the maps tangentially intersect the ω_n = ω_n+1 line. The intermittency thus generated could elicit a period of slow walking (as ω_n and ω_n+1 are less than -0.5 rad. s⁻¹) as observed in PD. The period 3 orbits are generated upon a further increase in τ_l. As can be seen from the maps in Figures 10–12B, a higher initial value of ω_n (e.g., ω_n > 0.45 rad. s⁻¹ for τ_l ≈ 3 N m) results in a further increase in ω_n+1 and gets attracted toward the periodic orbit of higher absolute value of angular velocity. This explains how swaying back and forth helps the PD patients in getting out of a freeze. Increasing τ_r results in almost opposite behavior as that of τ_l (Figure 12A). Varying ϕ can result in chaotic behavior as shown in Figure 9, and, Figure 11 indicates the variation in the maps which leads to this behavior. The neural control of the activity of plantar flexors is not explicitly modeled here. However, coming out of freeze could be the result of an increase of τ_r or decrease τ_l or increased initial absolute angular velocity generated by swaying. A low absolute value of angular velocity (voluntary or involuntary) or decrease of τ_r or increase of τ_l results in freezing (angular velocity moving to the region where ω_n = 0 rad. s⁻¹). This explains the higher chances of freezing episodes even when the subject reduces the velocity (voluntarily/involuntarily) near narrow passages. Increase in the step length or |θ_reset| results in freezing at comparatively higher absolute angular velocities (Figure 12B). But it may be noted that, typically, an increased step length is also associated with an increased absolute angular velocity due to inertia and therefore could be beneficial. There is likely an optimum step length for every subject as there is a trade-off between fatigue and initial angular velocity, which warrants further study.

4.4. Bifurcations of the One Dimensional System for the Inverted Pendulum Model

Even though for most of the regions, the slope of the map in relation to the ω_n = ω_n+1 can be identified visually, the stability of the system is not explicitly studied in the previous section. The contour of ${\tilde{f}}_{\omega }^n(x,{\tau }_l,{\tau }_r,\varphi )=x$ for n=1 and 3 are plotted for variation in parameters in Figures 13, 14, respectively. The stability is computed by taking the derivatives (numerically) for the maps described in Figures 10–12B. These contours show how the points of intersection with ω_n = ω_n+1 line for the maps shown in Figures 10–12B change upon variation in parameters.

FIGURE 14

Figure 14. Constant parameters used are ϕ = −π/2 rad., ω₀ = −0.4 rad. s⁻¹, and θ_reset = −0.1 rad. The stable period 3 region is shown in green. The intersection of the period 3 orbits (in blue) and the stable regions form the region of stable period three orbits. The presence of period 3 orbits implies orbits of all other periods and therefore chaos.

Period one orbits are the normal walking cycles. The existence of these orbits in both low and high angular velocity conditions and different parameter variations are shown in Figure 13. In Figure 13A, two fixed points comes closer to each other and completely vanish for high values of τ_l resulting in a complete lack of periodic solutions. Typically walking could be ascribed to the stable region for periodic orbits, but, when ω₀ is lower, and τ_l is non-zero, another periodic point emerges in the low-velocity regimes. This, therefore, results in slow-walking regions which under perturbations could lead to freezing. Also, at low-velocity regimes, the region is discontinuous and unstable for small perturbations of the parameter values or initial conditions. The stable periodic orbit moves to lower absolute value of angular velocities as τ_l is increased and eventually disappears.

The behavior observed while decreasing τ_r is analogous to increase in τ_l. Figure 13B illustrates how changing τ_r and ω₀ results in creation/destruction of the periodic orbits. It can be seen that at a sufficiently low value of τ_r the periodic orbit disappears. A higher value of τ_r results in the separation of the periodic orbits resulting in higher stable walking angular speeds. A similar behavior could be observed while decreasing τ_l in Figure 13A.

Initial angular velocity plays a major role in the behavior of the system. The effect of neural control parameters τ_l and τ_r in generating periodic behavior has been illustrated for lower and higher absolute angular velocity conditions in Figures 13C,D, respectively. In Figure 13C, the periodic orbit appears stable only for a tiny fraction of the parameters space. This is due to the highly discontinuous map shown previously. Conversely, at higher initial angular speeds the period one orbit is stable as shown in Figure 13D. It can be seen that an increase in τ_l moves the periodic orbit into an unstable region resulting in the possibility of a freeze. The presence of these orbits could only be seen in the low-velocity regions of the maps. Orbits of minimal period three indicate chaos and the presence of every other periodic orbits (Glendinning, 1994). The period 3 orbits for the variation of the parameters τ_l and τ_r is shown in the Figure 14. The Period 3 orbit is shown in blue and the period one in yellow.

5. Discussion, Summary, and Future Work

Freezing of gait results from a complex set of interacting physiological systems which consist of the brain, spinal cord, musculoskeletal system and external disturbances (Nutt et al., 2011). The model explains how a lack of coordination between central pattern generators of the plantar flexors of the leading leg and trailing leg (Nieuwboer et al., 2004) could lead to freezing and variability of walking.

A model of the torques generated by the plantar flexors acting on the stance leg has been proposed, and its effect on a biped and a reduced inverted pendulum model has been studied. The pattern of freezing observed in the model matches well with the behavior observed experimentally⁷ in Nutt et al. (2011) and Figure 5A). The equilibrium point description (Feldman, 1986; Sainburg, 2015) of the control of the muscles is avoided here and instead, we opted for an explicit control signal. However, variabilities in the “torque-length-characteristics” (Feldman, 1986) for a particular (set of) equilibrium point (points) can generate torques required for motion. Therefore, a parallel between the equilibrium point hypothesis of postural balance and our model can be drawn if the torques prescribed in the model are assumed to be the result of variabilities in the “torque-length-characteristics.”

Chaotic regions are observed to be closer to those regions where freezing ensues. In the inverted pendulum model, these regions show up only at low absolute angular velocity initial conditions. This may explain why freezing is a “rarely” occurring intermittent condition. This also may explain why freezing happens near obstacles or narrow paths where the subject voluntarily slows down to reduce the probability of collision. Obstacles could be either perceived or real. Hence, even though the pattern of freeze remains the same the causes could be varied. It might even be possible that the control of τ_l is driven by perceived obstacles or anxiety about the consequence of freezing (Ehgoetz Martens et al., 2014; Martens et al., 2016). Increase in τ_l, therefore, could be thought to be indirectly influenced by anxiety and perceived obstacles. However, this hypothesis warrants further experimentation.

Varying the parameter step length which controls the stride length is observed to affect the maps and therefore the freezing regions. The results indicated that keeping the steps closer to each other such that |θ_reset| is minimized, is safer for the PD patient. The stability of the period 2 & 3 orbits are highly sensitive to small variations of parameters (ϕ, τ_r, τ_l) which are proposed to be the reason for sporadic variabilities in gait seen in PD subjects. We also hypothesize that stable low absolute angular velocity regions of the state space for some parameter values form a “cantor set” and necessitates further study.

It can be speculated that a reason for the observed help of auditory/sensory cues (Rochester et al., 2005; Young et al., 2014; Amini et al., 2019) in reducing instances of freezing, is by indirectly forcing PD patients to make shorter steps with lesser variability, thus reducing the possibility of moving into the freezing region of walking. Variability in the walking times observed in the Inverted pendulum model translates to variability in step lengths in the biped model. Biped model shows a more complicated dependence on the parameters to eventual freezing (Figure 6A). This dependence is also a function of the initial conditions and could be investigated further in future work along with detailed bifurcation analysis.

The CPG activity is controlled by feedback mechanisms with delays, noise and input from the brain (which in turn is affected by different factors, including emotional state). The ground and other environmental conditions also play a role in walking. These variabilities are not accounted for in our model, which represents a limitation of the study. Like any other studies which are based on a mathematical model and numerical simulations, our results and conclusions also might not necessarily represent the entire spectrum of patients. Further extensive patient-based studies are to be performed prior to use of these ideas for the treatment of PD gait. As future work, a more detailed model is planned to include these variabilities. The future models will be compared with the simpler versions to understand the minimum set of variables generating the abnormal walking behavior. The key aspects explained using the proposed model can be summarized as follows

1. The higher variability in PD patients could be the result of parameters being closer to the point of chaos. A further change of the parameters can result in freezing. Therefore, increased variability should be looked at with caution (clinically) and should be treated to reduce it. The difficulty in the prediction of freezing also owes to the horseshoe near the freezing regions.

2. The pattern of reducing the step-sizes before freezing has been shown to be the result of slowing down (Figure 5A). Voluntary/involuntary reduction in angular velocity (in absolute terms) near the obstacles makes the subject more susceptible to freezing and highly irregular walking.

3. One plausible reason why sensory cues such as auditory or visual cues help in freezing is by reducing step lengths. The proposed model shows that the reduction in step length helps in reducing freezing episodes at lower absolute value of angular velocity conditions as it moves the patient away from the freezing region. Further experimental study is needed to understand the clinical applicability.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation, to any qualified researcher.

Author Contributions

MP and PM have contributed in design, implementation, analysis and editing of the study and manuscript. KT-A has contributed to the design and analysis of the study and editing of the manuscript. MW has contributed to the design of the study and editing of the manuscript. All authors contributed to the article and approved the submitted version.

Funding

MP and KT-A have got funding support from EPSRC (Engineering and Physical Science Research Council) grant No EP/N014391/1 funding the Centre for Predictive Modelling in Healthcare at the University of Exeter. The EPSRC Centre for Predictive Modelling in Healthcare (CPMH) brings together a world-leading team of mathematicians, statisticians and clinicians with a range of industrial partners, patients and other stakeholders to focus on the development of new methods for managing and treating chronic health conditions using predictive mathematical models. This unique approach is underpinned by the expertise and breadth of experience of the Centre's team and innovative approaches to both the research and translational aspects. As such the research presented in this manuscript has benefitted by the environment offered by the CPMH. MP has spent considerable time as part of the CPMH participating in various seminars, conferences and other activities contributed toward PhD. KT-A is the deputy director of CPMH.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The authors would like to express their appreciation to Prof. Peter Ashwin for useful inputs and valuable comments. KT-A and MP gratefully acknowledge the financial support of the EPSRC via grant EP/N014391/1.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fbioe.2020.552635/full#supplementary-material

Footnotes

1. ^Here a ↦ b indicates state a “maps to” state b after some time t where t ≥ 0,

2. ^${\left[x_1^{-},x_2^{-}\right]}^T$ and ${\left[x_1,x_2\right]}^{T-}$ are used interchangeably, where, x₁ and x₂ are components of some vector

3. ^p⁻ and p⁺ refers to any parameter p before and after collision, respectively, in a particular step cycle.

4. ^θ(t) and ω(t) are multi-valued functions at the point of collision.

5. ^Even though the symbol θ(T) is used at the point of collision, it may be noted that this is a one to many mapping and therefore is not a single valued function in the traditional sense of the word.

6. ^The solution to this optimization problem may not always exist.

7. ^Here we are referring to the Figure 1 in Nutt et al. (2011). Source of the figure :https://pubmed.ncbi.nlm.nih.gov/21777828/#&gid=article-figures&pid=figure-1-uid-0

References

Albin, R. L., Young, A. B., and Penney, J. B. (1989). The functional anatomy of basal ganglia disorders. Trends Neurosci. 12, 366–375. doi: 10.1016/0166-2236(89)90074-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Alexander, G. E., and Crutcher, M. D. (1990). Functional architecture of basal ganglia circuits: neural substrates of parallel processing. Trends Neurosci. 13, 266–271. doi: 10.1016/0166-2236(90)90107-L

PubMed Abstract | CrossRef Full Text | Google Scholar

Amini, A., Banitsas, K., and Young, W. R. (2019). Kinect4fog: monitoring and improving mobility in people with Parkinson's using a novel system incorporating the Microsoft Kinect V2. Disabil. Rehabil. Assist. Technol. 14, 566–573. doi: 10.1080/17483107.2018.1467975

PubMed Abstract | CrossRef Full Text | Google Scholar

Aoi, S., Ohashi, T., Bamba, R., Fujiki, S., Tamura, D., Funato, T., et al. (2019). Neuromusculoskeletal model that walks and runs across a speed range with a few motor control parameter changes based on the muscle synergy hypothesis. Sci. Rep. 9, 1–13. doi: 10.1038/s41598-018-37460-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Dai, H., and Tedrake, R. (2013). “L2-gain optimization for robust bipedal walking on unknown terrain,” in 2013 IEEE International Conference on Robotics and Automation (Karlsruhe: IEEE), 3116–3123. doi: 10.1109/ICRA.2013.6631010

CrossRef Full Text | Google Scholar

Davie, C. A. (2008). A review of Parkinson's disease. Brit. Med. Bull. 86:109–127. doi: 10.1093/bmb/ldn013

PubMed Abstract | CrossRef Full Text | Google Scholar

Delp, S. L., Anderson, F. C., Arnold, A. S., Loan, P., Habib, A., John, C. T., et al. (2007). Opensim: open-source software to create and analyze dynamic simulations of movement. IEEE Trans. Biomed. Eng. 54, 1940–1950. doi: 10.1109/TBME.2007.901024

PubMed Abstract | CrossRef Full Text | Google Scholar

Duan, X., Allen, R., and Sun, J. (1997). A stiffness-varying model of human gait. Med. Eng. Phys. 19, 518–524. doi: 10.1016/S1350-4533(97)00022-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Ehgoetz Martens, K. A., Ellard, C. G., and Almeida, Q. J. (2014). Does anxiety cause freezing of gait in Parkinson's disease? PLoS ONE 9:e106561. doi: 10.1371/journal.pone.0106561

PubMed Abstract | CrossRef Full Text | Google Scholar

Fathizadeh, M., Mohammadi, H., and Taghvaei, S. (2019). A modified passive walking biped model with two feasible switching patterns of motion to resemble multi-pattern human walking. Chaos Solitons Fractals 127, 83–95. doi: 10.1016/j.chaos.2019.06.018

CrossRef Full Text | Google Scholar

Fathizadeh, M., Taghvaei, S., and Mohammadi, H. (2018). Analyzing bifurcation, stability and chaos for a passive walking biped model with a sole foot. Int. J. Bifurcat. Chaos 28, 1850113. doi: 10.1142/S0218127418501134

CrossRef Full Text | Google Scholar

Feldman, A. G. (1986). Once more on the equilibrium-point hypothesis (lambda model) for motor control. J. Mot. Behav. 18, 17–54. doi: 10.1080/00222895.1986.10735369

PubMed Abstract | CrossRef Full Text | Google Scholar

Flash, T., and Hogan, N. (1985). The coordination of arm movements: an experimentally confirmed mathematical model. J. Neurosci. 5, 1688–1703. doi: 10.1523/JNEUROSCI.05-07-01688.1985

PubMed Abstract | CrossRef Full Text | Google Scholar

Gede, G., Peterson, D. L., Nanjangud, A. S., Moore, J. K., and Hubbard, M. (2013). “Constrained multibody dynamics with python: from symbolic equation generation to publication,” in International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, Vol. 55973 (Portland: American Society of Mechanical Engineers). doi: 10.1115/DETC2013-13470

CrossRef Full Text | Google Scholar

Genadry, W., Kearney, R., and Hunter, I. (1988). Dynamic relationship between EMG and torque at the human ankle: variation with contraction level and modulation. Med. Biol. Eng. Comput. 26, 489–496. doi: 10.1007/BF02441916

PubMed Abstract | CrossRef Full Text | Google Scholar

Giladi, N., McDermott, M., Fahn, S., Przedborski, S., Jankovic, J., Stern, M., et al. (2001). Freezing of gait in PD prospective assessment in the datatop cohort. Neurology 56, 1712–1721. doi: 10.1212/WNL.56.12.1712

PubMed Abstract | CrossRef Full Text | Google Scholar

Glendinning, P. (1994). Stability, Instability and Chaos: An Introduction to the Theory of Nonlinear Differential Equations, Vol. 11. Cambridge University Press. doi: 10.1017/CBO9780511626296

CrossRef Full Text | Google Scholar

Goswami, A., Thuilot, B., and Espiau, B. (1996). Compass-Like Biped Robot Part I: Stability and Bifurcation of Passive Gaits. Research Report RR-2996, INRIA.

Google Scholar

Grizzle, J. W., Abba, G., and Plestan, F. (2001). Asymptotically stable walking for biped robots: analysis via systems with impulse effects. IEEE Trans. Automat. Contr. 46, 51–64. doi: 10.1109/9.898695

CrossRef Full Text | Google Scholar

Heremans, E., Nieuwboer, A., and Vercruysse, S. (2013). Freezing of gait in Parkinson's disease: where are we now? Curr. Neurol. Neurosci. Rep. 13:350. doi: 10.1007/s11910-013-0350-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Hof, A., and van den Berg, J. (1977). Linearity between the weighted sum of the EMGs of the human triceps surae and the total torque. J. Biomech. 10, 529–539. doi: 10.1016/0021-9290(77)90033-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Iqbal, S., Zang, X., Zhu, Y., and Zhao, J. (2014). Bifurcations and chaos in passive dynamic walking: a review. Robot Auton. Syst. 62, 889–909. doi: 10.1016/j.robot.2014.01.006

CrossRef Full Text | Google Scholar

Kane, T. R., and Levinson, D. A. (1985). Dynamics, Theory and Applications. McGraw Hill. Available online at: https://hdl.handle.net/1813/638

Kharb, A., Saini, V., Jain, Y., and Dhiman, S. (2011). A review of gait cycle and its parameters. Int. J. Comput. Eng. Manage. 13, 78–83. Available online at: https://www.ijcem.org/papers72011/72011_14.pdf

Google Scholar

Latash, M. L. (2010). Motor synergies and the equilibrium-point hypothesis. Motor Control 14, 294–322. doi: 10.1123/mcj.14.3.294

CrossRef Full Text | Google Scholar

Latash, M. L., Scholz, J. P., and Schoner, G. (2002). Motor control strategies revealed in the structure of motor variability. Exerc. Sport Sci. Rev. 30, 26–31. doi: 10.1097/00003677-200201000-00006

PubMed Abstract | CrossRef Full Text | Google Scholar

Lozi, R. (2013). “Can we trust in numerical computations of chaotic solutions of dynamical systems?” in Topology and Dynamics of Chaos: In Celebration of Robert Gilmore's 70th Birthday (World Scientific), 63–98. doi: 10.1142/9789814434867_0004

CrossRef Full Text | Google Scholar

Lunze, J., and Lamnabhi-Lagarrigue, F. (2009). Handbook of Hybrid Systems Control: Theory, Tools, Applications. Cambridge: Cambridge University Press. doi: 10.1017/CBO9780511807930

CrossRef Full Text | Google Scholar

Mahmoodi, P., Ransing, R., and Friswell, M. (2013). Modelling the effect of ‘heel to toe’ roll-over contact on the walking dynamics of passive biped robots. Appl. Math. Model. 37, 7352–7373. doi: 10.1016/j.apm.2013.02.048

CrossRef Full Text | Google Scholar

Manchester, I. R., Tobenkin, M. M., Levashov, M., and Tedrake, R. (2011). Regions of attraction for hybrid limit cycles of walking robots. IFAC Proc. 44, 5801–5806. doi: 10.3182/20110828-6-IT-1002.03069

CrossRef Full Text | Google Scholar

Mandali, A., Rengaswamy, M., Chakravarthy, V. S., and Moustafa, A. A. (2015). A spiking basal ganglia model of synchrony, exploration and decision making. Front. Neurosci. 9:191. doi: 10.3389/fnins.2015.00191

PubMed Abstract | CrossRef Full Text | Google Scholar

Martens, K. E., Hall, J. M., Gilat, M., Georgiades, M. J., Walton, C. C., and Lewis, S. (2016). Anxiety is associated with freezing of gait and attentional set-shifting in Parkinson's disease: a new perspective for early intervention. Gait Post. 49, 431–436. doi: 10.1016/j.gaitpost.2016.07.182

PubMed Abstract | CrossRef Full Text | Google Scholar

McGeer, T. (1990). Passive dynamic walking. Int. J. Robot. Res. 9, 62–82. doi: 10.1177/027836499000900206

CrossRef Full Text | Google Scholar

Montazeri Moghadam, S., Sadeghi Talarposhti, M., Niaty, A., Towhidkhah, F., and Jafari, S. (2018). The simple chaotic model of passive dynamic walking. Nonlinear Dyn. 93, 1183–1199. doi: 10.1007/s11071-018-4252-8

CrossRef Full Text | Google Scholar

Muralidharan, V., Balasubramani, P., Chakravarthy, S., Lewis, S. J. G., and Moustafa, A. (2014). A computational model of altered gait patterns in Parkinson's disease patients negotiating narrow doorways. Front. Comput. Neurosci. 7:190. doi: 10.3389/fncom.2013.00190

PubMed Abstract | CrossRef Full Text | Google Scholar

Nieuwboer, A., Dom, R., De Weerdt, W., Desloovere, K., Janssens, L., and Stijn, V. (2004). Electromyographic profiles of gait prior to onset of freezing episodes in patients with Parkinson's disease. Brain 127, 1650–1660. doi: 10.1093/brain/awh189

PubMed Abstract | CrossRef Full Text | Google Scholar

Nutt, J. G., Bloem, B. R., Giladi, N., Hallett, M., Horak, F. B., and Nieuwboer, A. (2011). Freezing of gait: moving forward on a mysterious clinical phenomenon. Lancet Neurol. 10, 734–744. doi: 10.1016/S1474-4422(11)70143-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Olszewski, W., and Sandroni, A. (2011). Falsifiability. Am. Econ. Rev. 101, 788–818. doi: 10.1257/aer.101.2.788

CrossRef Full Text | Google Scholar

Parakkal Unni, M., Sinha, A., Chakravarty, K., Chatterjee, D., and Das, A. (2017). Neuro-mechanical cost functionals governing motor control for early screening of motor disorders. Front. Bioeng. Biotechnol. 5:78. doi: 10.3389/fbioe.2017.00078

PubMed Abstract | CrossRef Full Text | Google Scholar

Pekarek, D., Ames, A. D., and Marsden, J. E. (2007). “Discrete mechanics and optimal control applied to the compass gait biped,” in 2007 46th IEEE Conference on Decision and Control (New Orleans, LA: IEEE), 5376–5382. doi: 10.1109/CDC.2007.4434296

CrossRef Full Text | Google Scholar

Polese, J. C., Teixeira-Salmela, L. F., Nascimento, L. R., Faria, C. D. M., Kirkwood, R. N., Laurentino, G. C., et al. (2012). The effects of walking sticks on gait kinematics and kinetics with chronic stroke survivors. Clin. Biomech. 27, 131–137. doi: 10.1016/j.clinbiomech.2011.08.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Popper, K. (2002). Popper: The Logic of Scientific Discovery. New York, NY: Routledge Classics.

PubMed Abstract | Google Scholar

Ranavolo, A., Conte, C., Iavicoli, S., Serrao, M., Silvetti, A., Sandrini, G., et al. (2011). Walking strategies of visually impaired people on trapezoidal-and sinusoidal-section tactile groundsurface indicators. Ergonomics 54, 246–256. doi: 10.1080/00140139.2010.548533

PubMed Abstract | CrossRef Full Text | Google Scholar

Rochester, L., Hetherington, V., Jones, D., Nieuwboer, A., Willems, A.-M., Kwakkel, G., et al. (2005). The effect of external rhythmic cues (auditory and visual) on walking during a functional task in homes of people with Parkinson's disease. Arch. Phys. Med. Rehabil. 86, 999–1006. doi: 10.1016/j.apmr.2004.10.040

PubMed Abstract | CrossRef Full Text | Google Scholar

Ros, J., Font-Llagunes, J. M., Plaza, A., and Kövecses, J. (2015). Dynamic considerations of heel-strike impact in human gait. Multibody Sys. Dyn. 35, 215–232. doi: 10.1007/s11044-015-9460-0

CrossRef Full Text | Google Scholar

Sadeghian, H., and Barkhordari, M. (2020). Orbital analysis of passive dynamic bipeds; the effect of model parameters and stabilizing arm. Int. J. Mech. Sci. 178:105616. doi: 10.1016/j.ijmecsci.2020.105616

CrossRef Full Text | Google Scholar

Sainburg, R. L. (2015). Should the Equilibrium Point Hypothesis (EPH) be considered a scientific theory? Motor Control 19, 142–148. doi: 10.1123/mcj.2014-0056

PubMed Abstract | CrossRef Full Text | Google Scholar

Sinnet, R. W., Powell, M. J., Shah, R. P., and Ames, A. D. (2011). A human-inspired hybrid control approach to bipedal robotic walking. IFAC Proc. 44, 6904–6911. doi: 10.3182/20110828-6-IT-1002.03802

CrossRef Full Text | Google Scholar

Snijders, A. H., Nijkrake, M. J., Bakker, M., Munneke, M., Wind, C., and Bloem, B. R. (2008). Clinimetrics of freezing of gait. Mov. Disord. 23, S468–S474. doi: 10.1002/mds.22144

PubMed Abstract | CrossRef Full Text | Google Scholar

Taga, G. (1995). A model of the neuro-musculo-skeletal system for human locomotion. Biol. Cybern. 73, 97–111. doi: 10.1007/s004220050166

CrossRef Full Text | Google Scholar

Tamura, D., Aoi, S., Funato, T., Fujiki, S., Senda, K., and Tsuchiya, K. (2020). Contribution of phase resetting to adaptive rhythm control in human walking based on the phase response curves of a neuromusculoskeletal model. Front. Neurosci. 14:17. doi: 10.3389/fnins.2020.00017

PubMed Abstract | CrossRef Full Text | Google Scholar

Tresch, M. C., and Jarc, A. (2009). The case for and against muscle synergies. Curr. Opin. Neurobiol. 19, 601–607. doi: 10.1016/j.conb.2009.09.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Usherwood, J. R. (2005). Why not walk faster? Biol. Lett. 1, 338–341. doi: 10.1098/rsbl.2005.0312

PubMed Abstract | CrossRef Full Text | Google Scholar

Wahde, M., and Pettersson, J. (2002). “A brief review of bipedal robotics research,” in Proceedings of the 8th UK Mechatronics Forum International Conference (Mechatronics 2002), 480–488. Available online at: http://www.me.chalmers.se/~mwahde/AdaptiveSystems/Publications/WahdePetterssonMech2002.pdf

Google Scholar

Wolfram Research Inc. (2019). Mathematica, Version 12.0.

Wu, A. R., Simpson, C. S., van Asseldonk, E. H., van der Kooij, H., and Ijspeert, A. J. (2019). Mechanics of very slow human walking. Sci. Rep. 9, 1–10. doi: 10.1038/s41598-019-54271-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Young, W. R., Rodger, M. W., and Craig, C. M. (2014). Auditory observation of stepping actions can cue both spatial and temporal components of gait in Parkinson's disease patients. Neuropsychologia 57, 140–153. doi: 10.1016/j.neuropsychologia.2014.03.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Zelik, K. E., and Adamczyk, P. G. (2016). A unified perspective on ankle push-off in human walking. J. Exp. Biol. 219, 3676–3683. doi: 10.1242/jeb.140376

PubMed Abstract | CrossRef Full Text | Google Scholar

Znegui, W., Gritli, H., and Belghith, S. (2020). Design of an explicit expression of the Poincaré map for the passive dynamic walking of the compass-gait biped model. Chaos Solitons Fractals 130, 109436. doi: 10.1016/j.chaos.2019.109436

CrossRef Full Text | Google Scholar