Analytical CPG model driven by limb velocity input generates accurate temporal locomotor dynamics

CPG structure and function

The observations of neural activity in the absence of descending signals or sensory feedback led T.G. Brown to formulate the principle of intrinsic rhythmogenesis of spinal networks, the half-center oscillator hypothesis (Brown, 1911). Brown posited that “…the centres are paired, and that each pair consists of antagonistic opposites.” The intrinsic rhythmogenesis opposed the established view that the locomotor pattern is generated and shaped only by supraspinal and sensory feedback pathways. The bilateral CPG model in Fig. 1 was developed from a numerical model of a single-limb oscillator to describe phase dominance in fictive cat locomotion, which is a type of experimental behavior with diminished sensory contribution (Yakovenko et al., 2005). This model controlling two limbs consisted of two dedicated oscillators made of two reciprocally coupled half-center elements (gray area in Fig. 1). It can generate bilateral rhythm using the interactions within and between the half-center elements. Only the rhythm generating mechanism is captured by this feedforward rate model with time-varying inputs. The pattern formation mechanism responsible for the generation of motoneuronal input signals can be computationally decoupled from the temporal dynamics of rhythm generation (McCrea & Rybak, 2008).

The process of controlling locomotor phase durations is based on the ability of the network to integrate inputs until reaching a critical threshold causing a phase resetting within the CPG network, Fig. 2. We have previously developed the bilateral model (Yakovenko, 2011; Sobinov & Yakovenko, 2018) and describe it in brief here. The model was expressed as the system of differential equations consisting of two parts in Eq. (1): (i) the largely extrinsic signals (right side) and (ii) the intrinsic interactions (left side). The offset term (x ₀) could combine both intrinsic and extrinsic influences on the background excitability of spinal cord. The bilateral CPG model consists of a system of differential equations for four intrinsic states ($x={\left(x_1,x_2,x_3,x_4\right)}^T$) that represent flexor and extensor locomotor phases for each limb: (1)${\displaystyle \dot{x}-G_x^{UL}x=x_0+G_{u}u}$ where Gu matrix represents gains of input signals u, x ₀ are constant offset values, $G_x^{UL}$ matrix represents the strength of unilateral connections between the CPG half-centers (shown as arrows with weights r_ij in Fig. 1, the connections across the midline were removed). $G_x^{UL}$ matrix has the following form: (2)${\displaystyle G_x^{UL}=I\ast r_{\mathrm{leak}}}$ where I is the identity matrix, r_leak is the constant that determines intrinsic state-dependent feedback. This parametric transection effectively decouples the control of left and right limb. However, the bilateral coupling in this model could still be achieved through the common descending drive (as demonstrated for the heading direction control, below).

The internal states are limited to positive values with the switching threshold set to 1. Only one state from a pair, 1–2 (Fig. 2), is set to be active x ∈ (0, 1] to impose the reciprocal relationship between half-centers. This implementation assumes robust reciprocity between antagonistic states and enforces zero overlap between them. The single limb CPG would consist only of two reciprocal states ($x={\left(x_1,x_2\right)}^T$).

Even this simple model had many parameters that were largely undefined. Using an error-driven search algorithm in our previous study (Yakovenko, 2011) we found a set of optimal parameters (Table A1). These parameters were resolved by the minimization of the objective function with terms related to the errors in simulating swing and stance phases and the rate of their modulation for different overground speeds (Goslow, Reinking & Stuart, 1973; Halbertsma, 1983).

Solution using constant rate assumption

First, let us express explicitly all the term in Eq. (1) to describe only flexor and extensor states controlling a single limb. Here, x ₁ and x ₂ are the reciprocal state variables as shown in Fig.  1. The system of equations can then be stated as: (3)${\displaystyle \left\{\begin{array}{c}{\dot{x}}_1=x_{01}+g_{u1}u+r_{\mathrm{leak}}{x}_1\hfill \\ {\dot{x}}_2=x_{02}+g_{u2}u+r_{\mathrm{leak}}{x}_2.\hfill \end{array}\right.}$ Since r_leak is a small negative number (Table A1) the rate of state ($\dot{x}$) can be further approximated without this term using phase duration quantities as the difference of states for a given phase duration, i.e., the inverse of phase duration. Even for the time-variable input (u), the rate of state for a full phase duration can be simplified as: (4)${\displaystyle \dot{x}=\frac{max-min}{\tau }=\frac{1}{\tau }.}$ Figure 2 shows an example for this formulation based on Eq. (1) for a single limb with the assumption of the constant rate of integration (g_u1 and g_u2 are scalars, as in our previous studies). Each state (x ₁ and x ₂) integrates an input (u) only when active. The integration rate per step cycle can then be stated as in Eq. (4). This formulation is possible due to the removal of midline crossing connections (green in Fig. 1) between CPG states that complicate the relationship. Then the expression for cycle duration can be described as a sum of the antagonistic phases in the simplified system, Eq. (5): (5)${\displaystyle T_c={\tau }_1+{\tau }_2=\frac{1}{{\dot{x}}_1}+\frac{1}{{\dot{x}}_2}=\frac{{\dot{x}}_1+{\dot{x}}_2}{{\dot{x}}_1{\dot{x}}_2}.}$ Since the cycle duration, Tc, is a constant for a given constant input (u), the only time-varying variables are the states of the system, x ₁ and x ₂. In phase transition points, at t = τ₁ or t = τ₁ + τ₂, x ₁ and x ₂ are zero or a small value close to zero. We can further expand this equation with Eq. (3) and simplify it to all the known terms: (6)${\displaystyle T_c=\frac{x_{01}+x_{02}+\left(g_{u1}+g_{u2}\right)u}{\left(x_{01}+g_{u1}u\right)\left(x_{02}+g_{u2}u\right)}=\frac{a+bu}{\tilde{a}+\tilde{b}u+\tilde{c}{u}^2}}$ where the step cycle duration is expressed as a function of input (u) and all parameters a, b, $\tilde{a}$, $\tilde{b}$, $\tilde{c}$ are constants determined by the coefficients in the system of equations Eq. (3).

Solution using Taylor series

The same solution Eq. (6) was found by integrating the differential equations (3) between 0 and t. For this, Eq. (3) can be rewritten with the assumption of independent limb control: (7)${\displaystyle \dot{x}-rx=x_0+G_{u}u}$ where variables are as defined for Eq. (1), and r = r_leak. Note that the right-hand side can be assumed to be time-independent for constant input (u) and this type of equations has a general solution of the form e^kx. The left side of the above equation can be expressed as (8)${\displaystyle {\left(x{e}^{-rt}\right)}^{\prime }=\dot{x}{e}^{-rt}-rx{e}^{-rt}=\left(\dot{x}-rx\right)e^{-rt}.}$ Hence, Eq. (7) can be integrated and evaluated between 0 and t using (9)${\displaystyle {\left.\left(x{e}^{-rt}\right)\right|}_0^t={\int }_0^t\left(x_0+G_{u}u\right)e^{-rt}dt}$ (10)${\displaystyle x\left(t\right)e^{-rt}-0=\frac{x_0+G_{u}u}{-r}\left(e^{-rt}-1\right)}$ (11)${\displaystyle x\left(t\right)=\frac{x_0+G_{u}u}{r}\left(e^{rt}-1\right).}$ The exponential function can be further expanded with Taylor series and some components can be dropped since r is a number close to zero, so that rt ≈ 0 in the expansion: (12)${\displaystyle x\left(t\right)\approx \frac{x_0+G_{u}u}{r}\left(1+rt+\dots -1\right)\approx \left(x_0+G_{u}u\right)t.}$

Then, the full phase of each integrated state is (13)${\displaystyle t=\frac{1}{x_0+G_{u}u}.}$

Finally, the full cycle duration consisting of two reciprocal phases (t₁ + t₂) has the same form as Eq. (6) (14)${\displaystyle T_c=t_1+t_2=\frac{a+bu}{\tilde{a}+\tilde{b}u+\tilde{c}{u}^2}}$ where a, b, $\tilde{a}$, $\tilde{b}$, $\tilde{c}$ are constants defined by the examination of algebraic terms from Eq. (13).

Validation

Both methods converged on the same form, Eqs. (6) and (14), supporting the consistency of solutions with different assumptions. The relationship between step cycle duration and CPG input (Tc and u) is of the form T_c = au^−b. This simple analytical solution has a similar form to the phenomenological relationship between cycle duration and the velocity of overground forward progression T_c=0.5445V^−0.5925 (Goslow, Reinking & Stuart, 1973). Figure 3 shows the comparison of solutions with our analytical and the previous phenomenological model for the step cycle duration and velocity values. The simulated T_c data values were calculated with Eq. (7) using optimal parameters and u values selected with the regression equation u = (V + 0.1272)∕0.2357 (from Fig. 4 in our previous work (Yakovenko, 2011)) and plotted in Fig. 3C. The analytical solution (red) for leg speed was closely related to the empirical curve (black) calculated with the phenomenological functions that were calculated as the best-fit expressions for the experimental measurements (Goslow, Reinking & Stuart, 1973; Halbertsma, 1983) (Fig. 3A). Since the swing duration remains nearly constant as a function of either step cycle duration or the velocity of forward progression in cats (Halbertsma, 1983; Frigon et al., 2014), it can be approximated as a constant (≈0.25 s). Then, the stance duration is the same as in Eqs. (6) and (14) with the negative constant. This relationship may not be preserved for gaits with the large differences in limb speeds as those used in split-belt experiments (D’Angelo et al., 2014) and may require the consideration of bilateral inputs as in our previous study (Sobinov & Yakovenko, 2018). The analytical and empirical step cycle durations were highly correlated (Fig. 3B) for the linear relationship between CPG inputs (u) representing scaled forward velocity values (Fig. 3C).

The implementation of CPG with the limb speed inputs is expected to generate spatiotemporal step modulation appropriate for the locomotion on a curved path. We used the following equation to compute the relationship between the heading direction (γ) and limb speeds (V_R, V_L) : γ = Tc(V_R−V_L)∕W, where W is the interlimb step width (≈0.15 m in cat) and Tc is the step duration from Eq. (6). This equation was previously derived for the experimental and theoretical kinematics of locomotion along a curved path (Courtine, Papaxanthis & Schieppati, 2006; Sobinov & Yakovenko, 2018). Figure 4A shows the monotonic relationship with extremes occurring for large interlimb speed difference (∂V) at slow speeds. This range may be consistent with the “spin turning” when body spins around a supporting limb (Hase & Stein, 1999). Examples of simulated kinematics for three speeds (a, b, c = 0.5, 1, 2 m/s) with increasing interlimb speed differences in each consecutive step are shown in Fig. 4B. The green vector indicating the heading direction demonstrates the dependency not only on the interlimb speed difference (∂V), but also the overall magnitude of body’s velocity (V). The parsimonious analytical CPG model that includes computations for two limbs can generate steering.