Emergent Properties in a V1-Inspired Network of Hodgkin-Huxley Neurons

This article is devoted to the theoretical and numerical analysis of a network of excitatory and inhibitory neurons of Hodgkin-Huxley (HH) type, for which the topology is inspired by that of a single local layer of visual cortex V1. Our model is related to the recent CSY model and therefore differs from other classical models in the field. It combines a driven stochastic drive -- which may be interpreted as an ambient drive for each neuron -- with recurrent inputs resulting from the network activity. After a review of the dynamics of a single HH equation for both the deterministic and the stochastically driven case, we proceed to an analysis of the network. This analysis reveals emergent properties of the system such as partial synchronization and synchronization (defined here as a state of the network for which all the neurons spike within a short interval of time), correlation between excitatory and inhibitory conductances, and oscillations in the gamma-band frequency. The collective behavior enumerated herein is observed when the input-amplitude parameter $S^{EE}$ measuring excitatory-to-excitatory coupling (recurrent excitation) increases to within a certain range. Of note, our work indicates a distinct mechanism for obtaining the emergent properties, some of which have been classically observed. As a consequence our article contributes to the understanding of how assemblies of inhibitory and excitatory cells interact together to produce rhythms in the network. It also aims to bring problems from neuroscience to the realm of mathematics, where they can be analyzed rigorously.


Background and motivation
This paper deals with an analysis of the dynamics of a network of Hodgkin-Huxley (HH) ordinary differential equations (ODEs).We briefly recall that the HH equations were introduced in 1952, see [30], in order to describe the formation and propagation of action potentials along the giant squid nerve axon.They have served as a basis for numerous other mathematical models in neuroscience, e.g., the FitzgHugh-Nagumo (FHN) [27,38], the Moris-Lecar [37], and the Hindmarsh-Rose model [29], to cite only a few.For general textbooks on mathematical models related to the HH model, we refer to [15,21,25,26,31].All of the above-mentioned models (used to describe the electrical activity of a single neuron) have subsequently been implemented in various studies of neural ensembles, by way of networks of ODE's; related meanfield models have also been introduced [15,16,20,22,25,28,32,50,53].
Networks of ODEs have attracted a lot of attention from the applied mathematics and physics communities.We refer the reader to [11,12,39] for various standard concepts developed generally for the theory of networks with real-world topologies.Topics of interest regarding dynamical networks classically include synchrnonization, the existence of an attractor, pattern formation, bifurcation phenomena, wave propagations, etc... Our own contributions to this subject reside both in the context of networks of ODEs as well as in Reaction-Diffusion PDEs inspired by Neuroscience [1,2,4,5].The present paper represents a special thread in this line of work, as it is devoted specifically to studying emergent phenomena such as complete synchronization, partial synchronization, and rhythms in a random network with visual cortex V1-inspired topology.Biology abounds with examples of complete synchronization (all network elements coevolve with identical states), such as the synchronization of the flashes of fireflies or heart pacemaker cells.Further examples are given in [8,40,42,47] as well as the references cited therein.Another type of synchronization of specific interest in Neuroscience is the so-called partial synchronization (also sometimes referred to as cluster synchronization), observed when different subgroups of the total neural population synchronize their activity differently, see for example [13,23,43].From a physiological standpoint, the modeling of an assembly of neurons in the brain with networks of ODEs may reveal how neurons in primary sensory areas produce their responses, how cortical maps process information, etc...In particular, network models of ODEs can be used to address one of the typical questions prevalent in the physiological literature: To what extent are the responses of neuronal ensembles (and the subsequent perception of information) shaped by feed-forward vs. recurrent circuitry?This question is a current challenge and a subject of debate in neuroscience; a typical example is provided by the neurons in layer IV of the primary visual cortex, which receive sparse afferent inputs from the lateral geniculate nucleus (LGN) concurrent with abundant inputs from other cortical neurons; see [32].
In the present paper we focus further on identifying numerically a region of parameters for which a network of ODEs exhibits dynamics relevant to mathematical neuroscience.Accordingly, our setting features an ensemble of HH systems for which the coupling relies on the following assumptions: the feed-forward inputs are represented by a stochastic drive; the recurrent inputs are incorporated by way of excitatory and inhibitory synaptic coupling terms; and we are focused on three recognizable states of global activity -totally random activity, partial synchrony, and complete synchronization.
As discussed above, the basic idea of leveraging the interplay between feed-forward and recurrent connection structures is critical to the relevance of such models to neuroscience applications; among other information coming from real physical data, this has been used successfully in a series of recent papers aiming to describe the electrical activity of the visual cortex V1, see [17][18][19].See also [44,54].
The structure of the network model investigated herein is inspired by the studies carried out in these latter three articles.An important contribution of these papers is the introduction of the Recurrent Excitation Inhibition (REI) mechanism as way to obtain emergent properties observed in V1, such as the Gamma rhythm, correlation between g E and g I conductances, and partial synchronization.In addition to these contributions, this line of work features the discovery that when the parameters of the model are first tuned to the background regime, orientation selectivity arises by way of only a slight variation of the inputs, which is to say that the total energy contained in the inputs is very low in comparison to the energy coming from recurrent connections -and yet, these slight changes in the inputs have the outsize effect of inducing orientation selectivity.The authors also distinguished their contributions from a large body of work devoted to the so-called PING mechanism (Pyramidal-Interneuronal Network Gamma) [14,24,52]; see also [51], as well as [15] for a broad overview.
It is important to stress that the CSY model incorporates a variety of parameters intended to take into account physiological features whose investigation lies beyond the scope of our study.The model considered herein is separate and distinct from CSY and features HH neurons, which possess a richer structure than Leaky Integrate and Fire (LIF) models but are substantially more challenging for both theoretical and numerical work.For LIF models the reset and resting delay are set manually whereas in HH they are intrinsic properties of the dynamics.There are also further distinctions to our work: we employ Dirac impulses to model both the drive to the network and the signals between interacting units -a feature incorporated in [17] but subsequently modified with the introduction of CSY; moreover, our analysis leverages the network outputs in order to extract a single parameter (S EE ) that plays a crucial role in the onset of emergent behaviors.The ensuing parameter study indicates a path from the homogeneous stochastic state to a partial synchronization regime, global synchronization, the emergence of a Gamma rhythm, and g E -g I correlation.N.B.: in the global synchronized regime, we observe that I-population events completely encapsulate those by the E-population, i.e., the I-population spiking starts before and ends at the same time or shortly after the end of the E-spiking.This is a novel phenomenology that has to date not been observed in PING or REI.The work herein carries the potential to serve as a bridge to further rigorous mathematical analysis.Most results in this setting, including most of what we discuss in the present paper are of a computational nature.Nevertheless, to move further in the direction of rigor, our final discussion proposes three 2-dimensional models deduced from our work that can be studied mathematically.
With this background in mind, our goal is to proceed to the analysis of a model of a few hundred HH neurons driven individually by a stochastic input and connected with a topology inspired by V1 according to [17][18][19] but with the aforementioned distinctions.Proceeding hierarchically, we review and revisit some known facts on the dynamical characteristics of a single HH ODE model in Section 2. We focus specifically on the analysis for a range of the parameter I in which HH transitions from steady state to oscillatory behavior.This detail is critical to us going forward: since the current balance in individual cells resulting from inhibitory and excitatory inputs is achieved in this parameter regime, this balance replaces the parameter I in the ensuing treatment.In Section 3 we replace the parameter I by a stochastic Poisson input corresponding to a feed-forward drive, and we analyze the switch between quiescent and excitatory activity as the drive intensity is increased.Finally, we investigate the dynamics of the entire network in Section 4; we illustrate therein some numerical simulations of the network and discuss the emergence of organized behavior subject to the variation of appropriate parameters.We offer some concluding remarks and perspective regarding future work in Section 5. To set up an appropriate emphasis for our treatment, we first complete our introduction in Section 1.2 by way of describing in greater detail the network of HH equations to be considered in Section 4.

The HH ODE network model
The single unit of our network is the classical HH ODE system: In equation (1.1), V represents the membrane potential and n, m, h are gating variables modeling the opening and the closing of ionic channels; n is related to potassium fluxes, and m and h are related to sodium fluxes.Hodgkin and Huxley ( [30]) established these equations by a series of voltage measurements thanks to the Voltage Clamp Technique, which allowed them to maintain the membrane potential at constant value.Thanks to this approach, they were able to use data from such experiments in order to fit the functional parameters of their model.In brief, the model is obtained by considering the cell as an electrical circuit, and writing the Kirchoff law between internal and external currents.Thereafter the model assumes that the membrane acts as a capacitor.Ionic currents result from ionic channels acting as variable voltage-dependent resistances.The model takes into account three ionic currents: potassium (K + ), sodium (N a + ) and leakage (mainly chlorine, Cl − ).In equation (1.1), subscript t stands for the derivative d dt , I is the external membrane current, C is the membrane capacitance, g i , E i , i ∈ {K, N a, L} are the maximal conductances and the Nernst equilibrium potentials, respectively.The value of parameters is taken from [22,48].We refer to [15,21,25,31] for textbooks presenting this classical model.
The main goal of this article is to study networks of HH equations, in which each neuron receives excitatory and inhibitory inputs from its presynaptic neurons.Hence, adapting [17], we consider a network 1 ≤ i ≤ N of N = 500 HH-neurons with additional entries for excitatory and inhibitory presynaptic inputs.This leads to the following equation for the network: where Q is equal to E for E-neurons and I for I-neurons.This means that for each node in the network, we add two equations to the original HH ODE system.Those are the equations which contain coupling terms inputs coming from: 1. the network (presynaptic E and I-neurons) 2. a stochastic input drive, only for g E .These two variables stand here for gating variables and are added as such to the first equation.We now make explicit the further specifics regarding the structure of the final two equations, discuss details regarding the network topology of interest together with the relevant notation, and set up the values of fixed parameters for the study.

Network inputs
When the variable V j of a given neuron j crosses a threshold T r upward, a kick is generated in all of its postsynaptic neurons.In equation (1.3), the kicks coming from E neurons are represented mathematically by the Dirac term δ(t − s) in the g Ei equation.Accordingly, in equation (1.3), Γ E (i) denotes the set of E neurons that are presynaptic to neuron i.The notation N (j) refers to the set of times at which the neuron j crosses the threshold T r upwards.Similar notations hold for the g Ii equation.Each E neuron receives kicks with coupling strength S EE from presynaptic E neurons, and coupling strength S EI from presynaptic I neurons.Each I neuron receives kicks with coupling strength S IE from E-neurons, and coupling strength S II from I-neurons, see Figure 1.We also set Table 1.This table summarizes the value of fixed parameters used in numerical simulations of the network equation (1.3).
which makes kicks coming from presynaptic E neurons have a depolarizing effect on the membrane potential V i , while kicks coming from presynaptic I neurons to have hyperpolarizing effect.The value of the threshold T r is set to −10.

Poissonian Input
Driving conductances with Poissonian inputs are common in the literature [17,18,28,33].Accordingly, we have added a drive to the g Ei evolution equation by way of Poisson inputs of parameter ρ E = 0.9 for E-neurons and parameter ρ I = 2.7 for I neurons.Assuming that the time unit is ms, this implies that the mean value between two stochastically driven inputs is about 1 0.9 ms for E-neurons and 1 2.7 ms for I-neurons.Analogous notations to those used for the kicks coming from the network are used for the stochastic input drive: D(i) refers to the set of times at which the neuron i receives kicks corresponding to the input drive.Finally, note that the kicks have an amplitude given by the parameter S dr τ E = 0.04/2 = 0.02.

Network Topology
We consider a network of N = 500 neurons with N e = 375 E neurons and N i = 125 I neurons.The network is constructed as follows: each neuron is randomly assigned a subcollection of the network as presynaptic, according to the constraint that N •e are drawn from the E population, and N •i come from the I population.The • should be replaced with E or I, in accordance with the nature of the neuron in question.Once the selections are made for each network constituent, they remain fixed throughout the study.The values of these constants for the specific network are N ee = 50, N ei = 25, N ie = 190, and N ii = 25.The network so constructed stands as an example of an inhomogeneous random oriented graph, wherein the probability of network connections depends on cell type -I neurons have proportionally more presynaptic E neurons, as compared to the other three connection numbers.Figure 2 provides an illustration of the network.Note that this topology is inspired by the ratio implemented in [18].

Value of parameters
Table 1 summarizes the values of the fixed parameters S II , S EE , S IE , S EI and S dr .The dynamical effects of these parameters will be discussed in the sections to follow.

Analysis of the HH equation
In this section, we provide a brief expository overview of the single-cell model with an emphasis on properties of solutions to equation (1.1) that are critical for the network analysis to follow (c.f., [4,15,21,25,31]) with the declared values for the parameters.A simple analysis shows that the following theorem holds.
Proposition 2.1.There exists V m , V M ∈ R such that the compact set is positively invariant for system (1.1).
Proof.First note that the α • and β • functions are positive.For instance, this guarantees that the time-derivative of n is negative above 1 and positive below zero, indicating that the unit interval is a trapping region; this observation translates to m and h as well, altogether implying that [0, 1] 3 is positively invariant for n, m and h.As a consequence, from the first equation, V t becomes negative for V large enough.It becomes positive for −V large enough.This implies the result.
The next proposition follows from simple computations.
Proposition 2.2.The stationary solutions of (1.1) satisfy: with and Numerical simulations provide evidence that f is decreasing.Figure 3 illustrates the case I = 0, which gives a unique stationary solution with V ≃ −65.
In the network model, the terms contributing to the voltage equation give rise to effects similar to variation of the injected current I in the single-cell model.It is known, from numerical simulations [10,[34][35][36]45] and references therein, that for an interval of values for I, the system (1.1) exhibits a cascade of bifurcations giving rise to unstable and stable limit cycles, with persistence of the stable stationary point.At a certain critical value within this interval, the stationary point eventually becomes unstable trough a subcritical Hopf bifurcation.(Fig. 4, taken from [10] summarizes these facts.)The network model does not incorporate an explicit current injection; however, when in-network, the cells find themselves with contributions to the voltage equation that mimic a current injection in the appropriate range and may undergo some of the bifurcation transitions observed in the single cell.
To this point, we provide a qualitative analysis of the single-cell dynamics in a specific narrow range of the parameter I that is both dynamically interesting and relevant for its interpretation by way of the coupling contributions in the network paradigm.Figures 5-11 serve as graphical illustrations of the phenomenology we review in the rest of this section, for a pair of initial conditions and various values of the parameter I.
In Figure 5, we have drawn the bifurcation diagram for two initial conditions.It shows the coexistence of a stable limit cycle and a stable stationary solution for an interval containing [6.4, 7.5].In particular, Figure 6 illustrates this phenomenon with trajectories for I = 7.There is numerical evidence of coexistence of attractive limit cycle and stationary stable point for a region of I.
Figures 7 and 8 illustrate the temporal evolution of the variables and highlight that the dynamics is of slowfast type.Note that this property is remarkable especially because separation of timescales is not explicitly incorporated in the equation (no classical small parameter is present).
Figure 9 illustrates phase space dynamics together with nullclines.Together with Figure 8, the panels provide both dynamical and biological interpretation for the coevolution of the gating variables with V , over the same interval of V values [−80, 40]mV.We observe several distinct dynamic phases that alternate in rhythmic fashion and paint a collective picture for the evolution of the limit cycle -this picture is an important takeaway for understanding the dynamics at the network level.To start, one notes the proximity between the shapes of V and m as functions of time (see Fig. 8).The initial condition studied here leads to a very short transient prior to reaching the asymptotic regime, in this case, the stable limit cycle.For a more detailed description, we start at a time where V is at its maximum, for example t ≃ 18.This point corresponds roughly to a time at which m      (in green in Fig. 8) reaches also its maximum.On the other hand, a more attentive look at Figure 9, bottom-left (BL), shows that V reaches its maximum a little before m.Indeed, when V is at its maximum, m is increasing, n is increasing, and h is decreasing.Those dynamics are fast.The variable V decreases quickly.At t ≃ 20, V is at its minimum, and so is m (see Figs. 8 and 9 BL).Note the remarkable change in the dynamics of m and V for the period directly following this minimal point.At this time, one can distinguish the beginning of a slow phase during which the dynamics follows the m-V -projection of the nullcline.The variables V, m, h (h appears in blue in Fig. 8) increase (precipitating an increase to the sodium conductance), while n (in red in Fig. 8) decreases (and this leads to a decrease in the potassium conductance).At that point, the dynamics is still slow.Around t = 27, the dynamics of n and h changes again: n starts to decrease while h begins an increase.The dynamics are still slow until t ≃ 34, at which point we observe that V enters a fast dynamical phase that translates to the other variables as well.We observe a drastic increase in V (which means that the sodium flux dominates, recall that the E N a = 50).The variables m and n increase while h decreases.During this additional (increasing) fast phase, V reaches its maximum anew, and the loop is complete.It is worth noting that V has a fast increase followed by a fast decrease, which means that V t crosses 0 downwards in an abrupt fashion.Denoting by F (V, n, m, h) the first component of the right-hand side of (1.1), this turning point corresponds to ) is taken at the instant for which V is at its maximum and with δn > 0, δh < 0. The potassium channel plays an important role in repolarization.Analogously, m m , h m ) at the instant for which V is at its mimimum, and with δn < 0, δh > 0. It is worth noting how the trajectory in the (V, m)-plane follows the manifold n = α(V ) α(V )+β(V ) (Fig. 9 BL).Biological interpretation.According to the model, the variable n regulates the flux of potassium, while the variables m and h regulate the sodium flux.The spike in V occurs when it is pushed trough the sodium gradient induced by the value (E N a = 50), and corresponds to the depolarization of the membrane.Repolarization and hyperpolarization result from the potassium gradient (corresponding to E K = −77).According to the original paper by Hodgkin and Huxley, high values for m and n correspond to the activation of potassium and sodium gates, respectively.A low value of h corresponds to the inactivation of the sodium gate.
Action Potential and Excitability.The original HH paper was successful in part because it proposed a physiologically based mechanism for the generation of action potentials / spikes.From a dynamical modeling point of view, an action potential corresponds to a large and fast excursion in the phase space, especially in the V variable and away from the stationary stable point.This ability of the system (1.1) to go from rest into a spiking regime is known as excitability; we illustrate it for instance in Figures 7 and 8. Excitability is of fundamental importance for subsequent sections herein, since perturbations above threshold will induce a spike.
In this section, we have revisited the dynamics of the HH equation.The key point for subsequent sections is that by varying the parameter I, or alternatively, by changing the initial conditions, the HH system is able to produce spikes.This ends up being critically important for the network analysis because the N -dimensional temporal spiking sequence is the crude observable by which the behavior of the network is tracked.

One driven neuron
In this section, we focus on the dynamics of one driven neuron.We assume that I = 0.As discussed in paragraph 1, this implies that if there are no drive inputs the system evolves toward the steady state.The equation in this case writes as where D refers to the set of times at which the neuron receives kicks from the stochastic drive.In comparison with system (1.1), we have added one equation to account for the external drive contribution.Mathematically, we assume that the temporal evolution of g E incorporates jumps of amplitude S dr τ E that occur according to a Poisson process of rate λ.We order the event times in the Poisson process according to increasing value and assume the notation Recall that the Poisson sequence amounts to a realization of exponential laws with parameter λ, which provides the widths (s i+1 − s i ) of the interspike intervals.
Biological interpretation.A central question of neuroscience focuses on the mechanisms that shape the neural response in the presence of stimuli from the external environment.The cell spiking paradigm chosen in the present model reflects the interplay between external drive incorporated via Poisson kicks and recurrent inputs coming from the network coupling.We see here a direct relationship between the mathematical study (complex network science -issues of network topology, external drive, and the local node dynamics) and this fundamental biological question.The present section is aimed at understanding the effects of the external drive on a single cell.
The driven single-cell model is mathematically well defined.It features deterministic dynamics interrupted every so often by the events of a stochastic point process.The subdivision {t 0 , t 1 , ...} of the interval [0, T ] plays a critical role in the mathematical formulation.Note that the derivative g Et has to be taken in the sense of distributions.A classical computation leads to on the time interval [t i , t i+1 ), then at time t i+1 a kick arrives and, Therefore g E (t) is a discontinuous function with jumps, and is locally bounded.
In the discussion to follow, given the input drive, we study the behavior resulting from an increase in the parameter S dr .To start, we discuss theoretical results that bear relation to the cell's response to changes in this parameter.We follow this with some numerical investigations of the effects of such changes.

Some analytical results
Proposition 3.1 stands as an analogue to Proposition 2.1 and clarifies the mathematical well-definition of equation 3.1.
Proposition 3.1.The trajectories of the stochastic process generated by equation (3.1) are defined on (0, +∞) and piecewise-C 1 ; Furthermore, V is continuous, n, m, h are C 1 , g E has jumps.Furthermore, the set Proof.Existence and uniqueness on each interval [t i , t i+1 ) follows from the Cauchy theorem.At time t i+1 a jump occurs, which determines the value of g E in the next time interval.It follows that the solution g E is C ∞ in each interval (t i , t i+1 ).Next, we deal with the boundedness of trajectories.The proof of Proposition 2.1 remains valid, with the specific assumption that I = 0.There are jumps on V t but V is continuous.The derivative V t is negative if V is above E N a and positive if V is below E K .This implies the result.Remark 3.2.Note that the value of g E after the jump is therefore given by the following recurrence equation.
The following proposition follows from the above recurrence equation.
Proposition 3.3.We assume g E (0) = 0.Then, for t ∈ [t i , t i+1 ), g E is given by the following expression: Since the t k+1 − t k are governed by independent exponential laws, we can compute the value of the mean E[g E (t i )] just before a kick input.Computations lead to the following proposition.Proposition 3.4.Under the above assumptions, the following expression holds: And, Biological interpretation.Note that equation (3.2) gives simple quantitative information about the drive.Roughly speaking, it says that the mean value of g E after kicks and exponential decay is the product of the amplitude S dr and the input frequency λ (measured per ms).

Varying S dr
Next, we set values of λ = 0.9 and τ E = 2 and vary S dr .When increasing S dr , we reach a threshold at which the driven neuron starts to spike.Increasing S dr further leads to an increase in the spiking rate.This is illustrated in Figure 12.This numerical analysis serves as a basis for the study of the network for which we set ρ E = 0.9, ρ I = 2.7, τ E = 2, and S dr = 0.04.For theses values, one single driven neuron exhibits multiple spikes.It is worth noting that for the last two columns, neurons exhibit a frequencies of approximately 60 and 85Hz, respectively.These are precisely the values of parameters that will be set for the E and I-neurons; they characterize the stochastic drive, which is one of the three main elements of the network dynamics.We now turn to studying dynamics at the network level that reflect the balance between the stochastic drive and the effects of recurring inputs coming from the network coupling.

Emergent properties in a stochastically driven network
The following theorem provides a theoretical framework prior to a detailed numerical analysis.The proof is analogous to that of Proposition 3.1 and is therefore omitted.
Theorem 4.1.The trajectories of the stochastic process generated by equation (1.3) are defined on (0, +∞) and piecewise-C 1 ; for every i ∈ {1, ..., N }, V i is continuous, n i , m i , h i are C 1 ; g Ei , g Ii have jump discontinuities.Furthermore, the set (E K , E N a ) × (0, 1) 3 is positively invariant for (V i , n i , m i , h i ).
In this section, our aim is to illustrate how variation of the model parameters leads to emergent properties in the network.We present results from numerical simulations that explore the parameter space around where we take turns varying each of these parameters separately within the range [0.002; 0.03].Among these four parameters, S EE appears to be the most effective for establishing a path from stochastic homogeneity to synchronization in the network.We have structured the presentation of the results so as to highlight the following emergent phenomena: a path from homogeneity to partial synchronization and synchronization; the increase in correlation between g E and g I ; emergence of the Gamma rhythm -at some point the network has its own rhythm of oscillation consistent with the so called gamma frequency, which may be different from individual neuronal rhythms; variations in the mean spiking rates of E and I neurons that result from changes to the parameters.3.1) for S dr ∈ {0.006, 0.007, 0.008, 0.04} are illustrated; different columns correspond to different values of S dr in increasing order from left to right, except that the final two columns involve the highest value of this parameter but for two distinct values of λ.We observe a bifurcation between the spike-free and spiking regime.The first row represents V (t), the second row −g E V (t) which corresponds to I(t) for (1.1).The last row illustrates the projection in the n − V plane.The first spiking regime occurs for S dr = 0.007 (mean of −V g E ≃ 1.5).The value for λ is 0.9.For S dr = 0.006 (mean of−V g E ≃ 1.3), there are no spikes, while for S dr = 0.008, we obtain 6 spikes per second.Note that each spike for V occurs after an increase signal in −g E V .For S dr = 0.04 (mean of −V g E ≃ 7.4), there is 60 spikes per second (for comparison, for I = 7 in the deterministic case, i .e. for equation (1.1), we had 60 spikes per second).The last column corresponds to S dr = 0.04 and λ = 2.7.We observe a frequency of 84 spikes per second (mean of −V g E ≃ 22).The values of the parameters for the two last columns are those set to the input drive for E and I-neurons.

A path from random homogeneity toward partial synchronization and synchronization
We fix S EI = S IE = S II = 0.01 and vary S EE ∈ [0.01, 0.03].We then illustrate the dynamical behavior of the network for three distinct values of S EE at which the three typical states are observed.The main tool used to characterize these states is the rasterplot: for each time, we plot the neurons that are in a spiking state.The rasterplots relevant for this section are presented in the first row of Figure 17.Additionally, the reader should refer to Figures 13 and 14 indicating the statistics of interspike intervals as S EE is increased, as indicative of the transition toward organized behavior for the network.The latter will be discussed further in Section 4.6.

Random homogeneity
We start with the parameter value S EE = 0.01.The network exhibits a behavior for which no specific pattern seems to emerge; to the naked eye, it looks as though the picture could equivalently be generated by a uniform random process, see top left panel of Figure 17.For this reason, we refer to this state as randomly homogeneous.This observation is further enhanced by a side-by-side comparison of the network behavior to the simulation of a random process constructed by way of choosing HH-shaped spikes randomly and independently.Note that the mean value of spikes per second is around 11.49 for E-neurons and 48.48 for I neurons, see Table 3.A simple and relevant indicator of the total excitability of the network is given by the mean value of V across the network over time.In the left panel of Figure 16, we have plotted the mean value of V over all neurons as a function of   In this case, one can observe a clear periodicity, for which each period contains a refractory phase (when the peak of the distribution is moved to the left, this corresponds to t = 420 in dashed green), a resting state (when the peak of the distribution is around -65 mV, this corresponds to t = 429 in red) and excitation when the distribution stretches clearly to the right (this corresponds to t = 433 in dashed blue).3) for SII = SEI = SIE = SEE = 0.01.This plot shows the evolution of the mean value over (V i ) i∈{1,..,N } as a function of time.Right panel: this comparison signal was generated randomly with a mean of 10.375 spikes per ms; spikes are "generated" as follows: when appropriate, a temporal signal roughly mimicking the HH action potential over 10 ms is incorporated into voltage dynamics.This gives in mean, 10375 spikes over one second.This choice was made to approximate the 11.5 × 375 + 48.5 × 125 = 10375 spikes occurring in the network, see Table 3.Then, the signal was divided by 500 (to obtain a mean per neuron), and plotted.time.The right panel of this figure features the following content: at each time step (here the time step is set as dt = 0.01 ms), a spike shaped by analogy with the typical HH-V spike is generated with a probability of p.We choose p so that 1 dt × 100 × p ≃ 11.48 × 375 + 48.48 × 125, which means that the mean number of spikes generated by the simple stochastic iterative process is equal to the average number of spikes in the network.The signal is then divided by N = 500, so as to obtain the mean value per neuron.Note that the two plots are qualitatively similar; in particular, and especially relevant to the  present discussion, note that the scale of variation is the same for the two plots.An even more precise measure of asynchrony is the distribution of V across the network, as shown in panel A of Figure 15 and as discussed in Section 4.6.The statistics of the interspike intervals in this regime further support the characterization of random homogeneity; these measures are discussed in detail in Section 4.6 together with accompanying visual representations shown in Figures 13 and 14.

Biological interpretation
From a Neuroscience point of view, this behavior is consistent with the concept of background activity see [18,19] and references therein cited.and Note that the coupling has dramatically decreased the number of spikes per neuron, in comparison with simulations done in the previous section with neurons stimulated only by Poissonian drives.This emphasizes the inhibitory effect of I-neurons in the network dynamics: for the parameters considered here, the recurrent inputs have a global inhibitory effect.

Partial synchronization
As S EE is increased the synchronization phenomenon emerges.In this context, synchronization refers to a state of the network where all of the neurons fire within short time intervals of each other; see rasterplot in Figure 17, first row and last column.Between the random-homogeneous state and synchronization, a partial synchronization is observable, i.e., a state in which only a portion of the population will fire as part of an identifiable event.A prototypical example of this observed state occurs for the parameter values S EI = S IE = S II = 0.01 and S EE = 0.017.
The rasterplot for these parameters is reported in Figure 17, in the first row, second column.An examination of this rasterplot should be combined with a side-by-side review of all accompanying figures for these values of the parameters.
In Figure 18, second row, second column, we report the number of E and I spikes occurring in the time interval [425,450], which corresponds to an identified event.We observe that only around 190 spikes from E−neurons have been recorded during this interval.Around 190 spikes from I−neurons have also been recorded during the interval in question (some I−neurons have therefore spiked several times).Furthermore, the panel in the second row, second column of Figure 17, shows the time evolution of the potential V of a given E-neuron.More specifically, it shows 2 spikes over the interval [400,500], even though 4 events are identified in the network dynamics.The panel in the second row, second column of Figure 18, shows a specific I−neuron, which spikes 6 times during the same interval.

Synchronization
Synchronization occurs when S EE is further increased.We refer to columns 3 and 4 of Figures 17 and 18, for observation of this state.These columns correspond respectively to the following values of the parameters: Rasterplots are the most illustrative representation for the aforementioned activity and are reported in the first row of Figure 17.The first row of Figure 18, illustrates that in this case all of the neurons will spike during a given event.See also the second row of Figure 17, which illustrates a E-neuron (the specific cell indexed by 1) that spikes at each event.Synchronization in this context can be compared to excitation waves spreading over the whole network in short time intervals.The values of g E and g I clearly appear to be correlated.This can be observed in rows 3 of both Figures 17  and 18. Recall that, according to the equations, for a given neuron, g E results from the number of presynaptic E-spikes received, while g I results from the number of such I-spikes received.For a specific neuron, temporal correlations in g E and g I reflect concomitant augmentations to the spiking rates of its presynaptic E and I-neurons.This sort of correlation is typical during synchronization.It is also observed during partial synchronization, see Figure 18, row 3, column 2. Table 2 illustrates the evolution of the Pearson coefficient correlation as S EE increases.It shows a clear increase from 0.15 to 0.75.In our network, this suggests a strategy for detecting partial synchronization.Note that this type of correlation has been used in [7] to build a two-dimensional model that is able to produce various typical wandering brain rhythms.Such evidence of g E and g I correlation has been investigated thoroughly in experiments, see for example [9,41,46,49].

Gamma rhythm and neuronal frequencies
When partial synchronization and synchronization occur, events arise in the network at a frequency of 40 Hz; see the first row in Figure 17, which documents a typical oscillation in the Gamma regime.Note that in the state of partial synchronization, the frequency of the network is different from the frequency of individual E and I neurons.Note also that adjustments to S EE have a strong effect on the frequency of individual E-neurons, but limited effect on the frequency of I-neurons; see Table 3 and the third rows of Figures 17 and 18.We must also recall here that in isolation, neurons would have higher frequencies: 60 Hz for E-neurons and 84 Hz for I-neurons.For the case considered here, the network activity downshifts the spiking activity of each neuron; in some regimes, a Gamma rhythm oscillation emerges within the network.

Waves of excitation
During partial or complete synchronization, spikes do not occur simultaneously but instead spread through the network in a spatio-temporal fashion governed by the inter-neuronal connections.It would be of dynamical interest to follow these waves of excitation through the network, for the specific topology considered here.For example, see [3,6] in order to compare with other topologies or continuous equations.

I(t) and spikes
As discussed in Section 2, the spiking activity of a specific neuron depends on the value of the current I.In the network model (1.3), the parameter I corresponds to the currents induced by excitation and inhibition fluxes.This means that the dynamics of a specific neuron in the network are the same as the dynamics of a single neuron modeled by equation (1.1) with a corresponding I(t) equal to the sum of excitatory and inhibitory fluxes.Analogously, we denote by I(t) the expression This quantity is plotted in Figures 17 and 18, row 5. Some qualitative properties of individual neurons are remarkable and result from the corresponding I(t).For example, we observe that individual neurons may exhibit a kind of mixed mode oscillations (see Figs. 17 and 18, rows 3 and 6).We note also that if I(t) is varied while the neuron has already started to spike, it has a limited effect on the dynamics; see column 3 of Figure 17 -just after the second spike of the E-neuron, the current I(t) is high, but this does not generate another spike.Lastly, note that every spike has roughly the same shape, suggesting that a study of the attractor of the non-autonomous HH single equation would be of interest.

Statistics of spikes
We first report on the activity of the network by way of tracking the mean values of E and I spikes per second per neuron.We denote these outputs respectively by E ss and I ss .We have documented them in Tables 3-6.We observe that for the network and range of parameters considered here, increasing S EE has a strong effect on E ss but little effect on I ss .Next, we consider the statistics of interspike interval lengths for the four states described in Section 4.1.For E neurons, those statistics indicate that as S EE is increased, the distribution of the excitatory interspike interval lengths narrows (see Fig. 13).The mean and standard deviation of these distributions are reported in Table 7; note in particular the drastic decrease in the mean and the even more drastic decrease in the standard deviation.For I neurons, we observe also a remarkable change in the distribution of the inhibitory interspike interval lengths (see Fig. 14) as S EE is increased.There is a narrowing and a displacement of the distribution to the left.The mean and standard deviation of these distributions are reported in Table 8; we observe a slight decrease in the mean and a significant decrease in the standard deviation.
Finally, we consider the evolution of the distribution of E neurons in the voltage variable range.This is a natural analogy with meanfield and Boltzmann equations in statistical physics.This evolution provides an efficient characterization of the four states described above.In Figure 15, we represent the distribution for three specific times.In panel A, we consider S EE = 0.01 and in panel B, S EE = 0.03.This figure illustrates that in stochastic homogeneity, the distribution has very limited variations.In contrast, for S EE = 0.03, one can observe a clear periodicity, for which each period contains a refractory phase (when the peak of the distribution is moved to the left), a resting state (when the peak of the distribution is around -65) and excitation when the distribution extends significantly to the right.Remark 4.2.The main goal of this article is to provide an analysis of the network.All the numerics presented here correspond to a unique specific random given topology.It is legitimate to wonder if the results presented persist under another realizations of the random network.To answer this question we have simulated 50 different realizations of the network.The four states were observed for the same values of parameters of S EE : the results presented persist with different realizations of the network.

Conclusion
In this article, we have considered a network of inhibitory and excitatory neurons, where each cell was modeled by the HH equations.The topology chosen for the network was inspired by prior work on the visual cortex V 1.After reviewing and revisiting the dynamics of a single HH model, we have considered a single HH cell driven by an external Poissonian input.Our theoretical results and preliminary numerical calculations suggested a framework for studying the network by way of appropriate values for the parameters in the model.Continuing with this approach, we analyzed numerically the network by varying the coupling parameters S ab , a, b ∈ {I, E}.From there, we identified S EE as the most effective parameter (in the range of parameters considered) for reaching synchronization, and illuminated the transition from random homogeneity towards partial synchronization and synchronization.We also illustrated emergent phenomena such as: g E and g I correlation and Gamma-oscillations.A striking point of our study is that rhythms emerge as a property of the network activity itself: for example, in the synchronization regime, the network is oscillating at a Gamma rhythm of 40 hz, even though each individual cell features a natural oscillation frequency of 60 (for E-neurons) and 80 Hz (for I) in isolation.Finally, we would like to discuss a few two-dimensional models of interest with respect to some aspects treated in the paper.In the homogeneous stochastic regime, every neuron within its class, plays the same role and any synchronized effect in time emerges.Therefore, the dynamics could be described by two units: one E neuron and one I neuron.To fit the network outputs, these two neurons should be fed with stochastic inputs corresponding to the one observed in the network and the outputs should match the inputs: this is a kind of fixed point problem.In the synchronized regime, the E and I conductances are highly correlated.In this case the idea would be to construct coupled equations which match inputs, outputs with strong interaction between the E and I conductances (in contrast with the previous case).We refer to [7], for a two dimensional model of conductances (without membrane potential and ionic fluxes) related to this case.The partial synchronization regime would be modeled as the latter, with different input/outputs to match.

Figure 1 .
Figure 1.Schematic representation of the coupling in the network.Each E neuron receives kicks with coupling strength S EE from E neurons, coupling strength S EI from I neurons.Each I neuron receives kicks with coupling strength S IE from E neurons, coupling strength S II from I neurons.

Figure 2 .
Figure 2. Illustration of a network with fifteen E neurons and five I neurons with a similar topology than the network considered in the article (the network considered in the article is 25 times larger than the network represented here).

Figure 4 .
Figure 4. Bifurcation Diagrams as the parameter I is varied, reproduced as in [10], with permission of the authors.The figure at left indicates two bifurcating regions.One around I=10 and another around I=150.The latter one corresponds to a super-critical Hopf Bifurcation.The bifurcation around 10 is more complicated.A zoom is illustrated in the right figure.One can observe a sub-critical Hopf bifurcation around I=9.5, and bifurcation of cycles downwards.In particular, one can observe coexistence of a stationary stable point and a stable limit cycle for I approximately in the interval [6.3,9.5].

Figure 7 .
Figure 7. Temporal evolution for I = 7 and IC (V, n, m, h)(0) = (−65, 0.1, 0.1, 0.1).Left: n, m, h as functions of time, respectively in red, green and blue.Right: V as a function of time.These figures represent an action potential or spike.It is followed by a return to the stationary state.

Figure 9 .
Figure 9. Dynamics for I = 7 and IC (V, n, m, h)(0) = (−50, 0.5, 0.5, 0.5), which leads to evolution towards a limit cycle.Top left panel: figure provides the nullclines of n (red), m (green) and h (blue) as a function of V .One may rely on them to explain the dynamics.The three others pictures, illustrate respectively the dynamics of (V, n), (V, m) and (V, h), with their nullclines.These pictures highlight the dynamics of the HH model.

Figure 11 .
Figure 11.This figure illustrates the basin of attraction for I = 7, of the stationary solution and the limit-cycle.It is a projection in the V, n, h plan.In red, we plot IC which evolve to the stationary point.In blue, we plot IC which evolve to the limit cycle.

Figures 17 and 18
Figures 17 and 18 capture a significant proportion of the information content discussed in the rest of this section.

Figure 12 .
Figure12.Simulations of equation (3.1) for S dr ∈ {0.006, 0.007, 0.008, 0.04} are illustrated; different columns correspond to different values of S dr in increasing order from left to right, except that the final two columns involve the highest value of this parameter but for two distinct values of λ.We observe a bifurcation between the spike-free and spiking regime.The first row represents V (t), the second row −g E V (t) which corresponds to I(t) for (1.1).The last row illustrates the projection in the n − V plane.The first spiking regime occurs for S dr = 0.007 (mean of −V g E ≃ 1.5).The value for λ is 0.9.For S dr = 0.006 (mean of−V g E ≃ 1.3), there are no spikes, while for S dr = 0.008, we obtain 6 spikes per second.Note that each spike for V occurs after an increase signal in −g E V .For S dr = 0.04 (mean of −V g E ≃ 7.4), there is 60 spikes per second (for comparison, for I = 7 in the deterministic case, i .e. for equation (1.1), we had 60 spikes per second).The last column corresponds to S dr = 0.04 and λ = 2.7.We observe a frequency of 84 spikes per second (mean of −V g E ≃ 22).The values of the parameters for the two last columns are those set to the input drive for E and I-neurons.

Figure 15 .
Figure 15.Time Evolution of the Voltage Distribution of E neurons.Panel A: S EE = 0.01 (random homogeneity).In this case, the distribution has very limited variations as this is illustrated by the almost superposition of the three distributions corresponding to times t = 420, 429, and 433.Panel B: S EE = 0.03 (synchronization).In this case, one can observe a clear periodicity, for which each period contains a refractory phase (when the peak of the distribution is moved to the left, this corresponds to t = 420 in dashed green), a resting state (when the peak of the distribution is around -65 mV, this corresponds to t = 429 in red) and excitation when the distribution stretches clearly to the right (this corresponds to t = 433 in dashed blue).

Figure 16 .
Figure 16.Left panel: simulation of equation (1.3) for SII = SEI = SIE = SEE = 0.01.This plot shows the evolution of the mean value over (V i ) i∈{1,..,N } as a function of time.Right panel: this comparison signal was generated randomly with a mean of 10.375 spikes per ms; spikes are "generated" as follows: when appropriate, a temporal signal roughly mimicking the HH action potential over 10 ms is incorporated into voltage dynamics.This gives in mean, 10375 spikes over one second.This choice was made to approximate the 11.5 × 375 + 48.5 × 125 = 10375 spikes occurring in the network, see Table3.Then, the signal was divided by 500 (to obtain a mean per neuron), and plotted.

Figure 17 .
Figure 17.Simulation of system (1.3).This figure illustrates a path from random homogeneity to synchronization as the parameter S EE is increased.In this picture, the parameters SII = SEI = SIE = 0.01 are set and each column from left to right corresponds to a specific value of S EE .Respectively: SEE = 0.01, 0.017, 0.02 and 0.03.The first row represent the rasterplot: at each time the spiking neurons are represented by a point.On the top of the figure, in green, the I-neurons are plotted.E-neurons are plotted in red below.At left, the rasterplot illustrates a state where any event appears to be distinguishable.We call it random homogeneous activity.Increasing S EE induces synchronization.The second row represents the potential V 1 of neuron #1 as a function of a time.The third row, the E-conductance g E in red and the I-conductance g I in blue for the same neuron.The fourth row, the E-current denoted by I E in red and the I-current denoted by I I in blue.The fifth row, the sum of E and I currents which plays the role of I(t) in a single HH equation.The last row illustrates the projection of the trajectory of neuron #1 in the (V, n) phase space.

Figure 18 .
Figure 18.Simulation of system (1.3).This figure illustrates a path from random homogeneity to synchronization as the parameter S EE is increased.In this picture, the parameters SII = SEI = SIE = 0.01 are fixed and each column from left to right corresponds to a specific value of S EE .Respectively: SEE = 0.01, 0.017, 0.02 and 0.03.The first row represents the number of E and I-spikes occurring during an identified event.The rows 2 to 6 are analog to those of Figure 14, but for a I-neuron.
01 and S EE = 0.02 and S EI = S IE = S II = 0.01 and S EE = 0.03.

Table 2 .
Variation of S EE and its effect on the Pearson correlation coefficient between of g E and g I for one given excitatory neuron.

Table 3 .
Variation of S EE and its effect on the mean value of E and I-spikes per second per neuron.

Table 4 .
Variation of S IE .

Table 5 .
Variation of S EI .

Table 6 .
Variation of S II .

Table 7 .
Mean and standard deviation of EIS as S EE is increased.

Table 8 .
Mean and standard deviation of IIS as S EE is increased.EI has a notable effects both on I ss and E ss .increasing S II has a strong effect on I ss but little effect on E ss .increasing S IE has a strong effect on ss but little effect on E ss .