Bifurcation analysis of Wilson-Cowan model with singular impulses

TL;DR

This paper investigates the Wilson-Cowan neural model incorporating singular impulses, introducing a novel analytical technique to study solution existence and bifurcations, revealing complex attractor structures like medusas and rings.

Contribution

It introduces the concept of singular impulses in the Wilson-Cowan model and develops a new analysis method for bifurcation and solution existence.

Findings

Bifurcations caused by impulses and singularity affect neural dynamics.

Existence of solutions is established under new conditions.

Complex attractors such as medusas and rings are identified.

Abstract

The paper concerns with Wilson-Cowan neural model with impulses. The main novelty of the study is that besides the traditional singularity of the model, we consider singular impulses. A new technique of analysis of the phenomenon is suggested. This allows to consider the existence of solutions of the model and bifurcation in ultimate neural behavior is observed through numerical simulations. The bifurcations are reasoned by impulses and singularity in the model and they concern the structure of attractors, which consist of newly introduced sets in the phase space such that medusas and rings.

Equations

$$\dot{x}=-x+S(\rho +cx),$$

$$RC\frac{dV}{dt}=-V+IR,$$

$$V(t)=IR(1-e^{-\frac{t}{RC}})$$

$$\begin{array}{rl}{\mu }_e\frac{dE}{dt}& =-E+(k_e-r_{e}E)S_e(c_{1}E-c_{2}I+P),\\ {\mu }_i\frac{dI}{dt}& =-I+(k_i-r_{i}I)S_i(c_{3}E-c_{4}I+Q),\end{array}$$

$$S_e(x)=\frac{1}{1+\exp [-a_e(x-{\theta }_e)]}-\frac{1}{1+\exp (a_e{\theta }_e)},$$

$$\begin{array}{rl}& \Delta E{|}_{t={\theta }_i}=\bar{K}(E,I),\\ & \Delta I{|}_{t={\theta }_i}=\bar{J}(E,I),\end{array}$$

$$\begin{array}{r}{\mu }_e\Delta E{|}_{t={\eta }_j}=K(E,I,{\mu }_e),\\ {\mu }_i\Delta I{|}_{t={\eta }_j}=J(E,I,{\mu }_i),\end{array}$$

$$\begin{array}{rl}& {\mu }_e\frac{dE}{dt}=-E+(k_e-r_{e}E)S_e(c_{1}E-c_{2}I+P),\\ & {\mu }_i\frac{dI}{dt}=-I+(k_i-r_{i}I)S_i(c_{3}E-c_{4}I+Q),\\ & \Delta E{|}_{t={\theta }_i}=\bar{K}(E,I),\\ & \Delta I{|}_{t={\theta }_i}=\bar{J}(E,I),\\ & {\mu }_e\Delta E{|}_{t={\eta }_j}=K(E,I,{\mu }_e),\\ & {\mu }_i\Delta I{|}_{t={\eta }_j}=J(E,I,{\mu }_i),\end{array}$$

$$\begin{array}{rl}& 0=K(E,I,0),\\ & 0=J(E,I,0).\end{array}$$

$$(E_1,I_1),(E_2,I_2),...(E_k,I_k),(E_{k+1},I_{k+1}),\dots ,(E_l,I_l)$$

$$\left(\begin{array}{c}{E}_j\\ I_j\end{array}\right)+\left(\begin{array}{c}\bar{K}(E_j,I_j)\\ \bar{J}(E_j,I_j)\end{array}\right)=\left(\begin{array}{c}{E}_i\\ I_i\end{array}\right).$$

$$\underset{\begin{array}{c}(E,I)\to (E_j,I_j)\\ {\mu }_{e,i}\to 0\end{array}}{\lim }\left(\begin{array}{c}\frac{K(E,I,{\mu }_e)}{{\mu }_e}\\ \frac{J(E,I,{\mu }_i)}{{\mu }_i}\end{array}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right),j=1,2,\dots ,k.$$

$$\underset{{\mu }_{e,i}\to 0}{\lim }(E(t),I(t))=z_j(t)\text{for}0<t\le T,$$

$$\begin{array}{rl}\mu \frac{dE}{dt}& =-E+(0.97-E)S_e(12E-4I),\\ \mu \frac{dI}{dt}& =-I+(0.98-I)S_i(13E-11I).\end{array}$$

$$\left(\begin{array}{c}E\\ I\end{array}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right),\left(\begin{array}{c}0.44234\\ 0.22751\end{array}\right),\left(\begin{array}{c}0.18816\\ 0.067243\end{array}\right).$$

$$\left(\begin{array}{cc}-0.5468& -0.1511\\ 0.2250& -1.1904\end{array}\right),\left(\begin{array}{cc}-0.9895& -0.2829\\ 2.1299& -31045\end{array}\right),$$

$$\left(\begin{array}{cc}1.0001& -0.7469\\ 0.9879& -1.9096\end{array}\right),$$

$$\begin{array}{rl}& \Delta E{|}_{t={\theta }_i}=-2E+0.44234,\\ & \Delta I{|}_{t={\theta }_i}=-2I+\mathrm{0.22751.}\end{array}$$

$$\begin{array}{rl}& \mu \Delta E{|}_{t={\eta }_i}=-\mu E^{1/2}(E-0.44234{)}^2-\sin ({\mu }^2)I,\\ & \mu \Delta I{|}_{t={\eta }_i}=-\mu I^{1/3}(I-0.22751{)}^3-\sin ({\mu }^2)E,\end{array}$$

$$\left(\begin{array}{c}0\\ 0\end{array}\right)+\left(\begin{array}{c}0.44234\\ 0.22751\end{array}\right)=\left(\begin{array}{c}0.44234\\ 0.22751\end{array}\right)$$

$$\left(\begin{array}{c}0.44234\\ 0.22751\end{array}\right)+\left(\begin{array}{c}-0.88468+0.44234\\ -0.45502+0.22751\end{array}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right).$$

$$\underset{\begin{array}{c}(E,I)\to (E_j,I_j)\\ \mu \to 0\end{array}}{\lim }\left(\begin{array}{c}-E^{1/2}(E-0.44234{)}^2-\frac{1}{\mu }\sin ({\mu }^2)I\\ -I^{1/3}(I-0.22751{)}^3-\frac{1}{\mu }\sin ({\mu }^2)E\end{array}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right),j=1,2.$$

$$\underset{\mu \to 0}{\lim }(E(t,\mu ),I(t,\mu ))=\{\begin{array}{l}(0,0)\text{ if }t\in (0,{\theta }_1]\cup ({\theta }_2,{\theta }_3]\cup \dots \\ (0.44234,0.22751)\text{ if }t\in ({\theta }_1,{\theta }_2]\cup ({\theta }_3,{\theta }_4]\cup \dots \end{array},$$

$$\underset{\mu \to 0}{\lim }(E(t,\mu ),I(t,\mu ))=\{\begin{array}{l}(0.44234,0.22751)\text{ if }t\in (0,{\theta }_1]\cup ({\theta }_2,{\theta }_3]\cup \dots \\ (0,0)\text{ if }t\in ({\theta }_1,{\theta }_2]\cup ({\theta }_3,{\theta }_4]\cup \dots \end{array}.$$

$$\begin{array}{rl}& {\mu }_e\frac{dE}{dt}=-E+(0.97-E)S_e(13E-4I),\\ & {\mu }_i\frac{dI}{dt}=-I+(0.98-I)S_i(22E-2I),\\ & \Delta E{|}_{t={\theta }_i}=6.741E^2-3.58612E+0.45064,\\ & \Delta I{|}_{t={\theta }_i}=6.6087I^2-3.85682I+0.49,\\ & {\mu }_e\Delta E{|}_{t={\eta }_i}=-{\mu }_e{E}^{1/2}(E-0.44234{)}^2-\sin ({\mu }_e^2)I,\\ & {\mu }_i\Delta I{|}_{t={\eta }_i}=-{\mu }_i{I}^{1/3}(I-0.22751{)}^3-\sin ({\mu }_i^2)E,\end{array}$$

$$\begin{array}{rl}& S_e(x)=\frac{1}{1+\exp [-1.5(x-2.5)]}-\frac{1}{1+\exp (3.75)},\\ & S_i(x)=\frac{1}{1+\exp [-6(x-4.3)]}-\frac{1}{1+\exp (25.8)}.\end{array}$$

$$\left(\begin{array}{c}0\\ 0\end{array}\right),\left(\begin{array}{c}0.20353\\ 0.18691\end{array}\right),\left(\begin{array}{c}0.45064\\ 0.49\end{array}\right)$$

$$\left(\begin{array}{c}0.096205\\ 0\end{array}\right),\left(\begin{array}{c}0.37647\\ 0.49\end{array}\right).$$

$$\begin{array}{rl}& \frac{de}{dt}=-e+(0.97-e)\widetilde{S_e}(16e-12i+1.25),\\ & \frac{di}{dt}=-i+(0.98-i)\widetilde{S_i}(15e-3i),\end{array}$$

$$\begin{array}{rl}& \widetilde{S_e}(x)=\frac{1}{1+\exp [-1.3(x-4)]}-\frac{1}{1+\exp (5.2)},\\ & \widetilde{S_i}(x)=\frac{1}{1+\exp [-2(x-3.7)]}-\frac{1}{1+\exp (7.4)}.\end{array}$$

Full text

1 Introduction

Wilson and Cowan [1] proposed a model for describing the dynamics of localized populations of excitatory and inhibitory neurons. This model is a coarse-grained description of the overall activity of a large-scale neuronal network, employing just two differential equations [2]. It is used in the developing of multi-scale mathematical model of cortical electric activity with realistic mesoscopic connectivity [3]. On the other hand, sudden changes and the instantaneous perturbations in a neural network at a certain time, which are identified by external elements, are examples of impulsive phenomena which may influence the evolutionary process of the neural network [4]. In fact, the existence of impulse is often a source of richness for a model. That is to say, the impulsive neural networks will be an appropriate description of symptoms of sudden dynamic changes. Therefore, the models considered in this paper have impulsive moments.

The singularly perturbed problems depend on a small positive parameter, which is in front of the derivative, such that the solution varies rapidly in some regions and varies slowly in other regions. They arise in the various processes and phenomena such as chemical kinetics, mathematical biology, neural networks, fluid dynamics and in a variety models for control theory [5, 6, 7, 8, 9, 10, 11]. In this article, we will investigate the Wilson-Cowan model with singular impulsive function in which singular perturbation method has been used to analyze the dynamics of neuronal models.

Local bifurcations are ubiquitous in mathematical biology [12] and mathematical neuroscience [13, 14, 15], because they provide a framework for understanding behavior of the biological networks modeled as dynamical systems. Moreover, a local bifurcation can affect the global dynamic behavior of a neuron [13]. There are many neuronal models to consider the bifurcation analysis, for instance, the bifurcation for Wilson-Cowan model is discussed in the book of Hoppensteadt and Izhikevich [13] in which they consider the model of the following type

[TABLE]

where $x\in\mathbb{R}$ is the activity of the neuron, $\rho\in\mathbb{R}$ is the external input to the neuron, the feedback parameter $c\in\mathbb{R}$ characterizes the non-linearity of the system, and $S$ is a sigma shaped function. This system consists only one neuron or one population of neurons. When the bifurcation parameter $\rho$ changes, the saddle-node bifurcation occurs. In our paper, we will discuss two and four of population of neurons. These systems have impulses of prescribed moments of time. We will observe the local bifurcation in these models.

The attractors observed in our simulations do not resemble any attractors which have already been observed in the literature. This is why, we need to introduce a new terminology to describe an ultimate behavior of motion in the model. We call the recently introduced components of constructed attractors as medusas and rings. This “zoological” approach to dynamics is not unique in differential equations. For example, canards are cycles of singularly perturbed differential equations [16, 17, 18]. They were discovered in the van der Pol oscillator by Benoit et al [19]. This phenomenon explains the very fast transition upon variation of a parameter from a small amplitude limit cycle to a relaxation oscillation [16]. The fast transition is called canard explosion and happens within an exponentially small range of the control parameter. Because this phenomenon is hard to detect it was nicknamed a canard, after the French newspaper slang word for hoax. Furthermore, the shape of these periodic orbits in phase space resemble a duck; hence the name “canard,” the French word for duck. So the notion of a canard cycle was born and the chase after these creatures began [20]. It is important to note that both canards and medusas appear in the singularly perturbed systems.

Bifurcation occurred in this paper cannot be reduced to the existing local bifurcations in the literature, namely, saddle-node, pitchfork, Hopf bifurcations, etc. First of all, we are talking about the change of an attractor set in the four subpopulations of neurons of Wilson-Cowan model with impulses depending on the change of the small parameter. This time the bifurcation parameter is also the parameter of the singularity. Moreover, it is a parameter of the singularity not only in the differential equations of the model, but also in the impulsive part of it. Thus, the cause of bifurcation is not the change of eigenvalues, but it relates to the singular compartment and the impulsive dynamics of the model. This is why, theoretical approvement of the observed bifurcations has not been done in the paper. However, we see that the abrupt changes in the phase portrait through simulations. Additionally, we notice that in the numerical study attractors of the model can be described through the new picture’s elements which we call as medusa, medusa without ring and rings, which, in general, may not be considered invariant for solutions of the model despite that the elements are introduced for the first time. We are confident that they are very generic for differential equations with impulses and they will give a big benefit in the next investigations of discontinuous neural networks.

We will start by defining the membrane time constant since it will be used as the parameter of singularity and bifurcation.

2 Membrane Time Constant

The role of the membrane time constant is important in Wilson-Cowan models. In these models the frequency of the oscillation is determined primarily by the membrane time constants [21]. Let us define the membrane time constant $\mu$ for a simple circuit. Suppose that the membrane is characterized by a single membrane capacitance $C$ in series with a single voltage-independent membrane resistance $R,$ see Fig. 1.

Then, by Ohm’s law the dynamics of the potential $V$ across this circuit in response to a current injection $I$ changes as

[TABLE]

which has the solution

[TABLE]

The membrane time constant, here, is defined by the product of the membrane resistance and membrane capacitance $\mu=RC.$ The potential $V(t)$ is governed by exponential decay toward the steady-state $V=IR$ as $\mu\to 0.$ The membrane time constant is used to understand how quickly a neuron’s voltage level changes after it receives an input signal.

3 Singular Model with Singular Impulsive Function

The dynamics of excitatory and inhibitory neurons are described as follows [1]

[TABLE]

where $E(t)$ and $I(t)$ are the proportion of excitatory and inhibitory cells firing per unit time at time $t$, respectively, $c_{1}$ and $c_{2}$ are the connectivity coefficients, which are both positive, represent the average number of excitatory and inhibitory synaptic inputs per cell, $P(t)$ represents the external input to the excitatory subpopulation, the quantities $c_{3},c_{4}$ and $Q(t)$ are defined similarly for the inhibitory subpopulation. The nonzero quantities $\mu_{e}$ and $\mu_{i}$ represent the membrane time constants while $k_{e},k_{i},r_{e}$ and $r_{i}$ are associated with the refractory terms. Moreover, $S_{e}(x)$ is the sigmoid function of the following form

[TABLE]

where $\theta_{e}$ is the position of the maximum slope of $S_{e}(x)$ and $\max[\dot{S}_{e}(x)]=a_{e}/4,$ and $S_{i}$ is defined similarly.

Since the external inputs influence the neurons activities, $E(t)$ and $I(t)$ can change abruptly. It is natural to consider the previous continuous dynamics in the way that the membrane time constants proceed to be involved in the electrical processes and the impulsive equations have the form

[TABLE]

where the impulse moments $\theta_{i}$s are distinct, $\theta_{i}\in(0,T)$ and the equality $\Delta E|_{t=\theta_{i}}=E(\theta+)-E(\theta-)$ denotes the jump operator in which $t=\theta$ is the time when the external input influence $E(t)$, $E(\theta-)$ is the pre-impulse value and $E(\theta+)$ is the post-impulse value. Moreover, if one considers the impulsive equations as the limit cases of the differential equations, then at some moments impulsive changes of the activities can depend on the membrane time constants, similar to the ones for the system (1). More precisely, we will also study the equations of the form

[TABLE]

where the moments $\eta_{j}$s are, in general, different from $\theta_{i}$s. Finally, gathering all the dynamics details formulated above, our single Wilson-Cowan model with impulses has the following form

[TABLE]

with the initial activity $(E(0),I(0))=(E_{0},I_{0}).$

Define the function $F(E,I)=\begin{pmatrix}-E+(k_{e}-r_{e}E)S_{e}(c_{1}E-c_{2}I+P)\\ -I+(k_{i}-r_{i}I)S_{i}(c_{3}E-c_{4}I+Q)\end{pmatrix}.$

Suppose that $E,I\in\mathbb{R}$, $t\in[0,T],$ $F(E,I)$ is continuously differentiable on $D$, $K(E,I,\mu_{e}),J(E,I,\mu_{i})$ are continuous on $D\times[0,1]$ and $\bar{K}(E,I),\bar{J}(E,I)$ are continuous on $D$, $D$ is the domain $D=\{0\leq t\leq T,|E|<d,|I|<d\},$ $\theta_{i},i=1,2,\dots,p,$ and $\eta_{j},j=1,2,\dots,\bar{p},$ are distinct discontinuity moments in $(0,T).$

Substituting $\mu_{e}=\mu_{i}=0$ in (1) and (4), we obtain $F(E,I)=0$ and

[TABLE]

Assume that equations $F(E,I)=0$ and (6) have the steady states

[TABLE]

such that all of them are real and isolated in $\bar{D}$. They are considered to be states of low level background activities since such activities seem ubiquitous in neural tissue. $E(t)$ and $I(t)$ will be used to refer the activities in the respective subpopulations.

The following condition are required for system (1).

• (C1)

Jacobian matrices of $F(E,I)$ at the points $(E_{1},I_{1}),(E_{2},I_{2}),...,(E_{k},I_{k})$ are Hurwitz matrices (they have eigenvalues whose real parts are negative).

This condition implies that the states $(E_{1},I_{1}),(E_{2},I_{2}),...,(E_{k},I_{k})$ are stable steady states of the differential equation (1). Moreover, for the impulsive functions we need the following conditions.

• (C2)

For each $j\in\{1,2,\dots,k\}$ there exists $i\in\{1,2,\dots,k\}$ such that

[TABLE]

That is, after the each impulse moment $\theta_{j}$ the activity $(E(t),I(t))$ will be close to another stable steady state $\begin{pmatrix}E_{i}\\ I_{i}\end{pmatrix}.$

• (C3)

[TABLE]

In the denominator of the limit we have small parameters $\mu_{e}$ and $\mu_{i}$ which decay to zero. In order to avoid a blow up we need the last condition. In addition, the zero value of the limit gives us the privilege that the activities stay in the domain of attractions of the stable steady states.

Denote $D_{j}$ as the domain of attraction of stable steady state $(E_{j},I_{j}),j=1,2,\dots,k,$ such that $D_{i}\cap D_{j}=\emptyset$ if $i\neq j$ and $D_{j}\subset D,j=1,2,\dots,k.$ Also, $z_{j}(t)$ will be used for denoting the solution of $F(E,I)=$ and (6) such that if the initial value $(E_{0},I_{0})\in D_{j},$ then $z_{j}(t)=(E_{j},I_{j})$ for $t\in(0,\theta_{1}]$ and it alternates to the other stable steady states by condition (C2) for the next intervals $(\theta_{i},\theta_{i+1}],i=1,2,\dots,p-1.$

Theorem 3.1

Suppose that conditions (C1)-(C3) are true. If the initial value $(E_{0},I_{0})$ is located in the domain of attraction $D_{j}$ of the steady state $(E_{j},I_{j}),j=1,2,\dots,k,$ then the solution $(E(t),I(t))$ of (5) with $(E_{0},I_{0})$ exists on $[0,T]$ and it is satisfies the limit

[TABLE]

where $j=1,2,\dots,k(k-1)^{p}.$

The proof follows from the proof in [22].

**Example. ** Now, let us take the external forces $P(t)=Q(t)=0$, $\mu_{e}=\mu_{i}=\mu,$ and other coefficients in (1) as follows: $c_{1}=12,c_{2}=4,c_{3}=13,c_{4}=11,a_{e}=1.2,a_{i}=1,\theta_{e}=2.8,\theta_{i}=4,r_{e}=1,r_{i}=1,k_{e}=0.97,k_{i}=0.98.$ Then, one obtains

[TABLE]

Taking $\mu=0,$ one has the three equilibria (see Fig. 2), namely

[TABLE]

We have $F(E,I)=\begin{pmatrix}-E+(0.97-E)S_{e}(12E-4I)\\ -I+(0.98-I)S_{i}(13E-11I)\end{pmatrix}.$ Then, the Jacobian matrices of $F(E,I)$ on the steady states are

[TABLE]

[TABLE]

respectively. All eigenvalues of the first two matrices are negative, but last one has a positive eigenvalue. Therefore, the first two steady states are stable.

We extend model (8) with the following impulse functions

[TABLE]

[TABLE]

where $\theta_{i}=\frac{2i}{3},\eta_{i}=\frac{2i-1}{3},i=1,2,\dots,20.$ Let us check the conditions of Theorem 3.1. We have shown that the states $\begin{pmatrix}0\\ 0\end{pmatrix},\begin{pmatrix}0.44234\\ 0.22751\end{pmatrix}$ are stable. Moreover, they satisfy the equations (10) if $\mu=0.$ Condition (C2) holds since

[TABLE]

and

[TABLE]

Lastly, let us check the condition (C3):

[TABLE]

Clearly, all conditions are satisfied. Therefore, if the initial value $(E_{0},I_{0})$ is in the domain of attraction of the steady state $(0,0)$ then the activities $(E(t),I(t))$ approaches to the steady states as $\mu\to 0$, that is to say,

[TABLE]

and if it is in the domain of attraction of the steady state $(0.44234,0.22751),$ then

[TABLE]

To demonstrate the results via simulation, we take $(E_{0},I_{0})=(0.25,0)$ which is in the domain of attraction of $(0.44234,0.22751).$ Obviously, the results of the theorem can be seen in Fig. 3.

4 Bifurcation of New Attractor Composed of Medusa

In discontinuous dynamics, we show that a new type of attractor consisting medusa, medusa without ring, and rings exist. For this purpose, we study a pair of coupled Wilson-Cowan models. Here, we have four subpopulation in which each system has an excitatory and an inhibitory subpopulation. The first system admits three stable steady states and it has singular impulses. The second one has a limit cycle and it does not have any impulse effects. Also, in the latter system, the membrane time constants are equals to 1.

The first Wilson-Cowan model with impulsive singularity is of the following form:

[TABLE]

where the sigmoid functions are

[TABLE]

and impulse moments are $\theta_{i}=2i+4.95,\eta_{i}=2i-1+4.95,i=1,2,\dots,50.$ The differential equations in (11) have three stable states

[TABLE]

and two unstable steady states

[TABLE]

The second model, which has a limit cycle , is of the form

[TABLE]

where

[TABLE]

We couple system (11) and (12) as follows

[TABLE]

It is already known that differential equations in (11) has three stable steady states. Suppose that the membrane time constants in (13) are equal such that $\mu_{e}=\mu_{i}=\mu$ and the initial condition is $(0.4656,0.1101,0.1101,0.04766).$ Clearly, in Fig. 4, one can observe that a medusa exist for the value of parameter $\mu=0.05.$ Note that this is a single trajectory and its form looks like a medusa.

Fig. 4 is formed as follows. The (E,e,i)-coordinate which start at the given initial value approaches to the cycle. It moves around the cycle until the impulse moment $\eta_{1}.$ When the time reaches $t=\eta_{1},$ because of the impulse function the coordinate jump to $(E(\eta_{1}+,\mu),e(\eta_{1}+,\mu),i(\eta_{1}+,\mu)).$ Again it will approach the cycle and move until the impulse moment $t=\theta_{1}.$ Then the coordinate jumps to $(E(\theta_{1}+,\mu),e(\theta_{1}+,\mu),i(\theta_{1}+,\mu))$ and it will approach to the cycle. The (E,e,i)-coordinate moves in this pattern and finally the medusa in Fig. 4 is observed. The pattern is visualized in Fig. 5.

In the following figures, we will see that for different values of the parameter $\mu$ and for the different values of the initial conditions, we will obtain different medusas and rings. First of all, consider the system (13) with the initial values $(-0.01,0,0.17,0.25),$ $(0.21,0.20,0.20,0.15),$ and with the parameter $\mu=1$ to get Fig. 6. In this figure, there are two medusas without ring and a cycle. Indeed, they are a single trajectory, which is disconnected in the geometrical sense, but it is connected in the dynamics sense.

Next, we change the parameter to $\mu=0.2$ and use different initial activations $(-0.01,0,0.17,0.25),$ $(0.21,0.20,0.20,0.15),$ $(0.5,0.5,0.3,0.3).$

In Fig. 7, one medusa and two different rings are emerged. Geometrically, the attractor is disconnected. However, it is connected in the dynamics sense since it is a single attractor with three parts. There does not exist any limit cycle. The cycles which look like limits cycles are just parts of the whole trajectory.

Let us consider Fig. 8. In this figure initial activations are same as in 7.

The parameter is fixed and $\mu=0.1.$ Although the initial values are different, the trajectories eventually obtain the shape of the same medusa. There is an alone red trajectory. It is a part of the whole red trajectory. Therefore, it is neither a limit cycle nor a ring since the trajectory never comes to the neighborhood of it.

Finally, fix the parameter $\mu=0.05.$ In Fig. 9, any trajectory from the different initial values blue$(-0.01,0,0.17,0.25),$ red$(0.21,0.20,0.20,0.15),$ magenta$(0.5,0.5,0.3,0.3),$ ultimately gets the form of red or blue medusa. The blue and the magenta trajectories converges to the same medusa. This is why, we will say that the attractor consists of two disjoint medusas. They are disjoint since there is not a single trajectory which makes two medusas.

Note that as the parameter decreases the form of the trajectory becomes horizontal through the E-coordinate.

In conclusion, we see that the neuron’s dynamics have the following properties: in the case $\mu=1$ two medusas without ring and a cycle are obtained. When $\mu=0.2$ one medusa and two rings emerge and when $\mu=0.1$ one medusa emerges. Finally, if $\mu=0.05$ two medusas emerge. These results demonstrate that for different values of $\mu$ the qualitative changes in the behavior of trajectories of (13) occur ultimately. Therefore, we have a bifurcation. It is important to note that this bifurcation occurs because of the singularity and impulses. This is why, one cannot explain the bifurcations in this paper through the traditional types of bifurcations, saddle-node, pitchfork, Hopf bifurcation, etc. For example, the change of the numbers of medusas and rings in the local phase portrait depend on the impulsive jumps sizes. Bifurcation, here, also depends on the positions of cycles for unperturbed system.

5 Conclusion

It is the first time that only a single small parameter $\mu$ causes not only to the singularity, but also to the bifurcation. The singularity in this paper is a new kind such that it emerges both from the differential equation part and in the impulsive function. It is also important that the small parameter $\mu$ is a natural parameter which comes from the membrane time constant in Wilson-Cowan neuron model. We have shown the existence of bifurcation through the simulations. Theoretical proofs are not given since it is difficult to analyze the discontinuous dynamics of the model in which a single parameter causes both singularity and bifurcation. Therefore, bifurcation is not occurred by the change of eigenvalues, but it relates to the singular compartment and the impulsive dynamics of the model.

New type of attractor, which consists medusa, medusa without ring and rings, is defined. The name comes from the similarity of the form of trajectory and medusa.