Computational model of neuron-astrocyte interactions during focal seizure generation

Empirical research in the last decade revealed that astrocytes can respond to neurotransmitters with Ca2+ elevations and generate feedback signals to neurons which modulate synaptic transmission and neuronal excitability. This discovery changed our basic understanding of brain function and provided new perspectives for how astrocytes can participate not only to information processing, but also to the genesis of brain disorders, such as epilepsy. Epilepsy is a neurological disorder characterized by recurrent seizures that can arise focally at restricted areas and propagate throughout the brain. Studies in brain slice models suggest that astrocytes contribute to epileptiform activity by increasing neuronal excitability through a Ca2+-dependent release of glutamate. The underlying mechanism remains, however, unclear. In this study, we implemented a parsimonious network model of neurons and astrocytes. The model consists of excitatory and inhibitory neurons described by Izhikevich's neuron dynamics. The experimentally observed Ca2+ change in astrocytes in response to neuronal activity was modeled with linear equations. We considered that glutamate is released from astrocytes above certain intracellular Ca2+ concentrations thus providing a non-linear positive feedback signal to neurons. Propagating seizure-like ictal discharges (IDs) were reliably evoked in our computational model by repeatedly exciting a small area of the network, which replicates experimental results in a slice model of focal ID in entorhinal cortex. We found that the threshold of focal ID generation was lowered when an excitatory feedback-loop between astrocytes and neurons was included. Simulations show that astrocytes can contribute to ID generation by directly affecting the excitatory/inhibitory balance of the neuronal network. Our model can be used to obtain mechanistic insights into the distinct contributions of the different signaling pathways to the generation and propagation of focal IDs.


INTRODUCTION
The intracellular Ca 2+ elevations occurring in cultured astrocytes in response to a glutamate challenge (Cornell-Bell et al., 1990) was the initial observation that hinted at the existence of a form of excitability in astrocytes based on cytosolic Ca 2+ concentration changes. A few years later, Ca 2+ elevations in astrocytes from both cell cultures (Parpura et al., 1994) and brain slices (Pasti et al., 1997) were observed to result in Ca 2+ increases in nearby neurons mediated by astrocytic glutamate. Considering that astrocytes occupy non-overlapping spatial territories (Bushong et al., 2002;Halassa et al., 2007b) and that the processes of a single astrocyte can contact hundreds of synapses (Ventura and Harris, 1999), it was suggested that astrocyte-to-neuron communication may play a fundamental functional role in the brain. It was also found that astrocytes establish extensive contacts with cerebral blood vessels (Simard et al., 2003), which added further complexity to the functional role of neuron-to-astrocyte signaling. This neuronastrocyte-blood-vessel signaling pathway was revealed to be central in neurovascular coupling, the process by which episodes of intense neuronal activity at restricted brain regions trigger local increases in cerebral blood flow to satisfy the energy demand of active neurons (Zonta et al., 2003;Mulligan and MacVicar, 2004;Gordon et al., 2008).
These pioneering results lead to the idea that astrocytes and neurons establish a bidirectional communication in the brain which may play fundamental roles in the modulation of synaptic transmission and plasticity (Carmignoto, 2000;Haydon, 2001).
Over the last decade numerous studies provided evidence for the ability of astrocytes to listen and talk to the synapse by exerting both excitatory and inhibitory actions on neurons (Araque et al., 1999;Brockhaus and Deitmer, 2002;Zhang et al., 2003;Pascual et al., 2005;Panatier et al., 2006;Serrano et al., 2006;Jourdain et al., 2007;Perea and Araque, 2007). These studies revolutionized our view of how the brain works. The processing of sensory information in the brain, which has been for many years considered to be based exclusively on neuronal communication, is now viewed as a product of the dynamic signals that neurons and astrocytes constantly exchange in the brain network.
Such a bidirectional communication between neurons and astrocytes was conceptualized in the tripartite synapses in which the astrocyte composes with the pre-synaptic terminal and the postsynaptic target neuron, a third functional element of the synapse (Araque et al., 1999;Carmignoto, 2000;Halassa et al., 2007a;Perea et al., 2009).
The discovery that astrocytes are crucially involved in normal brain function raised the intriguing possibility that these cells may be involved also in brain dysfunction. The observation that glutamate released by astrocytes evokes episodes of synchronous activity in small groups of nearby neurons (Fellin et al., 2004(Fellin et al., , 2006, was the first clue suggesting that gliotransmission represent a relevant non-neuronal mechanism for neuronal synchrony, which may ultimately favor the generation of focal epileptiform activity (Kang et al., 2005;Tian et al., 2005). A new experimental protocol was recently developed by our group in rat enthorinal cortex (EC) slices in order to reproduce the spatial and temporal features of focal epileptiform discharges Losi et al., 2010). In this model, a pharmacological stimulation of neurons from a restricted cortical region induces a propagating seizure-like ictal discharge (ID). The ability to emulate an epileptogenic focus allows us to study the early cellular events that take place during the generation of epileptiform activity as it arises at a focal site and propagates to the surrounding brain tissue. By using this experimental protocol we recently provided evidence that neurons engage astrocytes into an excitatory loop that pushes the neuronal network toward the ID generation threshold .
There are currently many computational models of seizures generation, development and cessation (Pitkänen et al., 2006). The level of description ranges from mean field models (Wendling et al., 2002;Suffczynski et al., 2004) to biophysically detailed models (Destexhe, 1998;Bazhenov et al., 2004;Traub et al., 2005). We used here a simplified approach in the description of the dynamics of single neurons and astrocytes. With this simplified dynamics we implemented a computational network model that allowed us to investigate the network mechanisms of focal ID generation and the role of astrocytes at the onset of the ID.
We found that the positive feedback provided by the astrocytes influences the dynamics of the system and favors the generation of epileptiform activities. The computational model quantitatively reproduces the spatial and temporal features of ID generation and propagation and provides mechanistic insights into the astrocyte contribution.

NEURON MODEL
The computational model aims to reproduce the behavior of a brain network that in response to NMDA pulse stimulation generates a focal ID . The network consists of 320 excitatory and 80 inhibitory neurons randomly disposed and synaptically connected in a 2D configuration. As in our previous work (Reato et al., 2010), we used Izhikevich's model (Izhikevich, 2003) to describe the dynamics of single neurons. Briefly, the voltage dynamics of single neurons is characterized by four parameters: a, b, c, d as follows: With a reset of the dynamic variables u, v when a spike is generated: The choice of values for the four parameters defines different spiking behaviors. The parameters were chosen to reproduce the behavior of a regular spiking neuron for excitatory neurons (a = 0.02, b = 0.2, c = −65, d = 10) and of a fast spiking neuron for inhibitory neurons (a = 0.2, b = 0.26, c = −65, d = 0.5). The variable I represents the sum of the synaptic current and the external stimulation. The synaptic currents mimic AMPA, NMDA, GABA A , and GABA B receptor activation following (Izhikevich and Edelman, 2008). Briefly, the synaptic conductances are described by firstorder linear kinetics, dg x dt = − g x τ x + j s j δ t − t f (where x = AMPA, NMDA, GABA A , GABA B ) with, τ AMPA = 1 ms, τ NMDA = 2000 ms, τ GABA A = 6 ms, τ GABA B = 150 ms. Every time a presynaptic neuron m fires an action potential the conductance of the post-synaptic neurons increases instantaneously by s = s exc and s = s inh for pre-synaptic excitatory or inhibitory neurons respectively (s exc = 0.001, s inh = 0.01). The ratio of NMDA to AMPA receptors was set to be uniform at a value of 2, while GABA B to GABA A equal to 0.3 (s exc = 0.002 for NMDA and s inh = 0.003 for GABA B ). The synaptic current of a post-synaptic neuron is then given by: Where v (function of time, Equation 1) is the voltage of the post-synaptic neuron and v exc and v inh are the reversal potentials for excitatory and inhibitory synapses. Here we chose v inh = −90 mV, v exc = 0 mV. Each neuron receives excitatory inputs from a square of maximum 48 neighbors, while inhibition from maximum eight neurons. Using these parameters a single excitatory pre-synaptic spike induces a depolarization of maximum ∼0.1 mV, while an inhibitory pre-synaptic spike leads to maximum ∼0.5 mV hyperpolarization. All the main parameters of the simulations (the a, b, c, d parameters describing the dynamics of single neurons for both excitatory and inhibitory neurons and the s parameters for synaptic connections) were selected from a normal distribution with standard deviation equal to 1% of the average value. To mimic the onset of an ID, a few parameters of the network were chosen in order to place the network in a hyperexcitable state. The Frontiers in Computational Neuroscience www.frontiersin.org October 2012 | Volume 6 | Article 81 | 2 excitability of excitatory neurons was slightly increased by injecting depolarizing currents (amplitude equal to 2), that could mimic the effects of 4-AP (a K + channel blocker) used in the slice preparation. The high values chosen for both the conductance and the time constant of NMDA currents aim to reproduce the low Mg 2+ experimental conditions. Without stimulation, both excitatory and inhibitory neurons are completely silent.
The NMDA stimulation that in experimental slice preparations evoked an ID was simulated in the model by depolarizing a set of neurons within a 7 × 7 square area above threshold for 500 ms (49 neurons). We refer to this as a simulated pulse (SimP). Alternatively, the NMDA pulses could have been simulated by activating NMDA currents. However, since we are interested in analyzing the effects on NMDA currents during the ID onset, this would have resulted in "stimulation artifacts" (the NMDA current induced by the pulse). Since we were also interested in studying the mechanisms leading to ID generation, the intensity of the stimulation was set to a value that not necessarily induced an ID in all the simulations (see Figure 3D).
In each simulation, nine SimPs were applied. In unsuccessful simulations, the average firing rate in the network increases during each SimP, but it rapidly comes back to zero between successive SimPs. An ID was considered to be successfully generated when the firing rate in neurons remains sustained above 1 Hz. The ID onset was then defined as the number of SimPs which starts this process.
The cessation of the ID was obtained by a modulation of the parameter b in a firing specific way. More specifically, we assumed that an elevated spiking activity decreases the excitability of single neurons. Possible physiological correlates of this event are the inactivation of Na + channels (Bazhenov et al., 2004), the activation of Ca 2+ -or Na + -dependent K + channels (Alger and Nicoll, 1980;Schwindt et al., 1989;Bazhenov et al., 2004;Timofeev et al., 2004) or the exhaustion of metabolic support (Yamada et al., 2001;Kirchner et al., 2006).
The equation used is: Where R(t) is the spike train of a single neuron low-pass filtered (time constant equal to 150 s), m is the coupling constant between the spiking activity and b (chosen here to be 15), and b s the value of b in resting conditions (no spiking activity). The second term in the equation can be thought as a driving force to recover the normal neuronal functionality of the neuron, for example the metabolic support.
Because of the hyperexcitability of the network, i.e., neurons are firing intensively at ID onset, we had to integrate Izhikevich's equations with the method proposed in Izhikevich (2010) assuming the time step to be 1 ms: Where E(t) = g i (t)E i /g(t) with g(t) = g i (t) (the total sum of conductances) and E i = v exc , v inh for excitatory and inhibitory connections, respectively. This method is efficient and stable even for large synaptic currents.
"Excitation" refers to the sum of excitatory currents (AMPA and NMDA) averaged across neurons, and similarly "inhibition" refers to the average summed inhibitory currents (GABA A and GABA B ). Excitatory and inhibitory firing rate indicate the firing rate of excitatory and inhibitory neurons, respectively. Where otherwise indicated, excitation, inhibition and firing rates of single simulations were always filtered with a moving average filter using a 50 ms time window for better visualization. Postictal refractory period was estimated as the time between the end of the seizure (average firing rate back to zero) and the time at which the b variable recovers to the 95% of the initial value.
Under the conditions described above, our computational model is able to generate a neuronal network activity which resembles several characteristics of experimental focal IDs (see later in the "Results"): (i) the simulated ID originates from a small number of neurons in the network and propagates outside the focal area with a delay (Traub and Wong, 1982;Avoli et al., 2002); (ii) the simulated ID arises from an unbalance between inhibitory and excitatory activity at the focal area (Bradford, 1995;Ben-Ari, 2002); (iii) the simulated IDs have a cessation and a similar average duration (Jefferys, 1990;Traub et al., 1993;Pinto et al., 2005); (iv) the network enters into a period of postictal refractoriness (Jefferys, 1990); (v) the peak in the firing rate of the excitatory and inhibitory neurons during simulated IDs is compatible with that measured in the in vitro experimental models.
Our model failed to reproduce the bursting behavior which characterizes the firing discharges in individual neurons and the two main phases in ID development, i.e., the initial tonic and the delayed clonic activity. However, the main focus in this computational model was to include astrocytes in the neuronal network and gain insights into how these non-neuronal cells can affect the equilibrium between excitation and inhibition in the network.

ASTROCYTE MODEL
We introduce here a simple representation of astrocytes interacting with a neuronal network. The parameters related to the ability of astrocytes to respond to neuronal activity with cytosolic Ca 2+ elevations were captured from results obtained in experiments performed both in brain slices (Pasti et al., 1997;Porter and McCarthy, 1997) and in the living brain (Hirase et al., 2004;Wang et al., 2006;Kuga et al., 2011). To simulate the Ca 2+ dynamics of a single astrocyte we used a framework similar to the Izhikevich neuron model. The equations represent a dynamical system of two variables [(Ca 2+ ) and ϕ], without non-linear action potential or reset. The set of equation describing the Ca 2+ concentration has the following form: where ϕ is a recovery variable and σ j is assumed here to be the neuronal input when an action potential is generated by the neuron j, since astrocytes respond to neuronal releases of glutamate (Pasti et al., 1997;Porter and McCarthy, 1997). Although the equations are dimensionless, the values σ j where chosen to reproduce the pattern and amplitude of the Ca 2+ elevations that are experimentally observed in astrocytes in response to neuronal activity (Porter and McCarthy, 1996;Pasti et al., 1997). The values of σ j was chosen as been normally distributed with mean 0.00083 and standard deviation equal to 1% of the mean. Ca 2+ concentration was restricted to be non-negative. Similarly to the dynamics described by the Izhikevich's single neuron, different values of α and β determine different behaviors (time constant of changes and coupling with the recovery variable). Here we chose α = 0.001 and β = 0.01. This choice was made to reproduce the slow time course of Ca 2+ changes in astrocytes (Kawabata et al., 1996). When astrocytes were included in the whole network, these values were chosen to be normally distributed with a standard deviation equal to 1% of the mean.
To describe the release of astrocytic glutamate triggered by Ca 2+ elevations, we considered a first order dynamics (low pass filters), with a release of glutamate that can occur only when Ca 2+ reach a threshold (Pasti et al., 1997;Parpura and Haydon, 2000;Pasti et al., 2001): Where [Ca 2+ ] th = 0.0018 mM is the threshold for glutamate release, κ = 200 describes the coupling between the glutamate concentration [glu] and the recovery variable λ. Glutamate concentration was imposed to be non-negative. The time constants for the two variables were μ = 0.5 s and η = 10 s. The value of [Ca 2+ ] th was set based on available data showing that an increase in astrocytic Ca 2+ of a few hundreds of nM was able to trigger glutamate release (Parpura and Haydon, 2000). Assuming a value of 200 nM for a single synapse (Nadkarni and Jung, 2003) and considering that astrocytes in our model receive inputs from a maximum of nine neurons, the threshold value can be determined by multiplying the value for the single synapse by the number of inputs, as considered in other studies (Wade et al., 2011). The set of parameters used for a single astrocyte reproduces basic features of Ca 2+ dynamics and glutamate release in astrocytes. Increasing the input to an astrocyte, simulated as Poisson-distributed spike trains of increasing frequencies, leads to increasing Ca 2+ concentrations (Figures 1A1-A4). The Ca 2+ increase due to a single spike is less than 100 nM and lasts for about half a second (inset in Figure 1A1). These results are compatible with experimental evidences (Pasti et al., 1997;James et al., 2011) and previous computational models (Jefferys, 2003;Jung, 2004, 2007;Wade et al., 2011). The linear dependence of Ca 2+ increases as a function of simulated firing rate is reported in Figure 1B. Increasing the number of inputs by summing up Poisson-distributed spike trains (color scale from red to blue) also elevated Ca 2+ concentrations. Since the release of glutamate due to the Ca 2+ increases occurs only when Ca 2+ is above a threshold, only strong activation can drive the release. As an example, nine spike trains at 10 Hz induced transient releases of glutamate (Figures 1C1,C2, Ca 2+ threshold in red). Figure 1D summarizes the dependence of glutamate released by the astrocyte as a function of the firing rate and the number of inputs. For very low firing rates, there is no astrocytic glutamate release independently on the number of inputs. In the case of high firing rates, the release is linearly dependent on both the number of inputs and the firing rate. To further validate the parameters that we choose, we stimulated single astrocytes with a spike train from nine neurons from a simulated ID (see later in the "Results"). The neuronal activity leads to Ca 2+ increases in the astrocyte ( Figure 1E1) that caused a glutamate release ( Figure 1E2) only after the second pulse (see also below). Interestingly, when the ID was fully developed, Ca 2+ elevations reached a steady state value and glutamate was no longer released. It is known from experiments in cultures and in brain slices (Pasti et al., 1997;James et al., 2011) that upon intense stimulation the Ca 2+ level in astrocytes increases rapidly and remains at an elevated steadystate value for tens of seconds ( Figure 1E1). A single episode of glutamate release is experimentally observed only after the initial Ca 2+ rise.
Astrocytes were included in the network with a 1:1 ratio with neurons. The ratio of glia to neurons increases in phylogenesis and is 1.65 in the human frontal cortex (Oberheim et al., 2006;Sherwood et al., 2006). Given that astrocytes account for about 50% of the total number of glial cells, a 1:1 ratio seems to represent an acceptable approximation. The input from neuronal activity, σ for each m astrocyte, was considered as the excitatory input from neurons firing s exc (m) (the excitatory component) in a 3 × 3 square (inputs from nine excitatory neurons). This choice was made considering that the feedback of astrocytes on neuronal activity is thought to be local with four to eight neuronal somata enveloped by a single astrocyte (Halassa et al., 2007b), but with the processes from a single astrocyte associated with up to 600 dendrites and many thousands of synapses (Bushong et al., 2002;Oberheim et al., 2006). The glutamate released by astrocytes was used as input to the same neurons to which the astrocyte is exposed. This glutamate generated NMDA currents in these neurons by activating the NMDA channel (in the same way than s exc ). In some simulations (see "Results") we considered the effects of inhibitory inputs from astrocytes. This was done by considering that the astrocytic response activate GABA A receptors instead of NMDA (so simulating the effect of GABA release). All the simulations were performed using MATLAB (Mathworks), and the code is available at www.neuralengr. com/code.

FOCAL ID GENERATION IN ENTORHINAL CORTEX SLICES
As we previously reported , an episode of neuronal hyperactivity can generate a focal ID in EC slices perfused with the K + channel blocker 4-aminopyridine (4-AP) and low Mg 2+ . Figure 2 illustrates a typical ID that was generated in cortical layer V-VI by a double brief pressure pulse applied to an NMDA-containing glass pipette (Figures 2A,B).
Dual patch-clamp recordings revealed that the firing in neurons located within the focus (Figure 2C, neuron 1) evolved into a focal ID with some delay after the NMDA double pulse. Following the ID generation at the focus (Figure 2A, gray circle), neurons outside the focus (<400 μm from the NMDA pipette tip) were also recruited and exhibited a similar pattern of action potential firing ( Figure 2C, neuron 2). Given that the somatic Ca 2+ change in neurons reflects faithfully the action potential firing in these cells, in slices loaded with the Ca 2+ sensitive dye Oregon Green BAPTA1-AM (OG-B1-AM) we could monitor the activity of tens of neurons and follow how a focal ID is generated in the neuronal network. These experiments revealed that the NMDA stimulation The timing that corresponds to the different images (t 0 , t 1 , and t 2 ) is reported.
evoked a rapid Ca 2+ elevation in neurons located within the focal area, while neurons from the surrounding network were recruited into the ID only after a delay of 10.9 ± 0.8 s (30 IDs from 15 slices).

FOCAL ID GENERATION IN THE NEURONAL NETWORK MODEL
In the model we first examined how the neuronal network responds to a sequence of simulated NMDA pulses in the absence of astrocytes. To mimic the NMDA pulses at the focus, a depolarizing current pulse was injected for 500 ms in an area of 7 × 7 neurons (see "Methods"). The first SimP evoked robust spiking activity that remained restricted to neurons of the focus ( Figure 3A1). Upon successive SimPs the firing activity spread from the focus to the surrounding neurons approximately 10 s after the SimP onset (Figures 3A1-A4,B1). The neuronal firing discharge remained high thereafter for tens of seconds (61 ± 2 s) before a sudden cessation ( Figure 3B1). A postictal refractory period was observed with an average duration of 266 ± 1 s (see "Methods"). This pattern of activity resembles the focally evoked ID in slice preparations (see Figure 2D). A raster plot of the activity in a subpopulation of excitatory and inhibitory neurons within and outside the focus revealed that inhibitory neurons fire more intensively as compared to excitatory neurons during the SimPs, while the spiking activity in excitatory neurons increases with successive SimPs (Figure 3B2). The peak of the activity in the whole network was reached during the ID ( Figure 3B3) and its spectrogram clearly revealed two main components corresponding to the different activity in excitatory and inhibitory neurons that fire at about 15 and 60 Hz, respectively ( Figure 3B4). This pattern of activity in the two neuronal populations is consistent with experimental observations (Ziburkus et al., 2006). While both excitatory and inhibitory neurons at the focus were activated upon the initial stimulation (representative traces in Figures 3C1,C2), neurons outside the focus were recruited into the propagating ID with some delay (representative traces in Figures 3C3,C4).

ID GENERATION THRESHOLD
We consider as ID threshold the number of NMDA pulses that are needed to evoke an ID. This value is constant for a given slice , but it can vary for different slices.
Simulations with different parameters (see "Methods") showed that an ID could be generated in average by five SimPs and in more than 25% of cases no IDs could be evoked regardless the number of applied SimPs (Figure 3D, n = 250 simulations). The successive SimPs induced excitatory responses at the focus with increasing amplitude (Figure 3E1). The excitatory and inhibitory neurons outside the focus were not directly activated by the SimPs and increased their firing activity simultaneously, but with a marked delay (Figures 3E1,E2).

DYNAMICS OF EXCITATION AND INHIBITION AT THE FOCUS EXPLAIN ID GENERATION
We next investigated the interplay between excitation and inhibition in the genesis of the ID. We compared the simulations which successfully evoked an ID with those that failed to evoke an ID (in the different simulations excitatory and inhibitory neurons were randomly located within or outside the focus while maintaining their total number). For the cases in which an ID was successfully generated, we find that the ratios between the number of excitatory and inhibitory neurons, the average inhibitory and excitatory currents during the first SimP and the firing rate of inhibitory and excitatory neurons in the same time interval were lower compared to cases where an ID was not successfully generated (Figures 4A1-A3). We next tested whether a different strength in either excitation or inhibition at the focus changed the efficacy of the SimP in generating an ID. We analyzed the time course of excitation and inhibition at the focus and the average firing rate in the whole network. We examined three sets of network parameters chosen at random, but leading to different ID thresholds, i.e., no ID generation, high ID threshold (five SimPs) and low ID threshold (three SimPs) (Figures 4B-D respectively). In all cases, excitatory and inhibitory drive increased the firing rate ( Figures 4B2,B3,C2,C3,D2,D3). A detailed analysis of the dynamics revealed that after each SimP both excitation and inhibition were strongly but transiently activated (Figures 4B1,C1,D1). An important difference is that, in contrast to inhibition, excitation failed to recover the initial basal conditions, including the simulations in which no ID is generated ( Figure 4B). An additional striking difference between the three examples is the maximal inhibition level provided by the inhibitory neurons (the dynamic range). Inhibition reached its highest value in the high ID threshold condition. These results support the view that inhibition strength is a critical factor for focal ID onset. Notably, excitation rose faster than inhibition (slope > 1) driving the growth in firing rate forward. However, the ID occurred only after inhibition reached its maximal value (all inhibitory neurons were active). Therefore the ratio of excitatory versus inhibitory drive and the limiting dynamic range of inhibition are the two critical factors in ID generation. As a summary of results obtained, we report the distribution of points in the excitation-inhibition plane at the focus during the first seven SimPs in 250 Monte-Carlo simulations for the cases that evoked or failed to evoke an ID (normalized by the total area; Figures 4E1,E2). When inhibition at the focus reached high values, no IDs were generated and the ratio between excitation and inhibition remained low. This stands in contrast to the cases which lead to IDs, further supporting the notion that the relationship between excitation and inhibition determines not only the threshold for ID generation, but also whether or not an ID could be evoked. Data obtained from 250 runs also showed a clear correlation between the ID threshold and the average excitation and inhibition in the network during the first SimP (Figures 4F1,F2). This indicates that the overall response of the network, in terms of excitation and inhibition levels, is a good predictor of ID threshold: an increased excitation results in the lowering of the ID threshold and an opposite relationship holds for inhibition.

ASTROCYTE-TO-NEURON SIGNALING DECREASES THE ID THRESHOLD
The model of the single astrocyte (see "Methods") was incorporated into the network to test how astrocytes may affect ID threshold. Specifically, 400 astrocytes were added to the network model in a parallel 2D sheet of cells (see "Methods"). Astrocytes provide an excitatory feedback to neuronal activity in a Ca 2+dependent way ( Figure 5A). As illustrated in Figure 5B, in the presence of the astrocyte feedback signal, the ID was evoked by two SimPs, while in its absence a more intense stimulation of neurons was necessary. As illustrated in Figures 5F1,F2 the Ca 2+ change from a representative astrocyte at the focus was observed to follow rapidly the spiking activity in neurons (Figures 5C1,F1), and astrocytic glutamate release occurred upon the second SimP ( Figure 5C2). The average astrocytic Ca 2+ follow the neuronal activity (example in Figure 5D1) while the average glutamate release occurs transiently ( Figure 5D2). The Ca 2+ change and the release of glutamate from astrocytes outside the focus failed to affect focal ID threshold. The results from 250 Monte-Carlo runs show that the ID threshold was lowered after including the astrocytic feedback signal to neurons (Figure 5E). Once the ID was fully evolved, both the Ca 2+ change and the release of glutamate from astrocytes within and outside of the focus did not differ significantly (Figures 5F1,F2). However, the initially dominant activity of astrocytes at the focus was followed by an activity of the astrocytes outside the focus that became dominant immediately after the ID onset. These results are consistent with those from slice experiments which showed that when Ca 2+ elevation in astrocytes from the focus were selectively blocked (by the Ca 2+ chelator BAPTA) or stimulated (by TFLLR, a peptide agonist of thrombin PAR-1 receptors), the threshold of ID generation increased or decreased, respectively .

DOES AN ASTROCYTE INHIBITORY FEEDBACK SIGNAL TO NEURONS AFFECT ID THRESHOLD?
The ID threshold is mainly affected by the interplay between excitation and inhibition. Indeed, as we reported above, ID threshold can be increased by increasing the overall value of the inhibitory activity. The bar graph in Figure 6A reports the results from 250 simulations in different simulation settings, with and without an astrocytic contribution (mean and errors represent the results of a Poisson fit to the ID threshold distributions), while Figure 6B is the cumulative sum of the ID threshold distributions corresponding to the analyzed cases. Blue and red bars show that the ID threshold can be increased by increasing the overall strength of inhibitory connections (in this case from 0.01 to 0.015) in a purely neuronal network (no astrocytes). With higher inhibition, the simulated stimulation failed to induce an ID in 40% of simulations ( Figure 6B). As already shown, the introduction of an astrocytic excitatory feedback lowers the ID threshold (green bar) and decreases the number of failures to about 10%. In slice experiments we observed that the inhibition of Ca 2+ signals in astrocytes at the focus, but not outside the focus, increased the threshold of ID generation. These observations were fully reproduced in the computational model (dark blue and magenta bars) without further manipulations of the model over the results from the previous section. Astrocytic excitatory feedback in a network with stronger inhibitory connections (0.015 as for the red bar) brings back to baseline the ID threshold (yellow).
To explore other factors that may affect ID threshold, we considered the possibility that the activation of astrocytes, or of a subpopulation of astrocytes, results in a release of GABA that can lead to an overall increase of the inhibition strength in the neuronal network. Astrocytes can, indeed, release GABA ( Kozlov et al., 2006;Lee et al., 2011;Le Meur et al., 2012). As expected, when we included an astrocytic GABA release in the model, the threshold for ID generation increased (Figure 6, dark green). Note that this was an inhibitory-only feedback involving only GABA release and no glutamate. Surprisingly, however, the threshold for ID generation did not rise over the baseline condition with no astrocytic feedback (blue). A possible explanation for this could be the synchronizing action of an inhibitory GABA signal. Alternatively, the inhibitory feedback signal from astrocytes could be more effective in suppressing inhibitory than excitatory neurons. This action may ultimately generate a new level of complexity in the mechanism that governs the inhibition/excitation balance in the neuron-astrocyte network.

DISCUSSION
Increasing experimental evidence highlights the physiological significance of the tripartite synapse in which the astrocyte senses neurotransmitter release and, in turn, releases through a Ca 2+dependent mechanism gliotransmitters that have feedback modulatory actions on neurons. A number of recent studies in vitro and in vivo showed that the release of glutamate from astrocytes can, indeed, control both the basal excitability of neurons and some forms of long-term potentiation of synaptic strength (Serrano et al., 2006;Jourdain et al., 2007;Navarrete and Araque, 2010;Santello et al., 2011;Min and Nevian, 2012) and long-term depression (Zhang et al., 2003;Serrano et al., 2006;Han et al., 2012;Min and Nevian, 2012). The contribution of gliotransmission to brain dysfunctions remains, however, poorly understood. A model composed of a network of interacting neurons and astrocytes represents a useful tool in which the spatial-temporal features of focal seizure generation observed in slice models can be replicated and new mechanistic hypotheses can be tested. Over the last 10 years, different models have been advanced to describe the Ca 2+ dynamics of astrocytes in response to neuronal signals Jung, 2003, 2007;Silchenko and Tass, 2008;Di Garbo, 2009). These biophysical approaches described not only the astrocytic Ca 2+ response (Li and Rinzel, 1994), but also the possible feedback to neurons. More recently, the possible contribution of astrocytes in events related to the plasticity of synaptic transmission were also included in models (Nadkarni et al., 2008;De Pittà et al., 2011;Wade et al., 2011). While biophysical models are very useful to simulate basic units, like the tripartite synapse, they are hardly suitable for large scale simulations. In contrast, simplified models that include only the basic features of neuron-astrocyte interactions (Postnov et al., 2009) appear more appropriate to describe network dynamics and to investigate the role of astrocytes in epilepsy (Amiri et al., 2012). In our model we did not include any distinct biophysical features that characterize the physiological actions of either neurons or astrocytes. We rather describe the activity of a single astrocyte in terms of the specific input-output signals with which astrocyte and neurons interact. This simplified astrocyte model was then embed in a neuronal network model of IDs based on the Izhikevich's single neuron model.
The slow kinetics of the astrocyte Ca 2+ response to neuronal activity in the model reflect those of the mGluR-mediated Ca 2+ elevations that were evoked in astrocytes by axonal afferent stimulation in young rat hippocampal slices (Porter and McCarthy, 1995;Pasti et al., 1997;Perea and Araque, 2005). Indeed, the intracellular Ca 2+ variations in astrocytes depend primarily, although not exclusively, on activation of metabotropic receptors, phospholipase C-dependent inositol(1,4,5)-trisphosphate (IP3) production and, finally, stimulation of Ca 2+ release from IP3sensitive internal Ca 2+ stores (Kawabata et al., 1996). Glutamate release at the synapse triggers a Ca 2+ response in astrocytes that increases in both amplitude and frequency of oscillations according to increased levels of the neuronal activity (Pasti et al., 1997). These Ca 2+ changes trigger a SNARE-dependent exocytosis of glutamate that signals back to affect the excitability of neurons (Araque et al., 2000;Pasti et al., 2001;Parpura et al., 2004). Accordingly, in the model we reproduced the most essential features of glutamate release in response to Ca 2+ elevations in astrocytes. The release of glutamate is pulsatile and depends on the frequency of Ca 2+ oscillations (Pasti et al., 2001), while its efficacy is controlled by the amplitude of the Ca 2+ increase (Parpura and Haydon, 2000). In addition, a steady state Ca 2+ elevation may trigger only a single episode of release (Pasti et al., 2001).
Some approximations were applied to describe two features that characterize astrocyte signaling in our model. Firstly, we restricted our analysis to somatic Ca 2+ signals. These Ca 2+ increases can not be intended to fully represent the synapse-toastrocyte signaling occurring fundamentally at the proximal and the distal processes that are in contact with the synapse. Indeed, somatic Ca 2+ increases exhibit a lower frequency and slower kinetics with respect to those at the processes (Di Castro et al., 2011;Panatier et al., 2011). While these recent studies also showed that Ca 2+ elevations at the astrocytic processes can have a distinct functional role, it is noteworthy that the Ca 2+ elevation at the soma may represent a response that integrates the signals from the processes where astrocytes sense neurotransmitter release. Accordingly, Ca 2+ signals at the soma may adequately reflect the overall firing activity of surrounding neurons. Amplitude, frequency and general pattern of somatic Ca 2+ changes are, indeed, observed to vary according to different levels of neuronal activity (Pasti et al., 1997;Porter and McCarthy, 1997). Secondly, while glutamate (Parpura et al., 1994;Pasti et al., 1997;Bezzi et al., 1998), D-serine (Mothet et al., 2005;Henneberger et al., 2010), ATP (Arcuino et al., 2002;Serrano et al., 2006;Bowser and Khakh, 2007), and GABA (Kozlov et al., 2006;Lee et al., 2011) [for a review see Haydon and Carmignoto (2006)] are the main gliotransmitters mediating astrocyte-to-neuron signaling, in our model we fundamentally focused on glutamate because a large body of information is available about its modulatory action on both basal synaptic transmission (Fellin et al., 2004;Di Castro et al., 2011;Panatier et al., 2011) and long-term plasticity (Zhang et al., 2003;Panatier et al., 2006;Serrano et al., 2006;Jourdain et al., 2007;Navarrete and Araque, 2010;Santello et al., 2011;Han et al., 2012;Min and Nevian, 2012). In addition, the contribution of astrocytic glutamate in some forms of long-term potentiation of synaptic transmission has been recently confirmed in in vivo experiments (Takata et al., 2011;Navarrete et al., 2012). The potential role in focal seizure generation of ATP, D-serine and GABA will be the subject of future investigations. It is worth mentioning here that D-serine, and not glycine, is most likely the physiological co-agonist of the synaptic NMDA receptor in the brain (Mothet et al., 2000;Panatier et al., 2006;Fossat et al., 2012;Papouin et al., 2012). Given that D-serine is mainly, although not exclusively (Ding et al., 2011), synthesized in astrocytes and released through a Ca 2+ -dependent mechanism (Wolosker et al., 1999;Wolosker, 2011), astrocytic D-serine may cooperate with glutamate to enhance NMDA receptor openings and through this action favor neuronal excitability ultimately promoting epileptic discharges.
The pathological increase in brain network excitability that eventually leads to focal seizure generation is believed to derive from the activity of excitatory and inhibitory neurons as well as of astrocytes. The cellular events that favor or oppose seizure initiation and propagation remain, however, poorly defined. Our model offers the opportunity to study ID generation in simulated networks composed by either only neurons or interactive astrocytes and neurons. The results that we obtained are summarized in Figure 6 and can be, in our opinion, useful to understand how distinct signaling pathways may govern focal ID generation. Figure 6 plots the average threshold for ID generation in the different conditions (mean ± SD, Figure 6A) and the cumulative sum of the threshold histograms ( Figure 6B) showing failures. We found that in a network composed exclusively of neurons an ID can be generated by applying an intense stimulation of a group of neurons. The introduction of astrocytes into the network lowered ID threshold, while the inhibition of astrocyte signaling to neurons within, but not outside the focus, increased ID threshold. These results are fully consistent with those obtained in slice experiments  and demonstrate that focal IDs can be faithfully reproduced in our computational model. Accordingly, our model can be used to make predictions on the distinct contribution of different signaling pathways to ID generation. We present here some results regarding inhibitory signaling pathways. The ID threshold was observed to increase upon procedures that increase the strength of inhibition onto the principal neurons. This was achieved by either increasing the strength of the inhibitory transmission or by including in the model a distinct inhibitory feedback signal from astrocytes to neurons via GABA A receptors. These observations will be useful when addressing in future slice experiments the role of inhibitory signaling in ID generation.