Multivariable harmonic balance analysis of the neuronal oscillator for leech swimming

Abstract

Biological systems, and particularly neuronal circuits, embody a very high level of complexity. Mathematical modeling is therefore essential for understanding how large sets of neurons with complex multiple interconnections work as a functional system. With the increase in computing power, it is now possible to numerically integrate a model with many variables to simulate behavior. However, such analysis can be time-consuming and may not reveal the mechanisms underlying the observed phenomena. An alternative, complementary approach is mathematical analysis, which can demonstrate direct and explicit relationships between a property of interest and system parameters. This paper introduces a mathematical tool for analyzing neuronal oscillator circuits based on multivariable harmonic balance (MHB). The tool is applied to a model of the central pattern generator (CPG) for leech swimming, which comprises a chain of weakly coupled segmental oscillators. The results demonstrate the effectiveness of the MHB method and provide analytical explanations for some CPG properties. In particular, the intersegmental phase lag is estimated to be the sum of a nominal value and a perturbation, where the former depends on the structure and span of the neuronal connections and the latter is roughly proportional to the period gradient, communication delay, and the reciprocal of the intersegmental coupling strength.

Appendices

Appendix

1 Notation

The imaginary unit is denoted by ${\text{j}}$. Let e be a column vector with all entries being 1, and I is an identity matrix. A scalar function f acts elementwise when its argument is a matrix or vector. For instance, when X is a matrix, Y = f(X) means that Y _kl  = f(X _kl ) for all k and l where Y _kl and X _kl are the (k,l) entries of the matrices. The set of n dimensional real (complex) vectors is denoted by ℝⁿ (ℂⁿ). For a complex number x ∈ ℂ, the imaginary part is denoted by $\Im(x)$, and $\angle x$ is the polar angle of x in the complex plane, i.e., if $x=a+{\text{j}} b$, then $\angle x=\theta$ with tan(θ) = b/a. For a vector x ∈ ℂⁿ, we define y = |x| and $y=\angle x$ by y _k  = |x _k | and $y_k=\angle x_k$, respectively, for k = 1,...,n, and ${\mbox{diag}}(x)$ is the diagonal matrix with entries x _k . For a scalar x ∈ ℂ and a vector y ∈ ℂⁿ, we write z = x + y to mean z _k  = x + y _k for k = 1,...,n. The Kronecker product of vectors (or matrices) x and y is denoted by x ⊗ y.

2 Formulae for describing functions

The following are formulae for describing functions of the threshold nonlinearity ϕ(x): =  max (x,0):

$$\begin{array}{lll} \kappa_1(\mathfrak{a},\mathfrak{b})&:=&\frac{1}{\mathfrak{a}\pi}\int_{-\pi}^\pi \phi(\mathfrak{a}[\sin t+\mathfrak{b}]) \sin t dt \\ &=&\left\{\begin{array}{ll} (1/\pi)[\mathfrak{b}\sqrt{1-\mathfrak{b}^2}+\pi/2+\sin^{-1}\mathfrak{b}], & (|\mathfrak{b}|\leq 1)\\ 1, & (\mathfrak{b}>1) \\ 0, & (\mathfrak{b}<-1) \end{array}\right. \\ \kappa_2(\mathfrak{a},\mathfrak{b})&:=& \frac{1}{2\mathfrak{a}\pi}\int_{-\pi}^\pi \phi(\mathfrak{a}[\sin t+\mathfrak{b}]) dt \\ &=&\left\{\begin{array}{ll} (1/\pi)[\mathfrak{b}(\pi/2+\sin^{-1}\mathfrak{b})+\sqrt{1-\mathfrak{b}^2}], & (|\mathfrak{b}|\leq 1)\\ \mathfrak{b}, & (\mathfrak{b}>1) \\ 0, & (\mathfrak{b}<-1) \end{array}\right. \\ \end{array} $$

3 The leech CPG model

A systems-level model for leech swimming CPG has been developed in Zheng et al. (2007). The following is a slightly modified version of the model, where pure delays are used for the intersegmental coupling dynamics instead of first-order time lags:

$$\begin{array}{rcl} v_k&=&{\beta}+\sum\limits_{l=1}^m{\text{M}}_{kl}(s)\phi(v_l), {} (k=1,\ldots,m), \\ {\text{M}}_{kl}(s)&=&\left\{\begin{array}{ll} \mu f(\tau_k s) M & (l=k), \\[4pt] {\epsilon} e^{-|k-l|\tau_ds}M_{kl} & (l\neq k), \end{array}\right. \\ f(s)&=&\frac{1-r}{1+(1-r)s}, {} \phi(x):=\max(x,0), \\ M&=&\left[ \begin{array}{rrr} 0 & -1 & 0 \\ 0 & 0 & -1 \\ -1 & 0 & 0 \end{array}\right], {} M_A=\left[ \begin{array}{rrr} 0 & -1 & 0 \\ 0 & 0 & -1 \\ 0 & 0 & 0 \end{array} \right], {}\\ M_D&=&\left[ \begin{array}{rrr} 2 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{array}\right], \\ M_{kl}&=&\left \{\begin{array}{llll} M_A (k<l\leq k+q), \\[4pt] M_D (k-q\leq l<k), \\[4pt] 0 (\mbox{otherwise}), \end{array} \right. \end{array}$$

$$\begin{array}{l}{m = 17,\,\,\,\,\tau _o = 200\,\text{ms,}\,\,\,\,r = 0.3,} \\ \,\,{q = 5,\,\,\,\,\,\,\,\tau _d = 15\,\text{ms,}\,\,\,\,\,\,\,\mu = 6,} \\ {\qquad\qquad\,\,\,\beta = 30r\,\text{mV}\,\,\,\sigma: = {\varepsilon \mathord{\left/ {\vphantom {\varepsilon {\mu = 0.015.}}} \right. \kern-\nulldelimiterspace} {\mu = 0.015.}}} \end{array}$$

The nominal value τ _o can be chosen equal to either one of the τ _k values or their average value. Let Δτ _o be a nominal size of the variations in time constants; for instance,

$$ {\Delta}\tau_o:=\left[\max\limits_k(\tau_k)-\min\limits_k(\tau_k)\right]/(m-1). $$

We assume that the variation is small, i.e., ε: = Δτ _o /τ _o  ≪ 1. To the first-order approximation in ε, we have

$$ f(\tau_k s)\approx f(\tau_o s)+{\varepsilon} g_k(s), g_k(s):=\frac{\tau_k-\tau_o}{{\Delta}\tau_o}\tau_o sf'(\tau_o s). $$

(20)

As a result, the leech swimming CPG can be put in the form of Eq. (4) with

$$ M_o(s)=f(\tau_o s)M,{} \Delta_{kl}(s)=\left\{\begin{array}{ll} ({\varepsilon}/\sigma) g_k( s) M & (l=k), \\[4PT] e^{-|k-l|\tau_ds} M_{kl} & (l\neq k), \end{array}\right. $$

In Eq. (20), note that τ _k  − τ _o is determined by the time constant gradient Δτ ∈ ℝ^m − 1 defined by its kth entry

$$ {\Delta}\tau_k:=\tau_{k+1}-\tau_k $$

for k = 1,...,m − 1. Thus, the shape of the time constant variation along the body is encoded into g _k (s), while the magnitude of the variation is embedded in ε. For instance, if the nominal value τ _o is chosen equal to τ _p , then

$$ g_k(s)=\left(\sum\limits_{i=p}^{k-1}\zeta_i-\sum\limits_{i=k}^{p-1}\zeta_i\right) \tau_o sf'(\tau_o s) $$

where the summation from i = i ₁ to i ₂ is considered to be zero if i ₁ > i ₂, and ζ: = Δτ/Δτ _o represents the shape of the time constant variation. If, in particular, we consider the uniformly increasing time constant gradient, we have ζ = e.

4 Derivation of the MHB analysis condition for the nominal segmental oscillator

Let (ω, a,b,ψ) be the oscillation profile of the nominal segmental oscillator Eq. (3) with ε = 0. We define $\hbar:=Ae^{{\text{j}}\psi}$ the phasor of v(t) with $A:={\mbox{diag}}(a)$. We have the so-called harmonic balance equations, with $B:={\mbox{diag}}(b)$, shown as follows:

$$ 0 = \left[I- \mu M_o({\text{j}}\omega) \kappa_1(B)\right]\hbar, {} \hbar=Ae^{{\text{j}}\psi}, $$

(21)

$${\beta} = \left[ B- \mu M_o(0)\kappa_2(B)\right]a. $$

(22)

We see that the matrix $\mu M_o({\text{j}}\omega) \kappa_1(B)$ has an eigenvalue 1 and an associated eigenvector $\hbar$. For further analysis later, we denote the normalized left eigenvector by ℓ ∈ ℂⁿ:

$$ 0 = \ell^*\left[I- \mu M_o({\text{j}}\omega) \kappa_1(B)\right], {} \ell^*\hbar=1. $$

Under the uniformity assumptions $a=\mathfrak{a} e$ and $b= \mathfrak{b} e$, Eqs. (21) and (22) reduce Eqs. (8) and (9), respectively, and ℓ becomes a left eigenvector of M.

5 Derivation of the MHB analysis condition for weakly coupled oscillators

In this section, we will derive condition Eq. (11). For the uncoupled oscillators, we have their MHB equations as follows:

$$\begin{array}{lll} 0&=&[I-\mu{{\cal M}}_o({\text{j}}\omega)\kappa_1({{\cal B}}_o)]{\text{h}}_o, {} {\text{h}}_o:=\alpha\otimes\hbar, {} |{\text{h}}_o|={\text{a}}_o \\ {\beta} &=&[{{\cal B}}_o-\mu{{\cal M}}_o(0)\kappa_2({{\cal B}}_o)]{\text{a}}_o, {} {{\cal B}}_o:={\rm diag}(B,\ldots,B) \end{array}$$

for an arbitrary vector α ∈ ℂ^m satisfying |α| = e. The MHB equations for the coupled oscillators are given by

$$\begin{array}{lll} 0&=&[I-{{\cal M}}({\text{j}}\varpi)\kappa_1({{\cal B}})]{\text{h}}, {} |{\text{h}}|={\text{a}} \\ {\beta} &=&[{{\cal B}}-{{\cal M}}(0)\kappa_2({{\cal B}})]{\text{a}}. \end{array}$$

Let the variables after the weak coupling be approximated by

$$ \varpi\approx\omega, {} {\text{h}}\approx {\text{h}}_o+{\epsilon}\tilde{{\text{h}}},{} {{\cal B}}\approx {{\cal B}}_o+{\epsilon}\tilde{{{\cal B}}} $$

where we assume that the frequency perturbation is small and can be neglected.

Then, noting that

$$ \kappa_i({{\cal B}})\approx\kappa_i({{\cal B}}_o)+{\epsilon}\kappa_i'({{\cal B}}_o)\tilde{{\cal B}}, $$

the first MHB equation for the coupled oscillator becomes

$$ \begin{array}{lll} &&{\kern-6pt}[I-\mu{{\cal M}}_o({\text{j}}\omega)\kappa_1({{\cal B}}_o)]\tilde{\text{h}}\\ &&{\kern4pt} = [\mu{{\cal M}}_o({\text{j}}\omega)\kappa_1'({{\cal B}}_o)\tilde{{\cal B}}+ \Delta({\text{j}}\omega)\kappa_1({{\cal B}}_o)] (\alpha\otimes\hbar), \end{array} $$

(23)

where we used the MHB equations for the uncoupled oscillators and neglected the O(ε ²) terms. By the definitions of L and H, we have

$$ L^*(I-\mu{{\cal M}}_o({\text{j}}\omega)\kappa_1({{\cal B}}_o))=0, {} L^*H=I. $$

Multiplying Eq. (23) by L ^* from the left,

$$ (L^*\Upsilon H-\nabla)\alpha=0, $$

where

$$ \Upsilon:=-\kappa_1({{\cal B}}_o)^{-1}\kappa_1'({{\cal B}}_o)\tilde{{\cal B}}, {} \nabla:=L^*\Delta({\text{j}}\omega)\kappa_1({{\cal B}}_o)H. $$

We consider a chain of oscillators where the intersegmental coupling is uniform over the chain. That is, the coupling matrix from the (k − l)th oscillator to the kth is the same as that from the kth oscillator to the (k + l)th, so that the overall intersegmental coupling matrix has the block-Toeplitz structure. In this case, the perturbation vectors for all oscillators may be assumed approximately parallel to each other, i.e.,

$$ \tilde{{{\cal B}}} \approx \Gamma\otimes\tilde{B}, $$

where Γ is a real diagonal matrix, and we have

$$ L^*\Upsilon H=\rho\Gamma, {} \rho:=-\ell^*\kappa_1(B)^{-1}\kappa_1'(B)\tilde B\hbar. $$

Thus, the MHB equation implies

$$ (\rho\Gamma-\nabla)\alpha=0. $$

6 Numerical algorithm for solving the MHB condition

We propose the following algorithm for solving Eq. (11):

• 0. Let ε ₁, ε ₂ be small positive numbers and initialize k = 1 and $\gamma_k=[1\;\cdots\; 1]{^{\mbox{\tt T}}}/\sqrt{m}\in{\mathbb{R}}^m$.

• 1. Let $\Gamma_k=\mbox{diag}(\gamma_k)$. Let (ρ _k ,q _k ) be the maximal eigenvalue/eigenvector pair of $\nabla\Gamma_k^{-1}$ with ||q _k || = 1.

• 2. Update $\gamma_{k+1}:={\bar\gamma}_k/\|{\bar\gamma}_k\|$ where ${\bar\gamma}_k:=\gamma_k+{\epsilon}_1|q_k|$.

• 3. If ||γ _k + 1 − γ _k || ≤ ε ₂, then stop; otherwise, increment k and go to step 1 to iterate.

From $(\rho_kI-\nabla\Gamma_k^{-1})q_k=0$, we have $(\rho_k\Gamma_k-\nabla)\Gamma_k^{-1}q_k=0$. Thus, the first equation in Eq. (11) holds for Γ = Γ _k , ρ = ρ _k , and $\alpha=\Gamma_k^{-1}q_k$ at every step k. If the algorithm stops for a sufficiently small ε ₂, we have γ _k + 1 ≈ γ _k , implying that γ _k  ≈ |q _k |. Hence, we have |α| ≈ e on exit, satisfying the second equation in Eq. (11) approximately.

7 MHB analysis of the intersegmental phase lag

We shall derive the formula Eq. (18) for the average value of the intersegmental phase lags for a general chain of weakly coupled oscillators, where the coupling structures M, M _A , and M _D are arbitrary. We first consider the case where the segmental oscillators are identical (no period gradient), the biases within each oscillator are uniform ($b=\mathfrak{b} e$), and the intersegmental coupling is uniform over the chain. Then, the basic coupling matrix $\nabla$ is given by

$$ \nabla={\text{k}}\left[\begin{array}{cccccc} 0 & \nabla_1 & \cdots & \nabla_{{q_{\mbox{$A$}}}} && 0 \\[-1pt] \nabla_{-1} & 0 & \nabla_1 & \ddots & \ddots & \\[-1pt] \vdots & \nabla_{-1} & 0 & \nabla_1 & \ddots & \nabla_{{q_{\mbox{$A$}}}} \\[-1pt] \nabla_{-{q_{\mbox{$D$}}}} & \ddots & \nabla_{-1} & 0 & \ddots & \vdots \\[-1pt] & \ddots & \ddots & \ddots & \ddots & \nabla_1 \\[-1pt] 0 && \nabla_{-{q_{\mbox{$D$}}}} & \cdots & \nabla_{-1} & 0 \end{array}\right] $$

where ${q_{\mbox{$A$}}}$ and ${q_{\mbox{$D$}}}$ are the intersegmental projection spans in the ascending and descending directions, respectively; ${\text{k}}$ is defined by ${\text{k}}:=\kappa_1(\mathfrak{b})$; and

$$ \nabla_{ k}:={r_{\mbox{$A$}}} e^{{\text{j}}(\eta_A-k\omega\tau_d)}, {} {r_{\mbox{$A$}}} e^{{\text{j}}\eta_A}:=\ell^*M_A\hbar,\\ \nabla_{-k}:={r_{\mbox{$D$}}} e^{{\text{j}}(\eta_D-k\omega\tau_d)}, {} {r_{\mbox{$D$}}} e^{{\text{j}}\eta_D}:=\ell^*M_D\hbar, $$

for integers k ≥ 1. For the leech CPG, we have ${r_{\mbox{$A$}}}={r_{\mbox{$D$}}}=2/3$, η _A  = π/3, and η _D  = 0. Note that $\nabla$ is m×m, where m is the number of segmental oscillators. To estimate the average phase lag, let us consider the limiting case where the oscillator chain is infinitely long (m = ∞ ). In this case, a generic row of the MHB equation $(\rho I-\nabla)\alpha=0$ takes the following form:

$$ \sum\limits_{k=-{q_{\mbox{$D$}}}}^{{q_{\mbox{$A$}}}}\nabla_k\alpha_k=0, {} \nabla_0:=-\rho/{\text{k}}. $$

(24)

Due to the uniformity of the intersegmental connections over the chain (that makes $\nabla$ a Toeplitz matrix), the eigenvector α has the structure such that α _k + 1/α _k is constant over k. Consequently, the intersegmental phase lag is uniform over the chain, and we may let

$$ \alpha_k=re^{-{\text{j}} k\eta_o} $$

(25)

where η _o is the intersegmental phase lag per segment. Substituting Eq. (25) into Eq. (24) and solving for ρ, we have

$$ \rho={\text{k}}{r_{\mbox{$A$}}}\sum\limits_{k=1}^{{q_{\mbox{$A$}}}}e^{{\text{j}}(\eta_A-k(\eta_o+\omega\tau_d))} +{\text{k}}{r_{\mbox{$D$}}}\sum\limits_{k=1}^{{q_{\mbox{$D$}}}}e^{{\text{j}}(\eta_D+k(\eta_o-\omega\tau_d))}.\\[-6pt] $$

(26)

This is a parametrization of the set of infinitely many eigenvalues of the infinite matrix $\nabla$ in terms of η _o  ∈ ℝ. The profile of a stable oscillation can be estimated from the maximal eigenvalue, and hence, we are interested in finding η _o such that $f(\eta_o):=\Re[\rho]$ takes its maximum value. Taking the derivative and setting it to zero, we have

$$\begin{array}{lll} f'(\eta_o)/{\text{k}} & = & {r_{\mbox{$A$}}}\sum\limits_{k=1}^{{q_{\mbox{$A$}}}}k\sin(\eta_A-k(\eta_o+\omega\tau_d))\\[-2pt] &&- {r_{\mbox{$D$}}}\sum\limits_{k=1}^{{q_{\mbox{$D$}}}}k\sin(\eta_D+k(\eta_o-\omega\tau_d)) \\[-2pt] & \approx & \left[ {r_{\mbox{$A$}}}\sum\limits_{k=1}^{{q_{\mbox{$A$}}}}(k\eta_A-k^2\omega\tau_d)- {r_{\mbox{$D$}}}\sum\limits_{k=1}^{{q_{\mbox{$D$}}}}(k\eta_D-k^2\omega\tau_d) \right]\\[-2pt] &&- \left[ {r_{\mbox{$A$}}}\sum\limits_{k=1}^{{q_{\mbox{$A$}}}}k^2+ {r_{\mbox{$D$}}}\sum\limits_{k=1}^{{q_{\mbox{$D$}}}}k^2 \right]\eta_o=0 \end{array}$$

where we used the approximation sin(x) ≈ x for small x. Solving this equation for η _o , we obtain the formula in Eq. (18). It can be verified that the η _o thus found approximates the global maximizer of f(η _o ).

Next, we consider the case where the intrinsic periods are not uniform. Assume that each segmental oscillator has the nominal time constant τ _o except that the pth segment has the time constant τ _p  ≠ τ _o . With this change, the basic coupling matrix, denoted by $\bar \nabla$, is slightly modified from $\nabla$ by

$$ \bar\nabla=\nabla+\xi_pe_pe_p{^{\mbox{\tt T}}} $$

where e _p is the vector having zero entries except for the pth entry, which is one, and ξ _p is defined using g _k (s) in Eq. (20) as

$$ \xi_p:={\text{k}}({\varepsilon}/\sigma)(\ell^*M\hbar)g_p({\text{j}}\omega). $$

(27)

Note that $\ell^*M\hbar=e^{{\text{j}}\pi/3}$ for the leech CPG. Now, the MHB equation becomes $(\bar\rho I-\bar\nabla)\bar\alpha=0$. We assume that the effect of the time constant perturbation is local in the sense that $\bar\alpha\approx\alpha+e_p\delta_p$ for some δ _p  ∈ ℂ, where α is the maximal eigenvector of $\nabla$ satisfying $\nabla\alpha=\rho\alpha$ for the eigenvalue ρ defined in Eq. (26). We further assume that the perturbation in the eigenvalue is small; $\bar\rho\approx\rho$. With these approximations, we have

$$ (\rho I-\nabla-\xi_pe_pe_p{^{\mbox{\tt T}}})(\alpha+e_p\delta_p)\approx0. $$

(28)

Multiplying the equation by $e_p{^{\mbox{\tt T}}}$ from the left and using $\nabla\alpha=\rho\alpha$, we obtain

$$ \delta_p\approx\varsigma_p\alpha_p, {} \varsigma_p:=\xi_p/\rho $$

where |ς _p | is assumed small. As a result, we have

$$ \theta_p=-p\eta_o+\Im(\varsigma_p), {} \theta_k=-k\eta_o \mbox{~~for~}k\neq p. $$

(29)

To derive a formula for the phase lags in the case of an arbitrary period gradient, we consider the time constant perturbation in a single oscillator as before and start with Eq. (28). Taking the summation of Eq. (28) over p = 1,..., ∞, we have

$$ (\rho I-\nabla-\Xi)\delta-\Xi\alpha\approx0 $$

where Ξ is the (infinite) diagonal matrix having entries ξ _p with p = 1,... ∞, and δ and α are (infinite) vectors with entries δ _p and α _p , respectively. This equation implies

$$ (\rho I-\nabla-\Xi)(\alpha+\delta)\approx(\rho I-\nabla)\alpha=0. $$

Hence, if the coupling matrix $\nabla$ is perturbed to $\nabla+\Xi$ due to the period gradient, the resulting phases are perturbed from $\angle\alpha$ to $\angle(\alpha+\delta)$ or $\angle\alpha+\Im(\varsigma)$. In summary, the intersegmental phases for general period gradients, perturbed from a nominal τ _o , are estimated as

$$ \theta_k=-k\eta_o+\Im(\varsigma_k), {} \varsigma_k:=\xi_k/\rho, $$

for all k where η _o is given by Eq. (16) and ρ and ξ _k are defined by Eqs. (26) and (27), respectively. Using the definition of g _k (s), the intersegmental phase lags are given by

$$ \begin{array}{lll} \eta_k&=&\eta_o+({\varepsilon}/\sigma)({\Delta}\tau_k/{\Delta}\tau_o)\varphi, \\ \varphi&:=&-\Im[{\text{k}}(\ell^*M\hbar){\text{j}}\tau_o\omega f'({\text{j}}\tau_o\omega)/\rho] \end{array} $$

(30)

where φ depends on the communication delay τ _d through ρ.