A Linear Analysis of Coupled Wilson-Cowan Neuronal Populations

Let a neuronal population be composed of an excitatory group interconnected to an inhibitory group. In the Wilson-Cowan model, the activity of each group of neurons is described by a first-order nonlinear differential equation. The source of the nonlinearity is the interaction between these two groups, which is represented by a sigmoidal function. Such a nonlinearity makes difficult theoretical works. Here, we analytically investigate the dynamics of a pair of coupled populations described by the Wilson-Cowan model by using a linear approximation. The analytical results are compared to numerical simulations, which show that the trajectories of this fourth-order dynamical system can converge to an equilibrium point, a limit cycle, a two-dimensional torus, or a chaotic attractor. The relevance of this study is discussed from a biological perspective.


Introduction
In the last decades, the features of action potentials propagating along axonal membranes have been related to the oscillations found in electroencephalogram (EEG) records [1][2][3]. From a theoretical perspective, analogies between EEG and nonlinear dynamics [4] were formulated prior to knowledge of the seminal numerical simulation by Lorenz, published in 1963, about a chaotic system [5]. All these studies aim to disclose the neural code; that is, to understand how information is represented, carried, and transformed by neurons [6]. The ambition is to explain how cognitive functions emerge from neuron spikes. In this context, here we analytically investigate the activity of a simple neuronal assembly, in order to examine how its parameter values influence its dynamical behavior.
Consider a population formed by two groups of interacting neurons, so that the synapses of the first group are excitatory and the synapses of the second group are inhibitory. Consider also that these two groups are connected. In 1972, Wilson and Cowan proposed a mathematical model to describe the time evolution of the activity of this neuronal population [7]. In this model, the activity ( ) of each group obeys a nonlinear first-order ordinary differential equation with the following form: in which denotes the time, is a positive constant related to the natural decay of ( ), is a positive constant proportional to the refractory period, and ( ) represents the input to this group, written as a linear combination of the activities of both groups plus the influence of external stimuli. Resting activity corresponds to = 0; thus, a negative value of expresses a depression from this background level. The function ( ) must be sigmoidal. However, as stated by Wilson and Cowan [7], "no particular significance is to be attached to the choice of" ( ). Usual choices are 1/(1 + − ) − 1/2 [8], tanh( ) [9], and arctan( ) [10]. In several studies, the properties of a single population were explored via numerical simulations [8,9,11], since these usual sigmoidal functions pose difficulties to 2 Computational Intelligence and Neuroscience analytical works. The dynamics of the two-population model was also already numerically investigated [10,[12][13][14][15]. Monteiro et al. suggested to use ( ) = / √ 1 + 2 and some analytical results were derived for the Wilson-Cowan model [16]. Here, we also take this particular form of ( ) to analytically investigate the dynamics of two coupled Wilson-Cowan populations (composed of four groups). We also assume that (1) the isolated populations are identical (in the sense that they are characterized by the same parameter values, as considered in other studies [10,[12][13][14][15]); (2) the populations are coupled by links starting from their excitatory groups (because most synapses are excitatory; for instance, they are 84% in the cat cortex [17]); (3) the refractory period is negligible (which corresponds to taking = 0 in (1), as supposed in several works [9,10,12,14,[18][19][20]).
Some authors took ( ) as a piecewise linear function [14,[19][20][21]. In our analyses, we assume that the parameter values related to this neuronal assembly allow us to take ( ) ≃ ; thus, the model becomes linear. From this approximation, the occurrence of a Hopf bifurcation in the nonlinear model is analytically inferred. Other bifurcations are numerically found. Recall that bifurcation corresponds to a qualitative change in the dynamical behavior of a system caused by variation of parameter value(s).
This manuscript about the dynamics of a pair of coupled Wilson-Cowan neuronal populations is structured as follows. In Section 2, the model is presented and analyzed by taking into account the linear approximation just described. In Section 3, the analytical results are compared to numerical simulations performed by considering the nonlinear version of the model. We found that, as the time passes, the variables of the model can converge to an equilibrium point, a limit cycle, a two-dimensional torus, or a chaotic attractor. In Section 4, the results of this work are discussed from a biological point of view.

Analytical Results
Let two identical populations be coupled by links that can be symmetrical or asymmetrical. Assume that the variables ( ) and ( ) denote the activity of the excitatory group and the activity of the inhibitory group, respectively, of th population, with = 1, 2. According to Wilson and Cowan, the dynamics of these coupled populations can be described by The parameters and represent the natural (exponential) decay. The parameters , , , and are the strengths of the connections between the excitatory and inhibitory groups of an isolated population, as shown in Figure 1: is the inhibitory connection from to , is the excitatory connection from to , is the self-excitatory connection (from to ), and is the self-inhibitory connection (from to ). The strengths of the connections between the populations are denoted by the parameters and : 1 is the connection from 2 to 1 , 2 is the connection from 1 to 2 , 1 is the connection from 2 to 1 , and 2 is the connection from 1 to 2 . All these ten parameters are positive constants. The constants 1 , 1 , 2 , and 2 stand for stimuli from external sources reaching 1 , 1 , 2 , and 2 , respectively. Thus, the dimension of the parameter space is 14.
A steady state is a stationary solution, corresponding to an equilibrium point ( * 1 , * 1 , * 2 , * 2 ) in the four-dimensional state space 1 × 1 × 2 × 2 . The constants * 1 , * 1 , * 2 , and * 2 are obtained from 1 / = 0, with ≡ − and ≡ + . Note that if the populations are isolated, that is, if = = 0, then * = ( − )/( − ) and * = ( − )/( − ), for = 1, 2. Obviously, if 1 = 2 and 1 = 2 , then * 1 = * 2 and * 1 = * 2 . The local stability of an equilibrium point can be determined from the eigenvalues of the Jacobian matrix obtained from the system of (2) linearized around this point [22][23][24]. By the Hartman-Grobman theorem, such a point is locally asymptotically stable if all eigenvalues have negative real parts [22][23][24]. For this fourth-order system, the eigenvalues are the roots of the polynomial: Routh-Hurwitz criterion states that all eigenvalues have negative real parts if 1 > 0, 2 > 0, 3 > 0, 4 > 0, . For instance, in the case of isolated populations, local stability of the steady-state solution is assured if > and > . Observe that the stability conditions do not depend on the values of and . Thus, constant external stimuli do not modify the stability of the stationary solution (in the linear approximation). Observe also that a necessary condition for stability is > ; that is, < ≡ + + or > ≡ −( + ). Therefore, the equilibrium point is locally asymptotically stable only if the strength of the self-excitatory connection is below the critical value and/or if the strength of the selfinhibitory connection is above the critical value .
If 1 = 2 = 0, then 4 diminishes as 1 and/or 2 increase; if 1 = 2 = 0, then 2 , 3 , and 4 diminish as 1 and/or 2 increase. Thus, both kinds of connections between the populations, characterized by and , can singly destabilize the equilibrium point, if their strengths exceed critical numbers.
As shown in the next section, other bifurcations can take place by varying the parameter values of this dynamical system.

Numerical Simulations
In the simulations presented in this section, the values of , , , , , 1 , 2 , 1 , and 2 are kept fixed and the values of , 1 , 2 , 1 , and 2 are varied. Thus, we investigate how the dynamics is influenced by the strengths of the self-excitatory loops and of the links between the populations.
In (a) and (b), the populations are isolated. Therefore, the attractor can be either an equilibrium point or a limit cycle, because (2) split into two decoupled second-order systems. In fact, these are the attractors that can exist in the state space of second-order autonomous nonlinear systems [22][23][24].
In (a), the conditions for the asymptotical stability of the steady-state solution are satisfied (because = 10.01 > = 7.99 and = 200 > ≃ 80.0). In this case, 1 ( ) converges to 0.167, which is the same value obtained from the analytical expression * 1 = ( 1 − 1 )/( − ). Thus, the linear approximation of the sigmoidal function ( ) can be considered as valid to determine this equilibrium point and its stability. For = = 10.02, the system experiences a Hopf bifurcation. Hence, in (b), as = 12 > , 1 ( ) tends to a periodic solution.
It is easy to check that the steady state loses its stability when > , that is, when > . In (f), with = 12 > , state space converge to a two-dimensional torus. In fact, the Lyapunov exponents are 1 ≃ 2 ≃ 0.00 and 3 ≃ 4 ≃ −0.55. Therefore, by increasing , the system experiences a bifurcation concerning the birth of a toroidal attractor. In (g), by taking 1 = 2 = 2, the trajectories tend to a limit cycle, with 1 ≃ 0.00, 2,3 ≃ −0.54, and 4 ≃ −0.58. Thus, by increasing 1 and 2 , the transition from torus towards cycle occurs.