{'fraction_non_alphanumeric': 0.04513319911980898, 'fraction_numerical': 0.012968771946252166, 'mean_word_length': 4.0235183443085605, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Spontaneous Dynamics of Asymmetric Random Recurrent Spiking Neural Networks

We study in this paper the effect of an unique initial stimulation on random recurrent networks of leaky integrate and fire neurons. Indeed given a stochastic connectivity this so-called spontaneous mode exhibits various non trivial dynamics. This study brings forward a mathematical formalism that allows us to examine the variability of the afterward dynamics according to the parameters of the weight distribution. Provided independence hypothesis (e.g. in the case of very large networks) we are able to compute the average number of neurons that fire at a given time -- the spiking activity. In accordance with numerical simulations, we prove that this spiking activity reaches a steady-state, we characterize this steady-state and explore the transients.


Introduction
Many modern neurobiological problems have to confront the behaviors of large recurrent spiking neuron networks. Indeed, it is assumed that these observable behaviors are a result of a collective dynamics of interacting neurons. So the question that rises becomes : given a connectivity of the network and a single neuron property, what are the possible kinds of dynamics ?
In the case of homogeneous nets (same connectivity inside the network), some authors found sufficient conditions for phase synchronization (locking) or stability (Chow, 1998;Gerstner, 2001). Another author (Coombes, 1999) calculated Lyapunov exponents in a given symmetric connectivity map and showed that some neurons were "chaotic" (i.e. the highest exponent was positive). In the very general case (Golomb, 1994;Vreeswijk & Sompolinsky, 1996;Meyer & Vreeswijk, 2002) it has been shown that the dynamics can show a broad variety of aspects.
In the particular case of integrate and fire neurons, (Amit & Brunel, 1997a, 1997b used consistency techniques on nets of irregular firing neurons. This technique allowed them to derive a self-sustaining criterion. Using Fokker-Planck diffusion, the same kind of methods was used in the case of linear I&F neurons in (Fusi & Mattia, 1999;Mattia & Giudice, 2000), for stochastic networks dynamics with noisy input current (Giudice & Mattia, 2003) and in the case of sparse weight connectivity (Brunel, 2000).
However, stochastic recurrent spiking neurons networks are rarely studied in their spontaneous functioning. Indeed most the time the dynamics is driven by an external current -whether meaningful or noisy. Nevertheless, without this external current, the resulting dynamics is often conjectured. Our experimental results showed that large random recurrent networks exhibit interesting functioning modes. Depending on a coupling parameter between neurons (e.g. the variance of the distribution of weights) the network is able to follow a wide spectrum of spontaneous behavior -from the neural death (the initial stimulation is not enough to produce any further spiking activity) to the limit of maximal locking (some neurons fire all the time and the others never). In the intermediate states the average spiking activity grows. However we found out that the average spiking activity reaches a steady-state.
Thus, following basically the same ideas as (Amit & Brunel, 1997a;Fusi & Mattia, 1999), we try to predict these behaviors when using large random networks. In this case, we make an independence hypothesis and use mean field techniques -note that this socalled mean field hypothesis has been rigorously proved in a different neuronal network model (Moynot & Samuelides, 2002). More precisely, in our case, the connectivity weights will follow an independent identically distributed law and the neurons firing activities are supposed independent.
After introducing the spiking neural model, we propose a mathematical formalism which allows us to determine (with some approximations) the probability law of the spiking activity. Since no other hypothesis than independence are used, a reinjection of the dynamics is needed. It leads expectedly to a massive use of recursive equations. However non intuitive, these equations are a solid ground upon which many conclusions can rigorously be drawn.
Fortunately, the solutions of these equations are as expected (and often taken for granted), that is the average spiking activity (and as a consequence the average frequency) reaches a steady state very quickly. Moreover, this steady state depends only on the pa-rameters of the weight distribution. To keep the arguments simple, we detail the process for a weight matrix following a centered normal law. Extensions are proposed afterward for a non-zero mean and a sparse connectivity. All these results corroborate accurately with simulated neural networks data.

The neural model
The following series of equations describe the discrete leaky integrate and fire (I&F) model we use throughout this paper (Tuckwell, 1988). Our network consists of N all-to-all coupled neurons. Each time a given neuron fires, a synaptic pulse is transmitted to all the other neurons. This firing occurs whenever the neuron potential V crosses a threshold θ from below. Just after the firing occurred, the potential of the neuron is reset to 0. Between a reset and a spike, the dynamics of the potential is given by the following (discrete) temporal equation : The first part of the right hand side of the equation describes the leak currents (0 ≤ γ ≤ 1 γ is the leak). Obviously, a value of 0 for γ indicates that the neuron has no short term memory. On the other hand, γ = 1 describes a linear integrator.
The W ij are the synaptic influences (weights) and I i is an external input. δ(t) = 1 whenever t = 0 and 0 otherwise (Kronecker symbol). The T j are the times of firing of a neuron j. d ij is an axonal delay and (as well as T j ) is a multiple of the sample discretization time. The times of firing are defined formally for all neurons i as : and so we define T n i (the date of the n-th firing date) recursively as follow : We set that T 0 i = −∞. Moreover, once it has fired, the neuron's potential is reset to zero.
∀n, ∀i, lim In this discrete description of an I&F neuron, equations (2-2) and (2-4) can be disturbing. They only mean that when computing V i (T n i + 1) we set V i (T n i ) = 0 in equation (2-1).
In our study, we examine only spontaneous activity. It means that except at t = 0 there is no external input : ∀ i, ∀ t > 0, I i (t) = 0. Moreover, given the above description, it is possible to allow great heterogeneity for various parameters. In order to simplify, we restrict ourselves in this paper to a synaptic weight heterogeneity leaving the other parameters constant. So ∀ i, ∀ j, d ij = 0 (the axonal delay is 0). Similarly, ∀ i, θ i = θ (same threshold for all neurons). Finally, suppose that the weights follow a centered normal law N (0, σ 2 ) and let φ = σ √ N .

General Study
In this section we give a very general formulation of the distribution of the spiking activity defined as X t the numbers of firing neurons at a time step t for an event. We partition X t according to the instantaneous period of the neurons. That is we write : is the number of neurons who have fired at t and t − k but not in between.
Suppose that all neurons had a potential 0 at the starting process and only X 0 neurons were excited in order to make them fire, so V i (1) = X 0 j=1 W ij . Thus using equation (2-3) we have: where χ {V i (1)>θ} = 1 whenever V i (1) > θ and 0 otherwise. Keeping on we have : Indeed, the number of firing neurons at time step 2 are those that have fired twice (at t = 1 and t = 2 that is X (1) (2)>θ} since the reset potential is 0). We need to add those that have not fired at t = 1 (X Thus for t, taking into account the initial step, we have recursively : Now assume that neurons dynamics are independent, we can calculate the expectation of X t (by the so-called first Wald's identity) : The P (k, t) are the expectation of a Bernoulli distribution. It leads that : We are now able to retrieve the variance of X t (second Wald's identity) : More generally the moment generating function can be recursively computed :

Average Number Calculation
The equations (3-6) to (3-8) are useful when we can estimate the P (k, t) coefficients. We Thus the potential is a random sum of independent identically distributed normal laws (the weights). Unfortunately, this random sum is not, in general, a normal law itself. Nevertheless, as it is proved in appendix A, (see equation A-28), when the number N of neurons is large enough, and for a general class of distributions of the random variable X t , we can write, for a neuron that has its potential V i (t) = 0 : From now on we suppose that γ = 0. We need to know which are the neurons that have their potential V i (t) = 0. When all neurons receive a charge from a non-zero number of neurons, the probability to have the potential equal to zero vanishes. Thus only the neurons that fired at time t have this property (except for t = 0 when it is true for all the neurons). But taking into account a non-zero leak, we need to make here another independence assumption concerning the previous charges. Indeed the charge received by a neuron between t and t + k comes from γ k X t + . . . + γX t+k−1 + X t+k "neurons". But in order to proceed further we assume that theses charges are independent. Note that when γ is equal to zero this extra hypothesis is not needed. It leads to : where we set In order to lighten notations, we put x t = E(Xt) N , and : We recall that we set φ = σ √ N . So equation (4-10) becomes : Using the same notation we get the recursive computations of E(X t ): Equations (4-12) and (4-13) are our main result. We can see that there is no more reference of the number of neurons. However, the algorithm grows exponentially with the time.
Recalling a more classical spiking neurons formalism, we remark that equation (3-6) can be viewed as an integral over past charges of the form : This is exactly Gerstner's formula (Gerstner, 2001) to compute spiking activity defined (using our notations) as : In the case γ = 0 we can deduce the result independently from the above equations.
Indeed, since it means that at each time step the potential is reset to zero, we can directly write the probability of spiking of one neuron according to the number of neurons that have previously fired. Thus, with the same notations and hypothesis : But it is a special case of equation (4-12) when γ → 0. See appendix B for a detailed proof.

Analysis
The equation (4-12) is difficult to estimate. However some important things can be claimed.
Following from the definition of the P (k, t) we have for large enough t : It leads that x t is bounded and for large enough t monotonic. Therefore, x t converges towards x ⋆ when t → ∞. x ⋆ = 0 (neural death) is an obvious solution. For high enough φ another fixed point exists bounded by 1/2. Moreover, in this case, the P (k, t) are close to a geometric distribution of parameter p φ (x ⋆ ). Due to the definition of P (k, t), it leads to a geometric distribution of Inter-Spike Intervals that is : It enables us to define a network frequency f ⋆ defined as : But if we define the network average frequencyF (t) of a network over a period T and at a given time t by :F where δ i (t) = 1 if the neuron i has fired at time t and δ i (t) = 0 otherwise (in other words . Switching the sum symbol gives : for this realization of the distribution. Taking the expectation leads to : It means that the average frequency (on a time window T ) is the average of the spiking activity (over a period T ). Thus, when t → ∞ it leads to : Due to discrete timing, we generally don't have f ⋆ =f . Instead we have f ⋆ <f . It gives : The inequality becomes an equality if γ = 0. This is the case we now study.

Simple case
If we consider the case γ = 0, we recall that A solution x ⋆ = 0 exists when p φ (x) crosses the line y = x and is stable if and only if

Previous charges independence
In the case γ = 0, we now need the independence of charges hypothesis. However in the general case, we allowed the potential to have strong negative values. More than to be biologically wrong, it dramatically impedes the independence hypothesis. Indeed, some neurons with very low potential will never fire whatever happens.
In order to take this into account, we make a (biologically plausible) assumption : the potential is not allowed to decrease under a minimal value v min . This leads us to reconsider the charge function : We recall that under the hypothesis of independence and since we have a normal law for the weights it leads to : wherep φ is the new probability function we need to compute.
Let's assume that a neuron i has taken a charge C in the previous times and is subject to a charge x t at time t. The probability that the total charge exceeds the threshold must be split in two cases : the potential occurred by C is below v min or not. In the first case (below v min ), the charge will be x t on a potential v min with probability p = P (Φ < v min ), where Φ ∼ N (0, Cσ 2 ). In the other case it will be C + x t (with probability 1 − p). So the resulting probability will be : We can find the P (k, t) recursively noting that the probability p depends on previous x t .
To simplify, we suppose that v min = 0 (i.e the neuron cannot have negative values) it is even simpler. In this case, whatever the charge is, the probability p is always 1/2. So equation (5-26) becomes : It acts as if the decay rate was divided by two. Thus into equation (4-12) it leads to : where ∀m > 0,x m = x m andx 0 = 1.

Results and comparison
We conducted extensive numerical simulations to confront our formulas. The results we present here were made with 500 random networks of 1000 I&F neurons. The starting stimulated neurons were chosen randomly for each network tested -each neuron was scanned and set to fire with a probability x 0 . x 0 was chosen between 0.05, 0.30, 0.5, 0.8 and 1.0.
To be in concordance with equation (5-27) we set v min = 0.

Asymptotic Spiking Activity Comparison
We recorded the limit of x t according to different values of φ and for γ = 0.0, 0.5, 1.0. The results of the theoretical and experimental simulations are displayed in figure 1. The limit of x t is an increasing function of φ. Moreover lim φ→+∞ (lim t→+∞ x t ) = 1/2. As easily seen, the higher φ and lower γ give an extremely good prediction compared to experimental data. A high γ always gives a slight overestimation in the case of high φ. Indeed, the independence of charge never holds. On the other hand, when this hypothesis is rigorously true (i.e γ = 0) the prediction is strikingly accurate. For small values of φ, the prediction was correct also. However it exists an interval of intermediate values of φ (that is here between 1.5 and 2.0) where the prediction completely fails. Besides we shall see in the next section that this is also the case for the transient period. Indeed, for small enough coupling, the independence of the dynamics of neurons is not true. We expect that the length of this interval tends towards zero when N grows infinitely.

Transient period
We recorded the evolution of the spiking activity for the transient period before reaching the asymptotic limit for various values of φ and γ during the first 50 time steps. The results are displayed on figure 2. As in the previous section, a higher γ gives a slight overestimation but is completely correct when γ draws near zero. Moreover, we retrieve the same interval of φ where the neurons dynamics independence fails (and as a consequence the transient).
These results show that the steady-state is reached very quickly (less than 10 time steps).

Generalization
We supposed to make the arguments clear that the weight distribution followed a very simple law -a centered normal law. However, we do not have to restrict ourselves to this limit. Indeed, we look at two possible extensions of spiking activity equations. We show that it is possible to use a non-zero mean for the distribution as well as a sparse weight matrix in a certain sense. The results obtained showed the same characteristics.

Non-zero mean
In fact taking a non-zero mean is very easy. Indeed, provided that all the hypotheses remain valid, we can insert the mean into the definition of p φ which becomes p φ,µ defined as : . In order to relate this to a very simple case, we can look at what happens when φ → 0 -that is an homogeneous network. It leads to p φ,µ (y) = 0 when µy < θ, p φ,µ (y) = 1 when µy > θ and finally p φ,µ (y) = 1/2 when µy = θ . All depends on the value µx 0 θ .
Case µx 0 θ < 1 In this case, x 1 = p φ,µ (x 0 ) = 0 and it converges to neural death that is Case µx 0 θ > 1 Here, x 1 = p φ,µ (x 0 ) = 1 leads to ∀ t > 0 x t = 1. All neurons fire at the same time and with the highest frequency (one spike per time step).
The general behavior does not change with non-zero mean. The convergence is very fast and the transient are well predicted as shown in figure 3. The only difference with µ > 0 is that x ⋆ can be well above 1/2 and ultimately reach 1.

Sparse connectivity
A pure normal distribution is unlikely to occur in biological neurons. Most of the time, neurons are linked to a certain amount of neurons in the network. It can be viewed using a sparse weight matrix (a matrix with zero coefficient). In order to take this into account we compute sparse matrix as follows. A weight w has a probability p to be zero and a probability 1 − p to follow a normal law N (µ/N, φ/ √ N ). As in the previous section it leads to a newp φ,µ function. When calculating the charge, it came from "X" neurons leading to a sum of X normal laws. In the case of sparse matrix, it reduce to (1 − p)X neurons. So our new functionp φ,µ becomes:p φ,µ (y) = p φ,µ ((1 − p)y) It remains to insert this new function into equation (5-27). As shown in figure 4, the general behavior is not changed when introducing sparsity.

Discussion
Independence hypothesis can be a very powerful way to approach random networks of spiking neurons. So far, we proved that the "coupling factor" (φ = σ √ N) can characterize the average spiking activity. For instance whatever the initial stimulation (provided it is strong enough and independent of the neurons), the networks average spiking activity reaches the same steady-state. Since the death locking is also a possible steady-state, in dynamical systems terminology, depending on the value of φ, we exhibited a bifurcation.
In addition, this steady-state only depends on the coupling parameter (for a fixed type of neurons i.e. same leak). This dependence (if not obvious) is increasing and converges to 1/2 when φ grows larger. That is a limit state where half the neurons that are active fire all the time and the remaining stay silent. That is an extreme locking. Indeed, all neurons are periodic (either 1 or ∞) and those that fire do it synchronously (that is all the time).
The transient phase was also well determined -very quick convergence to the steadystate. We were also able to derive mathematically a sufficient condition for death locking (in a very simple case though).
However, as showed for a certain range of coupling factor, independence is deemed to fail and consequently the prediction -both in transient and asymptotically. Indeed, intuitively when a coupling value is intermediate, neurons can no longer behave as if it did not matter who fired. It seems somehow that this range of "failed independence" grows thin with increasing N (at the thermodynamical limit).
This shows that a spontaneous regime can be self-sufficient. It means that networks can afford discontinuous inputs without losing their internal activity. In other words, it is theoretically possible to observe dynamic cerebral areas without any permanent external stimulation.

A Random sum of random variables
Let us first prove a general result, which will have an application for the neuronal potential, written as a random sum of i.i.d. random variables.
Let f be a function so that lim k→∞ f (k) = α ∈ R, (p N k ) (k,N )∈N 2 a sequence satisfying ∀k ∈ N, lim N →∞ p N k = 0 and N k=1 p N k = 1. Let now define g(N ) = N k=1 p N k f (k), and prove that N 0 being fixed, it just remains, to complete the proof, to get a rank N 1 so that Application : let X (N ) be a sequence of random variables on [1, N ] so that

B Simple Case
We prove here that substituting γ = 0 in equation (4-12) give the equation : It yields : For t = 0, it gives x 1 = p φ (x 0 ). Using the recurrence hypothesis ∀m = 1 . . . t, x m = p φ (x m−1 ) we find that : Is is now enough to note that We showed by recurrence that we get equation (4-15).

C Sufficient condition for death
Let's go back to for y ∈ [0, 1]. Taking the derivative over y we have : We see immediately that p ′ φ (y) ≥ 0 so p φ is increasing. So, stable non-zero fixed points should appear for a value of y that crosses the line y = x from above. Then a sufficient condition is that ∀y : Let z = 1 y , τ = θ √ 2φ and g(z) = p ′ φ ( 1 z ) then : Taking the derivative of g(z) yields : So g ′ (z) has the same sign as 3 2τ 2 − z and since g ′ (0) = g ′ (+∞) = 0, it leads that the maximum of g(z) is obtain for z = 3 2τ 2 . But z = 1 y thus z ∈ [1, +∞[. It leads that So a sufficient condition for zero to be the only fixed point becomes :