Phase diagram of spiking neural networks

In computer simulations of spiking neural networks, often it is assumed that every two neurons of the network are connected by a probability of 2%, 20% of neurons are inhibitory and 80% are excitatory. These common values are based on experiments, observations, and trials and errors, but here, I take a different perspective, inspired by evolution, I systematically simulate many networks, each with a different set of parameters, and then I try to figure out what makes the common values desirable. I stimulate networks with pulses and then measure their: dynamic range, dominant frequency of population activities, total duration of activities, maximum rate of population and the occurrence time of maximum rate. The results are organized in phase diagram. This phase diagram gives an insight into the space of parameters – excitatory to inhibitory ratio, sparseness of connections and synaptic weights. This phase diagram can be used to decide the parameters of a model. The phase diagrams show that networks which are configured according to the common values, have a good dynamic range in response to an impulse and their dynamic range is robust in respect to synaptic weights, and for some synaptic weights they oscillates in α or β frequencies, independent of external stimuli.


INTRODUCTION
A simulation of neural network consists of model neurons that interact via network parameters; The model is a set of differential equations that describes the behaviors of neurons and synapses, captured by years of electro physiological studies (Nicholls et al., 2012). A model can be as simple as Integrate-and-fire or as realistic as Hodgkin-Huxley (Sterratt et al., 2011). Realistic models generate reliable results, but they are computationally expensive. On the other hand, simple models consume less computer resources, at the cost of loosing details and biological plausibility; therefore, there is a delicate balance between computational feasibility and biological plausibility, Fortunately, there are many reviews and books that compare different models (Izhikevich, 2004;Herz et al., 2006;Sterratt et al., 2011) and neural simulators (Brette et al., 2007). They help us to find a suitable model of a single neuron, set its parameters, and get reliable results in a reasonable time.
No matter how much a neural model is detailed and efficient, A carelessly constructed networks of perfect neural models, is misleading. It only produces wrong biologically plausible results, in a reasonable time. Unfortunately, network parameters are a lot more ambiguous.
There are four essential parameters in a neural network simulation that I investigated in this work -the ratio of excitatory neurons to inhibitory neurons, the topology and sparseness of connections, inhibitory and excitatory synaptic weights. To make it more general and realistic, one may consider the values follow a distribution. This distribution function may change the behavior.
The actual values in computer simulations usually are less diverse, In many computational studies, 20% of the neurons are inhibitory, the network topology is a random graph with sparseness of 2% and there are only a couple of fixed synaptic weights for inhibitory and excitatory neurons. Despite their simplicity, they observe many interesting phenomena like oscillations, synchrony, modulations (Song et al., 2000;Brette et al., 2007;Goodman and Brette, 2008;Buice and Cowan, 2009;Akam and Kullmann, 2010;Mejias and Longtin, 2012;Augustin et al., 2013) and even psychological disorders (Bakhtiari et al., 2012). There are also many works that show the results are not actually sensitive to the parameters (Prinz et al., 2004;Marder et al., 2007;Marder, 2011;Marder and Taylor, 2011;Gutierrez et al., 2013).
Why do this parameters appear in nature? What are their evolutionary advantages? How much is the result of a simulation robust in respect to the parameters? How does one decide about the network parameters in a simulation? Which other parameters might give the same results? These are few question that come to mind, while dealing with properties of neural networks, and they could be easily answered by looking at a phase diagram of neural networks, that shows the behavior of a network in respect to its parameters. Like the phase diagram of water, tells us about the properties of water in different pressures and temperatures. But such a phase diagram does not exist for spiking neural networks, yet. First, because of the large amount of computation needed and second, because of the ambiguity of the concept of phase or state in neural networks.
To answer these question, I could use sophisticated evolutionary algorithms, like the study of evolution of inteligence by McNally et al. (2012), however, I used a simple brute force search in the space of parameters, so I can also generate a phase diagram, similar to the work of Roy et al. (2013).

METHOD
I did a brute force search in the space of parameters, by simulating a population of 1000 neurons, let say a part of a single layer of cortex, that receives excitatory signals from previous layer of cortex or an external stimuli, and then sends excitatory signals to the next layer. There are lateral connections among all types of neurons (excitatory and inhibitory) in a single layer, but only excitatory neurons interconnects different layers. I stimulated this network and watched its response. This was repeated over a wide range of parameters. At the end I calculated fitness of each set of parameters and then investigated the overlap of the best parameters according to my simulation with the nature's choices (or computational neuroscientists' choices).
To do any sort of evolutionary analysis, first we need to agree on something to optimize -a fitness function, it can be the network capacity to learn and store information (Barbour et al., 2007), or it can be the speed and the quality of transmitted signals (Chklovskii et al., 2002). But here, to define the fitness function, I remind you of three obvious evolutionary facts: 1. High Dynamic Range is evolutionary favored, an animal that can see during days and nights, has a higher FIGURE 1 | This figure demonstrates a typical tuning curve and its dynamic range. Two consecutive stimuli are different by a factor of 2, also known as one stop.
survival chance than an animal which can see only during days.

Small Just-Noticeable Difference is evolutionary favored, an
animal which can discriminate more gray levels or colors, has a better chance to find its preys and foods. 3. Low Energy Consumption is also favored, an energy efficient animal, consumes less food and therefore has a better chance to survive (Niven and Laughlin, 2008).
Dynamic Range is the ratio of the largest possible value of stimuli to its smallest perceivable value. Just-Noticeable Difference  is the smallest detectable difference between two sensory stimuli, which is often roughly proportional to the magnitude of the stimulus. The ratio of Just-Noticeable Difference to the magnitude of the stimulus is known as Weber constant (Dayan and Abbott, 2001). Both Dynamic Range and Just-Noticeable are well presented in logarithmic scale, but instead of log 10 and dB, from now on, I will use log 2 and stop which are common in Optics and Photography (Figure 1).
Here we are left with a neural network and three concepts to play with-Dynamic Range, Just-Noticeable Different and Energy Consumption. To keep it simple, I only measure Dynamic Range while there are upper bounds for the other two. I require Just-Noticeable Different to be smaller than one stop, thus, the calculated Dynamic Range is a lower bound for the real Dynamic Range, and I require the Energy Consumption to be zero in absence of stimulus, which puts an upper bound for the Energy Consumption. Moreover, energy consumption can be a determinant factor when all other conditions are equal, when there are two sets of parameters with about the same Dynamic Range, the one which consumed less energy will be favored, and here, the measure of energy consumption is the number generated spikes (Attwell and Laughlin, 2001).
I stimulated the network with a set of stimuli that were equally spaced in the logarithmic scale, separated by 1 stop (×2) intervals, similar to Figure 1. I ran the simulation for a given time 1024 ms and then I calculated Dynamic Range over the range that Just-Noticeable Difference was smaller than 1 stop.
To be precise, I used a sequence of n stimuli plus the resting state (no stimulus), so the sequence includes n + 1 members s ∈ (s 0 , s 1 , s 2 , . . . , s n ). Except s 0 , the resting state, the difference between i-th and (i + 1)-th stimulus is one stop: s i+1 = 2 * s i . Each stimulus s i evokes a response r(s i ) ∈ R, where R = (r(s 0 ), r(s 1 ), r(s 2 ), . . . , r(s n )). Now, the Dynamic Range is simply the number of elements in the largest subsequence of R ⊂ R that fulfills the expectations, 1. For any r(s i ), r(s j ) ∈ R , we should have s i < s j ⇒ r(s i ) < r(s j ).
It means that R is strictly monotonic and it satisfies the constrain on the Just-Noticeable difference. 2. For any s i ∈ R , the network must return to resting state in acceptable time. This satisfy the Energy Consumption constrains.
For example, if stimuli (0, 1, 2, 4, 8, 16) leads to the responses (0, 0, 2, 3, 4, 4) then the subsequence (0, 2, 3, 4) ⊂ R is the largest subsequence of R which meets all the conditions; Therefore, the dynamic range is 4. In this way the minumum dynamic range can be 1 and its maximum can reach the total number of stimuli. In this study I assumed response r(s) as the maximum firing rate of population of excitatory neurons in response to the stimulus s. But the paradigm of this study is general and can be applied to any neural coding schemes.
I stimulated some of the neurons with an impulse, I forced s i excitatory neurons to fire once at time t = 0, then I monitored the propagation of the impulse. Here the number of excited neurons represents the magnitude of the stimuli s i . This is very similar to Optogenetics stimulation, in which, a laser pulse excites a specific population of neurons (Peron and Svoboda, 2011;Toettcher et al., 2011;Lim et al., 2012;Liu et al., 2012).
I used a simple compartment Izhikevic model (Izhikevich, 2003). It can imitate the behavior of Hodgkin-Huxley model, and its computational costs is comparable to the Integrate-and-fire model (Izhikevich, 2004). It is based on a system two-dimensional differential equations: Here white means there was not any simulation data at that bin, all the data were unreliable (they fired at the maximum possible population firing rate), or the network did not show any activity after 100 ms of simulation. The information are arranged according to the template of My network has a random graph topology (Erds and Rényi, 1960). In a random graph, every two neurons are connected by probability p, called sparseness. This is one of the simplest topology for a network or graph. I simulated networks of N = 1000 neurons with N I ∈ (100, 200, 300, 400, 500) inhibitory neurons. The neurons were connected by a random graph topology, with the sparseness p ∈ (0.01, 0.02, 0.04, 0.08, 0.16).
On the other hands, the synaptic weights are chosen randomly, but in a way to cover the region of interest in the w e , w i plane. In this way, it is always possible to add new points later to improve results. The new points could be just another pair of connection weights(w e , w i ) on the same random netwrok, or they could be synaptic weights of a whole new realiziation. I stimulated each network 10 times, 1,2,4,8,16,32,64,128,256) so in each trial, s neurons fire at t = 0. Then I watched the network for 1024ms. At the end, I calculated Dynamic Range of the network. I repeated the whole process for few realization of random networks and and then the median of all data in each hexagonal is calculated and displayed (Figure 2). The whole process needs a great amount of computational power. That is the reason I used NeMo (http://nemosim. sourceforge.net/), a neural simulator software that runs on GPU (Fidjeland et al., 2009;Fidjeland and Shanahan, 2010), but for the pilot study, I used Brian (http://briansimulator.org/) (Goodman and Brette, 2008). Both of them are open source and available on their websites.

RESULTS
I wanted to present the results as a function of 4 other parameters, for that I needed to map a 4-D space to a 2-D space. Here I used the template in Figure 3, to present the dynamic ranges in Figure 4, oscillations in Figure 5 , diuration of activities in FIGURE 6 | The time of last spike in each trial. Dark red means the simulation ended before the network returns to its resting state. Here white means there was not any simulation data at that bin, or all the data were unreliable. The information are arranged according to the template of Figure 3.

Frontiers in Computational Neuroscience
www.frontiersin.org March 2015 | Volume 9 | Article 19 | 5 Figure 6, maximum achived firing rate in each trial in Figure 7 and the time of achived maximum rate in Figure 8 First I made the plots of dynamic range in the space of synaptic weightsw e and w i . Then I put these plots next to each other according to sparseness of connections p and number of inhibitory neurons N i .

DISCUSSION
In Figure 4, top right corner of each figure is empty (white). This is where the excitatory neurons dominates and the network riches its maximum firing rate, this region is visible in Figure 7 as the darkest shade of red. The lesft side of each graph, where the inhibitory subnetwork dominates, is also empty, but this time because networks did not return to resting states. More detail about the length of activity demonstrated in Figure 6. A narrow band seperate this two regins, on this band, the inhibitory and excitatory network are balancing out each other. As we increase the number of connections, this band gets narrower, the delicate balance can be easily disturbed. The dynamic ranges on this band is low or moderate.
At the bottom of each graph, there is another region with valid data points. In this region, excitatory sub networks are not strong enough to saturate the system, so we have a good dynamic range regardless of the inhibitory synaptic weights. Nevertheles, stronger inhibitory sub networks improves the dynamic range.
For the calculation of dynamic range, only trials who returned to resting state have been used. But in many other trials the network sustained its activity upto the very last miliscond of simulation (Figure 6). These networks shows interesting oscillatory behaviures, that can be seen in Figures 5, 9. The oscillations are mostly α and β waves, with fewer case of γ waves. There are also δ and θ waves, but the simulation time of 1 s may not be enough to study these waves.
The networks are stimulated with only a pulse, and then they are left alone, even though there is not any spontaneous activity implemented in the model, yet, the networks show sustained oscillations.
The oscillation are mostly β-waves, and they happen in the region of dominating inhibitory sub networks. As we get to the balanced inhibitory-excitatory band, the oscillations become αwaves.
It seems that the results are not very sensitive to the parameters, We can have good dynamic range over some value of parameters, and we can have neural oscillators in different settings. These finding agrees with Prinz et al. (2004)  The results of this paper, specially Figures 4, 5, serves as a starting point to decide about parameters in a neural simulation, and they may help to find the networks with desired behavior.
This work by no mean is complete, I only studied the effect of four parameters, in just one model, over dynamic range and oscillations. Many other parameters could be included, like synaptic delay and distribution of weights or connections. Also it is interesting to know if my results holds for simpler models, like Integrate-and-fire.
At the end I suggest that any numerical simulation of neural network should be accompanied with such a phase diagram, it demonstrate the robustness of the results in respect to the parameters. One good example is the the work of Roy et al. (2013), where they made such diagram to compare experimental data of orientation selectivity index of mouse primary virtual cortex to their model.