Blow-up, steady states and long time behaviour of excitatory-inhibitory nonlinear neuron models

Excitatory and inhibitory nonlinear noisy leaky integrate and fire models are often used to describe neural networks. Recently, new mathematical results have provided a better understanding of them. It has been proved that a fully excitatory network can blow-up in finite time, while a fully inhibitory network has a global in time solution for any initial data. A general description of the steady states of a purely excitatory or inhibitory network has been also given. We extend this study to the system composed of an excitatory population and an inhibitory one. We prove that this system can also blow-up in finite time and analyse its steady states and long time behaviour. Besides, we illustrate our analytical description with some numerical results. The main tools used to reach our aims are: the control of an exponential moment for the blow-up results, a more complicate strategy than that considered in C\'aceres, Carrillo and Perthame (2011), for studying the number of steady states, entropy methods combined with Poincar\'e inequalities for the long time behaviour and, finally, high order numerical schemes together with parallel computation techniques in order to obtain our numerical results.


Introduction
One of the simplest self-contained mean field models for neural networks is the Network of Noisy Leaky Integrate and Fire (NNLIF) model.It is a mesoscopic description of the neuronal network, which gives, at time t, the probability of finding a neuron with membrane potential v.There is a great deal of works about Integrate and Fire neuron models, [4,27,31,2,29,3,4,30,14,16,32,19].
However, from a mathematical point of view the properties of the model are not well known yet.
In the last years, some works have pointed to this direction; in [9] some existence results have been proven: For a fully inhibitory network there is global in time existence, while for a purely excitatory network there is a global in time solution only if the firing rate is finite for every time.These results are consistent with the fact that in the excitatory case, solutions can blow-up in finite time if the value of the connectivity parameter of the network is large enough or if the initial datum is concentrated close enough to the threshold potential, [5,7,14,13,12].Moreover, in [5,7] the set of steady states for non coupled excitatory or inhibitory networks has been described, and the long time behaviour in the linear case has been studied.Later, in [8], for the nonlinear case with small connectivity parameters it has been proved exponential convergence to the unique steady state.
In the present work we extend these results to the excitatory-inhibitory coupled NNLIF model.This model was also studied in [3], where time delay and refractory period were included.Here we focus on other aspects: We prove that, although in a purely inhibitor network the solutions are global in time (see [9]), in the presence of excitatory neurons the system can blow-up in finite time.We also analyze the set of stationary states, which is more complicate than in the case of purely excitatory or inhibitory networks, and prove exponential convergence to the unique steady state when all the connectivity parameters are small.This exponencial convergence can be demonstrated, in terms of the entropy method, since for this case the system is a "small" perturbation of the linear one.Finally, the complexity of the coupled excitatory-inhibitory network is numerically described.
The study developed in [3], together with the results in the present paper, show that this simple model (excitatory-inhibitory coupled NNLIF) could describe phenomena well known in neurophysiology: Synchronous and asynchronous states.As in [3] we will call asynchronous the states in which the firing rate tends to be constant in time and synchronous every other state.Experimental and computational results exhibiting such phenomena can be found in [3] and references therein.Thus, on one hand, when there are asymptotically stable steady states there are asyncronous states, since the firing rate tends asymptotically to be constant in time.Moreover, the presence of several steady states could provide a rich behaviour of the network, since multi-stability phenomena could appear.
On the other hand, when the model does not have stable steady states, there are syncronous states.In this sense, the blow-up phenomenon could be understood as a synchronization of part of the network, because the firing rate diverges for a finite time.Possibly, this entails that a part of the network synchronizes, and thus fires at the same time.
Other PDE models describing the spiking neurons, which are related to the Fokker-Planck system considered in the present work, are based on Fokker-Planck equations including conductance variables, [26,6,25] (and references therein) and on the time elapsed models [22,23,24].In [15] the authors study the connection between the Fokker-Planck and the elapsed models.From the microscopic point of view, the first family of PDE models is obtained assuming that the spike trains follow a Poisson process, while recently the elapsed models have been derived as mean-field limits of Hawkes processes [11,10].
The structure of the paper is as follows: In Section 2, the mathematical model and the notion of solution considered are introduced.The finite time blow-up phenomenon is studied in Section 3, and the set of steady states as well as the long time behaviour of the system are analyzed in Section 4.1.
Finally, in Section 5 the analytical results are illustrated numerically.

The model and the definition of solution
We consider a neural network with n neurons (n E excitatory and n I inhibitory) described by the Integrate and Fire model, which depicts the activity of the membrane potential.The time evolution of the membrane potential V α (t) of an inhibitory neuron (α = I) or an excitatory one (α = E) is given by the following equation (see [3,4] for details) where C m is the capacitance of the membrane, g L is the leak conductance, V L is the leak reversal potential and I α (t) is the incoming synaptic current, which models all the interactions of the neuron with other neurons.In the absence of interactions with other neurons (I α (t) = 0), the membrane potential relaxes towards a resting value V L .However, the interaction with other neurons provokes the neuron to fire, that is, it emits an action potential (spike) when V α (t) reaches its threshold or firing value V F , and the membrane potential relaxes to a reset value V R .(Let us remark that These connections are assumed to be randomly choosen, and the network to be sparsely connected, namely, ǫ = C E n E = C I n I << 1, see [3].The synaptic current I α (t) takes the form of the following stochastic process where t i Ej and t i Ij are the times of the j th -spike coming from the i th -presynaptic neuron for excitatory and inhibitory neurons, respectively, CE = C E + C ext , and J α k , for α, k = E, I are the strengths of the synapses.The stochastic character is enclosed in the distribution of the spike times of the neurons.
The spike trains of all neurons in the network are supposed to be described by Poisson processes with a common instantaneous firing rate, ν α (t), α = E, I.These processes are supposed to be independent [3,5].By using these hypotheses, the mean value of the current, µ α C (t), and its variance, σ α2 C (t), take the form Many authors [3,4,20,21] then approximate the incoming synaptic current by a continuous in time stochastic process of Ornstein-Uhlenbeck type which has the same mean and variance as the Poissonian spike-train process.Specifically, I α (t) is approached by where B t is the standard Brownian motion.
Summing up, the approximation to the stochastic diferential equation model (2.1), taking the voltage and time units so that C m = g L = 1, finally yields with the jump process The firing rate or probability of firing per unit time of the Poissonian spike train, ν α (t), is calculated in [27] as where ν α,ext is the frequency of the external input and N α (t) is the mean firing rate of the population α.Also ν I,ext = 0 since the external connections are with excitatory neurons.
Going back to (2.4), a system of coupled partial differential equations for the evolution of the probability densities ρ α (v, t) can be written, where ρ α (v, t) denotes the probability of finding a neuron in the population α, with a voltage v ∈ (−∞, V F ] at a time t ≥ 0. In [3,4,20,21,28] a heuristic argument using Itô's rule gives a system of coupled Fokker-Planck or backward Kolmogorov equations with sources The right hand sides in (2.5) represent the fact that when neurons reach the threshold potential V F , they emit a spike over the network and reset their membrane potential to the reset value V R .The system (2.5) is completed with Dirichlet boundary conditions and an initial datum In order to simplify the notation, we denote and the variable v is translated with the factor ) The coupling of the system (2.5) is hidden in these two terms, since the mean firing rates N α obey to Moreover, (2.9) gives rise to the nonlinearity of the system (2.5), since firing rates are defined in terms of boundary conditions on distribution functions ρ α .On the other hand, since ρ E and ρ I represent probability densities, the total mass should be conserved: Before introducing the definition of solution considered in this work, let us summarize some notations.For 1 ≤ p < ∞, the space of functions such that f p is integrable in Ω is denoted by L p (Ω), L ∞ (Ω) is the space of essentially bounded functions in Ω, C ∞ (Ω) is the set of infinitely differentiable functions in Ω and L 1 loc,+ (Ω) denotes the set of non-negative functions that are locally integrable in Ω.
Definition 2.1 A weak solution of (2.5)-(2.9) is a quadruple of nonnegative functions (ρ E , ρ I , N E , N I ) Additionally, if test functions of the form ψ(t)φ(v) are considered, the formulation (2.10) is equivalent holds in the distributional sense for α = E, I. Checking that weak solutions conserve the mass of the initial data is a straightforward computation after choosing φ = 1 in (2.11),

Finite time blow-up
In [5] and [7] it was proved that weak solutions can blow-up in finite time for a purely excitatory network, when neurons are considered or not to remain at a refractory state for a time.However, for a purely inhibitory network it was shown in [17] that weak solutions are global in time.The following theorem claims that a network with excitatory and inhibitory neurons can blow up in finite time.
We remark that the theorem is formulated in a more general setting of drift terms h E than that considered in (2.7).The diffusion term (2.8) of the excitatory equation is assumed to not vanish at any time.For the inhibitory firing rate we assume (3.3), which is satisfied, for instance, if N I (t) is bounded for every time.This hypothesis should not be a strong constraint, because in [9] it was proved, in the case of only one population (in average excitatory or inhibitory), that if the firing rate is bounded for every time, then there exists a global solution in time.It could be natural to think that an analogous criterion should hold in a coupled network, although its proof seems much more complicated and remains as an open problem. ) ∀v ∈ (−∞, V F ] and ∀ N I , N E ≥ 0. Assume also that there exists some M > 0 such that Then, a weak solution to the system (2.5)-(2.9)cannot be global in time because one of the following reasons: and for b E E > 0 fixed.
and the multiplier φ(v) = e µv , a weak solution (ρ E (v, t), ρ I (v, t), N E (t), N I (t)) satisfies the following inequality where assumptions (3.1)-(3.2) and the fact that v ∈ (−∞, V F ) and N E (t)φ(V R ) > 0 were used.This inequality, Gronwall's lemma and the definition of µ provide the following inequality for the exponential Using the definition of µ and (3.3), we notice that and after some computations the right hand side of the previous inequality can be bounded by Finally, the following inequality holds We observe that if the initial state satisfies On the other hand, using the again definition of µ and (3.3), we observe that Thus, e µ t 0 f (s) ds ≥ e µV F t and consequently, considering (3.6), we obtain On the other hand, since ρ E (v, t) is a probability density and µ > 0, for all t ≥ 0: which leads to a contradiction if the weak solution is assumed to be global in time.Therefore, to conclude the proof there only remains to show inequality (3.5) in the two cases of the theorem.
1.For a fixed initial datum and b E E large enough, µ, M µ (0) and φ(V F ) are fixed, thus (3.5) holds.
2. For b E E > 0 fixed, if the initial data satisfy (3.4) then condition (3.5) holds immediately.Now, there only remains to show that such initial data exist.
For that purpose we can approximate an initial Dirac mass at V F by smooth probability densities, so that ρ 0 E ≃ δ(v − V F ).This gives the following condition . So, with our initial choice of µ we can ensure that the set of initial data we are looking for is not empty.
Remark 3.2 Hypothesis (3.3) could be relaxed by ds (where C > 0).Moreover, using a priori estimates (as done in [5]) it could be proved that s) ds, which seems not to be enough to reach the whole result.Precisely, it yields the blow-up for fixed initial data and large b E E , but not for fixed b E E and concentrated initial data.
Theorem 3.1 shows that blow-up occurs when the connectivity parameter of the excitatory to excitatory synapses, b E E , is large enough or if initially there are many excitatory neurons with a voltage value of their membrane potential very close to the threshold value, V F .From a biological perspective, and in view of the numerical results in Section 5 (Figs. 3, 4 and 5), in the first case, blow-up appears due to the strong influence of the excitatory population on the behavior of the network, poducing the incontrolled growth of the firing rate at finite time.In the second case, even with small connectivity parameter b E E , the abundance of excitatory neurons with membrane potential voltage values sufficiently close to V F causes the firing rate to diverge in finite time.For a microscopic description, at the level of individual neurons, we refer to [12] and [13] where the blow-up phenomenon is also analysed.
As mentioned above it was proved in [7] that fully excitatory networks can blow-up in finite time, when one includes the refractory state.This result can be extended to the case of excitatory-inhibitory networks by following the same ideas of Theorem 3.1.

Steady states and long time behavior 4.1 Steady states
For excitatory and inhibitory networks, the study of their steady states follows a similar strategy to that for the fully excitatory or inhibitory cases.However, for coupled networks the system solved by the stationary solutions is much more complicate, as we will see throughout this section.
As in [5], let us search for continuous stationary solutions (ρ E , ρ I ) of (2.5) such that ρ E , ρ I are C 1 regular, except possibly at V = V R where they are Lipschitz, and satisfy the following identity: in the sense of distributions, where H denotes the Heaviside function.Considering and the Dirichlet boundary conditions in (2.6), we find the following initial value problem for every whose solutions are of the form The expression (4.1) is not an explicit formula for ρ α , since the right hand side depends on N α , but provides a system of implicit equations for N α for which the conservation of mass ( ) has been used.Therefore, the stationary solutions (ρ E , ρ I ) have the profile (4.1),where (N E , N I ) are positive solutions of the implicit system (4.2), and there are as many as solutions of (4.2).Of course, in the particular case of a linear system, that is V α 0 (N E , N I ) = 0 and a α (N E , N I ) independent of the firing rates, there is a unique steady state.
Determining the number of solutions of the implicit system (4.2) is a difficult task.In this section we find some conditions on the parameters of the model in order to reach this goal.Firstly, we consider the following change of variables and notation: .
With these new variables, the system (4.2) is rewritten as 2 du dz, Now, using s = z−u 2 and s = z+u 2 , the functions I 1 and I 2 can be formulated as ext for all α = E, I. Then: 1.There is an even number of steady states or there are no steady states for (2.5)-(2.9)if If b E E is large enough in comparison with the rest of connectivity parameters, there are no steady states.Specifically, there are no steady states if both (4.4) and are fulfilled, where

.6)
2. There is an odd number of steady states for (2.5)-(2.9)if If b E E is small enough in comparison with the rest of connectivity parameters, there is a unique steady state.
Proof.The proof reduces to study the existence of solutions to (4.3).It is organized in several steps.
Step 1.To prove that 1 and following analogous ideas as in [5], it is easy to observe that there is a unique This fact is a consequence of the following properties of I 2 : 2. For every N E fixed, I 2 (N E , N I ) is an increasing strictly convex function on N I , since for all 2 s k−1 (e s wF − e s wR ) ds. (4.9) Thus, lim Thus, lim 4. Using expression (4.3) for I 2 , we have I 2 (N E , 0) < ∞, for every N E fixed, since  Step 2. Properties of N I (N E ).
Obtaining an analytical expression of N I (N E ) is not easy.However, some properties of N I (N E ) are enough to prove the theorem: with respect to N E .

lim
) is a positive increasing function, thus its limit could be infinite or a constant C > 0. The finite limit leads to a contradiction since , and then lim , while lim It is a consequence of the fact that where Expression (4.12) is obtained using (4.9), (4.10), (4.11) and the definition of wF and wR .

lim
, where lim This limit is obtained after some tedious calculations, using a Taylor expansion of e sV F / √ a I − e sV R / √ a I at s = 0 and computing the following limit: for every k ∈ R, b > 0 and n ∈ N.

As a consequence of the two previous properties, lim
The limit of N ′ I (N E ) ensures that the rate of increase for N I (N E ) (in point 3 of Step 2 it was proved that N I (N E ) tends to infinity) is constant at infinity, and can be controlled in terms of the parameters b I E , b I I , V R and V F .
Step 3. Properties of the function For every fixed Step 1).To conclude the proof there remains to determine the number of solutions to With this aim, we define Depending on the properties of F we can find a different number of solutions to (4.13).First, following analogous ideas as in [5], let us show some properties of I 1 : 2. For every N E ∈ [0, ∞) fixed, I 1 (N E , N I ) is an increasing strictly convex function on N I , since for all integers k ≥ 1 Thus, lim 3. For every Thus, lim , then I 1 (N E ) decreases for N E large enough and lim , then I 1 (N E ) increases for N E large enough and lim ) is an increasing function and lim These properties are proved using 2 (e sw F − e sw R )ds.
Therefore, we have that . Now, using properties 4. and 7. of N I (N E ) the monotonicity properties of I 1 (N E ) are immediate.
Taking into account the previous points, the following properties of F are obtained: 1. F(0) = 0.
increasing function for all N E ≥ 0.

Close to
The monotonicity of F is given by the sign of its derivative 2 (e sw F − e sw R ) ds.
We observe that lim 2 (e sw F − e sw R )ds , where . So, the monotonicity properties of F can be proven as follows.
This limit is calculated using the limit of N ′ I (N E ) and proceeding in a similar way than in [5].

If b
In this case the limit of F is a product of increasing functions with limit ∞.
Step 4. The proof of the theorem is a consequence of the previous steps.
The monotonicity of F and its limit, calculated in step 3, provide the number of steady states in terms of (4.4) and (4.7), since these conditions give the range of the parameter values for which the limit of F can be compared to 1.
Under assumptions (4.4) and (4.5), there are no steady states.The reason is that for these values On the other hand, in an analogous way as in [5][Theorem 3.1(iv)], we can prove Therefore, to conclude this part of the proof we show (4.14).Defining g(N E ) := the monotonicity of g is not known, there might be several values of N * E that solve the equation for N E ≥ NE and we obtain that w F < 0 for N E ≥ NE .The inequality and due to (4.5) and (4.6).
To conclude the proof we note that there is a unique steady state for parameters where lim And, for b E E small but such that the limit of F is finite, we deduce that there is a unique stationary solution in an analogous way as in [5][Theorem 3.1(i)] for a purely excitatory network.
As proved in [5] we can find conditions for the connectivity parameters in order to have at least two steady states.We explain it in the following remark.
Remark 4.2 If the parameters of the model satisfy (4.4) and there are at least two stationary solutions to (2.5)-(2.9).
The proof is as follows.As F(0) = 0 and lim )), we have to prove that holds, and for all N E ∈ D we have Thus, w R > 0 for N E ∈ D. Finally, in an analogous way as in case (ii) of Theorem 3.1 in [5], using (4.3) we obtain 2 du dz

Long time behaviour
In [5] it was proved exponential convergence to the steady state for the linear case.Later, these results were extended in [8] to the nonlinear case, for a purely excitatory or inhibitory network with small connectivity parameters.In both papers the main tools used were the relative entropy and Poincaré inequalities.These techniques can be adapted for a coupled excitatory-inhibitory network, where the diffusion terms are considered constant, a α , where α = E, I.This is the goal of this subsection.
In Theorem 4.1 it was proved that for small connectivity parameters and constant diffusion terms there is only one stationary solution of the system (2.5)-(2.9),( . Therefore, for any smooth convex function G : R + → R, we can define the relative entropy for α = E, I as dv. Following similar computations to those developed in [5,7,8] for fast-decaying solutions to system (2.5)-(2.9),i.e. weak solutions for which the weak formulation holds for all test functions growing algebraically in v, it can be obtained, for α = E, I: where Nα (t) = N α (t) − N ∞ α .Therefore, the entropy production of the relative entropy (its time derivative) for the excitatory (resp.inhibitory) population depends on the firing rate of the inhibitory (resp.excitatory) population.In other words, both entropy productions are linked by means of the firing rates.However, for the quadratic entropy, G(x) = (x − 1) 2 , a control on the sum of them, enough and initial data (ρ 0 E , ρ 0 I ) such that Then, for fast decaying solutions to (2.5)-(2.9)there is a constant µ > 0 such that, for all t ≥ 0, Proof.The proof follows analogous steps as in [8][Theorem 2.1], with the main difficulty that, in this case, the entropy productions for excitatory and inhibitory populations are linked.This is the reason why the total relative entropy, given by the functional (4.17), has to be considered.Thus, the entropy production is the sum of the entropy productions of each population.In this way, the terms of Nα can be gathered and bound properly, as it is shown below.
Using (4.16) with G(x) = (x − 1) 2 for each term of the entropy (4.17), its time derivative can be written as where the last two terms were obtained using that G(x) − xG ′ (x) = 1 − x2 , integrating by parts and taking into account the boundary conditions (2.6).
Considering the inequality (a + b) 2 ≥ ǫ a 2 − 2b 2 , for a, b ∈ R and 0 < ǫ < 1/2, the Sobolev injection of L ∞ (I) in H 1 (I) for a small neighborhood I of V R , where ρ ∞ α is bounded from below and the Poincaré inequality (see [5][Appendix] and [8] for details), the third and fourth terms in (4.19) satisfy: for some positive constant C α 0 .
The estimate of the last two terms of (4.19) is quite more involved.First, each term is split into four addends, so that Then, Cauchy-Schwarz's and Young's inequalities provide Getting these bounds together In this way, for b α k small enough such that b < 0, the first and second terms of the right hand side of (4.20) are negative, thus , the entropy production can be bounded as follows Now, using (4.18) and Gronwall's inequality it can be proved that E[t] decreases for all times and , for all t > 0. Thus, Applying the Poincaré inequality on each term, there exist γ, γ ′ > 0 such that Finally, Gronwall's inequality concludes the proof.

Numerical results
The analytical results proved in previous sections are shown numerically in the present section.The numerical scheme considered for this purpose is based on a fifth order conservative finite difference WENO (Weighted Essentially Non-Oscillatory) scheme for the advection term, standard second order centered finite differences for the diffusion term and an explicit third order TVD (Total Variation Diminishing) Runge-Kutta scheme for the time evolution.To reduce the computation time, parallel computation techniques for a two cores code is developed.Thus, the time evolution for both equations of the system is calculated simultaneously.Each core handles one of the equations.MPI (Message Passing Interface) communication between the cores has been included in the code, since the system is coupled by the firing rates.Therefore, at the end of each Runge-Kutta step, each core needs to know the value of the firing rate of the other core.
For the simulations, an uniform mesh for v ∈ [−V lef t , V F ] is considered.The value of −V lef t is chosen such that ρ α (−V lef t , t) ∼ 0 (since ρ α (−∞, t) = 0).The time step size is adapted dynamically during the simulations via a CFL (Courant-Friedrich-Levy) time step condition.Some parameter values are common to most simulations, V F = 2, V R = 1 and ν E,ext = 0 and the diffusion terms a α (N E , N I ) have been taken constant as a α = 1.In the simulations where these values are different, the considered values are indicated in their figures and explanations.In most cases, the initial condition is where k is a constant such that However, in order to analyze the stability of steady states, stationary profiles are taken as initial conditions with and where N α is an approximated value of the stationary firing rate.
To get an idea about the number of steady states, the system (4.3) is solved numerically.For every N E fixed, N I (N E ) is calculated as the root of N I (N E )I 2 (N E , N I (N E )) − 1 = 0 using the bisection method with tolerance 10 −8 , and then a numerical approximation of F(N E ) is developed by quadrature formulas to find the number of intersections with the function 1. Fig. 2 shows the graphics of (4.12)) and F(N E ) for parameter values for which the limit of F(N E ) is finite.
The numerical approximation of the limits is in agreement with their analytical limits.
Blow-up.Numerical results for three situations in which the solutions are not global-in-time are described.In Fig. 3 we depict a blow-up situation produced by a big value of b E E , where the initial datum is far from V F .However, for a value of the connectivity parameter b E E small, we show a blow-up situation originated by an initial condition quite concentrated around V F in Fig. 4. In both cases we observe that an extremely fast increase of the firing rate of the excitatory population causes the blow-up of the system.Furthermore, we notice that the firing rate of the inhibitory population also starts to grow sharply.Nevertheless, the values it takes are quite small in comparison with those from the excitatory population.
For a purely inhibitory network the global existence of its solutions was proved in [9].Therefore, one could think that high values of b I I could prevent the blow-up of the excitatory-inhibitory system.
However, Theorem 3.1 shows that this is not the case and Fig. 5 describes this situation; although the value of b I I is big, a high value of b E E causes the divergence of N E (t) and the blow-up of the system.
Steady states.The proof of Theorem 4.1 provides a strategy to find numerically the stationary firing rates, which consists of finding the intersection points between the functions F(N E ) and constant 1 (see (4.13)).With this idea, we have plotted both functions for different parameter values.
The first case of Theorem 4.1 (there is an even number of steady states or there are no steady states) is shown in Fig. 6.In this situation, the relation between the parameters implies lim The second case of Theorem 4.1 (there is an odd number of steady states), which implies lim 1, is depicted in Fig. 7.In the left plot, b E E is small enough such that F is an increasing function, and therefore there is a unique steady state.Also, in the center plot there is only one stationary solution, but in this case the values of connectivity parameters do not guarantee the monotonicity of F. Finally, the right plot shows values of connectivity parameters for which there are three steady states.
To conclude the analysis of the number of steady states, in Fig. 8  as described in [5].For small values of b E E there is a unique steady state.As b E E increases, it appears another steady state, that merges with the first one and then disappears, apparently yielding a saddle node bifurcation.In Fig. 10 we show the time evolution of the firing rates, for the case described in Fig. 6 (right), with two steady states.We use as initial condition profiles like the one presented in (5.2).The numerical results show that the larger steady state is unstable, while the lower one seems to be stable.Therefore, numerical stability analysis indicates that the unique steady state is stable, while, when there are two steady states the highest seems to be unstable.
In this direction, another interesting bifurcation analysis that can be obtained in terms of the parameter b E I is depicted in Fig. 9.For large values of b I E there is a unique stable steady.Then, as it decreases it appears another steady state that bifurcates and gives rise to two equilibria.The largest disappears, and the lower one approaches the smallest steady state.Afterwards, both disappear, probably through a saddle node bifurcation.Numerical stability analysis determines that the lowest steady state is always stable, while the other ones are all unstable.The Fig. 7 (right) depicts three steady states.The stability analysis of this situation is more complicated than the previous one (the case of two stationary solutions).Fig. 11 shows the time evolution of the solutions for different initial data; the steady state with less firing rate seems to be stable, while the other two steady states are unstable.
In conclusion, in the present work we have extended the known results for purely excitatory or inhibitory networks [5] to excitatory-inhibitory coupled networks.We have proved that the presence of an inhibitory population in a coupled network does not avoid the blow-up phenomenon, as happens in a purely inhibitory network.Besides, we have analysed the number of steady states of the system, which is a more complicate issue than in the case of uncoupled systems.For small connectivity parameter values, we have shown that solutions converge exponentially fast to the unique steady state.Our analytical and numerical results contribute to support the NNLIF system as an appropriate model to describe well known neurophysiological phenomena, as for example synchronization/asynchronization of the network, since the blow-up in finite time might depict a synchronization of a part of the network, while the presence of a unique asymptotically stable stationary solution represents an asynchronization of the network.In addition, the abundance in the number of steady states, in terms of the connectivity parameter values, that can be observed for this simplified model, probably will help us to characterize situations of multi-stability for more complete NNLIF models and also other models including conductance variables as in [6].In [7] it was shown that if a refractory period is included in the model, there are situations of multi-stability, with two stable and one unstable steady state.In [6] bi-stability phenomena were numerically described.Multi-stable networks are related, for instance, to the visual perception and the decision making [18,1].In these directions, some tasks which remain open for the NNLIF model are: The analytical study of the stability of the network when there is more than one steady state, the possible presence of periodic solutions and the existence of solutions for coupled networks, maybe following a similiar strategy as in [17].A more extensive bifurcation analysis could probably contribute to find answers to these questions.6), for different initial conditions: ρ 0 α − 1, 2 are given by the profile (5.2) with (N E , N I ) stationary values and ρ 0 α − 3 is a normalized Maxwellian with mean 0 and variance 0.25 (see (5.1)).For both firing rates, the lower steady state seems to be asymptotically stable whereas the higher one seems to be unstable.
neuron receives C ext connections from excitatory neurons outside the network, and C = C E + C I connections from neurons in the network; C E = ǫ n E from excitatory neurons and

s 2 2 sTheorem 4 . 1
(e s wF − e s wR ) ds.When b E I = b I E = 0 the equations are uncoupled and the number of steady states can be determined in the same way as in [5], depending on the values of b E E , since for the inhibitory equation there is always a unique steady state.The following theorem establishes a classification of the number of steady states in terms of the model parameters, in the case that the connectivity parameters b E I and b I E do not vanish, and in comparison with the values of the pure connectivity parameters b E E and b I I .Assume that the connectivity parameters b E I and b

Figure 1
Figure 1 depicts the function I 2 in terms of N I for different values of N E fixed.In this figure, the properties of I 2 enumerated above can be observed.

Figure 1 :
Figure 1: The function I 2 in terms of N I , for different values of N E fixed.I 2 is an increasing function on N I and decreasing on N E .

4. I 1 4 . 1 . 4 . 2 .
(N E ) := I 1 (N E , N I (N E )) has the following properties of monotonicity: If b I E b E I < b E E b I I , then I 1 (N E ) is a decreasing function with lim N E →∞ I 1 (N E ) = 0.If b I E b E I > b E E b I I 4.2.1.and b

From a biological point
of view, the previous analysis of the number of steady states shows the complexity of the set of stationary solutions in terms of the model parameters: the reset and threshold potentials, V R and V F , and the connectivity parameters, b E E , b I E , b E I and b I I , which describe the strengths of the synapses between excitatory and inhibitory neurons.For example, in function of the connectivity parameter b

Theorem 4 . 3
) can be obtained for small connectivity parameters, proving an analogous result as in [8][Theorem 2.1].Assume a α constant for α = E, I, the connectivity parameters b E E , b I E , b E I and b I I small N E →∞ F(N E ) < 1.In the left plot, the pure connectivity parameters (b E E and b I I ) are high in comparison with the connectivity parameters b I E and b E I , in such a way that there are no steady states.In the right plot, there are two steady states because the pure connectivity parameters are small, since in this case the maximum value of F is bigger than one.
it is a comparison between an uncoupled excitatory-inhibitory network (b I E = b E I = 0) and a coupled network with small b I E and b E I .In Fig. 8 (left) the number of steady states for the uncoupled system (b I E = b E I = 0) is analyzed, while Fig. 8 (right) investigates the number of steady states when the parameters b I E and b E I are small (b I E = b E I = 0.1).As expected, it can be observed that the number of steady states does not change if we choose small values for b I E and b E I .It depends only on the value of b E E in the same way

Figure 2 :
Figure 2: Study of the limits of F(N E ) and N 2 I (N E )I(N E ) (see (4.12)) when it is finite.Left figure: b E I = 0.75, b I E = 0.5, b I I = 0.25, a E = 1, a I = 1 for different values of b E E and V F .Right figure: b E E = 1.8, b E I = 0.75, b I E = 0.5, b I I = 0.25 a E = 1 for different values of a I and V F .

Figure 3 :Figure 4 :
Figure 3: Firing rates and probability densities for b E E = 3, b E I = 0.75, b I E = 0.5, b I I = 0.25, in case of a normalized Maxwellian initial condition with mean 0 and variance 0.5 (see (5.1)).N E blows-up because of the large value of b E E .

Figure 5 :Figure 6 :
Figure 5: Firing rates and probability densities for b E E = 3, b E I = 0.75, b I E = 0.5, b I I = 3, in case of a normalized Maxwellian initial condition with mean 0 and variance 0.5 (see (5.1)).The blow-up of N E cannot be avoided by a large value of b I I .

Figure 7 :
Figure 7: F(N E ) for different parameter values corresponding to the second case of Theorem 4.1: there is an odd number of steady states.Left figure: b E E = 0.5, b E I = 0.5, b I E = 3 and b I I = 0.5 (one steady state).Center figure: b E E = 3, b E I = 9, b I E = 0.5, b I I = 0.25 (one steady state).Right figure: b E E = 3, b E I = 7, b I E = 0.5 and b I I = 0.25 (three steady states).

Figure 8 :Figure 9 :Figure 10 :
Figure 8: Comparison between an uncoupled excitatory-inhibitory network (b I E = b E I = 0) and a coupled network with small b I E and b E I .The qualitative behavior is the same in both cases.Left figure: b E I = b I E = 0, b I I = 0.25, and different values for b E E .Right figure: b E I = b I E = 0.1, b I I = 0.25 and different values for b E E .

Figure 11 :
Figure 11: Stability analysis for the case of three steady states (right in Fig. 7).Left figure: firing rates for different initial conditions: ρ 0 α − 1, 2, 3 which are given by the profile (5.2) with (N E , N I ) stationary values.Only the lowest steady state seems to be asymptotically stable.Center and right figure: evolution of the probability densities.(Simulations were developed considering v ∈ [−6, 2]).