Stability and bifurcations of twodimensional systems with distributed delay and applications to a WilsonCowan model
Abstract
A generalization of the wellknown WilsonCowan 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 WilsonCowan model without timecoarse graining, the resulting stability domain is unbounded, while in the case of a discrete timedelay, 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 WilsonCowan model without time delays.
1 Introduction
1.1 Modeling background
Computational modeling of neuronal behavior covers a large range of spatiotemporal 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 HogkinHuxley model was in its original form limited to the two voltagedependent 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 widespread possibility has been creating reduced models as modifications of HodgkinHuxley 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 WilsonCowan 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:
(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 WilsonCowan model [5] was then obtained at this point applying timecoarse 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 coarsegraining 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 timedelays, 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:
(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
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 WilsonCowan 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:
where are bounded realvalued continuous functions defined on .
2 Main stability and bifurcation results
The equilibrium states of system (2) are the solutions of the following algebraic system:
(3) 
The linearized system at an equilibrium state is
(4) 
where and .
Applying the Laplace transform to the linearized system (4), we obtain:
(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:
(6) 
and hence, the characteristic equation associated to the equilibrium state is
(7) 
where
2.1 Delayindependent stability and instability results
The following delayindependent 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 (Delayindependent stability and instability).

In the nondelayed case, the equilibrium state of system (2) is locally asymptotically stable if and only if
(8) 
If the following inequality holds:
(9) then the equilibrium state of system (2) is locally asymptotically stable for any delay kernel .

If the following inequality holds:
(10) then the equilibrium state of system (2) is unstable for any delay kernel .
Proof.
1. In the nondelayed case (, for any ), the characteristic equation (7) becomes:
The necessary and sufficient condition (8) for the asymptotic stability of the equilibrium state of system (2) follows from the RouthHurwitz stability test.
2. If we assume that the characteristic equation (7) has a root in the right halfplane (), it follows that
and hence, from (7) we deduce:
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 .
2.2 Saddlenode bifurcation
Theorem 2.2 (Saddlenode bifurcation).
A saddlenode bifurcation takes place at the equilibrium state of system (2), regardless of the delay kernel , if and only if and
(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 .
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:
The Laplace transform of is
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:
with , , , , for any .
With the change of variable , the characteristic equation (7) becomes
(12) 
Denoting , it follows that the characteristic equation is equivalent to
(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:
(14) 
Indeed, we prove the following:
Lemma 2.3.
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
from Viete’s formulas, we obtain:
Therefore, we obtain the equivalence from the statement. ∎
It is worth noting that the curve intersects the saddlenode 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
(15) 
and belong to a line
(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
Assuming that such a root exists, let us denote it by . Hence:
As is a root of , we obtain:
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:
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:
At the boundary of the stability domain , the following bifurcation phenomena take place in a neighborhood of the equilibrium of system (2):

Saddlenode bifurcations take place along the open line segment ;

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

BogdanovTakens bifurcation at ;

DoubleHopf bifurcation at ;

ZeroHopf 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 saddlenode 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:
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
(17) 
Taking the derivative with respect to , we obtain:
and we express:
(18) 
Therefore, based on the parametric equations of the curve given by (14), we deduce:
Applying the real part, we finally obtain:
Taking the derivative with respect to in equation (17), we obtain:
and we express:
(19) 
Again, using the expressions of the parametric curve given by (14), we deduce:
Once again, applying the real part, we arrive at:
and thus, we obtain the gradient vector:
The tangent vector to the curve at the point is , where
Thus, we obtain the following tangent vector to the curve :
Therefore, fixing the orientation of the curve in the direction of increasing , a rightpointing normal vector to the curve is:
Hence, the directional derivative is:
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
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:
Using the fact that, , we deduce:
and therefore, using the fact that (as in the proof of Lemma 2.4), we obtain
Hence, it is easy to see that:
as . The following gradient vector is obtained:
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.