On Dynamics of Integrate-and-Fire Neural Networks with Conductance Based Synapses

We present a mathematical analysis of networks with integrate-and-fire (IF) neurons with conductance based synapses. Taking into account the realistic fact that the spike time is only known within some finite precision, we propose a model where spikes are effective at times multiple of a characteristic time scale δ, where δ can be arbitrary small (in particular, well beyond the numerical precision). We make a complete mathematical characterization of the model-dynamics and obtain the following results. The asymptotic dynamics is composed by finitely many stable periodic orbits, whose number and period can be arbitrary large and can diverge in a region of the synaptic weights space, traditionally called the “edge of chaos”, a notion mathematically well defined in the present paper. Furthermore, except at the edge of chaos, there is a one-to-one correspondence between the membrane potential trajectories and the raster plot. This shows that the neural code is entirely “in the spikes” in this case. As a key tool, we introduce an order parameter, easy to compute numerically, and closely related to a natural notion of entropy, providing a relevant characterization of the computational capabilities of the network. This allows us to compare the computational capabilities of leaky and IF models and conductance based models. The present study considers networks with constant input, and without time-dependent plasticity, but the framework has been designed for both extensions.

Neuronal networks have the capacity to treat incoming information, performing complex computational tasks (see (Rieke, Warland, Steveninck, & Bialek, 1996) for a deep review), including sensory-motor tasks.It is a crucial challenge to understand how this information is encoded and transformed.However, when considering in vivo neuronal networks, information treatment proceeds usually from the interaction of many different functional units having different structures and roles, and interacting in a complex way.As a result, many time and space scales are involved.Also, in vivo neuronal systems are not isolated objects and have strong interactions with the external world, that hinder the study of a specific mechanism (Frégnac, 2004).In vitro preparations are less subject to these restrictions, but it is still difficult to design specific neuronal structure in order to investigate the role of such systems regarding information treatment (Koch & Segev, 1998).In this context models are often proposed, sufficiently close from neuronal networks to keep essential biological features, but also sufficiently simplified to achieve a characterization of their dynamics, the most often numerically and, when possible, analytically (Gerstner & Kistler, 2002b;Dayan & Abbott, 2001).This is always a delicate compromise.At one extreme, one reproduces all known features of ionic channels, neurons, synapses ... and lose the hope to have any (mathematics and even numeric) control on what is going on.At the other extreme, over-simplified models can lose important biological features.Moreover, sharp simplifications may reveal exotic properties which are in fact induced by the model itself, but do not exist in the real system.This is a crucial aspect in theoretical neuroscience, where one must not forget that models are subject to hypothesis and have therefore intrinsic limits.
For example, it is widely believed that one of the major advantages of the Integrate-and-Fire (IF) model is its conceptual simplicity and analytical tractability that can be used to explore some general principles of neurodynamics and coding.However, though the first IF model was introduced in 1907 by Lapicque (Lapicque, 1907) and though many important analytical and rigorous results have been published, there are essential parts missing in the state of the art in theory concerning the dynamics of IF neurons (see e.g.(Mirollo & Strogatz, 1990;Ernst, Pawelzik, & Geisel, 1995;Senn & Urbanczik, 2001;Timme, Wolf, & Geisel, 2002;Memmesheimer & Timme, 2006;Gong & Leeuwen, 2007;Jahnke, Memmesheimer, & Timme, 2008) and references below for analytically solvable network models of spiking neurons).Moreover, while the analysis of an isolated neuron submitted to constant inputs is straightforward, the action of a periodic current on a neuron reveals already an astonishing complexity and the mathematical analysis requires elaborated methods from dynamical systems theory (Keener, Hoppensteadt, & Rinzel, 1981;Coombes, 1999b;Coombes & Bressloff, 1999).In the same way, the computation of the spike train probability distribution resulting from the action of a Brownian noise on an IF neuron is not a completely straightforward exercise (Knight, 1972;Gerstner & Kistler, 2002a;Brunel & Sergi, 1998;Brunel & Latham, 2003;Touboul & Faugeras, 2007) and may require rather elaborated mathematics.At the level of networks the situation is even worse, and the techniques used for the analysis of a single neuron are not easily extensible to the network case.For example, Bressloff and Coombes (Bressloff & Coombes, 2000b) have extended the analysis in (Keener et al., 1981;Coombes, 1999b;Coombes & Bressloff, 1999) to the dynamics of strongly coupled spiking neurons, but restricted to networks with specific architectures and under restrictive assumptions on the firing times.Chow and Kopell (Chow & Kopell, 2000) studied IF neurons coupled with gap junctions but the analysis for large networks assumes constant synaptic weights.Brunel and Hakim (Brunel & Hakim, 1999) extended the Fokker-Planck analysis combined to a mean-field approach to the case of a network with inhibitory synaptic couplings but under the assumptions that all synaptic weights are equal.However, synaptic weight variability plays a crucial role in the dynamics, as revealed e.g. using mean-field methods or numerical simulations (VanVreeswijk & Hansel, 1997;VanVreeswijk & Sompolinsky, 1998;VanVreeswijk, 2004).Mean-field methods allow the analysis of networks with random synaptic weights (Amari, 1972;Sompolinsky, Crisanti, & Sommers, 1988;Cessac, Doyon, Quoy, & Samuelides, 1994;Cessac, 1995;Hansel & Mato, 2003;Soula, Beslon, & Mazet, 2006;Samuelides & Cessac, 2007) but they require a "thermodynamic limit" where the number of neurons tends to infinity and finite-size corrections are rather difficult to obtain.Moreover, the rigorous derivation of the mean-field equations, that requires large-deviations techniques (BenArous & Guionnet, 1995), has not been yet done for the case of IF networks with continuous time dynamics (for the discrete time case see (Soula et al., 2006;Samuelides & Cessac, 2007)).
Therefore, the "analytical tractability" of IF models is far from being evident.In the same way, the "conceptual simplicity" hides real difficulties which are mainly due to the following reasons.IF models introduce a discontinuity in the dynamics whenever a neuron crosses a threshold: this discontinuity, that mimics a "spike", maps instantaneously the membrane potential from the threshold value to a reset value.The conjunction of continuous time dynamics and instantaneous reset leads to real conceptual and mathematical difficulties.For example, an IF neuron without refractory period (many authors have considered this situation), can, depending on parameters such as synaptic weights, fire uncountably many spikes within a finite time interval, leading to events which are not measurable (in the sense of probability theory).This prevents the use of standard methods in probability theory and notations such as ρ(t) = n i=1 δ(t − t i ) (spike response function ) simply lose their meaning1 .Note also that the information theory (e.g. the Shannon theorem, stating that the sampling period must be less than half the period corresponding to the highest signal frequency) is not applicable with unbounded frequencies.But IF models have an unbounded frequencies spectrum (corresponding to instantaneous reset).From the information theoretic point of view, it is a temptation to relate this spurious property to the erroneous fact that the neuronal network information is not bounded.These few examples illustrate that one should not be abused by the apparent simplicity of IF models and must be careful in pushing too much their validity in order to explore some general principles of neurodynamics and coding.
The situation is not necessarily better when considering numerical implementations of IF neurons.Indeed, it is known from a long time that the way the membrane potential is reset in a neuronal network simulation have significant consequences for the dynamics of the model.In particular, Hansel et al (Hansel, Mato, Meunier, & Neltner, 1998) showed that a naive implementation of integrate-and-fire dynamics on a discrete time grid introduces spurious effects and proposed an heuristic method to reduce the errors induced by time discretization.In parallel, many people have developed event based integration schemes (Brette et al., 2007), using the fact that the time of spike of a neuron receiving instantaneous spikes from other neurons can be computed analytically, thus reducing consequently the computation time and affording the simulation of very large networks.In addition, exact event based computational schemes are typically used for the above-mentioned analytically tractable model classes, see, e.g.(Mirollo & Strogatz, 1990;Timme et al., 2002).Unfortunately, this approach suffers two handicaps.If one considers more elaborated models than analytically tractable models, one is rapidly faced to the difficulty of finding an analytical expression for the next spike time (Rudolph & Destexhe, 2006).Moreover, any numerical implementations of a neural network model will necessarily introduce errors compared to the exact solution.The question is : how does this error behave when iterating the dynamics ?Is it amplified or damped ?In IF models, as set previously, these errors are due to the discontinuity in the membrane potential reset and to the time discretization.This has been nicely discussed by Hansel et al (Hansel et al., 1998) paper.These authors point out two important effects.When a neuron fires a spike between time t and t + ∆t a local error on the firing time is made when using time discretization.First, this leads to an error on the membrane potential and second this error is propagated to the other neurons via the synaptic interaction term.Unfortunately, this analysis, based on numerical simulations, was restricted to a specific architecture (identical excitatory neurons) and therefore the conclusions drawn by the authors cannot be extended as it is to arbitrary neural architectures.Indeed, as we show in the present paper, the small error induced by time discretization can be amplified or damped, depending on the synaptic weights value.This leads to the necessity of considering carefully (that is mathematically) the spurious effects induced by continuous time and instantaneous reset in IF models, as well as the effects of time discretization.This is one aspect discussed in the present paper.
More generally, this work contains several conclusions forming a logical chain.After a discussion on the characteristic times involved in real neurons and comparison to the assumptions used in Integrate and Fire models we argue that discrete time IF models with synchronous dynamics can be used to model real neurons as well, provided that the time scale discretization is sufficiently small.More precisely, we claim that IF equations are inappropriate if one sticks to much on the instantaneous reset and spike time, but that they provide a good and mathematically tractable model if one allows reset and spike to have some duration.We therefore modify the reset and spike definition (while keeping the differential equation for the dynamics of the membrane potential below the threshold).The goal is however NOT to propose yet another numerical scheme for the numerical integration of continuous time IF models.Instead, our aim is to analyze mathematically the main properties of the corresponding dynamical system, describing the evolution of a network with an arbitrary, finite, size (i.e.we don't use neither a mean-field approach nor a thermodynamic limit).We also consider an arbitrary architecture.Finally, in our analysis the time discretization step is arbitrary small, (thus possibly well below the numerical precision).For this, we use a dynamical system approach developed formerly in (Blanchard, Cessac, & Krueger, 2000;Cessac, Blanchard, Krueger, & Meunier, 2004).In particular, in (Cessac, 2008), a discrete time version of a leaky Integrate-and-Fire network, was studied.It was shown that the dynamics is generically periodic, but the periods can become arbitrary large (in particular, they can be larger than any accessible computational time) and in (non generic) regions of the synaptic weights space, dynamics is chaotic.In fact, a complete classification of the dynamical regimes exhibited by this class of IF models was proposed and a one-toone correspondence between membrane potential trajectories and raster plots was exhibited (for recent contributions that study periodic orbits in large networks of IF neurons, see (Jahnke et al., 2008;Gong & Leeuwen, 2007).).Beyond these mathematical results, this work warns one about some conclusions drawn from numerical simulations and emphasizes the necessity to have, when possible, a rigorous analysis of the dynamics.
The paper (Cessac, 2008) dealt however with a rather simple version of IF neurons (leaky Integrate and Fire) and one may wonder whether this analysis extend to models closer to biology.In the present paper we extend these results, and give a mathematical treatment of the dynamics of spikes generated in synaptic coupled integrate-and-fire networks where synaptic currents are modeled in a biophysically plausible way (conductance based synapses).As developed in the text, this extension is far from being straightforward and requires a careful definition of dynamics incorporating the integration on the spikes arising in the past.This requires a relatively technical construction but this provides a setting where a rigorous classification of dynamics arising in IF neural networks with conductance based synapse can be made, with possible further extended to more elaborated models.
The paper is organized as follows.In section 1 we give a short presentation of continuous time Integrate-and-Fire models.Then, a careful discussion about the natural time scales involved in biological neurons dynamics and how continuous time IF models violate these conditions is presented.From this discussion we propose the related discrete time model.Section 2 makes the mathematical analysis of the model and mathematical results characterizing its dynamics are presented.Moreover, we introduce an order parameter, called d(Ω, S), which measures how close to the threshold are neurons during their evolution.Dynamics is periodic whenever d(Ω, S) is positive, but the typical orbit period can diverge when it tends to zero.This parameter is therefore related to an effective entropy within a finite time horizon, and to the neural network capability of producing distinct spikes trains.In other words, this is a way to measure the ability of the system to emulate different input-output functions.See (Langton, 1990;Bertschinger & Natschläger, 2004) for a discussion on the link between the system dynamics and its related computational complexity 2 .The smaller d(Ω, S), the larger is the set of distinct spikes trains that the neural network is able to produce.This implies in particular a larger variability in the responses to stimuli.The vanishing of d(Ω, S) corresponds to a region in the parameters space, called "the edge of chaos", and defined here in mathematically precise way.In section 3 we perform numerical investigations of d(Ω, S) in different models from leaky Integrate-and-Fire to conductance based models.These simulations suggest that there is a wide region of synaptic conductances where conductance based models display a large effective entropy, while this region is thinner for leaky Integrate-and-Fire models.This provides a quantitative way to measuring how conductances based synapses and currents enhances the information capacity of integrate and fire models.Section 4 proposes a discussion on these results.
1.1 General structure of Integrate and Fire models.
We consider the (deterministic) evolution of a set of N neurons.Call V k (t) the membrane potential of neuron k ∈ {1 . . .N } at time t and let V(t) be the vector [V k (t)] N k=1 .We denote by V ≡ V(0) the initial condition and the (forward) trajectory of V by: where time can be either continuous or discrete.In the examples considered here the membrane potential of all neurons is uniformly bounded, from above and below, by some values V min , V max .Call M = [V min , V max ] N .This is the phase space of our dynamical system.We are focusing here on "Integrate and Fire models", which always incorporate two regimes.For the clarity of the subsequent developments we briefly review these regimes (in a reverse order).
Fix a real number θ ∈ [V min , V max ] called the firing threshold of the neuron 3 .Define the firing 2 It has been proposed that optimal computational capabilities are achieved by systems whose dynamics is neither chaotic nor ordered but somewhere in-between order and chaos.This has led to the idea of computation at "the edge of chaos".Early evidence for this hypothesis has been reported by Kauffman (Kauffman, 1969) and Langton (Langton, 1990) considering cellular automata behavior, and Packard (Packard, 1988) using a genetic algorithm.See (Bertschinger & Natschläger, 2004) for a review.In relation, with these works, theoretical results by Derrida and co-authors (Derrida & Pomeau, 1986;Derrida & Flyvbjerg, 1986) allow to characterize analytically the dynamics of random Boolean networks and for networks of threshold elements (Derrida, 1987).Recently (Bertschinger & Natschläger, 2004) have contributed to this question, considering numerical experiments in the context of real-time computation with recurrent neural networks.
3 We assume that all neurons have the same firing threshold.The notion of threshold is already an approximation which is not sharply defined in Hodgkin-Huxley (Hodgkin & Huxley, 1952) or Fitzhugh-Nagumo (FitzHugh, 1961;Nagumo, Arimoto, & Yoshizawa, 1962) models (more precisely it is not a constant but it depends on the dynamical variables).Recent experiments (Naundorf, Wolf, & Volgushev, 2006;McCormick, Shu, & Yu, 2007;Naundorf, Wolf, & Volgushev, 2007) even suggest that there may be no real potential threshold.
times of neuron k, for the trajectory4 V, by: where The firing of neuron k corresponds to the following procedure.If V k (t) ≥ θ then neuron membrane potential is reset instantaneously to some constant reset value V reset and a spike is emitted toward post-synaptic neurons.In mathematical terms firing reads5 : where is called the "reset potential".In the sequel we assume, without loss of generality, that V reset = 0.This reset has a dramatic effect.Changing the initial values of the membrane potential, one may expect some variability in the evolution.Now, fix a neuron k and assume that there is a time t > 0 and an interval [a, b] Then, after reset, this interval is mapped to the point V reset .Then, all trajectories born from [a, b] collapse on the same point and have obviously the same further evolution.Moreover, after reset, the membrane potential evolution does not depend on its past value.This induces an interesting property used in all the Integrate and Fire models that we know (see e.g.(Gerstner & Kistler, 2002b)).The dynamical evolution is essentially determined by the firing times of the neurons, instead of their membrane potential value.
Below the threshold, V k < θ, neuron k's dynamics is driven by an equation of form: where C is the membrane capacity of neuron k.Without loss of generality we normalize the quantities and fix C = 1.In its most general form, the neuron k's membrane conductance g k > 0 depends on V k (see e.g.Hodgkin-Huxley equations (Hodgkin & Huxley, 1952)) and time t, while the current i k can also depend on V, the membrane potential vector, on time t, and also on the collection of past firing times.The current i k can include various phenomenological terms.Note that (3) deals with neurons considered as points instead of spatially extended objects.
Let us give two examples investigated in this paper.
(i) The leaky integrate and fire model.
In its simplest form equation (3) reads: where g k is a constant, and τ k = g k C is the characteristic time for membrane potential decay when no current is present.This model has been introduced in (Lapicque, 1907).
(ii) Conductance based models with α profiles More generally, conductance and currents depend on V only via the previous firing times of the neurons (Rudolph & Destexhe, 2006).Namely, conductances (and currents) have the general form 6 , ) where t (n) j is the n-th firing time of neuron j and t is the list of firing times of all neurons up to time t.This corresponds to the fact that the occurrence of a post-synaptic potential on synapse j, at time t (n) j , results in a change of the conductance g k of neuron k.As an example, we consider models of form: where the first term in the r.h.s. is a leak term, and where the synaptic current reads: ) , where E ± are reversal potential (typically E + ≃ 0mV and E − ≃ −75mV ) and where: In this equation, M j (t, V) is the number 7 of times neuron j has fired at time t.G ± kj is the synaptic efficiency (or synaptic weight) of the synapse j → k.(It is zero if there is no synapse j → k), where + [−] expresses that synapse j → k is excitatory [inhibitory].The α function mimics the conductance time-course after the arrival of a post-synaptic potential.A possible choice is: 6 The rather cumbersome notation g k (t, t ) simply expresses that in conductance based models the conductance depends on the whole set (history) of (past) firing times.Note that membrane potentials are reset after neuron firing, but not neuron conductances.
7 Henceforth, one assumes that there are finitely many spikes within a finite time interval.For continuous time dynamics, this fact is not guaranteed when neglecting the refractory period.Note also that this number, as well as the list t (n) j t , depends on the initial condition V and a small change in the initial condition may induce a drastic change of M j (t, V) at time t, as discussed later.This effect is sometimes disregarded (Coombes, 1999b).This issue has also been discussed (for current based IF-like models) as "phase history functions" in (Ashwin & Timme, 2005;Broer, Efstathiou, & Subramanian, 2008) (we thank one of the reviewers for this remark).
with H the Heaviside function and τ ± being characteristic times.This synaptic profile, with α(0) = 0 while α(t) is maximal for t = τ , allows us to smoothly delay the spike action on the post-synaptic neuron.We are going to neglect other forms of delays in the sequel.
Then, we may write (5) in the form (3) with: and: 1.2 Discrete time dynamics.
1.2.1 Characteristic time scales in neurons dynamics.
Integrate and Fire models assume an instantaneous reset of the membrane potential corollary to an infinite precision for the spike time.We would like to discuss shortly this aspect.Looking at the spike shape reveals some natural time scales: the spike duration τ (a few ms); the refractory period r ≃ 1ms; and the spike time precision.Indeed, one can mathematically define the spike time as the time where the action potential reaches some value (a threshold, or the maximum of the membrane potential during the spike), but, on practical ground, spike time is not determined with an infinite precision.An immediate conclusion is that it is not correct, from an operational point of view, to speak about the "spike time", unless one precises that this time is known with a finite precision δτ .Thus the notion of list of firing time t (n) j t used in section 1.1, must be revisited, and a related question is "what is the effect of this indeterminacy on the dynamical evolution ?".Note that this (evident ?) fact is forgotten when modeling e.g.spike with Dirac distributions.This is harmless as soon as the characteristic time δτ is smaller than all other characteristic times involved in the neural network.This is essentially true in biological networks but it is not true in Integrate and Fire models.These time scales arise when considering experimental data on spikes.When dealing with models, where membrane potential dynamics is represented by ordinary differential equations usually derived from Hodgkin-Huxley model, other implicit times scales must be considered.Indeed, Hodgkin-Huxley formulation in term of ionic channel activity assumes an integration over a time scale dt which has to be (i) quite larger than the characteristic time scale τ P of opening/closing of the channels, ensuring that the notion of probability as a meaning; (ii) quite larger than the correlation time τ C between channel states ensuring that the Markov approximation used in the equations of the variable m, n, h is legal.This means that, although the mathematical definition of d dt assumes a limit dt → 0, there is a time scale below which the ordinary differential equations lose their meaning.Actually, the mere notion of "membrane potential" already assumes an average over microscopic time and space scales.Note that the same is true for all differential equations in physics !But this (evident ?) fact is sometimes forgotten when dealing with Integrate and Fire models.Indeed, to summarize, the range of validity of an ODE modeling membrane potential dynamics is max(τ C , τ P ) << dt << δτ << τ .But the notion of instantaneous reset implies τ = 0 and the mere notion of spike time implies that δτ = 0 !!There is a last time scale related to the notion of raster plot.It is widely admitted that the "neural code" is contained in the spike trains.Spike trains are represented by raster plots, namely bi-dimensional diagrams with time on abscissa and some neurons labeling on ordinate.If neuron k fires a spike "at time t k " one represents a vertical bar at the point (t k , k).Beyond the discussion above on the spike time precision, the physical measurement of a raster plot involves a time discretization corresponding to the time resolution δ A of the apparatus.When observing a set of neurons activity, this introduces an apparent synchronization, since neurons firing between t and t + δ A will be considered as firing simultaneously.This raises several deep questions.In such circumstances the "information" contained in the observed raster plot depends on the time resolution δ A (Golomb., Hertz, Panzeri, Treves, & Richmond, 1997;Panzeri & Treves, 1996) and it should increase as δ A decreases.But is there a limit time resolution below which this information does not grow anymore ?In Integrate and Fire models this limit is δ A = 0.This may lead to the conclusion that neural networks have an unbounded information capacity.But is this a property of real neurons or only of Integrate and Fire models ?
The observation of raster plots corresponds to switching from the continuous time dynamics of membrane potential to the discrete time and synchronous dynamics of spike trains.One obtains then, in some sense, a new dynamical system, of symbolic type, where variables are bits ('0' for no spike, and '1' otherwise).The main advantage of this new dynamical system is that it focuses on the relevant variables as far as information and neural coding is concerned i.e. one focuses on spikes dynamics instead of membrane potentials.In particular, membrane potentials may still depend continuously on time, but one is only interested in their values at the times corresponding to the time grid imposed by the raster plot measurement.In some sense this produces a stroboscopic dynamical system, with a frequency given by the time resolution δ A , producing a phenomenological representation of the underlying continuous time evolution.
This has several advantages.(i) this simplifies the mathematical analysis of the dynamics avoiding the use of delta distributions, left-right limits, etc ... appearing in the continuous version; (ii) provided that mathematical results do not depend on the finite time discretization scale, one can take it arbitrary small; (iii) it enhances the role of symbolic coding and raster plots.
Henceforth, from now on, we fix a positive time scale δ > 0 which can be mathematically arbitrary small, such that (i) a neuron can fire at most once between [t, t + δ[ (i.e.δ << r, the refractory period); (ii) dt << δ, so that we can keep the continuous time evolution of membrane potentials (3), taking into account time scales smaller than δ, and integrating membrane potential dynamics on the intervals [t, t + δ[; (iii) the spike time is known within a precision δ.Therefore, the terminology, "neuron k fires at time t" has to be replaced by "neuron k fires between t and t + δ"; (iv) the update of conductances is made at times multiples8 of δ.

Raster plot.
In this context, we introduce a notion of "raster plot" which is essentially the same as in biological measurements.A raster plot is a sequence and is 0 otherwise.Note however that for mathematical reasons, explained later on, a raster plot corresponds to the list of firing states {ω(t)} ∞ t=0 over an infinite time horizon, while on practical grounds one always considers bounded times.Now, for each k = 1 . . .N , one can decompose the interval This splitting induces a partition P of M, that we call the "natural partition".The elements of P have the following form.
This is a N dimensional vector with binary components 0, 1.We call such a vector a firing state.Then M = ω∈Λ M ω where: Therefore the partition P corresponds to classifying the membrane potential vectors according to their firing state.Indeed, to each point V(t) of the trajectory Ṽ corresponds a firing state ω(t) whose components are given by: where Z is defined by : where χ is the indicator function that will later on allows us to include the firing condition in the evolution equation of the membrane potential (see ( 20)).On a more fundamental ground, the introduction of raster plots leads to a switch from the dynamical description of neurons, in terms of their membrane potential evolution, to a description in terms of spike trains where ω provides a natural "neural code".From the dynamical systems point of view, it introduces formally a symbolic coding and symbolic sequences are easier to handle than continuous variables, in many aspects such as the computation of topological or measure theoretic quantities like topological or Kolmogorov-Sinai entropy (Katok & Hasselblatt, 1998).A natural related question is whether there is a one-to-one correspondence between the membrane potential trajectory and the raster plot (see theorem 2).Note that in the deterministic models that we consider here, the evolution, including the firing times of the neurons and the raster plot, is entirely determined by the initial conditions.Therefore, there is no need to introduce an exogenous process (e.g.stochastic) for the generation of spikes (see (Destexhe & Contreras, 2006) for a nice discussion on these aspects).
Furthermore, this definition has a fundamental consequence In the present context, current and conductances at time t become functions of the raster plot up to time t.Indeed, we may write (7) in the form: where [ω] t = {ω(s)} t s=0 is the raster plot up to time t and M j (t, ω) is the number of spikes emitted by neuron j up to time t, in the raster plot ω (i.e.M j (t, ω) = t s=1 ω j (s)).But now t (n) j is a multiple of δ.
Remark.In continuous time Integrate and Fire models ω can assume uncountably many values.This is because a neuron can fire at any time and because firing is instantaneous.Therefore, the same property holds also if one considers sequences of firing states over a bounded time horizon.This is still the case even if one introduces a refractory period, because even if spikes produced by a given neuron are separated by a time slot larger or equal than the refractory period, the first spike can occur at any time (with an infinite precision).If, on the opposite, we discretize time with a time scale δ small enough to ensure that each neuron can fire only once between t and t + δ, ω, truncated at time T δ can take at most 2 N T values.For these reasons, the "computational power" of Integrate and Fire models with continuous time is sometimes considered as infinitely larger than a system with clocked based discretization (Maass & Bishop, 2003).The question is however whether this computational power is something that real neurons have, or if we are dealing with a model-induced property.

Integrate regime.
For this regime, as we already mentioned, we keep the possibility to have a continuous time (dt << δ) evolution of membrane potential (3).This allows us to integrate V on time scales smaller than δ.But, since conductances and currents depends now on the raster plot ω, we may now write (3) in the form: When neuron k does not fire between t, t + δ one has, integrating the membrane potential on the interval t, t + δ (see appendix): where: and where: is the integrated current with: Remarks.
1.In the sequel, we assume that the external current (see ( 8)) is time-constant.Further developments with a time dependent current, i.e. in the framework of an input-output computation (Bertschinger & Natschläger, 2004), will be considered next.
2. We note the following property, central in the subsequent developments.Since 1.2.4 Firing regime.
Let us now consider the case where neuron k fires between t and t + δ.In classical IF models this means that there is some t . This corresponds to Fig. 1b.Doing this one makes some error compared to the real spike shape depicted in Fig. 1a.In our approach, one does not know exactly when firing occurs but we use the approximation that the spike is taken into account at time t + δ.This means that we integrate V k until t + δ then reset it.In this scheme V k can be larger than θ as well.This explains why Z(x) = χ(x ≥ θ).This procedures corresponds to Fig. 1c (alternative I).One can also use a slightly different procedure.We reset the membrane potential at t + δ but we add to its value the integrated current between [t, t + δ[.This corresponds to Fig. 1d (alternative II).We have therefore three types of approximation for the real spike in Fig. 1a.Another one was proposed by (Hansel et al., 1998), using a linear interpolation scheme.Other schemes could be proposed as well.One expects them to be all equivalent when δ → 0. For finite δ, the question whether the error induced by these approximations is crucial and discussed in section 2.6.
In this paper we shall concentrate on alternative II though mathematical results can be extended to alternative I in a straightforward way.This corresponds to the initial choice of the Beslon-Mazet-Soula Model motivating the paper (Soula et al., 2006) and the present work.
In this case, the reset corresponds to : (recall that V reset = 0).Integrate and Fire regime can now be included in a unique equation, using the function Z defined in (11): where we set δ = 1 from now on.
Figure 1: A. "Real spike" shape; the sampling window is represented at a scale corresponding to a "small" sampling rate to enhance the related bias.B. Spike shape for an integrate and fire model with instantaneous reset, the real shape is in blue.C. Spike shape when reset occurs at time t + δ (alternative I).D. Spike shape with reset at time t + δ plus addition of the integrate current (green curve) (alternative II).
Consider the leaky integrate and fire model, where conductances are constant.Set ) for excitatory (inhibitory) synapses.Then, replacing the α-profile by a Dirac distribution, (20) reduces to: This model has been proposed by G. Beslon, O. Mazet and H. Soula in (Soula et al., 2006).A mathematical analysis of its asymptotic dynamics has been done in (Cessac, 2008) and we extend these results to the more delicate case of conductance based IF models in the present paper.(Note that having constant conductances leads to a dynamics which is independent of the past firing times (raster plot).In fact, the dynamical system is essentially a cellular automaton but with a highly non trivial dynamics).

Alpha-profile conductances.
Consider now a conductance based model of form (3), leading to: where K is a constant: while, using the form (6) for α gives: One has therefore to handle an exponential of an exponential.This example illustrates one of the main problem in IF models.IF models have been introduced to simplify neurons description and to simplify numerical calculations (compared e.g. with Hodgkin-Huxley's model (Hodgkin & Huxley, 1952)).Indeed, their structure allows one to write an explicit expression for the next firing times of each neurons, knowing the membrane potential value.However, in case of α exponential profile, there is no simple form for the integral and, even in the case of one neuron, one has to use approximations with Γ functions (Rudolph & Destexhe, 2006) which reduce consequently the interest of IF models and event based integration schemes.
2 Theoretical analysis of the dynamics.
2.1 The general picture.
In this section we develop in words some important mathematical aspects of the dynamical system (20), mathematically proved in the sequel.
Singularity set.The first important property is that the dynamics (20) (and the dynamics of continuous time IF models as well) is not smooth, but has singularities, due to the sharp threshold definition in neurons firing.The singularity set is: This is the set of membrane potential vectors such that at least one of the neurons has a membrane potential exactly equal to the threshold 9 .This set has a simple structure: it is a finite union of N − 1 dimensional hyperplanes.Although S is a "small" set both from the topological (non residual set) and metric (zero Lebesgue measure) point of view, it has an important effect on the dynamics.
Local contraction.The other important aspect is that the dynamics is locally contracting, because γ k (t, [ω] t ) < 1 (see eq. ( 18)).This has the following consequence.Let us consider the trajectory of a point V ∈ M and perturbations with an amplitude < ǫ about V (this can be some fluctuation in the current, or some additional noise, but it can also be some error due to a numerical implementation).Equivalently, consider the evolution of the ǫ-ball B(V, ǫ).If B(V, ǫ) ∩ S = ∅ then we shall see in the next section that the image of B(V, ǫ) is a ball with a smaller diameter.This means, that, under the condition B(V, ǫ) ∩ S = ∅, a perturbation is damped.Now, if the images of the ball under the dynamics never intersect S, any ǫ-perturbation around V is exponentially damped and the perturbed trajectories about V become asymptotically indistinguishable from the trajectory of V. Actually, there is a more dramatic effect.If all neurons have fired after a finite time t then all perturbed trajectories collapse onto the trajectory of V after t + 1 iterations (see prop. 1 below).
Initial conditions sensitivity.On the opposite, assume that there is a time, t 0 , such that the image of the ball B(V, ǫ) intersects S. By definition, this means that there exists a subset of neurons {i 1 , . . ., i k } and For example, some neuron does not fire when not perturbed but the application of an ǫ-perturbation induces it to fire (possibly with a membrane potential strictly above the threshold).This requires obviously this neuron to be close enough to the threshold.Clearly, the evolution of the unperturbed and perturbed trajectory may then become drastically different (see Fig. 2).Indeed, even if only one neuron is lead to fire when perturbed, it may induce other neurons to fire at the next time step, etc . . ., inducing an avalanche phenomenon leading to unpredictability and initial condition sensitivity 10 .
It is tempting to call this behavior "chaos", but there is an important difference with the usual notion of chaos in differentiable systems.In the present case, due to the sharp condition defining the threshold, initial condition only occurs at sporadic instants, whenever some neuron is close enough to the threshold.Indeed, in certain periods of time the membrane potential typically is quite far below threshold, so that the neuron can fire only if it receives strong excitatory input over a short period of time.It shows then a behavior that is robust against fluctuations.On the other hand, when membrane potential is close to the threshold a small perturbation may induce drastic change in the evolution.
Stability with respect to small perturbations.Therefore, depending on parameters such as the synaptic weights, the external current, it may happen that, in the stationary regime, the typical trajectories stay away from the singularity set, say within a distance larger than ǫ > 0, which can be 9 A sufficient condition for a neuron i to fire at time t is V i (t) = θ hence V(t) ∈ S.But this is not a necessary condition.Indeed, as pointed in the footnote 5, there may exist discontinuous jumps in the dynamics, even if time is continuous, either due to noise, or α profiles with jumps (e.g.α(t) = 1 τ e − t τ , t ≥ 0).Thus neuron i can fire with V i (t) > θ and V(t) / ∈ S. In the present case, this situation arises because time is discrete and one can have V (t − δ) < θ and V (t) > θ.This holds as well even if one uses numerical schemes using interpolations to locate more precisely the spike time (Hansel et al., 1998).
arbitrary small, (but positive).Thus, a small perturbation (smaller than ǫ) does not produce any change in the evolution.At a computational level, this robustness leads to stable input-output transformations.In this case, as we shall see, the dynamics of ( 20) is asymptotically periodic (but there may exist a large number of possible orbits, with a large period).In this situation the system has a vanishing entropy11 .This statement is made rigorous in theorem 1 below.
On the other hand, if the distance between the set where the asymptotic dynamics lives12 and the singularity set is arbitrary small then the dynamics exhibit initial conditions sensitivity, and chaos.Thus, a natural question is: is chaos a generic situation ?How does this depend on the parameters ?A related question is: how does the numerical errors induced by a time discretization scheme evolve under dynamics (Hansel et al., 1998) ?
"Edge of chaos".It has been shown, in (Cessac, 2008) for the BMS model, that there is a sharp transition13 from fixed point to complex dynamics, when crossing a critical manifold usually called the "edge of chaos" in the literature.While this notion is usually not sharply defined in the Neural Network literature, we shall give a mathematical definition which is moreover tractable numerically.Strictly speaking chaos only occurs on this manifold, but from a practical point of view, the dynamics is indistinguishable14 from chaos, close to this manifold.When the distance to the edge of chaos further increases the dynamics is periodic with typical periods compatible with simulation times.This manifold can be characterized in the case where the synaptic weights are independent, identically distributed with a variance σ 2 N .In BMS model (e.g., time discretized gIF model with constant conductances) it can be proved that the chaotic situation is non generic (Cessac, 2008).We now develop the same lines of investigation and discuss how these result extend to the model (20).Especially, the "edge of chaos" is numerically computed and compared to the BMS situation.
Let us now switch to the related mathematical results.Proofs are given in the appendix.

Piecewise affine map.
Let us first return to the notion of raster plot developed in section 1.2.2.At time t, the firing state ω(t) ∈ Λ can take at most 2 N values.Thus, the list of firing states ω(0), . . ., ω(t) ∈ Λ t+1 can take at most 2 N (t+1) values.(In fact, as discussed below, only a subset of these possibilities is selected by the dynamics).This list is the raster plot up to time t and we have denoted by [ω] t .Once the raster plot up to time t has been fixed the coefficients γ k and the integrated currents J k in (20) are determined.
Fixing the raster plot up to time t amounts to construct branches for the discrete flow of the dynamics, corresponding to sub-domains of M constructed iteratively, via the natural partition (9), in the following way.
Fix t > 0 and [ω] t .Note: This is the set (possibly empty) of initial membrane potentials vectors V ≡ V(0) whose firing pattern at time s is ω(s), s = 0 . . .t. Consequently, ∀V ∈ M [ω] t , we have: as easily found by recursion on (20).We used the convention Then, define the map: with V k (t + 1) given by ( 25) and F 0 ω = Id.Note that F t+1 ω is affine.Finally define: such that the restriction of F t+1 to the domain M [ω] t is precisely F t+1 ω .This mapping is the flow of the model ( 20) where: A central property of this map is that it is piecewise affine and it has at most 2 N (t+1) branches F t+1 ω parametrized by the legal sequences [ω] t which parametrize the possible histories of the conductance/currents up to time t.
Let us give a bit of explanation of this construction.Take V ≡ V(0) ∈ M ω(0) .This amounts to fixing the firing pattern at time 0 with the relation , where γ k , J k do not depend on V(0) but only on the spike state of neurons at time 0. Therefore, the mapping is affine (and continuous on the interior of M ω(0) ).Since ω(0) is an hypercube, F 1 ω M ω(0) is a convex connected domain.This domain typically intersects several domains of the natural partition P.This corresponds to the following situation.Though the pattern of neuron firing at time 0 is fixed as soon as V(0) ∈ M ω(0) , the list of neurons firing at the next time depends on the value of the membrane potentials V(0), and not only on the spiking pattern at time 0. But, by definition, the domain: 1)ω(0) , the spiking pattern at time 0 is ω(0) and it is ω(1) at time 1.If the intersection is empty this simply means that one cannot find a membrane potential vector such that neurons fire according to the spiking pattern ω(0) at time t = 0 then fire according to the spiking pattern ω(1) at time t = 1.If the intersection is not empty we say that "the transition ω(0) → ω(1) is legal15 ".
Proceeding recursively as above one constructs a hierarchy of domains Figure 2: Schematic representation, for two neurons, of the natural partition P and the mapping discussed in the text.In this case, a firing state is a vector with components ω = ω1 ω 2 labeling the partition elements.A set of initial conditions, say a small (L ∞ ) ball in M ω , is contracted by leak (neuron 1 in the example) and reset (neuron 2 in the example), but its image can intersect the singularity set.This generates several branches of trajectories.Note that we have given some width to the projection of the image of the ball on direction 2 in order to see it on the picture.But since neuron 2 fires the width is in fact 0.
As stated before, M [ω] t is the set of membrane potential vectors V such that the firing patterns up to time t are ω(0), . . ., ω(t).If this set is non empty we say that the sequence ω(0), . . ., ω(t) is legal.Though there are at most 2 N (t+1) possible raster plots at time t the number of legal raster plots is typically smaller.This number can increase either exponentially with t or slower.We shall denote by Σ + Λ the set of all legal (infinite) raster plots (legal infinite sequences of firing states).Note that Σ + Λ is a topological space for the product topology generated by cylinder sets (Katok & Hasselblatt, 1998).The set [ω] t of raster plots having the same first t + 1 firing patterns is a cylinder set.

Phase space contraction.
Now, we have the following: is affine, with a Jacobian matrix and an affine constant depending on t, [ω] t .Moreover, the Jacobian matrix is diagonal with eigenvalues Proof The proof results directly from the definition ( 26) and ( 25) with γ k (s, [ω] s ) < 1, ∀s ≥ 0 (see ( 18) ).
Since the domains M ω of the natural partition are convex and connected, and since F is affine on each domain (therefore continuous on its interior), there is a straightforward corollary:

Corrollary 1 The domains M [ω] t are convex and connected.
There is a more important, but still straightforward consequence: Corrollary 2 F t+1 is a non uniform contraction on M where the contraction rate in direction k is Then, we have the following : That is, all orbits born from the domain M [ω] t converge to the same orbit in a finite time.

If
Proof Statement 1 holds since, under these hypothesis, all eigenvalues of F t+1 ω are zero.For 2, since DF t+1 ω is diagonal, the Lyapunov exponent in direction k is defined by λ k (ω) = lim t→∞ 1 t+1 t s=0 log(σ k (t, ω)) whenever the limit exists (it exists almost surely for any F invariant measure from Birkoff theorem).
Remark.An alternative definition of Lyapunov exponent has been introduced by Coombes in (Coombes, 1999b) and (Coombes, 1999a), for Integrate and Fire neurons.His definition, based on ideas developed for impact oscillators (Muller, 1995), takes care of the discontinuity in the trajectories arising when crossing S. Unfortunately, his explicit computation at the network level (with continuous time dynamics), makes several implicit assumptions (see eq. ( 6) in (Coombes, 1999a)): (i) there is a finite number of spikes within a bounded time interval; (ii) the number of spikes that have been fired up to time t, ∀t > 0, is the same for the mother trajectory and for a daughter trajectory, generated by a small perturbation of the mother trajectory at t = 0; (iii) call T k i , in Coombes' notations, the k-th spike time for neuron i in the mother trajectory, and T k i the k-th spike time for neuron i in the daughter trajectory.Then T k i = T k i + δ k i , where δ k i is assumed to become arbitrary small, ∀k ≥ 0, when the initial perturbation amplitude tends to zero.While assumption (i) can be easily fulfilled (e.g. by adding a refractory period) assumptions (ii) and (iii) are more delicate.
Transposing this computation to the present analysis, this requires that both trajectories are never separated by the singularity set.A sufficient condition is that the mother trajectory stays sufficiently far from the singularity set.In this case the Lyapunov exponent defined by Coombes coincides with our definition and it is negative.On the other hand, in the "chaotic" situation (see section 2.4.3),assumptions (ii) and (iii) can typically fail.For example, it is possible that neuron i stops firing after a certain time, in the daughter trajectory, while it was firing in the mother trajectory, and this can happen even if the perturbation is arbitrary small.This essentially means that the explicit formula for the Lyapunov exponent proposed in (Coombes, 1999a) cannot be applied as well in the "chaotic" regime.

Asymptotic dynamics.
2.4.1 Attracting set A and ω-limit set.
The main notion that we shall be interested in from now on concerns the invariant set where the asymptotic dynamics lives.
Definition 1 (From (Guckenheimer & Holmes, 1983;Katok & Hasselblatt, 1998)) A point y ∈ M is called an ω-limit point for a point x ∈ M if there exists a sequence of times {t k } +∞ k=0 , such that x(t k ) → y as t k → +∞.The ω-limit set of x, ω(x), is the set of all ω-limit points of x.The ω-limit set of M, denoted by Ω, is the set Ω = x∈M ω(x).
Equivalently, since M is closed and invariant, we have Ω = ∞ t=0 F t (M).Basically, Ω is the union of attractors.But for technical reasons, related to the case considered in section 2.4.3, it is more convenient to use the notion of ω-limit set.
The proof is given in the appendix A.2.
Note that conditions (1) and ( 2) are not disjoint.The meaning of these conditions is the following.(1) corresponds to imposing that a neuron has fired after a finite time T (uniformly bounded, i.e. independent of V and i).(2) amounts to requiring that after a certain time t 0 , the membrane potential stays below the threshold value and it cannot accumulate on θ.We essentially want to avoid a situation when a neuron can fire for the first time after an unbounded time (see section 2.4.3 for a discussion of this case).Thus assumptions (1),(2) look quite reasonable.Under these assumptions the asymptotic dynamics is periodic and one can predict the evolution after observing the system on a finite time horizon T , whatever the initial condition.Note however that T can be quite a bit large.
There is a remarkable corollary result, somehow hidden in the proof given in the appendix.The neurons that do not fire after a finite time are still driven by the dynamics of firing neurons.It results that, in the asymptotic regime, non firing neurons have a membrane potential which oscillates below the threshold.This exactly corresponds to what people call "sub-threshold oscillations".In particular, there are times where those membrane potentials are very close to the threshold, and a small perturbation can completely changes further evolution.This issue is developed in the next section.Possible biological interpretations are presented in the discussion section.

Ghost orbits.
The advantage of the previous theorem is that we define conditions where one can relatively easily controls dynamics.However, what happens if if we consider the complementary situation corresponding to the following definition ?
Definition 2 An orbit Ṽ is a ghost orbit if ∃i such that: Namely there exists at least one initial condition V and one neuron i such that one cannot uniformly bound the first time of firing, and V i (t) approaches arbitrary close the threshold.In other words subthreshold oscillations drive the neuron "dangerously close" to the threshold, though we are not able to predict when the neuron will fire.This may look a "strange" case from a practical point of view, but it has deep implications.This indeed means that we can observe the dynamics on arbitrary long times without being able to predict what will happen later on, because when this neuron eventually fire, it may drastically change the evolution.This case is exactly related to the chaotic or unpredictable regime of IF models.
One may wonder whether the existence of ghost orbits is "generic".To reply to this question one has first to give a definition of genericity.In the present context, it is natural to consider the dynamical system describing the time evolution of our neural network as a point in a space H of parameters.These parameters can be e.g., the synaptic weights, or parameters fixing the time scales, the reversal potentials, the external currents, etc... Varying these parameters (i.e., moving the point representing our dynamical system in H) can have two possible effects.Either there is no qualitative change in the dynamics and observable quantities such as e.g., firing rates, average inter-spikes interval, etc, are varying continuously.Or, a sharp change (bifurcation) occurs.This corresponds to the crossing of a critical or bifurcation manifold in H. Now, a behavior is generic if it is "typical".On mathematical grounds this can have two meanings.Either this behavior is obtained, with a positive probability, when drawing the parameters (the corresponding point in H) at random with some natural probability (e.g., Gaussian).In this case one speaks of "metric genericity".Or, this behavior holds in a dense subset of H.One speaks then of "topological genericity".The two notions usually do not coincide.
In the BMS model, ghost orbits are non generic in both senses (Cessac, 2008).The proof does not extend to more general models such as (20) because it heavily uses the fact that the synaptic current takes only finitely many values in the BMS model.As soon as we introduce a dependence in ω this is not the case anymore.We don't know yet how to extend this proof.

Edge of chaos.
On practical grounds ghost orbits involve a notion of limit t → +∞ which has no empirical meaning.Therefore the right question is: are there situations where a neuron can fire for the first time after a time which is well beyond the observation time ?One way to analyze this effect is to consider how close the neurons are to the threshold in their evolution.On mathematical grounds this is given by the distance from the singularity set to the ω-limit set: The advantage of this definition, is that it can easily be adapted to the plausible case where observation time is bounded (see section 3).Now, the following theorem holds.
1.If d(Ω, S) > 0 then Ω is composed by finitely many periodic orbits with a finite period.
2. There is a one-to-one correspondence between a trajectory on Ω and its raster plot.
The proof is exactly the same as in (Cessac, 2008) so we don't reproduce it here.It uses the fact that if d(Ω, S) > 0 then there is a finite time T , depending on d(Ω, S) and diverging as d(Ω, S) → 0, such that F T has a Markov partition (constituted by local stable manifolds since dynamics is contracting) where the elements of the partition are the domains M [ω] T .Note however that d(Ω, S) > 0 is a sufficient but not a necessary condition to have a periodic dynamics.In particular, according to theorem 1 one can have d(Ω, S) = 0 and still have a periodic dynamics, if at some finite time t, for some neuron i, V i (t) = θ.This strict equality is however not structurally stable, since a slight change e.g. in the parameters will remove it.The main role of the condition d(Ω, S) > 0 is therefore to avoid situations where the membrane potential of some neuron accumulates on θ from below (ghost orbits).See the discussion section for a possible biological interpretation on this.
But d(Ω, S) plays another important role concerning the effect of perturbations on the dynamics.Indeed, if d(Ω, S) > 0 then the contraction property (corollary 2) implies that any perturbation smaller than d(Ω, S) will be damped by dynamics.In particular, making a small error on the spike time, or any other type of error leading to an indeterminacy of V smaller than d(Ω, S) will be harmless.
Let us now discuss the second item of theorem 2. It expresses that the raster plot is a symbolic coding for the membrane potential trajectory.In other words there is no loss of information on the dynamics when switching from the membrane potential description to the raster plot description.This is not true anymore if d(Ω, S) = 0.
The first item tells us that dynamics is periodic, but period can be arbitrary long.Indeed, following (Cessac, 2008) an estimate for an upper bound on the orbits period is given by: where < γ > denotes the value of γ averaged over time and initial conditions16 (see appendix for details).Though this is only an upper bound this suggests that periods diverge when d(Ω, S) → 0. In BMS model, this is consistent with the fact that when d(Ω, S) is close to 0 dynamics "looks chaotic".There- fore, d(Ω, S) is what a physicist could call an "order parameter", quantifying somehow the dynamics complexity.The distance d(Ω, S) can be numerically estimated as done in ( 33) and (34), section 3. Before, we need the following list of (operational) definitions.

Definition 3 (Edge of chaos)
The edge of chaos is the set of points E in the parameter space H where d(Ω, S) = 0.
The topological structure of E can be quite complicated as we checked in the simplest examples (e.g., the BMS model with Laplacian coupling) and suggested by the papers (Bressloff & Coombes, 2000b, 2000a) (see discussion).There are good reasons to believe that this set coincides with the set of points where the entropy is positive (see (Kruglikov & Rypdal, 2006a, 2006b) and discussion below).The set of points where the entropy is positive can have a fractal structure even in the simplest examples of one dimensional maps (Mackay & Tresser, 1986;Gambaudo & Tresser, 1988).Therefore, there is no hope to characterize E rigorously in a next future.Instead, we use below a numerical characterization.
The edge of chaos is a non generic set in the BMS model, and the same could hold as well in model ( 20).Nevertheless, it has a strong influence on the dynamics, since crossing it leads to drastic dynamical changes.Moreover, close to E dynamics can be operationally indistinguishable from chaos.More precisely, let us now propose another definition.

Definition 4 (Effective entropy)
Fix T o called "the observational time".This is the largest accessible duration of an experiment.Call n(t) the number of (legal) truncated raster plots up to time t.Then, the effective entropy is; Note that in the cases where raster plots provide a symbolic coding for the dynamics then lim To→∞ h (ef f ) = h (top) , the topological entropy.
On practical grounds, this definition corresponds to the following notion.The larger the effective entropy, the more the system is able to produce distinct neural codes.This provides one way to measure the "complexity" of the dynamics.On more "neuronal" grounds this quantity measures the variability in the dynamical response of the neuronal network to a given stimulus (external current) or its ability to produce distinct "functions" (a function being the response to a stimulus in terms of a spikes train).Thinking of learning mechanisms (e.g., Hebbian) and synaptic plasticity (LTD,LTP,STDP) one may expect to having the largest learning capacities when this entropy is large.This aspect will be developed in a separated paper.(For the effect of Hebbian learning and entropy reduction in firing rate neural networks see (Siri, Berry, Cessac, Delord, & Quoy, 2008)).
Finally, a positive effective entropy means that the system essentially behaves like a chaotic system during the time of the experiment.Indeed, the entropy is closely related to the distance d(Ω, S), since, from (30), a rough estimate/bound of h (ef f ) is easily obtained from (30), (31): The notion of effective entropy provides some notion of "width" to the edge of chaos E. For a fixed T o the system behaves chaotically in a thick region E To ⊃ E in H such that h (ef f ) > 0. And from (32) one expects that this entropy gets larger when d(Ω, S) gets smaller.

Effects of time discretization.
Under the light of the previous results, let us reconsider the approximation where spikes are taken into account at multiple of the characteristic time scale δ, for the conductances update.Doing this, we make some error in the computation of the membrane potential, compared to the value obtained when using the "exact" spike time value.Now, the question is whether this error will be amplified by the dynamics, or damped.As we saw, dynamics (20) is contracting but the effect of a small error can have dramatic consequences when approaching the singularity set.The distance d(Ω, S) provides a criterion to define a "safe" region where a small error of amplitude ǫ > 0 in the membrane potential value is harmless, basically, if ǫ < d(Ω, S).On the other hand, if we are in a region of the parameters space where d(Ω, S) = 0 then a slight perturbation has an incidence on the further evolution.Since δ can be arbitrary small in our theorems we have a good control on the dynamics of the continuous time IF models except at the edge of chaos where d(Ω, S) = 0.This enhances the question of mathematically characterizing this region in the parameter space H.Note indeed that numerical investigations are of little help here since we are looking for a parameter region where the distance d(Ω, S) defined as an asymptotic limit, has to be exactly zero.The problem is that even sophisticated schemes (e.g., event based) are also submitted to round off errors.Therefore, as a conclusion, it might well be that all numerical schemes designed to approximate continuous time Integrate and Fire models display trajectory sensitivity to spike time errors, when approaching d(Ω, S) = 0.
3 A numerical characterization of the "edge of chaos".
As mentioned in the previous section an exact analytic computation of the edge of chaos in the general case seems to be out of reach.However, a "coarse grained" characterization can be performed at the numerical level and possibly some analytic approximation could be obtained.For this, we choose the synaptic weights (resp.the synaptic conductances) (and/or the external currents) randomly, with some probability P W (P i (ext) ), where W is the matrix of synaptic weights (W ij = E ± G ± ij ) and i (ext) the vector of external currents (recall that external currents are time constant in this paper).A natural starting point is the use of Gaussian independent, identically distributed variables, where one varies the parameters, mean and variance, defining the probability definition (we call them statistical parameters, see (Cessac & Samuelides, 2007;Samuelides & Cessac, 2007) for further developments on this approach).Doing these, one performs a kind of fuzzy sampling of the parameters space, and one somehow expects the behavior observed for a given value of the statistical parameters to be characteristic of the region of W, i (ext) that the probabilities P W , P i (ext) weight (more precisely, one expects to observe a "prevalent" behavior in the sense of (Hunt, Sauer, & Yorke, 1992)).
The idea is then to estimate numerically d(Ω, S) in order to characterize how it varies when changing the statistical parameters.As an example, in the present paper, we select conductances (resp.synaptic weights) randomly with a Gaussian probability with a fixed mean and a variance σ 2 N , and we study the behavior of d(Ω, S) when σ is varying.Note that the neural network is almost surely fully connected.We compute numerically an approximation of the distance d(Ω, S), where we fix a transient time T r and an observation time T o and average over several initial conditions V (n) , n = 1 . . .N CI for a given sample of synaptic weights.Then we perform an average over several synaptic weights samples W m , m = 1 . . .N W .In a more compact form, we compute: where: In this way, we obtain a rough and coarse grained location of the edge of chaos where the distance d(Ω, S) vanishes.
We have performed the following experiments with two control parameters.
• Variance of random synaptic weights.We randomly select the synaptic strength which modulates the synaptic conductance using a Gaussian distribution so that 80% of the synapses are excitatory and 20% inhibitory.The average standard-variation σ is varied.The value σ = 0.5 corresponds, in our units, to what is chosen in the literature when considering the biological dispersion in the cortex (e.g., (Rudolph & Destexhe, 2006)).Note however that, at the present stage of investigation, Dale's principle is not taken into account.
• Membrane leak time-constant.As an additional control parameter we vary the membrane leak around the usual τ L = 1 • • • 20ms values.This choice is two-fold.The value τ L = 20ms corresponds to in-vitro measurement, while τ L → 1ms allows to represent in-vivo conditions in the cortex.On the other hand, it acts directly on the average contraction rate < γ > which is a natural control parameter.
Each simulation randomly selects the initial potential values in a −70 . . .30mV range.For each condition the simulation is run for N CI = 100 initial conditions and N W = 10 random selection of the synaptic weights.With a sampling period of 0.1ms, the network is run during T r = 1000 steps in order to "skip" the transients17 and then observed during T o = 1000 steps.In order to explore the role of history dependent conductances on the dynamics we considered different models from the biologically plausible IF model to BMS model.More precisely we explore four modeling variants: 1. Model I, defined in (20).20) with a fixed γ.The evolution equation of membrane potentials is thus given by:

Model II (
where the average γ is the value observed during the numerical simulation of model I.Note that γ depends on the parameters σ, τ L .The goal of this simulation is to check the role of the fluctuations of γ(t, [ω] t ), controlling the instantaneous contraction rate, compared to a mean-field model where γ(t, [ω] t ) is replaced by its average.This corresponds to what is called "current based" synapses instead of "conductance based" synapses in the literature, (see e.g.(Brette et al., 2007)).20) approximation with a fixed γ and simplified synapses.The evolution equation of membrane potentials is given by:

Model III (
In addition to the previous simplification, we consider so-called "current jump" synapses where the synaptic input simply corresponds to an impulse, added to the membrane potential equation.
Here the magnitude of the impulse and its delay δ − = 2ms and δ + = 10ms in order to keep both characteristics as closed as possible to the previous condition.
4. Model IV ( 20) with a fixed γ and instantaneous simplified synapses.The evolution equation of membrane potentials is given by: where in addition to the previous simplification, the delay has been suppressed and where This last condition strictly correspond to the original BMS model (Soula et al., 2006).
The results are given below.We have first represented the average value γ for the model I in the range σ ∈ [0.01, 1], τ L ∈ [10, 40]ms (see Fig. 3).The quantity related to the contraction rate, is remarkably constant (with small variations within the range [0.965, 0.995]).Then, we have considered the average value of the current J k (t, [ω] t ), averaged over time, initial conditions and neurons and denoted by I to alleviate the notations (Fig. 4), the logarithm of the distance d(Ω, S) (Fig. 5), and the average Inter Spike Interval (ISI-Fig.6), for the four models.The main observations are the following.Average current and Inter Spike Intervals have essentially the same form for all models.This means that these quantities are not really relevant if one wants to discriminate the various models in their dynamical complexity.The observation of the distance d(Ω, S) is quite more interesting.First, in the four models, the distance becomes very small when crossing some "critical region" in the plane τ L , γ.This region has a regular structure for the BMS model, but its structure seems more complex for (20).Note however that the numerical investigations used here do not allow us to really conclude on this point.The most remarkable fact is that, in models III and IV, the distance increases when σ increases beyond this region, while it does not in models I and II.This corresponds to the following observation.When the d(Ω, S) is small, one observes a complex dynamics with no apparent period.One naturally concludes to a chaotic regime.As we saw, strictly speaking it is in fact periodic but since periods are well beyond observable times, the situation is virtually chaotic 18 .When the distance increases, the orbits period decreases.
18 Moreover, it is likely that the phase space structure has some analogies with spin-glasses.For example, if γ = 0 the dynamics is essentially equivalent to the Kauffman's cellular automaton (Kauffman, 1969).It has been shown by Derrida and coworkers (Derrida & Flyvbjerg, 1986;Derrida & Pomeau, 1986) that the Kauffman's model has a structure similar to the Sherrington-Kirckpatrick spin-glass model (Mézard, Parisi, & Virasoro, 1987).The situation is even more complex when γ = 0.It is likely that we have in fact a situation very similar to discrete time neural networks with firing rates where a similar analogy has been exhibited (Cessac, 1994(Cessac, , 1995)).29 Therefore, there is a range of σ values where period become smaller than observational time and one concludes that dynamics is periodic.
The situation is different for model I,II since the distance does not apparently increases with σ.This suggests that introducing conductance based synapses and currents enhances considerably the width of the edge of chaos.On practical grounds, this means that model I,II have the capacity to display a very large number of distinct codes for wide choices of parameters.This is somewhat expected since the opposite conclusion would mean that introducing spike dependent conductances and current does not increases the complexity and information capacity of the system.But it is one thing to guess some behavior and another thing to measure it.Our investigations on the distance d(Ω, S), a concept based on the previous mathematical analysis, makes a step forward in this direction.
One step further, we have represented examples of raster plots in Fig. 8 and 9 for models I and IV.The figure 8 essentially illustrates the discussion above on the relation between the distance d(Ω, S) and the dynamics; for σ = 0.05, where d(Ω, S) is "large", and dynamics is periodic; and for σ = 0.4, where d(Ω, S) is small, and dynamics looks more chaotic, for the two models.The difference between the two models becomes more accurate as σ increases.Fig. 9 represents raster plots for models I,IV, with σ = 10, where we study the effect of a small amount of noise, of amplitude 10 −4 × θ in the external current.This has no effect on model IV while it changes slightly the raster plot for model I, as expected.There is another remarkable difference.The code is sparser for model I than for model IV.This suggests that model I is in some sense optimal with respect to coding since it is able to detect very small changes in an input but the changes is not drastic and the neural code remains very sparse.

Discussion
We have thus an operational definition for the "edge of chaos" where an "order parameter", the distance of orbits to the singularity has been defined.This parameter has a deep meaning.It controls how much the system is sensitive to perturbations.Such perturbations can be noise, but they can also be a small variation in the external current, corresponding e.g. to an input.If the amplitude of this perturbation is smaller than d(Ω, S) then it has no effect on the long term dynamics, and the neural code (raster plot) is unchanged.On the other hand, when the distance is small, even a tiny perturbation has a dramatic effect on the raster plot: the system produces a different code.As a corollary, the effective entropy is maximal when the distance is minimal.On practical ground, having a positive distance with a large effective entropy corresponds to situations where the system is able to produce a large number of distinct codes within the observational time, while this code is nevertheless robust to small perturbations of the input.Thus, we have a good compromise between the variability of the responses to distinct inputs and robustness of the code when an input is subject to small variations.
Several questions are now open.A first one concerns the way how we measured this distance.We used a random sampling with independent synaptic weights.But these weights are, in reality, highly correlated, via synaptic plasticity mechanism.What is the effect of e.g.STPD or Hebbian learning on the effective entropy is a perspective for a future work.Recent results in (Soula, 2005) and (Siri et al., 2008;Siri, Berry, Cessac, Delord, & Quoy, 2007) suggest that synaptic plasticity reduces the entropy by diminishing the variability of raster plots and increasing the robustness of the response to an input.Some general (variational) mechanism could be at work here.This aspect is under investigation.
Another important issue is the effect of noise.It is usual in neural network modeling to add Brownian The control parameter is τ L = 20ms.As visible in Fig. 5, σ = 0.05 corresponds to a small order dynamics where the periodic behavior is clearly visible, and σ = 0.40 to the "edge of chaos".One blob width is 1msec noise to the deterministic dynamics.This noise accounts for different effects such as the diffusion of neurotransmitters involved in the synaptic transmission, the degrees of freedom neglected by the model, external perturbations, etc ... Though it is not evident that the "real noise" is Brownian, using this kind of perturbations has the advantage of providing a tractable model where standard theorems in the theory of stochastic processes (Touboul & Faugeras, 2007) or methods in non equilibrium statistical physics (e.g.Fokker-Planck equations (Brunel & Hakim, 1999)) can be applied.Though we do not treat explicitly this case in the present work, the formalism has been designed to handle noise effects as well.As a matter of fact, the effect of Brownian noise on the dynamics of our model can be analyzed with standard techniques in probability theory and stochastic perturbations of dynamical systems (Freidlin & Wentzell, 1998).In particular, the probability distribution of the membrane potential trajectory can be obtained by using The noise is added to the membrane potential and its magnitude is very small (10 −4 × θ).One blob width is 1msec.Cessac, 2007).Noise have several effects.Firstly, the stochastic trajectories stay around the unperturbed orbits until they jump to another attraction basin, the characteristic time depending on the noise intensity ("Arrhenius law").This has the effect of rendering the dynamics uniquely ergodic, which somehow simplifies the statistical analysis.The effect of noise will be essentially prominent in the region where d(Ω, S) is small, leading to an effective initial condition sensitivity and an effective positive Lyapunov exponent, that could be computed using mean-field approaches (Cessac, 1995).It is possible to estimate the probability that a trajectory approaches the singularity set S within a finite time T and a distance d by using Freidlin-Wentzell estimates (Freidlin & Wentzell, 1998).One can also construct a Markov chain for the transition between the attraction basin of the periodic orbits of the unperturbed dynamics.
The overall picture could be very similar (at least for BMS model) to what happens when stochastically perturbing Kauffman's model (Golinelli & Derrida, 1989), with possibly a phase space structure rem-iniscent of spin-glasses (where noise plays the role of the temperature).This study is under investigations.
Yet another important issue relates to the fact that spikes can also be lost.This aspect is not yet taken into account in the present formalism, but annihilation of spikes is a future issue to address.
A final issue is the relation of this work with possible biological observations.We would like in particular to come back to the abstract notion of ghost orbit.As said in the text, this notion corresponds to situation where the membrane potential of some "vicious" neuron fluctuates below the threshold, and approaches it arbitrary close, with no possible anticipation of its first firing time.This leads to an effective unpredictability in the network evolution, since when this neuron eventually fire, it may drastically change the dynamics of the other neurons, and therefore the observation of the past evolution does not allow one to anticipate what will be the future.In some sense, the system is in sort of a metastable state but it is not in a stationary state.Now, the biological intuition tends to consider that a neuron cannot suddenly fire after a very long time, unless its input changes.This suggests therefore that "vicious" neurons are biologically implausible.However, this argument, to be correct, must precisely define what is a "very long time".In fact, one has to compare the time scale of the experiment to the characteristic time where the vicious neurons will eventually fire.Note also that since only a very small portion of neurons can be observed e.g. in a given cortex area, some "vicious" neurons could be present (without being observed since not firing), with the important consequence discussed in this paper.The observation of "temporarily silent" neurons which firing induces a large dynamic change would be an interesting issue in this context.
As a final remark we would like to point out the remarkable work of Latham and collaborators discussing the effects induced by the addition or removal of a single spike in a raster plot.A central question is whether this "perturbation" (which is not necessarily "weak") will have a dramatic effect on the further evolution (see (Latham, Roth, Hausser, & London, 2006) and the talk of P. Latham available on line at http : //www.archive.org/details/RedwoodCenter 2006 09 25 Latham).Especially the questions and discussions formulated during the talk of P. Latham are particularly salient in view of the present work.As an additional remark note that a perturbation may have an effect on trajectories but not on the statistics build on these trajectories (e.g.frequency rates) (Cessac & Sepulchre, 2006).
This leads to: If neuron k does not fire between t and t+δ we have, integrating the previous equation for t 1 ∈ [t, t+δ[ and setting t 2 = t + δ : A.2 Proof of theorem 1.
The proof uses the following lemma.
Proof of lemma 1 Call Π F (resp.Π F ) the projection onto the subspace generated by the basis vectors e i , i ∈ F (resp. e j , j ∈ F ) and set V F = Π F V (V F = Π F V), F F = Π F F (F F = Π F F). Since each neuron j ∈ F is such that (25): for t sufficiently large, (larger than the last (finite) firing time t j ), these neurons do not act on the other neurons and their membrane potential is only a function of the synaptic current generated by the neurons ∈ F. Thus, the asymptotic dynamics is generated by the neurons i ∈ F.Then, ∀V ∈ Ω(Γ F ,T,ǫ ), V F (t + 1) = F F [V F (t)] and V F (t + 1) = F F [V F (t)].One can therefore focus the analysis of the ω limit set to its projection Ω F (Γ F ,T,ǫ ) = Π F Ω(Γ F ,T,ǫ ) (and infer the dynamics of the neurons j ∈ F via (38)).
Construct now the partition P (T ) , with convex elements given by M [ω] T , where T is the same as in the definition of Γ F ,T,ǫ .By construction, F T +1 is continuous on each element P (T ) and fixing M [ω] T amounts to fix the affinity constant of F T +1 .By definition of T , DF T +1 F V , the derivative of F T +1 F at V, has all its eigenvalues equal to 0 whenever V ∈ Ω F (Γ F ,T,ǫ ) (prop.1).Therefore the image of Ω F (Γ F ,T,ǫ ) under F T +1 F is a finite union of points belonging to M. Since, Ω F (Γ F ,T,ǫ ) is invariant, this a finite union of points, and thus a finite union of periodic orbits.
The dynamics of neurons ∈ F is driven by the periodic dynamics of firing neurons and it is easy to see that their trajectory is asymptotically periodic.Finally, since M = ∪ F Γ F ,T,ǫ the ω limit set of M is a finite union of periodic orbits.

A.3 Average of a function.
Since the dynamics is not uniquely ergodic (there are typically many periodic attractors), one has to be careful with the notion of average of a function φ.We have first to perform a time average for each attractor i, φ(i) = lim T →∞ T t=1 φ(V (i) (t)), where V (i) is an initial condition in the attraction basin of attractor i.Then, we have to average over all attractors, with a weight corresponding to the Lebesgue measure µ (i) of its attraction basin.This gives: where N is the number of attractors.

Figure 7 :
Figure 7: Average value of log(d(Ω, S)) for the model I. (left) σ ∈ [1, 10], τ ∈ [10, 40] ms.We present here a projection in the plane σ, log d(Ω, S), and the vertical bars correspond to the variations with τ L .It allows us to verify the stability of the previous result for higher variability of the synaptic weights.(middle) τ L ∈ [1, 1 . . .2]ms below the usual 20ms value, σ ∈ [1, 10].Such range corresponds to cortical neurons in high-conductance state.It allows to check the behavior of d(Ω, S) in this case.(right) Sampling period of 1ms, in order to verify the robustness of the numerical results with respect to the sampling rate.

Figure 8 :
Figure 8: Examples of raster plots for the conductance based model (Model I, top row) and the leaky integrate and fire model (Model IV, bottom row).A time window of 100 samples is shown in each case.The control parameter is τ L = 20ms.As visible in Fig.5, σ = 0.05 corresponds to a small order dynamics where the periodic behavior is clearly visible, and σ = 0.40 to the "edge of chaos".One blob width is 1msec

Figure 9 :
Figure 9: Raster plots for models I (upper row) and IV (lower row), with σ = 10.00 and the same condition as in Fig. 8. First column: model I and model without noise.Second column: same realization of synaptic weights and same initial conditions but with a small amount of noise in the external current.The noise is added to the membrane potential and its magnitude is very small (10 −4 × θ).One blob width is 1msec.