Dynamical consequences of sensory feedback in a half-center oscillator coupled to a simple motor system

Abstract

We investigate a simple model for motor pattern generation that combines central pattern generator (CPG) dynamics with a sensory feedback (FB) mechanism. Our CPG comprises a half-center oscillator with conductance-based Morris–Lecar model neurons. Output from the CPG drives a push–pull motor system with biomechanics based on experimental data. A sensory feedback conductance from the muscles allows modulation of the CPG activity. We consider parameters under which the isolated CPG system has either “escape” or “release” dynamics, and we study both inhibitory and excitatory feedback conductances. We find that increasing the FB conductance relative to the CPG conductance makes the system more robust against external perturbations, but more susceptible to internal noise. Conversely, increasing the CPG conductance relative to the FB conductance has the opposite effects. We find that the “closed-loop” system, with sensory feedback in place, exhibits a richer repertoire of behaviors than the “open-loop” system, with motion determined entirely by the CPG dynamics. Moreover, we find that purely feedback-driven motor patterns, analogous to a chain reflex, occur only in the inhibition-mediated system. Finally, for pattern generation systems with inhibition-mediated sensory feedback, we find that the distinction between escape- and release-mediated CPG mechanisms is diminished in the presence of internal noise. Our observations support an anti-reductionist view of neuromotor physiology: Understanding mechanisms of robust motor control requires studying not only the central pattern generator circuit in isolation, but the intact closed-loop system as a whole.

Appendices

Table of Parameter Values

Table 2 specifies the parameter values used for simulations

Table 2 Parameter Values

Model equations

For completeness, we list here the full equations of the model, as introduced in Sect.  2. The simulation codes are available from https://github.com/zhuojunyu-appliedmath/CPG-FB. Instructions for reproducing each figure and table in the paper are provided (see the README file at the github site).

For cell i ($\ne j$),

$$\begin{aligned} C\frac{dV_i}{dt}&=I_\text {ext}-g_\text {L}(V_i-E_\text {L})-g_\text {Ca}M_\infty (V_i)(V_i-E_\text {Ca})\\&\quad -g_\text {K}N_i(V_i-E_\text {K})\\&\quad -g_\text {syn}^\text {CPG}S_\infty ^\text {CPG}(V_j)\left( V_i-E_\text {syn}^\text {CPG}\right) \\&\quad -g_\text {syn}^\text {FB}S_\infty ^\text {FB}(L_j)(V_i-E_\text {syn}^\text {FB}),\\ \frac{dN_i}{dt}&=\lambda _N(V_i)(N_\infty (V_i)-N_i), \end{aligned}$$

and

$$\begin{aligned} M_\infty (V_i)&=\frac{1}{2}\left( 1+\tanh \left( \frac{V_i-E_1}{E_2}\right) \right) ,\\ N_\infty (V_i)&=\frac{1}{2}\left( 1+\tanh \left( \frac{V_i-E_3}{E_4}\right) \right) ,\\ \lambda _N(V_i)&=\phi _N\left( \cosh \left( \frac{V_i-E_3}{2E_4}\right) \right) ,\\ S_\infty ^\text {CPG}(V_j)&=\frac{1}{2}\left( 1+\tanh \left( \frac{V_j-E_\text {thresh}}{E_\text {slope}}\right) \right) ,\\ S_\infty ^\text {FB}(L_j)&=1-\frac{1}{2}\tanh \left( \frac{L_j-x_0\ell }{L_\text {slope}}\right) . \end{aligned}$$

For muscle i,

$$\begin{aligned} F_i(t)=F_0\cdot a_i(t)\cdot LT(L_i(t)). \end{aligned}$$

Here,

$$\begin{aligned}&LT(L_i(t))=-\left( 5.27\times 10^{-4}\right) (L_i(t)/\ell )^2\\&\quad +0.1054(L_i(t)/\ell )-4.27,\\&L_i(t)=(50\pm 0.8x(t))\ell ,\\&\quad L_j(t)=(50\mp 0.8x(t))\ell ,\quad -50<x<50, \end{aligned}$$

and

$$\begin{aligned} a_i(t)=g[A_i(t)-a_0],\quad \quad 0<a_i\le 1, \end{aligned}$$

where

$$\begin{aligned} dA_i/dt= & {} \tau ^{-1}\{U_i-[\beta +(1-\beta )U_i]A_i\},\\ U_i= & {} 1.03-4.31e^{(-0.198u_i)},\quad \quad 0\le U_i\le 1,\\ u_i= & {} (1/2)V_i. \end{aligned}$$

The position of the pendulum is determined by

$$\begin{aligned} \frac{dx}{dt}=\frac{1}{b\ell }(F_2-F_1). \end{aligned}$$

Flow-invariant set

Here we show that the state space of our system has a a compact flow-invariant domain, that is, a subset $\Omega \subset \mathbb {R}^7$ that is invariant under the flow generated by the model, forward in time. We describe the set as follows:

$$\begin{aligned}&(V_1,V_2,N_1,N_2,A_1,A_2,x)\in \Omega \nonumber \\&\quad = [V_\text {min},V_\text {max}]^2 \times [0,1]^2 \nonumber \\&\qquad \times [0,A_\text {max}]^2 \times [-x_\text {max},x_\text {max}]. \end{aligned}$$

(23)

For the specific parameters, we use $V_\text {min}\approx -76.63$ mV, $V_\text {max}= 110$ mV, $A_\text {max}\approx 1.021$, and $x_\text {max}\approx 8.05$.

For the range of $N_i$, since $N_\infty \in [0,1]$, $N_i\in [0,1]$ for all time as in Eq. (3). The voltage $V_i\in [V_\text {min},V_\text {max}]$ because if $V<V_\text {min}$, $\frac{dV}{dt}>0$ for all values of the other variables’ ranges, and if $V>V_\text {max}$, $\frac{dV}{dt}<0$ always holds. Considering the reversal potentials of all conductances, when $V\ge E_\text {Ca}$, the maximal possible inward current would occur if the potassium, CPG and FB conductances were shut off ($N=0, S_\infty ^\text {CPG}=S_\infty ^\text {FB}=0$), with extremal value of the calcium gate ($M_\infty$) to be determined. Thus,

$$\begin{aligned} \frac{dV_i}{dt}&\le I_\text {ext}-g_\text {L}(V_i-E_\text {L})-g_\text {Ca}M_{\infty }(V_i-E_\text {Ca})\\&=(I_\text {ext}+g_\text {L}E_\text {L}+g_\text {Ca}M_{\infty }E_\text {Ca})-(g_\text {L}+g_\text {Ca}M_\infty )V_i. \end{aligned}$$

To make $\frac{dV_i}{dt}<0$ hold $\forall \, V_i>V_\text {max}$, then,

$$\begin{aligned} V_\text {max}= & {} \max _{M_\infty \in [0,1]} \left\{ \frac{I_\text {ext}+g_\text {L}E_\text {L}+g_\text {Ca}M_{\infty }E_\text {Ca}}{g_\text {L}+g_\text {Ca}M_\infty }\right\} \\= & {} \max _{M_\infty \in [0,1]}\left\{ \frac{I_\text {ext}+g_\text {L}E_\text {L}+E_\text {Ca}(g_\text {L}+g_\text {Ca}M_\infty )-E_\text {Ca}g_\text {L}}{g_\text {L}+g_\text {Ca}M_\infty }\right\} \\= & {} E_\text {Ca}+\max _{M_\infty \in [0,1]}\left\{ \frac{I_\text {ext}+g_\text {L}E_\text {L}-E_\text {Ca}g_\text {L}}{g_\text {L}+g_\text {Ca}M_\infty }\right\} . \end{aligned}$$

As $M_\infty$ takes its minimum at 0, $V_\text {max}=110$ mV for parameters we choose in Table 2.

In order to find $V_\text {min}$, we turn off the $\text {Ca}^{2+}$ channel and allow for maximal inhibitory feedback:

$$\begin{aligned} \frac{dV_i}{dt}&\ge I_\text {ext}-g_\text {L}(V_i-E_\text {L})-g_\text {K}N_i(V_i-E_\text {K})\\&\quad -g_\text {syn}^\text {CPG}S_\infty ^\text {CPG}\left( V_i-E_\text {syn}^\text {CPG}\right) -g_\text {syn}^\text {FB}S_\infty ^\text {FB}(V_i-E_\text {syn}^\text {FB})\\&=(I_\text {ext}+g_\text {L}E_\text {L}+g_\text {K}N_iE_\text {K}\\&\quad +g_\text {syn}^\text {CPG}S_\infty ^\text {CPG}E_\text {syn}^\text {CPG}+g_\text {syn}^\text {FB} S_\infty ^\text {FB}E_\text {syn}^\text {FB})\\&\quad -(g_\text {L}+g_\text {K}N_i+g_\text {syn}^\text {CPG}S_\infty ^\text {CPG}+g_\text {syn}^\text {FB}S_\infty ^\text {FB})V_i. \end{aligned}$$

If $\frac{dV_i}{dt}>0$ always holds for arbitrary $V_i<V_\text {min}$, then

$$\begin{aligned} V_\text {min}&=\min _{N_i,S^\text {CPG}_\infty ,S_\infty ^\text {FB}\in [0,1]}\\&\quad \left\{ \frac{I_\text {ext}+g_\text {L}E_\text {L}+g_\text {K}N_iE_\text {K}+g_\text {syn}^\text {CPG}S_\infty ^\text {CPG}E_\text {syn}^\text {CPG} +g_\text {syn}^\text {FB}S_\infty ^\text {FB}E_\text {syn}^\text {FB}}{g_\text {L}+g_\text {K}N_i+g_\text {syn}^\text {CPG}S_\infty ^\text {CPG} +g_\text {syn}^\text {FB}S_\infty ^\text {FB}}\right\} \\&=\min _{N_i,S^\text {CPG}_\infty ,S_\infty ^\text {FB}\in [0,1]}\\&\quad \left\{ \frac{I_\text {ext}+g_\text {L}E_\text {L}-80(g_\text {L}+g_\text {K}N_i+g_\text {syn}^\text {CPG}S_\infty ^\text {CPG}+g_\text {syn}^\text {FB}S_\infty ^\text {FB}) +80g_\text {L}}{g_\text {L}+g_\text {K}N_i+g_\text {syn}^\text {CPG}S_\infty ^\text {CPG}+g_\text {syn}^\text {FB}S_\infty ^\text {FB}}\right\} \\&=-80+\min _{N_i,S^\text {CPG}_\infty ,S_\infty ^\text {FB}\in [0,1]}\left\{ \frac{I_\text {ext}+g_\text {L}E_\text {L}+80g_\text {L}}{g_\text {L}+g_\text {K}N_i+g_\text {syn}^\text {CPG}S_\infty ^\text {CPG}+ g_\text {syn}^\text {FB}S_\infty ^\text {FB}}\right\} . \end{aligned}$$

When $N_i$, $S_\infty ^\text {CPG}$ and $S_\infty ^\text {FB}$ simultaneously reach their maximum, $V_\text {min}\approx -76.63$ mV. For comparison, the voltage range in Fig. 2 roughly corresponds to $V_\text {min}<V<V_\text {max}$.

For the variables related to the biomechanics, as shown in Eq. (12), the range of $U_i$ is

$$\begin{aligned} {[}0,U_\text {max}]=[0,f(V_\text {max}/2)]\subset [0,1.03). \end{aligned}$$

The U-A relation, Eq. (13), gives

$$\begin{aligned} A_i\in [0,A_\text {max}]=\left[ 0,\frac{U_\text {max}}{\beta +(1-\beta )U_\text {max}}\right] \subset [0,1.021). \end{aligned}$$

Since $L_1=(50+0.8x)\ell$ and $L_2=(50-0.8x)\ell$, then as in Eq. (11),

$$\begin{aligned} LT(L_1)= & {} -(5.27\times 10^{-4})(50+0.8x)^2\\&+0.1054(50+0.8x)-4.27,\\ LT(L_2)= & {} -(5.27\times 10^{-4})(50-0.8x)^2\\&+0.1054(50-0.8x)-4.27. \end{aligned}$$

For convenience, we drop the length scale $\ell \,(=1\text {mm})$ since it does not affect the calculation result. Equation (14) defines the range for $a_i$: $a_i\in (0,a_\text {max}]=(0,g(A_\text {max}-a_0)]\subset (0,1.712)$. Because $F_i$ is determined by $a_i$ and x (Eq. (10)), and x is in turn determined by $F_i$ (Eq. (16)), we can compute $x_\text {min}$ and $x_\text {max}$ through the ODE $dx/dt=\frac{1}{b}(F_2-F_1)$, i.e.,

$$\begin{aligned} \frac{dx}{dt}&=\frac{1}{b}(F_2-F_1)\\&=\frac{1}{b}(F_0\cdot a_2\cdot LT(L_2)-F_0\cdot a_1\cdot LT(L_1))\\&=\frac{F_0}{b}\{a_2[-(5.27\times 10^{-4})(50-0.8x)^2\\&\quad +0.1054(50-0.8x)-4.27]\\&\quad -a_1[-(5.27\times 10^{-4})(50+0.8x)^2\\&\quad +0.1054(50+0.8x)-4.27]\}\\&=\frac{F_0}{b}[(3.3728\times 10^{-4})(a_1-a_2)x^2\\&\quad -0.04216(a_1+a_2)x+0.3175(a_1-a_2)]\\&:=\frac{F_0}{b}[W(a_1-a_2)x^2-Y(a_1+a_2)x+Z(a_1-a_2)]. \end{aligned}$$

On the one hand, if $x>x_\text {max}$, then $\frac{dx}{dt}<0$ for all other variables within their respective ranges. Thus, we aim to find $x_\text {max}\in (-50,50)$ such that for all $x_\text {max}<x<50$,  $W(a_1-a_2)x^2-Y(a_1+a_2)x+Z(a_1-a_2)<0$ for arbitrary $a_1,\, a_2\in (0,1.712)$. If $a_1\le a_2$, the inequality holds for all $x>0$. When $a_1>a_2$, $x=50$ always satisfies the inequality, indicating that $x_\text {max}(<50)$ exists and should be the maximal abscissa of the left intersection of the quadratic function and horizontal axis, denoted by $x_\text {L}$:

$$\begin{aligned} x_\text {L}&=\frac{Y(a_1+a_2)-\sqrt{Y^2(a_1+a_2)^2-4WZ(a_1-a_2)^2}}{2W(a_1-a_2)}\\&=\frac{Y}{2W}\left( \frac{1+\frac{a_2}{a_1}}{1-\frac{a_2}{a_1}}\right) -\frac{1}{2W}\sqrt{Y^2\left( \frac{1+\frac{a_2}{a_1}}{1-\frac{a_2}{a_1}}\right) ^2-4WZ}. \end{aligned}$$

Let $t=\frac{a_2}{a_1}\in (0,1)$ and let $r=\frac{1+t}{1-t}=1+\frac{2}{\frac{1}{t}-1}$. Then r is a monotonically increasing function from 1 to $\infty$. Hence,

$$\begin{aligned}x_\text {L}=\frac{Y}{2W}r-\frac{1}{2W}\sqrt{Y^2r^2-4WZ}.\end{aligned}$$

It is not hard to find that $x_\text {L}$ monotonically decreases by calculating its derivative, and thus,

$$\begin{aligned}x_\text {max}=(x_\text {L})_\text {max}=x_\text {L}\big |_{r=1}=\frac{Y-\sqrt{Y^2-4WZ}}{2W}\approx 8.05.\end{aligned}$$

On the other hand, if $x<x_\text {min}$, then $\frac{dx}{dt}>0$ for all other variables within their respective ranges. Given the symmetry of our model, $x_\text {min}=-x_\text {max}\approx -8.05$. Therefore, the flow-invariant set for the solutions of our HCO–muscle model has the form given by (23), with $V_\text {min}\approx -76.63$ mV, $V_\text {max}= 110$ mV, $A_\text {max}\approx 1.021$ and $x_\text {max}\approx 8.05$ for the specific parameters used in this paper.