Stability and bifurcations of two-dimensional systems with distributed delay and applications to a Wilson-Cowan model

# Stability and bifurcations of two-dimensional systems with distributed delay and applications to a Wilson-Cowan model

Eva Kaslik    Emanuel-Attila Kokovics    Anca Rădulescu
###### Abstract

A generalization of the well-known Wilson-Cowan model of excitatory and inhibitory interactions in localized neuronal populations is presented, by taking into consideration distributed time delays. A stability and bifurcation analysis is undertaken for the generalized model, with respect to two characteristic parameters of the system. The stability region in the characteristic parameter plane is determined and a comparison is given for several types of delay kernels. It is shown that if a weak Gamma delay kernel is considered, as in the original Wilson-Cowan model without time-coarse graining, the resulting stability domain is unbounded, while in the case of a discrete time-delay, the stability domain is bounded. This fact reveals an essential difference between the two scenarios, reflecting the importance of a careful choice of delay kernels in the mathematical model. Numerical simulations are presented to substantiate the theoretical results. Important differences are also highlighted by comparing the generalized model with the original Wilson-Cowan model without time delays.

## 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 cell. The Hogkin-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 popular option has been using state variables to characterize the meanfield spiking activity in a population of cells. 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.

Among meanfield models, the Wilson-Cowan model is perhaps the most popular. The model, derived in 1972 by Wilson and Cowan [5], describes the localized interactions in a pair of excitatory and inhibitory neuronal populations. At each time instant , the proportions of excitatory and inhibitory cells firing per unit of time are captured by the two state variables and . The original model considers the effect that an external input has on the E/I system, based not only on the coupling stengths between the two units, but also on the history of firing in each. More specifically, and are considered to be equal to the proportion of cells which are sensitive (i.e. not refractory) and which also receive at least threshold excitation at the moment of time . This leads to the following system of integral equations:

 ⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩E(t+τ)=(1−∫tt−rE(s)ds)⋅Se[∫t−∞h(t−s)(c1E(s)−c2I(s)+Pe(s))ds]I(t+τ′)=(1−∫tt−r′I(s)ds)⋅Si[∫t−∞h(t−s)(c3E(s)−c4I(s)+Pi(s))ds] (1)

In this system, the first factors in the right hand side represent the proportion of sensitive excitatory / inhibitory cells, where is the absolute refractory period (msc), the functions , are sigmoid threshold functions, their arguments denoting the mean field level of excitation / inhibition generated in an excitatory /inhibitory cell at time (assuming that individual cells sum their inputs and that the effect of the stimulation decays exponentially with a time course ). Moreover, are connectivity coefficients representing the average number of excitatory / inhibitory synapses per cell and denote external inputs.

The model historically known as the Wilson-Cowan model [5] was then obtained at this point applying time-coarse graining. This final form, consisting of a system of ordinary differential equations without time delays, is very convenient and was extensively analyzed and used in the modeling literature. However, the coarse-graining procedure obscures potentially vital timing information related to the readiness of a cell in either population to generate a spike. If this information is crucial to the neural function that the model is aiming to address, it would be useful to return to the original equations, prior to coarse graining. In this context, generalizations of the standard model, including discrete time-delays, have been investigated in several papers, often considering refractory periods being equal to zero.

### 1.2 Our model

Based on the integral terms appearing in the original model (1) as arguments of the threshold functions, the following model with distributed delays will be analyzed in this paper:

 ⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩˙u(t)=−u(t)+f[θu+∫t−∞h(t−s)(au(s)+bv(s))ds]˙v(t)=−v(t)+f[θv+∫t−∞h(t−s)(cu(s)+dv(s))ds] (2)

where and represent the synaptic activities of the two neuronal populations, are connection weights and , are background drives. The activation function is considered to be increasing and of class on the real line.

In system (2), the delay kernel is a probability density function representing the probability that a particular time delay occurs. It is assumed to be bounded, piecewise continuous and satisfy

 ∫∞0h(s)ds=1,with the average time delayτ=∫∞0sh(s)ds<∞.

The particular case of discrete time delays (Dirac kernels) has been discussed in [6]. However, there are other important classes of delay kernels often used in the literature, such as Gamma kernels or uniform distribution kernels. It is worth mentioning that in the original Wilson-Cowan model [5], a weak Gamma kernel has been used, so this case should be the original reference point. Analyzing mathematical models with particular classes of delay kernels (e.g. weak Gamma kernel or strong Gamma kernel ) may shed a light on how distributed delays affect the dynamics differently from discrete delays. However, in the modeling of real world phenomena, one usually does not have access to the exact distribution, and approaches using general kernels may be more appropriate [7, 8, 9, 10, 11, 12, 13, 14, 15].

Initial conditions associated with system (2) are of the form:

 u(s)=φ(s),v(s)=ψ(s),∀s∈(−∞,0],

where are bounded real-valued continuous functions defined on .

## 2 Main stability and bifurcation results

The equilibrium states of system (2) are the solutions of the following algebraic system:

 {u=f(θu+au+bv)v=f(θv+cu+dv) (3)

The linearized system at an equilibrium state is

 ⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩˙u=−u+ϕ1∫t−∞h(t−s)(au(s)+bv(s))ds˙v=−v+ϕ2∫t−∞h(t−s)(cu(s)+dv(s))ds (4)

where and .

Applying the Laplace transform to the linearized system (4), we obtain:

 {zU(z)−u(0)=−U(z)+ϕ1H(z)(aU(z)+bV(z))zV(z)−v(0)=−V(z)+ϕ2H(z)(cU(z)+dV(z)) (5)

where and represent the Laplace transforms of the state variables and respectively, while is the Laplace transform of the delay kernel .

System (5) is equivalent to:

 (z+1−aϕ1H(z)−bϕ1H(z)−cϕ2H(z)z+1−dϕ2H(z))(U(z)V(z))=(u(0)v(0)) (6)

and hence, the characteristic equation associated to the equilibrium state is

 Δ(z)=(z+1)2−αH(z)(z+1)+βH2(z)=0 (7)

where

### 2.1 Delay-independent stability and instability results

The following delay-independent stability and instability results are easily obtained, based on the properties of the Laplace transform and the particularities of the characteristic equation (7):

###### Theorem 2.1 (Delay-independent stability and instability).

1. In the non-delayed case, the equilibrium state of system (2) is locally asymptotically stable if and only if

 α
2. If the following inequality holds:

 |α|+|β|<1 (9)

then the equilibrium state of system (2) is locally asymptotically stable for any delay kernel .

3. If the following inequality holds:

 β<α−1 (10)

then the equilibrium state of system (2) is unstable for any delay kernel .

###### Proof.

1. In the non-delayed case (, for any ), the characteristic equation (7) becomes:

 Δ(z)=z2+(2−α)z+β−α+1=0

The necessary and sufficient condition (8) for the asymptotic stability of the equilibrium state of system (2) follows from the Routh-Hurwitz stability test.

2. If we assume that the characteristic equation (7) has a root in the right half-plane (), it follows that

 |H(z)|=∣∣∣∫∞0e−zth(t)dt∣∣∣≤∫∞0|e−zt|h(t)dt=∫∞0e−R(z)th(t)dt≤∫∞0h(t)dt=1,

and hence, from (7) we deduce:

 |z+1|2 =|αH(z)(z+1)−βH2(z)| ≤|α||z+1|+|β|

Considering the polynomial , from inequality (9) it follows that and , and hence for any . From the above inequality, we have , and hence, we deduce , which is absurd, since .

Therefore, all the roots of the characteristic equation (7) are in the left half-plane and the equilibrium state of system (2) is asymptotically stable, regardless of the delay kernel .

3. Condition (10) is equivalent to . We also have that as , and therefore, the characteristic equation (7) has at least one positive real root. Hence, the equilibrium state of system (2) is unstable, regardless of the delay kernel . ∎

A saddle-node bifurcation takes place at the equilibrium state of system (2), regardless of the delay kernel , if and only if and

 β=α−1 (11)
###### Proof.

Condition (11) is equivalent to (i.e. the characteristic equation has a zero root). Moreover, is a simple root of the characteristic equation if and only if .

Let us denote by the root of the characteristic equation (7) which satisfies . Taking the derivative with respect to in equation (7), we obtain:

 2(z+1)dzdβ−αH′(z)(z+1)dzdβ−αH(z)dzdβ+2βH′(z)H(z)dzdβ+H2(z)=0,

and hence:

 dzdβ=−H2(z)2(z+1)−αH′(z)(z+1)−αH(z)+2βH′(z)H(z).

As and , it follows that:

 dzdβ∣∣∣β=α−1=1(α−2)(τ+1)≠0.

This completes the proof. ∎

### 2.3 Hopf bifurcation

In the following, we will show that the average time delay of the delay kernel plays an important role in the Hopf bifurcation analysis.

Let , for any . The function is a probability density function with the mean value:

 ∫∞0t^h(t)dt=∫∞0tτh(τt)dt=τ−1∫∞0uh(u)du=τ−1τ=1.

The Laplace transform of is

 ^H(z)=∫∞0e−zt^h(t)dt=∫∞0e−ztτh(τt)dt=∫∞0e−zuτh(u)du=H(zτ)

We assume from now on that does not depend on the average time delay . In fact, for the most important classes of delay kernels we have:

• Dirac kernel: ;

• Gamma kernel: ;

We write in the polar form as:

 ^H(iω)=ρ(ω)e−iθ(ω)

with , , , , for any .

With the change of variable , the characteristic equation (7) becomes

 ^Δ(z)=τ2Δ(zτ)=(z+τ)2−ατ^H(z)(z+τ)+βτ2^H2(z)=0 (12)

Denoting , it follows that the characteristic equation is equivalent to

 Q2τ(z)−αQτ(z)+β=0 (13)

The characteristic equation (7) has a pair of complex conjugated roots on the imaginary axis if and only if there exists such that is a root of the polynomial .

Case 1: .

In this case, the polynomial has complex conjugated roots and , and hence, Hopf bifurcations may only take place at the equilibrium state of system (2) along the curve from the -plane, defined by the following parametric equations:

 (γτ):⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩α=ατ(ω)=2R(Qτ(iω))=2ρ(ω)[cosθ(ω)−ωτsinθ(ω)]β=βτ(ω)=|Qτ(iω)|2=1ρ(ω)2(1+ω2τ2), ω>0. (14)

Indeed, we prove the following:

###### Lemma 2.3.

If then the characteristic equation (7) has a pair of complex conjugated roots on the imaginary axis if and only if belong to the curve defined by (14).

###### Proof.

The characteristic equation (7) has a pair of complex conjugated roots on the imaginary axis if and only if there exists such that the polynomial has complex conjugated roots and . As

 Qτ(iω) =iω+ττ^H(iω)=iω+ττρ(ω)e−iθ(ω)=(iω+τ)(cosθ(ω)+isinθ(ω))τρ(ω)= =τcosθ(ω)−ωsinθ(ω)+i(ωcosθ(ω)+τsinθ(ω))τρ(ω)

from Viete’s formulas, we obtain:

 α =Qτ(iω)+¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯Qτ(iω)=2R(Qτ(iω))=2⋅τcosθ(ω)−ωsinθ(ω)τρ(ω)=2ρ(ω)[cosθ(ω)−ωτsinθ(ω)]; β =Qτ(iω)¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯Qτ(iω)=|Qτ(iω)|2=∣∣∣iω+ττρ(ω)e−iθ(ω)∣∣∣2=ω2+τ2τ2ρ(ω)2=1ρ(ω)2(1+ω2τ2).

Therefore, we obtain the equivalence from the statement. ∎

It is worth noting that the curve intersects the saddle-node bifurcation line at the point of coordinates , corresponding to .

Case 2: .

In this case, the following result holds:

###### Lemma 2.4.

If then the characteristic equation (7) has a pair of complex conjugated roots on the imaginary axis if and only if is a root of the equation

 ωcosθ(ω)+τsinθ(ω)=0 (15)

and belong to a line

 (l):β=αρ(ω)cosθ(ω)−1(ρ(ω)cosθ(ω))2, (16)
###### Proof.

The characteristic equation (7) has a pair of complex conjugated roots on the imaginary axis if and only if there exists such is a real root of the polynomial . Therefore, and hence, it follows that is a root of the equation

 ωcosθ(ω)+τsinθ(ω)=0.

Assuming that such a root exists, let us denote it by . Hence:

 Qτ(iω∗) =R(Qτ(iω∗))=τcosθ(ω∗)−ω∗sinθ(ω∗)τρ(ω∗)=τcosθ(ω∗)+(ω∗)2τcosθ(ω∗)τρ(ω∗)

As is a root of , we obtain:

 β=αQτ(iω⋆)−Q2τ(iω⋆)=αρ(ω∗)cosθ(ω∗)−1(ρ(ω∗)cos2θ(ω∗))2.

Therefore, the proof is complete. ∎

The existence of the roots of the equation (15), and hence, of the lines given by the previous Lemma, is a particularity of the delay kernel that is considered, and will be discussed separately.

###### Theorem 2.5.

Assuming that the equation (15) has at least one positive real root, let us denote:

 ωτ=min{ω>0 : τsinθ(ω)+ωcosθ(ω)=0}andμτ=(ρ(ωτ)cosθ(ωτ))−1.

The boundary of the stability region of the equilibrium state of system (2) is given by the union of the line segments and curve given below:

 (l0): β=α−1 , α∈[1+μτ,2]; (lτ): β=μτ(α−μτ) , α∈[2μτ,1+μτ]; (γτ): ⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩α=2ρ(ω)[cosθ(ω)−ωτsinθ(ω)]β=1ρ(ω)2(1+ω2τ2) , ω∈(0,ωτ).

At the boundary of the stability domain , the following bifurcation phenomena take place in a neighborhood of the equilibrium of system (2):

• Saddle-node bifurcations take place along the open line segment ;

• Hopf bifurcations take place along the open line segment and curve ;

• Bogdanov-Takens bifurcation at ;

• Double-Hopf bifurcation at ;

• Zero-Hopf bifurcation .

###### Proof.

It is easy to see that for , the characteristic equation has negative real roots, and therefore, the equilibrium state of system (2) is asymptotically stable. As the roots of the characteristic function (or ) continuously depend on the parameters and , the number of roots such that may change if and only if a root or a pair of pure imaginary roots appear, i.e. along the curve or the lines or in the -plane. A simple analysis shows that the line segments and curve segment given in the statement above enclose a connected region of the -plane containing the origin, i.e. the stability region of the equilibrium of system (2). Moreover, Theorem 2.2 shows that a saddle-node bifurcation takes place along the line segment , and hence, statement a. holds.

b. We will first show that Hopf bifurcations take place along the curve segment , with .

Lemma 2.3 provides that the characteristic equation (13) has a pair of pure imaginary roots if belong to the curve . Let us consider an arbitrary and denote by the root of the characteristic equation (13) satisfying where . Our aim is to prove the transversality condition:

 ∇¯¯¯uR(z)(α∗,β∗)>0,

for any outward pointing vector from the region , i.e. , where denotes the outward pointing normal vector to the curve .

From (13) it follows that is a solution of the equation

 Q2τ(z)−αQτ(z)+β=0. (17)

Taking the derivative with respect to , we obtain:

 2Qτ(z)Q′τ(z)∂z∂α−Qτ(z)−αQ′τ(z)∂z∂α=0,

and we express:

 ∂z∂α=Qτ(z)Q′τ(z)[2Qτ(z)−α]. (18)

Therefore, based on the parametric equations of the curve given by (14), we deduce:

 ∂z∂α∣∣∣(α∗,β∗)=Qτ(iω)Q′τ(iω)[2Qτ(iω)−α∗]=Qτ(iω)Q′τ(iω)[2Qτ(iω)−2R(Qτ(iω))]=Qτ(iω)2Q′τ(iω)⋅i⋅I(Qτ(iω)).

Applying the real part, we finally obtain:

 ∂R(z)∂α∣∣∣(α∗,β∗)=R[−i2⋅Qτ(iω)Q′τ(iω)I(Qτ(iω))]=12⋅I[Qτ(iω)Q′τ(iω)I(Qτ(iω))]=12I(Qτ(iω))⋅I[Qτ(iω)Q′τ(iω)].

Taking the derivative with respect to in equation (17), we obtain:

 2Qτ(z)Q′τ(z)∂z∂β−αQ′τ(z)∂z∂β+1=0,

and we express:

 ∂z∂β=1Q′τ(z)[α−2Qτ(z)]. (19)

Again, using the expressions of the parametric curve given by (14), we deduce:

 ∂z∂β∣∣∣(α∗,β∗)=12Q′τ(iω)[−iI(Qτ(iω))]=i2Q′τ(iω)I(Qτ(iω)).

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

 ∂R(z)∂β∣∣∣(α∗,β∗)=−12⋅I[1Q′τ(iω)I(Qτ(iω))]=−12I(Qτ(iω))⋅I[1Q′τ(iω)]

and thus, we obtain the gradient vector:

 ∇R(z)(α∗,β∗) =(∂R(z)∂α,∂R(z)∂β)∣∣∣(α∗,β∗)=12I(Qτ(iω))(I[Qτ(iω)Q′τ(iω)],−I[1Q′τ(iω)])= =12I(Qτ(iω))|Q′τ(iω)|2(I[Qτ(iω)¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯Q′τ(iω)],I(Q′τ(iω))).

The tangent vector to the curve at the point is , where

 α′(ω) =ddω(2R(Qτ(iω)))=2R[ddωQτ(iω)]=2R(i⋅Q′τ(iω))=−2I(Q′τ(iω)) β′(ω) =i⋅Q′τ(iω)¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯Qτ(iω)−i⋅Qτ(iω)¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯Q′τ(iω)=i⋅(Q′τ(iω)¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯Qτ(iω)−Qτ(iω)¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯Q′τ(iω))= =−2I[Q′τ(iω)¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯Qτ(iω)]=2I[Qτ(iω)¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯Q′τ(iω)].

Thus, we obtain the following tangent vector to the curve :

 (α′(ω),β′(ω))=2(−I(Q′τ(iω)),I[Qτ(iω)¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯Q′τ(iω)]).

Therefore, fixing the orientation of the curve in the direction of increasing , a right-pointing normal vector to the curve is:

 ¯¯¯n(ω)=(I[Qτ(iω)¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯Q′τ(iω)],I(Q′τ(iω))).

Hence, the directional derivative is:

 ∇¯¯¯uR(z)(α∗,β∗)=⟨∇R(z)(α∗,β∗),¯¯¯u⟩=⟨¯¯¯n2I(Qτ(iω))|Q′τ(iω)|2,¯¯¯u⟩=⟨¯¯¯n,¯¯¯u⟩2I(Qτ(iω))|Q′τ(iω)|2>0,

as and for any . Therefore, the transversality condition holds, and it follows that a Hopf bifurcation takes place along the curve segment which bounds the stability region .

In the following, we will show that Hopf bifurcations take place along the line segment , with . If belong to this line segment, Lemma 2.4 shows that the characteristic equation (13) has a pair of pure imaginary roots . Let us denote by the root of the characteristic equation (13) with the property where is arbitrarily fixed and . We will prove the transversality condition

 ∇¯¯¯uR(z)(α∗,β∗)>0,

for any vector pointing outward from the region , i.e. , where is an outward pointing normal vector to the line segment .

Similarly as in (18) and (19) we have:

 ∂z∂α=Qτ(z)Q′τ(z)[2Qτ(z)−α]and∂z∂β=1Q′τ(z)[α−2Qτ(z)].

Using the fact that, , we deduce:

 Q′τ(z)=^H(z)−(z+τ)^H′(z)τ^H2(z)=1τ^H(z)−z+ττ^H(z)⋅^H′(z)^H(z)=1−τQτ(z)^H′(z)τ^H(z)

and therefore, using the fact that (as in the proof of Lemma 2.4), we obtain

 Q′τ(iωτ)=1−τμτe−iθ(ωτ)[−iρ′(ωτ)−ρ(ωτ)θ′(ωτ)]τρ(ωτ)e−iθ(ωτ)=eiθ(ωτ)+τμτ[iρ′(ωτ)+ρ(ωτ)θ′(ωτ)]τρ(ωτ).

Hence, it is easy to see that:

 R[Q′τ(iωτ)]=cosθ(ωτ)+τμτρ(ωτ)θ′(ωτ)τρ(ωτ)=μτ(τ−1cos2θ(ωτ)+θ′(ωτ))<0,

as . The following gradient vector is obtained:

 ∇R(z)(α∗,β∗) =R[(Qτ(iωτ)Q′τ(iωτ)[2Qτ(iωτ)−α∗],1Q′τ(iωτ)[α∗−2Qτ(iωτ)])]= =1(2μτ−α∗)⋅R(1Q′τ(iωτ))⋅(μτ,−1)= =R(Q′τ(iωτ))(2μτ−α∗)|Q′τ(iω)|2⋅¯¯¯n.

As the scalar term appearing in front of is positive, it follows that the gradient vector is indeed a normal vector to the line pointing outward from the stability region . In conclusion, the transversality condition is now deduced as in the case of the curve segment , and it follows that a Hopf bifurcation takes place along the line segment .

The points c., d. and e. from the statement of the theorem follow easily taking into consideration the intersections between the lines , and the curve . ∎

The following result follows similarly in the case when equation (15) does not admit any positive roots. In this case, the stability domain will be unbounded.

###### Theorem 2.6.

If the equation (15) does not admit any positive real root, the boundary of the stability region of the equilibrium state of system (2) is given by the union of the half-line and curve given below:

 (l0): β=α−1 , α∈(−∞,2]; (γτ): ⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩α