Stability and bifurcations in Wilson-Cowan systems with distributed delays, and an application to basal ganglia interactions

TL;DR

This paper investigates how distributed delays affect the stability and bifurcations of Wilson-Cowan neuronal models, providing a comprehensive mathematical analysis and applying it to basal ganglia circuits related to Parkinson's Disease.

Contribution

It extends Wilson-Cowan models to include general delay distributions, offering new stability and bifurcation results with rigorous proofs and numerical validation.

Findings

Delay distributions significantly influence system stability.

Hopf bifurcations are characterized under general delay conditions.

Application to basal ganglia illustrates relevance to Parkinson's Disease.

Abstract

The traditional Wilson-Cowan model of excitatory and inhibitory mean field interactions in neuronal populations considers a weak Gamma distribution of time delays when processing inputs, and is obtained via a time-coarse graining technique that averages the population response. Previous analyses of the stability of the Wilson-Cowan model focused on more simplified cases, where the delays were either not present, constant or were of a specific type. Since these simplifications may significantly alter the behavior of the model, we focus on understanding the behavior of the system before time-course graining, and for a wider range of delay distributions. For these generalized delay equations, we perform stability and bifurcation analyses with respect to parameters that capture both the coupling profile, and the time delay. The investigation is done through the examination of the system's associated characteristic equation. Under mild assumptions, we give complete mathematical proofs of our theoretical results, for the model with general delay distributions and prove the transversality condition for the possible Hopf bifurcations, in a generalized context. The stability region in this parameter space is described theoretically for several types of delay kernels, and numerical simulations are presented to substantiate the theoretical results. We illustrate these theoretical principles in an application to a basal ganglia circuit, in which $\beta$-band oscillations have been associated with Parkinson's Disease.

Tables

Table 1. Table 1: Synaptic connection weight parameters and other parameters estimated in [ 12 ]

Parameter	Healthy state	Diseased state	Other Parameters	Value
$w_{SS}$	$0$	$0$	$Ctx$	$27$ spk/s
$w_{SG}$	$19.0$	$20.0$	$Str$	$2$ spk/s
$w_{GS}$	$1.12$	$10.7$	$M_S$	$300$ spk/s
$w_{GG}$	$6.60$	$12.3$	$B_S$	$17$ spk/s
$w_{CS}$	$2.42$	$9.2$	$M_G$	$400$ spk/s
$w_{XG}$	$15.1$	$139.4$	$B_G$	$75$ spk/s

Equations

$$\left\{\begin{array}[]{l}\displaystyle E(t+{\tau_{1}})=\left(1-\int_{t-r}^{t}E(s)ds\right)\cdot\mathcal{S}_{e}\left[\int_{-\infty}^{t}h(t-s)\left(c_{1}E(s)-c_{2}I(s)+P_{e}(s)\right)ds\right]\\ \displaystyle I(t+{\tau_{2}})=\left(1-\int_{t-r^{\prime}}^{t}I(s)ds\right)\cdot\mathcal{S}_{i}\left[\int_{-\infty}^{t}h(t-s)\left(c_{3}E(s)-c_{4}I(s)+P_{i}(s)\right)ds\right]\end{array}\right.$$

$$\left\{\begin{array}[]{l}\displaystyle\dot{E}(t)=-{\tau_{1}}E(t)+\mathcal{S}_{e}\left[\int_{-\infty}^{t}h(t-s)\left(c_{1}E(s)-c_{2}I(s)+P_{e}(s)\right)ds\right]\\ \displaystyle\dot{I}(t)=-{\tau_{2}}I(t)+\mathcal{S}_{i}\left[\int_{-\infty}^{t}h(t-s)\left(c_{3}E(s)-c_{4}I(s)+P_{i}(s)\right)ds\right]\end{array}\right.$$

$$\left\{\begin{array}[]{l}\displaystyle\dot{u}(t)=-u(t)+f_{1}\left[\theta_{u}+\int\limits_{-\infty}^{t}h(t-s)\left(au(s)+bv(s)\right)ds\right]\\ \displaystyle\dot{v}(t)=-v(t)+f_{2}\left[\theta_{v}+\int\limits_{-\infty}^{t}h(t-s)\left(cu(s)+dv(s)\right)ds\right]\\ \end{array}\right.$$

$$\underset{0}{\overset{\infty }{\int }}h(s)ds=1,\text{with the average time delay}\tau =\underset{0}{\overset{\infty }{\int }}sh(s)ds<\infty .$$

$$u(s)={\psi }_1(s),v(s)={\psi }_2(s),\forall s\in (-\infty ,0],$$

$$\Vert \psi {\Vert }_{\infty ,\mu }=\underset{t\in (-\infty ,0]}{\sup }{e}^{\mu t}|\psi (t)|.$$

$$\left\{\begin{array}[]{l}\displaystyle u=f_{1}(\theta_{u}+au+bv)\\ \displaystyle v=f_{2}(\theta_{v}+cu+dv)\\ \end{array}\right.$$

$$\left\{\begin{array}[]{l}\displaystyle\dot{u}=-u+\phi_{1}\int\limits_{-\infty}^{t}h(t-s)\left(au(s)+bv(s)\right)ds\\ \displaystyle\dot{v}=-v+\phi_{2}\int\limits_{-\infty}^{t}h(t-s)\left(cu(s)+dv(s)\right)ds\\ \end{array}\right.$$

$$\left\{\begin{array}[]{l}\displaystyle zU(z)-u(0)=-U(z)+\phi_{1}H(z)\left(aU(z)+bV(z)\right)\\ \displaystyle zV(z)-v(0)=-V(z)+\phi_{2}H(z)\left(cU(z)+dV(z)\right)\\ \end{array}\right.$$

$$\left(\begin{array}[]{cc}z+1-a\phi_{1}H(z)&-b\phi_{1}H(z)\\ -c\phi_{2}H(z)&z+1-d\phi_{2}H(z)\\ \end{array}\right)\left(\begin{array}[]{c}U(z)\\ V(z)\\ \end{array}\right)=\left(\begin{array}[]{c}u(0)\\ v(0)\\ \end{array}\right)$$

$$\Delta (z)=(z+1{)}^2-\alpha H(z)(z+1)+\beta H^2(z)=0$$

$$\alpha$$

$$\beta$$

$$\alpha <\min \{2,\beta +1\}$$

$$|\alpha |+|\beta |<1$$

$$\beta <\alpha -1$$

$$\Delta (z)=z^2+(2-\alpha )z+\beta -\alpha +1=0$$

$$|H(z)|=\left|\underset{0}{\overset{\infty }{\int }}{e}^{-zt}h(t)dt\right|\le \underset{0}{\overset{\infty }{\int }}|e^{-zt}|h(t)dt=\underset{0}{\overset{\infty }{\int }}{e}^{-\Re (z)t}h(t)dt.$$

$$|H(z)|\le \underset{0}{\overset{\infty }{\int }}h(t)dt=1.$$

$$|z+1{|}^2$$

$$\le |\alpha ||z+1|+|\beta |$$

$$H(0)=\underset{0}{\overset{\infty }{\int }}h(t)dt=1,\text{for any }t\ge 0.$$

$$\Delta (0)=1-\alpha H(0)+\beta H^2(0)=1-\alpha +\beta .$$

$$\beta =\alpha -1$$

$$2(z+1)\frac{dz}{d\beta }-\alpha H^{\prime }(z)(z+1)\frac{dz}{d\beta }-\alpha H(z)\frac{dz}{d\beta }+2\beta H^{\prime }(z)H(z)\frac{dz}{d\beta }+H^2(z)=0,$$

$$\frac{dz}{d\beta }=-\frac{H^2(z)}{2(z+1)-\alpha H^{\prime }(z)(z+1)-\alpha H(z)+2\beta H^{\prime }(z)H(z)}.$$

$${\frac{dz}{d\beta }|}_{\beta =\alpha -1}=\frac{1}{(\alpha -2)(\tau +1)}\ne 0.$$

$$\underset{0}{\overset{\infty }{\int }}t\hat{h}(t)dt=\underset{0}{\overset{\infty }{\int }}t\tau h(\tau t)dt={\tau }^{-1}\underset{0}{\overset{\infty }{\int }}uh(u)du={\tau }^{-1}\tau =1.$$

$$\hat{H}(z)=\underset{0}{\overset{\infty }{\int }}{e}^{-zt}\hat{h}(t)dt=\underset{0}{\overset{\infty }{\int }}{e}^{-zt}\tau h(\tau t)dt=\underset{0}{\overset{\infty }{\int }}{e}^{-z\frac{u}{\tau }}h(u)du=H\left(\frac{z}{\tau }\right)$$

$$\hat{H}(i\omega )=\rho (\omega )e^{-i\theta (\omega )}$$

$$(A_1):\rho \text{ and }\theta \text{ are of class C1 on R};$$

Full text

For these generalized delay equations, we perform stability and bifurcation analyses with respect to parameters that capture both the coupling profile, and the time delay. The investigation is done through the examination of the system’s associated characteristic equation. Under mild assumptions, we give complete mathematical proofs of our theoretical results, for the model with general delay distributions and prove the transversality condition for the possible Hopf bifurcations, in a generalized context. The stability region in this parameter space is described theoretically for several types of delay kernels, and numerical simulations are presented to substantiate the theoretical results.

We found that the stability regions and bifurcations differ significantly between different types of delay distributions: weak Gamma distributions promote stable firing rates, while strong Gamma distributions are associated with regular oscillations, and Dirac distributions appear to facilitate more complex , aperiodic patterns. This supports the unexplored possibility of different delay distributions being used as the substrate for different functional behaviors, and emphasizes the importance of a careful choice of the delay kernel in the mathematical model.

We illustrate these theoretical principles in an application to a basal ganglia circuit, in which $\beta$-band oscillations have been associated with Parkinson’s Disease.

1 Introduction

1.1 Modeling background

Computational modeling of neuronal behavior covers a large range of spatio-temporal scales, from detailed, membrane potential based models of single spiking neurons, to broad network models of interacting brain regions. At the lower scale, typically associated with electrode recordings from single cell in vitro preparations, modeling difficulties often arise from the inherent complexity and high dimensionality associated with considering detailed molecular mechanisms that govern ionic currents and spiking activity in a single cell. The Hodgkin-Huxley model was in its original form limited to the two voltage-dependent currents found in the squid giant axon, but it has to be extended to dozens of equations per neuron if it includes other ion channels involved in neuronal excitability [1, 2]. In the context of studying behavior in functional neuronal networks, which may involve thousands of neurons, the dimensionality of a network formed of single cells may become an obstacle to computational feasibility. At the opposite end, often associated with functional imaging data, modeling activity within entire brain regions as a whole lacks specificity, and prevents clear interpretations of what “activity” of a state variable may actually represent [3].

The middle range between the two ends consists of a rich variety of models. One wide-spread possibility has been creating reduced models as modifications of Hodgkin-Huxley equations, by modeling the behavior of multiple ion channels into one comprehensive variable [4]. Another practical option has been using state variables to characterize the meanfield spiking activity in a population of cells [5]. This type of model is still able to incorporate information on spiking mechanisms, and efficiently illustrates resulting firing patterns by using only one variable per population.

One of the most historically significant and best understood meanfield models, the Wilson-Cowan model [6] is perhaps the most popular. The model, derived in 1972, describes the localized interactions in a pair of excitatory and inhibitory neuronal populations. At each time instant $t$, the proportions of excitatory and inhibitory cells firing per unit of time are captured by the two state variables $E(t)$ and $I(t)$. The original model considers the effect that an external input $P$ has on the E/I system, based not only on the coupling strengths between the two units, but also on the history of firing in each. More specifically, $E(t+{\tau_{1}})$ and $I(t+{\tau_{2}})$ represent the proportion of cells which are sensitive (i.e. not refractory) and which also receive at least threshold excitation at the moment of time $t$. This leads to the following system of integral equations:

[TABLE]

Here, the first factor on each right hand side represents the proportion of sensitive excitatory / inhibitory cells, where $r$ and $r^{\prime}$ are the absolute refractory periods (msc). The functions $\mathcal{S}_{e}$, $\mathcal{S}_{i}$ are sigmoid threshold functions, their arguments denoting the mean field level of excitation / inhibition generated in an excitatory /inhibitory cell at time $t$. The summation coefficients $c_{i}>0$ are connectivity weights representing the average number of excitatory / inhibitory synapses per cell and $P_{e},P_{i}$ denote external inputs. It is additionally assumed that a cell’s weighted summation of inputs is modulated longitudinally by a decreasing kernel $h(t)$. The original paper [6] suggests the use of a weak Gamma kernel $h(t)$, before proceeding with further mathematical simplifications, for convenience of the analysis.

Indeed, with the additional simplifying assumption of negligible refractory periods ($r=r^{\prime}=0$) and with a degree one Taylor approximation of the left-hand terms for small ${\tau_{1}}$ and ${\tau_{2}}$, the Equations (1) can be easily rewritten as the following system of intergo-differential equations:

[TABLE]

For simplicity, the equations are often considered to have the same time scale ${\tau_{1}=\tau_{2}=1}$. The model historically used as the Wilson-Cowan model [6] was then obtained from (2) by applying further time-coarse graining. This final form, consisting of a system of non-delayed ordinary differential equations, is very convenient and was extensively analyzed and used in the modeling literature [7]. However, this final course-graining step obscures potentially vital information on the the shape of the temporal integration of synaptic inputs in both populations (by assuming, for example, that $h(t)$ is close to one for $t<r,r^{\prime}$ respectively, and that is decays rapidly to zero for $t>r,r^{\prime}$). This information may be crucial to the neural function that the model is aiming to address, hence a few generalizations of the standard model found it useful to further explore the original equations (prior to coarse graining), imposing loser conditions on temporal integration, such as discrete time-delays, and often considering null refractory periods $r,r^{\prime}$ [8, 9, 10, 11].

However, constant delays do not constitute the most biologically realistic possibility. Since $h(t-s)$ represents the effect of presynaptic input at time s on the membrane potential at time t, it can be thought to include two aspects: the synaptic transmission delay, and the time-course of the post-synaptic potential (PSP). In this context, constant delays imply not only equal transmission delays between the neurons, but also instantaneous PSPs. While often delays due to PSP are ignored, and a constant kernel is chosen on the basis on experimentally observed synaptic transmission delays [12], such choices are rather made in the interest of simplicity, and do not incorporate the full picture, and other types of distributed delays may deliver a more realistic model of synaptic transmission delays.

In this paper, we consider and analyze the Wilson-Cowan model in its original integro-differential form, for a general delay distribution and with the only simplifying assumption (retained for simplicity) being that of zero refractory periods $r=r^{\prime}=0$. We then compare the results for several types of delay kernels. We analyze the long-term behavior in each case, then we discuss the main differences between these behaviors, and potential advantages of using delay equations versus the simpler and more main-stream course grain approximation.

We emphasize, that our main theorems have already been briefly presented (without proofs) in [13], accompanied by several numerical simulations, mostly in the case of a discrete time-delay. Thus, in this paper we illustrate rigorous mathematical proofs of our previously presented statements [13] and of our additional lemmas and propositions added here. Choosing different parameters similar numerical simulations are presented, as in [13, 8]. Finally a Parkinson’s Disease (PD) model of the basal ganglia is adapted after [12], followed by a series of simulations, treating both the case of healthy basal ganglia and a diseased/parkinsonian one, using different delay distributions. We then present detailed neuro-biological interpretations and argumentations of our simulations of the PD model.

1.2 Our model

To address the integral terms from the original equations  (2) we analyze in this paper the general model with distributed delays defined by

[TABLE]

where $u(t)$ and $v(t)$ represent the firing activity in the two neuronal populations, $a,b,c,d$ are the connection weights and $\theta_{u}$, $\theta_{v}$ are background drives. The activation functions $f_{1}$ and $f_{2}$ are smooth and increasing on the real line. The original model (2) can be recovered for equal time scales ${\tau_{1}=\tau_{2}}$, if the activation functions are particularly set to sigmoidals $f_{1}=\mathcal{S}_{e}$ and $f_{2}=\mathcal{S}_{i}$, with background drives $\displaystyle\theta_{u}=\int_{-\infty}^{t}h(t-s)P_{e}(s)ds$ and $\displaystyle\theta_{v}=\int_{-\infty}^{t}h(t-s)P_{i}(s)ds$, respectively.

The delay kernel $h:[0,\infty)\to[0,\infty)$ in system (3) is suggested to have an exponential form (weak Gamma kernel) in the original Wilson-Cowan reference [6]. We will consider the extensive case of a general probability density function, representing the probability that a particular time delay occurs. The delay kernel is considered to be bounded and piecewise continuous, satisfying

[TABLE]

Our analysis of the system (3) will focus around understanding the changes in its long term behavior that are triggered by changes in spatial and temporal summation of inputs (as reflected in the synaptic weights and in the integration kernel).

In recent years, significant work has been dedicated to understanding the effects of spatial integration of inputs [14, 15]. An entire direction in computational neuroscience is centered around investigating the effects of connectivity and synaptic weights on dynamics in coupled neural populations [16]. Seminal empirical work has shown that the synaptic weight profile in a network changes with location [17], with task-related circumstances (e.g., when learning and memory formation are involved [18]), but also with general circumstances (e.g., sleep and circadian rhythms [19, 20]). such modeling work is shedding increasingly more light on the mechanisms by which different synaptic landscapes modulate different outcomes in networks of neuronal populations.

Comparatively, less is known about temporal profiles in the integration of inputs. While physiological evidence supports that exponential decay may represent an appropriate model for certain type of neural contexts [21], other modulation patterns have been explored [22, 23]. It is generally agreed that effects of stimulation decay with the time course, but it is quite possible that the timeline and profile of this decay (the shape of the post-synaptic integration) depend on the properties of the neurons involved (receptor types, voltage-dependent currents in the soma and dendrites) [24].

The Wilson-Cowan type system (3) simplifies the behavior of a self-interacting E/I network into a mean-field, two-dimensional model without the geometric, spatial structure described in other, more detailed models. The distribution of synaptic weights contributing to spatial summation of inputs is then entirely encapsulated in the set of mean field coefficients $a,b,c,d$. The temporal summation of inputs is described by the kernel $h$. In this paper, we study, for different types of kernels, the types of asymptotic regimes accessible to the system. For each kernel type, we further study the conditions on the synaptic coefficients and distributed delay that produce transitions between these regimes.

The specific case of discrete time delays (Dirac kernels) has been previously analyzed in [8, 25]. However, as previously mentioned, the discrete time integration model does not have a strong physiological backup; other important types of delay kernels are frequently used in the literature, such as uniform distribution kernels [26, 27, 28] or Gamma kernels (recall that the original Wilson-Cowan model itself proposed a weak Gamma (exponential) kernel $h(t)=\tau^{-1}\exp(-t/\tau)$) [6]. Analyzing and comparing mathematical models which incorporate different types of delay kernels (e.g. weak Gamma kernel or strong Gamma kernel $h(t)=4\tau^{-2}t\exp(-2t/\tau)$) may shed a light on the difference between using different distributed delays, as well as on the differences between using continuous versus discrete delays on the system’s dynamics. Furthermore, when modeling natural systems, one usually does not have access to the exact delay distribution, and approaches using general kernels may prove to be more revealing [29, 28, 26, 30, 31, 27, 32, 33, 34].

The results in Section 2 are delivered within this general context, keeping additional conditions on the kernels to a minimum. In Sections 3.1 and 3.2, we further illustrate the behaviors for specific theoretical and biologically driven examples, using particular kernels (the weak and strong Gamma kernels and the Dirac kernel). In Section 4, we draw comparisons between the behaviors reported for different types of distributed delays, and discuss potential contributions of the delay scheme to dynamic control mechanisms in coupled neuronal populations.

2 Stability and bifurcation results

The initial conditions that are associated to system (3) are:

[TABLE]

where $\psi_{1}$ and $\psi_{2}$ belong to the Banach space $C_{0,\mu}(\mathbb{R}_{-},\mathbb{R})$ (where $\mu>0$) of continuous real valued functions defined on $(-\infty,0]$ such that $\lim\limits_{t\rightarrow-\infty}e^{\mu t}\psi(t)=0$, endowed with the norm:

[TABLE]

For existence and uniqueness results regarding the solution of the initial value problem associated to system (3) we cite [35, 36]. An important tool in the stability analysis of system (3) is the principle of linearized stability (see for example [37]), characterizing the local stability properties of an equilibrium state of (3) in terms of the stability of the null solution of the linearized system at that particular equilibrium. For the principle of linearized stability in the context of systems of differential equations with infinite time delay, we refer to [30].

The equilibrium states of system (3) are found by solving the following algebraic system:

[TABLE]

Linearizing at an equilibrium state $(u^{\star},v^{\star})$ leads to

[TABLE]

where $\phi_{1}=\phi_{1}(u^{\star},v^{\star})=f_{1}^{\prime}(\theta_{u}+au^{\star}+bv^{\star})>0$ and $\phi_{2}=\phi_{2}(u^{\star},v^{\star})=f_{2}^{\prime}(\theta_{v}+cu^{\star}+dv^{\star})>0$.

Either looking for solutions of the form $(u(t),v(t))=e^{zt}(u(0),v(0))$ or directly employing the Laplace transform technique to the linearized system (6), we obtain:

[TABLE]

where $U(z)$, $V(z)$ and $H(z)$ denote the Laplace transforms of state variables $u$ and $v$, and of the delay kernel $h$, respectively.

System (7) is equivalently written as:

[TABLE]

which leads to the following characteristic equation (associated to the equilibrium $(u^{\star},v^{\star})$):

[TABLE]

where

[TABLE]

2.1 Delay-independent stability and instability results

We first discuss several delay-independent stability and instability properties of system (3), emphasizing that these results are particularly useful when the exact delay distributions are unknown.

Theorem 2.1**.**

Assume that the delay kernel $h(t)$ in system (3) satisfies the properties (4).

In the non-delayed case (i.e. $h(t)=\delta(t)$), the equilibrium $(u^{\star},v^{\star})$ of system (3) is locally asymptotically stable if and only if the following inequality holds:

[TABLE] 2. 2.

If the following inequality holds:

[TABLE]

then the equilibrium $(u^{\star},v^{\star})$ of system (3) is locally asymptotically stable, regardless of the choice of the delay kernel $h(t)$. 3. 3.

If the following inequality holds:

[TABLE]

then the equilibrium state $(u^{\star},v^{\star})$ of system (3) is unstable for any delay kernel $h(t)$.

Proof.

1. In the non-delayed case ($H(z)=1$, for any $z\in\mathbb{C}$), the characteristic equation (9) becomes:

[TABLE]

The necessary and sufficient condition (10) for the asymptotic stability of the equilibrium state $(u^{\star},v^{\star})$ of system (3) follows from the Routh-Hurwitz stability test.

2. Using basic inequality techniques, we have:

[TABLE]

Now, if we assume that the characteristic equation (9) has a root $z$ in the right half-plane ($\Re(z)\geq 0$), it follows that $e^{-\Re(z)t}\leq 1$, for any $t\geq 0$ and hence:

[TABLE]

From (9) we deduce:

[TABLE]

Considering the polynomial $P(x)=x^{2}-|\alpha|x-|\beta|$, from inequality (11) it follows that $P(1)>0$ and $P^{\prime}(1)>0$, and hence $P(x)>0$ for any $x\geq 1$. From the above inequality, we have $P(|z+1|)\leq 0$, and hence, we deduce $|z+1|<1$, which is absurd, since $|z+1|^{2}=|z|^{2}+2\Re(z)+1\geq 1$.

Therefore, all the roots of the characteristic equation (9) are in the left half-plane and the equilibrium state $(u^{\star},v^{\star})$ of system (3) is asymptotically stable, regardless of the delay kernel $h(t)$.

Using the delay kernel properties (4) and the Laplace transform definition, we have

[TABLE]

Based on (9), it follows that

[TABLE]

Thus, condition (12), i.e. $\beta<\alpha-1$ is equivalent to $\Delta(0)<0$. Given that the characteristic function $\Delta(z)$ is continuous on $[0,\infty)$ and $\Delta(z)\rightarrow\infty$ as $z\rightarrow\infty$, the characteristic equation (9) has at least one positive real root. Hence, the equilibrium state $(u^{\star},v^{\star})$ of system (3) is unstable, regardless of the delay kernel $h(t)$. ∎

2.2 Saddle-node bifurcation

Theorem 2.2** (Saddle-node bifurcation).**

A saddle-node bifurcation takes place at the equilibrium state $(u^{\star},v^{\star})$ of system (3), regardless of the delay kernel $h(t)$, if and only if $\alpha\neq 2$ and

[TABLE]

Proof.

Via relationship (13), one can see that, condition (14) is equivalent to $\Delta(0)=0$ (i.e. the characteristic equation has a zero root). Moreover, $z=0$ is a simple root of the characteristic equation if and only if $\alpha\neq 2$.

Let us denote by $z(\beta)$ the root of the characteristic equation (9) which satisfies $z(\alpha-1)=0$. Taking the derivative with respect to $\beta$ in equation (9), we obtain:

[TABLE]

and hence:

[TABLE]

As $H(0)=1$ and $H^{\prime}(0)=-\tau$, it follows that:

[TABLE]

This completes the proof of transversality. While we have not verified the conditions for nondegeneracy (which are rather difficult to check analytically), one would expect the nonlinearity (hence the saddle-node bifurcation) to be generic. ∎

2.3 Hopf bifurcation

In the following, we will show that the average time delay $\tau$ of the delay kernel $h(t)$ plays an important role in the Hopf bifurcation analysis.

Let $\hat{h}(t)=\tau h(\tau t)$, for any $t\geq 0$. The function $\hat{h}$ is a probability density function with the mean value:

[TABLE]

The Laplace transform $\hat{H}(z)$ of $\hat{h}(t)$ is

[TABLE]

We assume from now on that $\hat{H}(z)$ is independent of the average time delay $\tau$. In fact, we note that this is true for several classes of delay kernels, such as:

• •

for the Dirac kernel (discrete time delay): $\displaystyle\hat{H}(z)=e^{-z}$;

• •

for the $p$-Gamma kernel: $\displaystyle\hat{H}(z)=\left(\frac{p}{p+z}\right)^{p}$.

Based on [38], we write $\hat{H}$ in the polar form as:

[TABLE]

with $\rho(0)=1$, $\theta(0)=0$. We further assume that:

[TABLE]

It is easy to verify that for the above mentioned probability densities, these assumptions are satisfied. It is also important to note that $\omega\mapsto\hat{H}(-i\omega)$ is in fact the characteristic function of a random variable with the mean equal to $1$.

With the change of variable $z\mapsto\displaystyle\frac{z}{\tau}$, the characteristic equation (9) becomes

[TABLE]

Denoting $\displaystyle Q_{\tau}(z)=\frac{z+\tau}{\tau\hat{H}(z)}$, it follows that the characteristic equation is equivalent to

[TABLE]

The characteristic equation (9) has a pair of complex conjugate roots on the imaginary axis if and only if there exists $\omega>0$ such that $Q(i\omega)$ is a root of the polynomial $P(\lambda)=\lambda^{2}-\alpha\lambda+\beta$.

Case 1: $\alpha^{2}-4\beta<0$.

In this case, the polynomial $P(\lambda)$ has complex conjugate roots $Q_{\tau}(i\omega)$ and $\overline{Q_{\tau}(i\omega)}$, and hence, Hopf bifurcation may only take place at the equilibrium state $(u^{\star},v^{\star})$ of system (3) along the curve $(\gamma_{\tau})$ from the $(\alpha,\beta)$-plane, defined by the following parametric equations:

[TABLE]

Indeed, we prove the following:

Lemma 2.3**.**

If $\alpha^{2}-4\beta<0$ then the characteristic equation (9) has a pair of complex conjugate roots on the imaginary axis if and only if $(\alpha,\beta)$ belong to the curve $(\gamma_{\tau})$ defined by (17).

Proof.

Since $\alpha^{2}-4\beta<0$, the characteristic equation (9) has a pair of complex conjugate roots on the imaginary axis if and only if there exists $\omega>0$ such that the polynomial $P(\lambda)$ has complex conjugate roots $Q_{\tau}(i\omega)$ and $\overline{Q_{\tau}(i\omega)}$. As

[TABLE]

from Vieta’s formulas, we obtain:

[TABLE]

Therefore, we obtain the equivalence from the statement. ∎

Case 2: $\alpha^{2}-4\beta\geq 0$.

In this case, the following preliminary lemma will be needed to establish further results:

Lemma 2.4**.**

Let $\tilde{\theta}=\lim\limits_{\omega\rightarrow\infty}\theta(\omega)$ and $K=\lceil\tilde{\theta}/\pi\rceil$ (if $\tilde{\theta}=\infty$ then $K=\infty$). The function

[TABLE]

has exactly $N$ roots $0=\omega_{0}(\tau)<\omega_{1}(\tau)<...<\omega_{N-1}(\tau)$, such that $\theta(\omega_{n}(\tau))\in\left(n\pi-\dfrac{\pi}{2},n\pi\right)$, for $n\in\{1,2,...,N-1\}$, where:

[TABLE]

Proof.

From $\theta(0)=0$ it follows that $\omega_{0}(\tau)=0$ is a root of the function $q_{\tau}(\omega)$. As $\theta$ is increasing, for any $n\in\{1,2,...,K\}$ it is easy to see that

[TABLE]

Therefore, for each $n\in\{1,2,...,K-1\}$, there exists a unique root $\omega_{n}(\tau)$ of the function $q_{\tau}(\omega)$ which satisfies $\theta(\omega_{n}(\tau))\in\left(n\pi-\dfrac{\pi}{2},n\pi\right)$. The $(K+1)$-th root exists if and only if $(-1)^{K}q_{\tau}(\theta^{-1}(\tilde{\theta}))>0$. ∎

The following result holds:

Lemma 2.5**.**

If $\alpha^{2}-4\beta\geq 0$ then the characteristic equation (9) has a pair of complex conjugate roots $\pm i\omega\tau^{-1}$ on the imaginary axis if and only if $\omega$ is a positive root of (18) and $(\alpha,\beta)$ belong to the line

[TABLE]

Proof.

Since $\alpha^{2}-4\beta\geq 0$, the characteristic equation (9) has a pair of complex conjugate roots on the imaginary axis if and only if there exists $\omega>0$ such that $Q_{\tau}(i\omega)$ is a real root of the polynomial $P(\lambda)$. Therefore, $\Im(Q(i\omega))=0$ and hence, it follows that $\omega$ is a root of the equation

[TABLE]

Assuming that such a root exists, let us denote it by $\omega^{*}$. Hence:

[TABLE]

As $Q_{\tau}(i\omega^{\star})$ is a root of $P(\lambda)$, we obtain:

[TABLE]

Therefore, the proof is complete. ∎

Lemma 2.6**.**

The characteristic equation (9) has a double pure imaginary root if and only if $\alpha^{2}=4\beta$.

Proof.

Assume that there exists $\omega>0$ such that $z=i\omega\tau^{-1}$ is a double root of (9), i.e.

[TABLE]

Assumption $(A_{2})$ guarantees that $H(z)\neq 0$. Eliminating $\beta$ from system (20) we obtain the equality:

[TABLE]

On one hand, if the first term in the above equation is zero, replacing into the first equation of system (20), it follows that $\alpha^{2}=4\beta$. On the other hand, if the second term of the above equation is zero, as $z=i\omega\tau^{-1}$ and $H(i\omega\tau^{-1})=\rho(\omega)e^{\displaystyle-i\theta(\omega)}$, we obtain:

[TABLE]

Taking the imaginary part in the previous equality, leads to $\tau\rho^{\prime}(\omega)=-\omega\rho(\omega)\theta^{\prime}(\omega)$. Next, taking the real part, we get

[TABLE]

From assumption $(A_{2})$, we have $\rho(\omega)>0$ and from $(A_{3})$, we have $\rho^{\prime}(\omega)\leq 0$. Therefore, the above equality is a contradiction. In conclusion, the characteristic equation (9) has a double pure imaginary root if and only if $\alpha^{2}=4\beta$. ∎

Remark 2.1**.**

From Lemma 2.6 it follows that the characteristic equation (9) has a double pure imaginary root if and only if $(\alpha,\beta)$ are at the intersection points of the curve $(\gamma_{\tau})$ with the lines provided by Lemma 2.5 (if any, depending on the existence of the roots of equation (18)). These intersection points lie on the parabola $\alpha^{2}=4\beta$ (as in Proposition 2.7).

2.4 Main results

Combining the previous four lemmas, we deduce that the lines defined for $n\in\{0,1,...,N-1\}$ by

[TABLE]

as well as the curve $(\gamma_{\tau})$ defined parametrically by (17) play an important role in the bifurcation analysis. Hence, some remarkable properties of these lines and curve are proved in the following:

Proposition 2.7**.**

For the lines $(l_{n})$ defined by (21) and the curve $(\gamma_{\tau})$ defined by (17) the following properties hold:

• i.

$|\mu(\omega_{n}(\tau))|$* is increasing with respect to $n$, with $\text{sign}(\mu(\omega_{n}(\tau)))=(-1)^{n}$ and $|\mu(\omega_{n}(\tau))|>1$;*

• ii.

the line $(l_{n})$ is included in the region $\alpha^{2}-4\beta\geq 0$ and is tangent to the parabola $\alpha^{2}=4\beta$ at the point ${P_{n}}\left(2\mu(\omega_{n}(\tau)),\mu(\omega_{n}(\tau))^{2}\right)$;

• iii.

the intersection of the line $(l_{n})$ with the vertical axis $\alpha=0$ is the point of $\beta$-coordinate $-\mu(\omega_{n}(\tau))^{2}$.

• iv.

the curve $(\gamma_{\tau})$ is a simple curve included in the region $\alpha^{2}-4\beta\leq 0$;

• v.

the line $(l_{n})$, the curve $(\gamma_{\tau})$ and the parabola $\alpha^{2}=4\beta$ intersect at the point ${P_{n}}\left(2\mu(\omega_{n}(\tau)),\mu(\omega_{n}(\tau))^{2}\right)$.

Proof.

For the proof of (i.) From equation (18) it easily follows that

[TABLE]

and hence:

[TABLE]

Therefore, as $\omega\mapsto\rho(\omega)$ is decreasing, we deduce that the function $n\mapsto|\mu(\omega_{n}(\tau))|$ is increasing. Moreover, from Lemma 2.4 it follows that $\theta(\omega_{n}(\tau))\in\left(n\pi-\dfrac{\pi}{2},n\pi\right)$ and hence $\text{sign}(\mu(\omega_{n}(\tau)))=\text{sign}(\cos(\theta(\omega_{n}(\tau))))=(-1)^{n}$.

Moreover, assumption $(A_{3})$ implies the function $\beta_{\tau}(\omega)$ is strictly increasing, hence, the curve $(\gamma_{\tau})$ does not cross itself, hence statement (iv.) follows.

The rest of the statements can be easily proved by elementary reasoning. ∎

We are now ready to prove the main theoretical results, which are based on the previously described root locus technique.

Theorem 2.8**.**

If $N\geq 2$, denoting

[TABLE]

the stability region $S_{\tau}(\alpha,\beta)$ of the equilibrium $(u^{\star},v^{\star})$ of system (3) is bounded by the line segments and curve given below:

[TABLE]

At the boundary of the stability region $S_{\tau}(\alpha,\beta)$, the following bifurcation phenomena take place in a neighborhood of the equilibrium $(u^{\star},v^{\star})$ of system (3):

• a.

Saddle-node bifurcations take place along the open line segment $(l_{0})$;

• b.

Hopf bifurcations take place along the open line segment $(l_{\tau})$ and curve $(\gamma_{\tau})$;

• c.

Bogdanov-Takens bifurcation at $(\alpha,\beta)=(2,1)$;

• d.

Double-Hopf bifurcation at $(\alpha,\beta)=\left(2\mu_{\tau},\mu_{\tau}^{2}\right)$;

• e.

Zero-Hopf bifurcation $(\alpha,\beta)=\left(1+\mu_{\tau},\mu_{\tau}\right)$.

Proof.

It is easy to see that for $\alpha=\beta=0$, the characteristic equation has negative real roots, and therefore, the equilibrium state $(u^{\star},v^{\star})$ of system (3) is asymptotically stable. As the roots of the characteristic function $\Delta(z)$ (or $\hat{\Delta}(z)$) continuously depend on the parameters $\alpha$ and $\beta$, the number of roots $z$ such that $\Re(z)>0$ may change if and only if a root $z=0$ or a pair of pure imaginary roots appear, i.e. along the curve $(\gamma_{\tau})$ or the lines $(l_{n})$, $n\in\{0,1,...,N-1\}$ in the $(\alpha,\beta)$-plane. A simple geometrical reasoning based on Proposition 2.7 shows that the line segments and curve segment given in the statement above enclose a connected region of the $(\alpha,\beta)$-plane containing the origin, i.e. the stability region $S_{\tau}(\alpha,\beta)$ of the equilibrium $(u^{\star},v^{\star})$ of system (3). Moreover, Theorem 2.2 shows that a saddle-node bifurcation takes place along the line segment $(l_{0})$, and hence, statement a. holds.

b. We will first show that Hopf bifurcations take place along the curve segment $(\gamma_{\tau})$, with $\omega\in(0,\omega_{\tau})$.

Lemma 2.3 provides that the characteristic equation (16) has a pair of pure imaginary roots $\pm i\omega$ if $(\alpha,\beta)$ belong to the curve $(\gamma_{\tau})$. Let us consider an arbitrary $\omega\in(0,\omega_{\tau})$ and denote by $z(\alpha,\beta)$ the root of the characteristic equation (16) satisfying $z(\alpha^{*},\beta^{*})=i\omega$ where $(\alpha^{*},\beta^{*})=(\alpha_{\tau}(\omega),\beta_{\tau}(\omega))\in(\gamma_{\tau})$. Our aim is to prove the transversality condition:

[TABLE]

for any outward pointing vector $\overline{u}$ from the region $S_{\tau}(\alpha,\beta)$, i.e. $\langle\overline{u},\overline{n}\rangle>0$, where $\overline{n}$ denotes the outward pointing normal vector to the curve $(\gamma_{\tau})$.

From (16) it follows that $z(\alpha,\beta)$ is a solution of the equation

[TABLE]

Taking the derivative with respect to $\alpha$, we obtain:

[TABLE]

and we express:

[TABLE]

Therefore, based on the parametric equations of the curve $(\gamma_{\tau})$ given by (17), we deduce:

[TABLE]

Applying the real part, we finally obtain:

[TABLE]

Taking the derivative with respect to $\beta$ in equation (22), we obtain:

[TABLE]

and we express:

[TABLE]

Again, using the expressions of the parametric curve $(\gamma_{\tau})$ given by (17), we deduce:

[TABLE]

Once again, applying the real part, we arrive at:

[TABLE]

and thus, we obtain the gradient vector:

[TABLE]

The tangent vector to the curve $(\gamma_{\tau})$ at the point $(\alpha^{*},\beta^{*})$ is $(\alpha^{\prime}(\omega),\beta^{\prime}(\omega))$, where

[TABLE]

Thus, we obtain the following tangent vector to the curve $(\gamma_{\tau})$:

[TABLE]

Therefore, fixing the orientation of the curve $(\gamma_{\tau})$ in the direction of increasing $\omega$, a right-pointing normal vector to the curve $(\gamma_{\tau})$ is:

[TABLE]

Hence, the directional derivative is:

[TABLE]

as $\rho(\omega)>0$, $\langle\overline{n}(\omega),\overline{u}\rangle>0$ and $q_{\tau}(\omega)>0$ for any $\omega\in(0,\omega_{\tau})$, from the proof of Lemma 2.4. Therefore, the transversality condition holds, and it follows that a Hopf bifurcation takes place along the curve segment $(\gamma_{\tau})$ which bounds the stability region $S_{\tau}(\alpha,\beta)$.

In general, along other segments of the curve $(\gamma_{\tau})$, it is important to emphasize that if $\omega\in(\omega_{n}(\tau),\omega_{n+1}(\tau))$, where $n\in\{0,1,...,N-1\}$ and $z(\alpha,\beta)$ denotes the root of the characteristic equation (16) satisfying $z(\alpha^{*},\beta^{*})=i\omega$ where $(\alpha^{*},\beta^{*})=(\alpha_{\tau}(\omega),\beta_{\tau}(\omega))\in(\gamma_{\tau})$, by a similar reasoning as above, we get

[TABLE]

where the orientation of the curve $(\gamma_{\tau})$ is fixed in the direction of increasing $\omega$, and $\overline{n}(\omega)$ is a right-pointing normal vector to the curve $(\gamma_{\tau})$ (e.g. see Figure 2). From Lemma 2.4, it follows that $\text{sign}(q_{\tau}(\omega))=(-1)^{n}$ and we can also notice that $\text{sign}(\langle\overline{n}(\omega),\overline{u}\rangle)=(-1)^{n}$, for any vector $\overline{u}$ pointing towards the region above the curve segment of $(\gamma_{\tau})$, with $\omega\in(\omega_{n}(\tau),\omega_{n+1}(\tau))$. Therefore, $\nabla_{\overline{u}}\Re(z)(\alpha^{*},\beta^{*})>0$, and hence, stability cannot be regained when crossing the curve segment from the connected region below the curve segment to the region above.

In the following, we will show that Hopf bifurcations take place along the line segment $(l_{\tau})$, with $\alpha\in\left[2\mu_{\tau},1+\mu_{\tau}\right]$. If $(\alpha,\beta)$ belong to this line segment, Lemma 2.5 shows that the characteristic equation (16) has a pair of pure imaginary roots $\pm i\omega_{\tau}$. Let us denote by $z(\alpha,\beta)$ the root of the characteristic equation (16) with the property $z(\alpha^{*},\beta^{*})=i\omega_{\tau}$ where $\displaystyle\alpha^{*}\in\left[2\mu_{\tau},1+\mu_{\tau}\right]$ is arbitrarily fixed and $\displaystyle\beta^{*}=\mu_{\tau}(\alpha-\mu_{\tau})$. We will prove the transversality condition

[TABLE]

for any vector $\overline{u}$ pointing outward from the region $S_{\tau}(\alpha,\beta)$, i.e. $\langle\overline{u},\overline{n}\rangle>0$, where $\displaystyle\overline{n}=\left(\mu_{\tau},-1\right)$ is an outward pointing normal vector to the line segment $(l_{\tau})$.

Similarly as in (23) and (24) we have:

[TABLE]

Using the fact that, $\displaystyle Q_{\tau}(z)=\frac{z+\tau}{\tau\hat{H}(z)}$, we deduce:

[TABLE]

and therefore, using the fact that $Q_{\tau}(i\omega_{\tau})=\mu_{\tau}$ (as in the proof of Lemma 2.5), we obtain

[TABLE]

Hence, it is easy to see that:

[TABLE]

as $\mu_{\tau}<0$. The following gradient vector is obtained:

[TABLE]

As the scalar term appearing in front of $\overline{n}$ is positive, it follows that the gradient vector is indeed a normal vector to the line $(l_{\tau})$ pointing outward from the stability region $S_{\tau}(\alpha,\beta)$. In conclusion, the transversality condition is now deduced as in the case of the curve segment $(\gamma_{\tau})$, and it follows that a Hopf bifurcation takes place along the line segment $(l_{\tau})$. We also emphasize that the transversality condition may be checked in a similar way as before along each line $(l_{n})$, $n\in\{1,2,...,N-1\}$.

The points c., d. and e. from the statement of the theorem follow easily taking into consideration the intersections between the lines $(l_{0})$, $(l_{\tau})$ and the curve $(\gamma_{\tau})$. ∎

Example 2.1**.**

If a Dirac kernel $h(t)=\delta(t-\tau)$ is considered in (3), we have $\rho(\omega)=1$ and $\theta(\omega)=\omega$ and it can be shown that $S_{\tau_{2}}(\alpha,\beta)\subset S_{\tau_{1}}(\alpha,\beta)$, for any $\tau_{1}<\tau_{2}$. Hence, the intersection of all stability regions in this case will be:

[TABLE]

which is shown in the first panel of Fig. 3 as the shaded triangle.

Example 2.2**.**

If a strong Gamma kernel $h(t)=4\tau^{-2}t\exp(-2t/\tau)$ is considered in system (3), as $\rho(\omega)=4/(4+\omega^{2})$ and $\theta(\omega)=2\arctan(\omega/2)$, the explicit equation of the curve $(\gamma_{\tau})$ is:

[TABLE]

Moreover, as $\mu_{\tau}=-(\tau+2)^{2}/\tau$, the equation of the line $(l_{\tau})$ is

[TABLE]

Therefore, the stability region given by Theorem 2.8 is

[TABLE]

Using basic inequality techniques, it can be shown that $S_{2}(\alpha,\beta)\subset S_{\tau}(\alpha,\beta)$, for any $\tau>0$, where:

[TABLE]

In conclusion, in the case of a strong Gamma kernel, $S_{2}(\alpha,\beta)$ (see third panel of Fig. 3) represents the intersection of all stability regions $S_{\tau}(\alpha,\beta)$, for any $\tau>0$.

The following result follows similarly in the case when equation (18) does not admit any positive roots. In this case, the stability region $S_{\tau}(\alpha,\beta)$ will be unbounded.

Theorem 2.9**.**

If equation (18) does not admit any strictly positive real root (i.e. $N=1$), the boundary of the stability region $S_{\tau}(\alpha,\beta)$ of the equilibrium state $(u^{\star},v^{\star})$ of system (3) is given by the union of the half-line and curve given below:

[TABLE]

At the boundary of the stability region $S_{\tau}(\alpha,\beta)$, the following bifurcation phenomena take place in a neighborhood of the equilibrium $(u^{\star},v^{\star})$ of system (3):

• a.

Saddle-node bifurcations take place along the open half-line $(l_{0})$;

• b.

Hopf bifurcations take place along the curve $(\gamma_{\tau})$;

• c.

Bogdanov-Takens bifurcation at $(\alpha,\beta)=(2,1)$.

Example 2.3**.**

If a weak Gamma kernel $h(t)=\tau^{-1}\exp(-t/\tau)$ is considered in system (3), as $\rho(\omega)=(1+\omega^{2})^{-1/2}$ and $\theta(\omega)=\arctan\omega$, it can be easily seen that equation (18) does not have any positive roots. The explicit equation of the curve $(\gamma_{\tau})$ (which is in fact an arc of a parabola) is:

[TABLE]

Therefore, the stability region given by Theorem 2.9 is

[TABLE]

As $\tau+\dfrac{1}{\tau}\geq 2$, it follows that $S_{1}(\alpha,\beta)\subset S_{\tau}(\alpha,\beta)$, for any $\tau>0$. Hence, in this case, $S_{1}(\alpha,\beta)$ (shown in the second panel of Fig. 3) represents the intersection of all stability regions $S_{\tau}(\alpha,\beta)$, for all $\tau>0$.

Comparable stability characterisations have been obtained for networks of identical neurons in [27], without specifically focusing on the occurrence of bifurcation phenomena. However, in this previous work, the stability region is characterized in terms of the eigenvalues of the network’s connection matrix. While our results from Theorems 2.8 and 2.9 can be easily translated in similar terms, for our two-dimensional Wilson-Cowan model it seems to be more convenient to use $\alpha$ and $\beta$ as characteristic parameters, allowing us to directly discuss the combined effect of the coupling coefficients and the average time delay on the stability regime of the considered model, as described below.

While Theorem 2.1 guarantees that there is a kernel-invariant stability region common to a large class of delay distributions (including all Dirac and Gamma kernels), Figure 3 shows that the size of the whole stability region varies widely around this region, depending on the kernel type. As an illustration of Examples 2.1, 2.2 and 2.3, the panels of Figure 3 show respectively the intersection of stability regions (over all $\tau>0$) found theoretically for the Dirac kernel (dark triangle $S_{\infty}(\alpha,\beta)$ in left panel), the weak Gamma kernel (shaded region $S_{1}(\alpha,\beta)$ in the center panel) and for the strong Gamma kernel (shaded region $S_{2}(\alpha,\beta)$ in the right panel). It is then easy to see that the Gamma kernels support a larger subset of parameters for which the system converges to a stable equilibrium than the Dirac kernel. The weak Gamma kernel in particular produces an unbounded stability region (according to Theorem 2.9). Since the subsequent analysis in the original Wilson-Cowan model was done for a weak Gamma kernel [6], one interesting avenue for future work would be to compare the behaviors between the model with Gamma kernel as it was adapted in the original paper (with coarse-grain approximation), and in its unaltered form (without the coarse-grain approximation).

Our results also reveal that, for each kernel type, there are multiple mechanisms (relying on changes to both the parameters $(\alpha,\beta)$ and to the delay $\tau$) that may be used by the system to move from a regime with a unique stable equilibrium (steady firing activity) to a regime of stable oscillations. A transition based solely on altering the integration parameters $(\alpha,\beta)$ is more easily accomplished for the Dirac kernel than for the strong Gamma kernel, since the latter has a significantly larger stability region. For example, smaller perturbation needed to be applied to $(\alpha,\beta)=(0,0)$ (which is in the kernel-invariant stability region) in order to produce oscillations for the Dirac kernel versus the strong Gamma kernel. Moreover, notice that in both the case of the Dirac and the strong Gamma kernel, the transition in and out of the stability region may be promoted by either an increase in $\beta$, or by a decrease in $\alpha$. In contrast, in the case of the weak Gamma kernel (for which the stability region is unbounded), changing $\beta$ does not alter the stability of the equilibrium (i.e., does not transition the system into and oscillatory regime) if $\alpha$ is negative and sufficiently large in absolute value.

Physiologically, these paths can be tied back to changes in the coupling weights $a,b,c,d$, since $\alpha=af_{1}^{\prime}(\theta_{u}+au^{\star}+bv^{\star})+df_{2}^{\prime}(\theta_{v}+cu^{\star}+dv^{\star})$ and $\beta=(ad-bc)f_{1}^{\prime}(\theta_{u}+au^{\star}+bv^{\star})f_{2}^{\prime}(\theta_{v}+cu^{\star}+dv^{\star})$. While the dependence of $(\alpha,\beta)$ on $(a,b,c,d)$ is too complex to draw any direct interpretations, increasing the cross-population coupling coefficients $b$ and $c$ may be way of lowering $\beta$ values, while self-coupling coefficients $a$ and $d$ seems to promote increases in both $\alpha$ and $\beta$. It can then be speculated that, in the case of Dirac and strong Gamma kernels, transitions in and out of oscillations can be accomplished by changing either self-coupling or inter-coupling strength (to cross either one of the Hopf curves delimiting stable state from oscillatory regimes). More explicit relationships can be worked out between the computational parameters $(\alpha,\beta)$ and the connectivity weights $(a,b,c,d)$ under specific modeling assumptions for the integrating functions $f_{1}$ and $f_{2}$, the type of kernel and distributed delay $\tau$. To illustrate this possibility, we explored numerically a particular theoretical example (in Section 3.1), and a neurophysiology-driven application to dynamic behavior in the basal ganglia (in Section 3.2).

More generally, our theoretical work suggests that different types of delay distributions may be optimal for different dynamic behaviors and transitions between them, and supports the relatively unexplored possibility that the delay kernel may act as a potential specialization mechanism of the cell to perform circuit-specific and task-specific functions. For example, functions which generally require constant firing rates may be easily performed by neural populations specialized in weak Gamma integration of inputs. Functions that require increased flexibility in and out of oscillations may be best accomplished by strong Gamma integrators. Functions characterized by transitions between different and more complex oscillatory patterns may require a Dirac cell specialization (as additionally illustrated by our simulations in the next section). This will be further discussed and contextualized in the Discussion section.

3 Applications and numerical simulations

In our previous numerical work on the model [13], we showed how periodic behaviors, but also quasi-periodic and chaotic patterns are also possible in the $(u,v)$-phase-plane. In this study, we aim to further illustrate how access to these behaviors, as well as the transitions between them, are governed not only by the coupling strengths and by the delay $\tau$, but also by the shape of the delay distribution. To this end, we performed numerical simulations, using our in house codes written in Wolfram Mathematica [39], and made available on github [40].

3.1 Simulations in a theoretical context

In this section, a brief numerical experiment aims to simulate the kinetic details in a series of behaviors that were classified directly by our theoretical results. To facilitate comparison with the original reference, we considered the same activation functions $f_{1}(x)=f_{2}(x)=f(x)=(1+\exp(-\delta x))^{-1}$ and input drives $\theta_{u}=0.1$, $\theta_{v}=0.2$. The coupling parameters in our current simulation are fixed to new values $a=d=-19$, $b=c=10$, and $\delta=10$, corresponding overall to new characteristic parameters $\alpha=-17.8796$ and $\beta=57.7268$.

As $\alpha^{2}>4\beta$ and $\beta>\alpha-1$ from the theoretical analysis presented in the previous section, it follows that a Hopf bifurcation occurs only in the conditions of Theorem 2.8, when $(\alpha,\beta)\in(l_{\tau})$. To find the critical value of the time delay $\tau$ for a particular delay kernel, we proceed as follows, according to Theorem 2.8:

• •

find the roots $\mu$ of the quadratic equation $\beta=\mu(\alpha-\mu)$ in the interval $(-\infty,-1)$;

• •

for each $\mu$ computed at the previous step, numerically compute the root $\omega>0$ of the equation $\rho(\omega)\cos\theta(\omega)=\mu^{-1}$ such that $\theta(\omega)\in(\pi/2,\pi)$, and retain the smallest such value $\omega^{*}$;

• •

the critical value of the average time delay is $\tau^{*}=-\omega^{*}\cot\theta(\omega^{*})$.

It is also important to note that the frequency of oscillations at the Hopf bifurcation which occurs at the critical value $\tau^{*}$ of the average time delay is computed according to the formula $\dfrac{\omega^{*}}{2\pi\tau^{*}}$, based on the fact that, according to our theoretical results, $\pm i\omega^{*}/\tau^{*}$ are the critical roots of the characteristic function $\Delta(z)$.

While for these fixed parameter values the system always has a steady state at $(u^{\star},v^{\star})=(0.0478985,0.0511112)$, the stability of this equilibrium changes when the time delay $\tau$ is increased past a supercritical Hopf bifurcation value $\tau=\tau^{*}$, with the creation of a stable limit cycle (that reflects onset of oscillations in the corresponding cells). The position of $\tau^{*}$, as well as the geometry, duty cycle and evolution of the oscillation born through this bifurcation further depend on the type of kernel used (as shown in Figure 4). For example, based on Theorem 2.8, the critical value is $\tau^{\star}_{0}=0.120766$ for the Dirac kernel, $\tau^{\star}_{2}=0.433992$ for the strong Gamma kernel. The frequencies of oscillations occurring due to the Hopf bifurcation at the critical values of the average time delay are $2.16675$ (for the Dirac kernel) and $0.87829$ (for the strong Gamma kernel), respectively (see Fig. 4). As for the values of $\alpha$ and $\beta$ given above, we have $\alpha^{2}>4\beta$ and $\beta>\alpha-1$, Theorem 2.9 guarantees that, for the specific case of a weak Gamma kernel, a Hopf bifurcation never occurs for any $\tau>0$ (the equilibrium $(u^{\star},v^{\star})$ stays locally asymptotically stable). Hence, no oscillations are expected to occur in a neighborhood of the equilibrium if a weak Gamma kernel is considered in the mathematical model.

This reflects the potential importance of using different delay kernel types when modeling coupled population activity, since different delay distribution profiles may promote different behaviors. For example, our results suggest that the Dirac and strong Gamma kernels are more conducive of oscillations in the mean field firing rates, while the weak Gamma kernel promotes steady firing rates. Moreover, our numerical simulations also revealed cascades of bifurcations between periodic and quasi-periodic oscillations in conjunction with a discrete time delay (see Fig. 5), which could not be observed in the case of strong Gamma kernels with the same system parameters. This suggests that the limiting case of discretely distributed delays, while likely the most delicate in terms of a physiological implementation, may also offer the optimal physiological substrate for the complex aperiodic oscillations associated to task-specific functions in certain brain areas. To further investigate this possibility, we considered a specific application of this theoretical framework to a simplified model of the basal ganglia circuit involved in Parkinson’s Disease.

3.2 Applications to the Parkinson’s Disease

We will now use this theoretical layout to contextualize an existing computational model (based upon anatomical and electrophysiological research). This was developed by Holgado et al. [12] to describe an application to the subthalamic nucleus (STN) - globus pallidus (GPe) network involved in Parkinson’s Disease:

[TABLE]

where $STN$ and $GP$ represent the firing rates of the subthalamic nucleus and globus pallidus, respectively; $Ctx$, $Str$ are cortex, and striatum constant inputs; $\tau_{S}$, $\tau_{G}$ are time constants associated with the $STN$, $GP$ populations; $\Delta t_{MN}$ are trasmission delays from population $M$ to population $N$ and $F_{S}(\cdot)$ and $F_{G}(\cdot)$ are activation functions for each population. The parameter values used in [12] are given in Table 1. The synaptic connection weight parameters ($w_{SS},w_{GS},w_{SG},w_{GG}$) have been estimated both in a healthy state and a parkinsonian/diseased state, and are also given in Table 1. Moreover, for simplicity, equal time constants and time delays are considered (i.e. $\tau_{S}=\tau_{G}=6$ (ms), $\Delta t=\Delta t_{GS}=\Delta t_{SG}=\Delta t_{GG}$ (ms)). In our further investigation of this application, we will also assume equal transmission delays, although our theoretical setup is otherwise more general (for example, it transcends the need to linearize the activation functions, which was a building step in the original analysis in the reference).

We consider a time rescaling $t\mapsto\tau_{S}t$, the rescaled state variables $u(t)=STN(\tau_{S}t)$, $v(t)=GP(\tau_{S}t)$, and the parameters $a=w_{SS}$, $b=-w_{GS}$, $c=w_{SG}$, $d=-w_{GG}$, $\theta_{u}=w_{CS}Ctx$, $\theta_{v}=-w_{XG}Str$. The activation functions are as in [12]:

[TABLE]

The state variables $u$ and $v$ verify a system of the type (3) with an average time delay $\tau=\Delta t/\tau_{S}$ (taking into account that in the original model from [12], discrete time delays are considered).

For the healthy state parameters given in Table 1, the characteristic parameter values are $(\alpha_{H},\beta_{H})=(-3.06805,2.24878)$. As in this case, $\alpha_{H}^{2}>4\beta_{H}$, a Hopf bifurcation may occur only in the conditions of Theorem 2.8, when $(\alpha_{H},\beta_{h})$ belong to the line segment $(l_{\tau})$. In the Dirac kernel case, the critical value of $\tau=\Delta t/\tau_{S}$ responsible for the occurrence of a Hopf bifurcation is computed as described in the previous section: $\tau_{0,H}^{*}=1.367$, and the frequency of oscillations occurring in system (28) near the bifurcation point is $\omega^{*}/(2\pi\tau_{S}\tau^{*})=41.5133~{}Hz$, which agrees with the numerical simulations presented in Fig. 6. On the other hand, in the case of a strong Gamma kernel, it may be easily checked that $(\alpha_{H},\beta_{H})$ belongs to the ”smallest” stability region $S_{2}(\alpha,\beta)$ indicated in Example 2.2, and therefore the equilibrium of system (28) is asymptotically stable for any $\tau$ (see Fig. 7). Likewise, in the case of a weak Gamma kernel, $(\alpha_{H},\beta_{H})$ belongs to the ”smallest” stability region $S_{1}(\alpha,\beta)$ indicated in Example 2.3, and hence the equilibrium of system (28) is asymptotically stable for any $\tau$, as shown in Fig. 8.

For the pathological state parameters given in Table 1, the characteristic parameter values are $(\alpha_{P},\beta_{P})=(-2.53928,11.2213)$. This time, as $\alpha_{P}^{2}<4\beta_{P}$, a Hopf bifurcation occurs only when $(\alpha_{P},\beta_{P})\in(\gamma_{\tau})$. To find the critical value of $\tau$ for a particular delay kernel, we numerically solve for smallest positive $\omega$ and $\tau$ the system which provides the parametric equations of the curve $(\gamma_{\tau})$, as in Theorem 2.8 or Theorem 2.9, and we obtain the following values for the considered kernels:

• •

Dirac: $\tau_{0,P}^{*}=0.216411$ and corresponding frequency of oscillations in (28): $84.8049~{}Hz$;

• •

weak Gamma: $\tau_{1,P}^{*}=0.619418$ and corresponding frequency of oscillations in (28): $50.7756~{}Hz$;

• •

strong Gamma: $\tau_{2,P}^{*}=0.283222$ and corresponding frequency of oscillations in (28): $72.5652~{}Hz$;

Again, the frequencies have been computed according to the formula $\omega^{*}/(2\pi\tau_{S}\tau^{*})$, due to the time rescaling which has been used. The theoretically predicted frequencies correspond to the numerically computed values, as seen in Figs. 6, 7 and 8.

We simulated the behaviour of the system (28), with the above chosen parameters and activation functions. The behaviors found illustrate in a clinical context our theoretical results obtained in this paper (e.g. Theorems 2.8 and 2.9) and expand on the modeling results in the reference [12].

As a first step, Figures 6, 7 and 8 illustrate the behavior of each of the two state variables $u$ and $v$ for different delay kernels: for discrete delays (Figure 6) for strong Gamma delays (Figure 7) and for weak Gamma delays (Figure 8). In each figure, the left panels represent the behavior of the system with coupling parameters within a healthy functioning regime, while the right panels illustrate the behavior for coupling parameters within the PD range (as shown in Table 1). The different rows of the figures illustrate different values of the delay $\Delta t/\tau_{S}\in[0,2]$, interval compatible with the physiological range (as per the table in the reference [12]). Notice that, in the PD case (right panels), sustained oscillations in the Beta range can be readily triggered by delay values within this range ($\Delta t/\tau_{S}\sim 0.21$ in the case of the discrete kernel, $\Delta t/\tau_{S}\sim 0.28$ in the case of the strong Gamma kernel, $\Delta t/\tau_{S}\sim 0.61$ in the case of the weak Gamma kernel). Interestingly, while the oscillation amplitudes increase with the delay (top to bottom panels), the frequency is higher for shorter delays.

In contrast, note that much higher delay values are required in order to trigger stable oscillations in the healthy regime (they appear at a critical value of $\Delta t/\tau_{S}\sim 1.36$ for the discrete kernel, and they never occur in the case of weak or strong Gamma kernel). This is in agreement with the idea that periodic oscillations in the Beta range are greatly enhanced in PD, as proposed in the original reference [12]. Empirical studies support the link between increased Beta oscillatory activity and PD symptoms (such as bradykinesia and rigidity) [41], as discussed in more detail in the last section.

The differential propensity of the system for oscillations in the PD versus the healthy regime is further illustrated in Figures 9, 11, 11, 13 and 13. These deliver a more comprehensive description of the subtle interplay (between coupling strengths, delay kernel and delay value) that orchestrates the triggering and maintenance of oscillatory behavior.

Figure 9 illustrates differences between the onset of oscillations in the PD system versus the healthy system, in the case of discrete delays. For each panel, the coupling parameters $w_{SS}$, $w_{GG}$, $w_{CS}$ and $w_{XG}$ were fixed (not shown): to the healthy state values in Table 1 (in the left panel) and to the PD values (in the right panel). Then, for each pair of cross-connectivity parameters $w_{SG}$ and $w_{GS}$ in the appropriate range, the critical value (i.e., the value corresponding to onset of oscillations) of $\Delta t/\tau_{S}$ was computed. The corresponding value was then plotted at each parameter point $(w_{SG},w_{GS})$ as an associated color in a blue to red color map (color bar shown on the right of each parameter plot).

An immediate comparison between the two figure panels reiterates (in this broader parameter frame) the idea that, for Dirac distributed delays, sustained oscillations are accessible in the PD coupling range at significantly shorter delay times (consistent with the physiological values) than in the healthy coupling range. A more in depth analysis of the color patterns may suggest candidate mechanisms, based on synaptic remodeling, for onset and cessation of oscillatory behavior in each case. For example, in healthy systems, a weaker GP to STN coupling is more likely to result in sustained oscillations (since it lowers the critical delay to values within the physiological range). In PD systems, operating with short delay times ($\Delta t/\tau_{S}=0.2-0.3$), strengthening of the STN to GP synaptic coupling may readily trigger sustained oscillations in PD patients, while weakening it can lead to ceasing the oscillatory behavior. This suggests the importance of having access to information not only on the synaptic coupling profile, but also on knowledge of the type of delay range and distribution involved in the neural integration, in the circuit in question.

Continuing along with this thought, the rest of this section revisits, in the context of this application, the importance of the subtle, but crucial differences in dynamics induced by different profiles in the delay distribution. Figures 11 and  11 illustrate and compare the same exact coupling ranges as in Figure 9 (for the healthy state on the left, and for the PD state on the right), for the case of the strong Gamma kernel. Notice that, in the pathological case, the critical values of $\Delta t/\tau_{S}$ are sightly higher throughout the parameter plane than the similar values in the case of the Dirac kernel. However, the candidate mechanisms proposed for the trigger and stop of oscillations in the case of the Dirac kernel may equally apply.

On the other hand, by contrast with the Dirac case, oscillations are no longer possible in the healthy range of coupling parameters in the case of the strong Gamma kernel. Theoretically, this is because the healthy coupling range (as shown in the left hand panel of Figure 9 and re-expressed as the red subset in the $(\alpha,\beta)$ domain, in Figure 11) is fully contained in this case in the stability region of the system (shown as a grey shaded area in the $(\alpha,\beta)$ domain in Figure 11), as computed in Theorem 2.7. Thus, previously shown in Figure 6a, any tendency towards oscillations is quickly suppressed, and the system converges to a steady firing regime.

In line with our previous analysis, stable oscillations are gated by even longer delay values when using weak Gamma distributed delays. Indeed, Figures 8, 13 and 13 illustrate the behavior of the system using a weak Gamma kernel, with the same fixed parameters as for the Dirac and strong Gamma kernels. The healthy coupling range remains a subset of the stability domain (as shown in Figure 13), hence stable oscillations are not possible as a long-term outcome in this scenario (Figure 8a). For the pathological coupling range, oscillations are possible for values of $\Delta t/\tau_{S}$ larger than $\sim 0.5$, but only if both the STN-GP and GP-STN synaptic pathways are strong enough (Figure 13 shows that oscillations are possible in this case only in the upper right corner of the $(w_{SG},w_{GS})$ parameter region).

A further comparison between oscillations in the case of Dirac, weak and strong Gamma distributed delays is shown in Figure 14. For the PD coupling range, each row of panels captures the emergence of oscillations when increasing the delay $\Delta t/\tau_{S}$ past the critical value: for the Dirac kernel (top row), for the weak Gamma kernel (middle row), and for the strong Gamma kernel (bottom row). As previously discussed, the critical delay for the Dirac distribution ($\Delta t/\tau_{S}\sim 0.21$) is slightly lower than that for the strong Gamma distribution ($\Delta t/\tau_{S}\sim 0.28$). In addition, the figure confirms that the critical value for the weak Gamma distribution is significantly larger than both ($\Delta t/\tau_{S}\sim 0.61$), as one would have expected from the theoretical results illustrated in Figure 3. The significance of these differences will be further discussed in the next section.

4 Discussion

We have accomplished a local stability and bifurcation analysis of a generalized version of the Wilson-Cowan model of excitatory and inhibitory interactions in localized neuronal populations, incorporating general distributed delays. Essential differences have been pointed out for different scenarios involving diverse delay kernels, (ranging from simple periodic to quasi-periodic and aperiodic behaviors). This emphasizes the potentially crucial importance of the kernel choice to the dynamic behavior of the system.

Differentiating between specific behaviors is very important, since they represent different physiological rhythms observed empirically in neuronal E/I circuits.

Recall that the system variables $u(t)$ and $v(t)$ represent meanfield firing activities in the excitatory and respectively inhibitory Wilson Cowan coupled populations. Then a stable equilibrium corresponds to convergence of the two populations’ firing to almost constant rates, while a stable cycle represents meanfield oscillations between higher and lower firing rates in the two populations.

Oscillatory behavior is ubiquitous in the brain, as a way for neurons to efficiently assemble in a synchronized fashion, optimal for receiving and sending information [42]. Oscillations are typically viewed as a result of balancing excitation and inhibition, with either gaining ground at different times within an oscillation cycle.

Brain rhythms cover frequency bands from lower than 1Hz to 600 Hz, with different rhythms associated with different brain areas, functions and states. Moreover, this correspondence is not one-to-one, in the sense that two different brain circuits may generate oscillations in the same frequency band (likely by using different underlying mechanisms). In turn, the function of a particular oscillatory pattern largely depends on the underlying brain structure. For example, Beta rhythms (13-30 Hz) in the motor system are associated with the absence of movement.

Rhythms alterations in a brain circuit outside of its healthy range have been found by empirical studies to be associated with mental or neurological disorders (for example, exacerbated synchronization is the mark of epilepsy, Gamma oscillations have been found to be diminished in schizophrenia, and Parkinson’s Disease has been tied to enhanced Beta oscillations in motor areas).

However, the relationship between pathology and oscillations is rather complex, going beyond simply locating increased amplitude oscillations within a specific frequency band. Transition paths between healthy and pathological rhythms are still under investigation for most disorders, and recent research has been striving to define and understand them. Along these lines, our bifurcation analyses may reveal potential mechanisms that govern transitions between oscillatory behaviors in PD.

With this in mind, our results show that a weak Gamma kernel may be a preferable mode of input integration in populations whose function is based primarily on sustaining an approximately constant rate of firing. A strong Gamma kernel may be used as an input integration profile when the basic function of the circuit requires transitions into simple, regular oscillations. Finally, the Dirac kernel (which can mathematically be viewed as a subtler, limit-like profile for distributed delays) may be the signature of finer oscillators, whose function is modulated by transitions between more complex firing patterns. For each type of kernel, our results suggest that more specific behaviors, and transitions between dynamic regimes can be further controlled by the connectivity strengths, and by the delay $\tau$.

The importance of the type and length of delay to functional brain rhythms, supported by our theoretical results and numerical simulations, was further illustrated in an application to a Wislon-Cowan type model with distributed delays of the basal ganglia, inspired by the work of Holgado et. al.[12], investigating control mechanisms of Beta rhythms in Parkinson’s Disease.

In the context of normal executive control, Beta oscillations have been associated with promoting existing movement, at the cost of initiating a new motor set [41]. In this context, it does not seem surprising that oscillations in the basal ganglia are significantly more pronounced in PD. Deep brain stimulation studies suggest in fact a nontrivial relationship, in which beta oscillatory behavior contributes both quantitatively (higher amplitudes) but also causally to the motor impairment in PD. Hence understanding the timing and conditions that favor onset of Beta activity is crucial, since understanding this causal relationship may hold potential for therapeutic interventions.

The Holgado model places these questions in the framework of a dynamical system formed of two coupled state variables, representing activation of the STN and GP, as they respond to outside stimulation, and to modulations from other brain areas (the cortex and the striatum). Based on empirical evidence, the authors defined in the reference different ranges of pairwise synaptic coupling strengths between these brain areas, one characteristic to healthy brain functioning, the other corresponding to the dynamic patterns found in Parkinson’s Disease. They found a prevalence for stable oscillations in the PD regime of the model, corresponding to physiological rhythms in the Beta-range in the corresponding brain areas of PD patients. We advanced the analysis of this model by investigating the importance of the coupling range to the system’s behavior, in conjunction with the size and the distribution of the delays. We found that, for all kernels, stable oscillations are accessible to the system in the PD coupling range at shorter delays than to the system in the healthy range (exactly how short, depends on the type of kernel used).

Since the shorter delays are compatible with physiological values, lower critical delays can be seen as a gate to promoting the Beta-range oscillations in PD patients. Comparably short delay values are significantly more likely to place a system with healthy synaptic coupling within the stability domain (as shown in Figure 15a). This corresponds to what is known about healthy basal ganglia function, which promotes desynchronization of Beta oscillations during executive responses. In contrast, in a system within the PD functioning range, Beta oscillations persist even as the execution is being carried out. This may contribute to the motor abnormalities found in PD (such as bradykinesia and rigidity, which have been empirically associated with enhanced Beta rhythms). Altogether, our analysis of Holgado’s model of basal ganglia suggests that it may be important to investigate neural mechanisms that have the potential to reshape a cell’s memory profile when processing inputs, as a possible means to mitigate executive function abnormalities in PD patients.

We additionally observed particular patterns for the critical delay (onset of oscillations) as well as for entering the Beta range activity, in the case of each kernel (Fingure 15b). We considered a system within the PD coupling regime, and with a distribution of delays within the physiological range, this combination placing it on a dynamic path of perpetual oscillations, even in absence of any additional stimulation. We explored the possibility that cessation of oscillations may be obtained by enhancing a target synaptic pathway. We noticed that, depending on the context, the perturbation may not necessarily have to be as large as to place the system in the healthy coupling regime. The adjustment would need to be just significant enough to move the system to a coupling profile where a shorter time delay is required in order to maintain oscillatory behavior; hence, with the current time delay, oscillations will be suppressed. Our future work aims to further study the potential physiological implementation of such mathematical solutions, and their potential towards development of clinical treatments.

A significant limitation of the application in its current form is that brought by using the original, two-dimensional Wilson-Cowan coupling scheme. A two-dimensional system is a great place to start, since it is substantially easier to analyze and understand than higher dimensional, more complex systems. However, the neuronal circuit that governs executive function (and that has traditionally been studied in conjunction with conditions associated with executive disturbances, like Parkinson’s Disease, or Obsessive Compulsive Disorder) crucially encompasses additional brain areas and connections, not represented in the Holgado model. For example, PD has been linked in both experimental and computational literature to modifications in the striatal-GP pathway [43, 44], suggesting that the striatum should be included and studied as an additional state variable in the coupled system.

This underlines the general importance of extending the analysis of the Wilson-Cowan model with distributed delays to higher dimensional networks. We expect theoretical approaches to become more difficult, being affected by the fact that dynamic complexity increases vastly past dimension two, by the increased size of the coupling parameter space, and by the additional contributions of the network architecture to the dynamics). However, the study of a more complete version of the importance of delays in the circuit controlling executive function could complement in novel and crucial ways existing work on modeling this circuit in absence of delays.