We build a random network of N = 1, 000 adaptive exponential integrate-and-fire neurons (Brette and Gerstner, 2005) with probability p for the formation of connections among them equal to 0.1. The network consists of 80% excitatory and 20% inhibitory neurons (Noback et al., 2005). The dynamics of each neuron i, i = 1, …, N in the network is given by the set of equations

The membrane potential V_i and adaptation current w_i represent the state of each neuron i. The capacitance membrane C_m is set to C_m = 200 pF, the leak conductance to g_L = 12 nS, the resting potential to E_L = −70 mV, the slope factor to Δ_T = 2.0 mV and the spike threshold to V_T = −50 mV. The adaptation current depends on the adaptation time constant τ_w = 300 ms and the level of subthreshold adaptation a_i that is randomly distributed in the interval [0.19, 0.21] nS. We consider the injection of current I_i to each neuron i in terms of the relative rheobase current r_i = I_i/I_rheobase (Naud et al., 2008). The rheobase is the minimum amplitude of the applied current to generate a single or successive firings. The application of this constant current allows neurons to change their potentials from resting potentials to spikes. The value of the rheobase depends on the neuron parameters. The external current arriving at neuron i is represented by Γ_i. We consider the external current according to a SCP with amplitude A_I and time duration T_I. The random connections in the network are described by the binary adjacency matrix M_ij with entries either equal to 1 when there is a connection from i to j or 0 in the absence of such a connection. g_i is the synaptic conductance, τ_s the synaptic time constant, and V_REV the synaptic reversal potential. We consider τ_s = 2.728 ms, V_REV = 0mV for excitatory synapses, and V_REV = −80 mV for inhibitory synapses. The synaptic conductance decays exponential with a synaptic time constant τ_s. When the membrane potential of neuron i is above the threshold V_i > V_thres (Naud et al., 2008), the state variable is updated by the rule

where g_s assumes the value of g_exc when neuron i is excitatory (i ≤ 0.8N) and g_inh when neuron i is inhibitory (i > 0.8N). In this work, we study the parameter space (g_exc, g_inh) and consider a relative inhibitory synaptic conductance g = g_inh/g_exc. We consider parameter values in which the individual uncoupled neurons perform spike activities. The initial values of V and w are randomly distributed in the interval [−70, −50] mV and [0, 70] pA, respectively. The initial g_i value is equal to 0.

The synchronous behavior in the network can be identified by means of the complex phase order parameter (Kuramoto, 1984)

where R(t) and Φ(t) are the amplitude and angle of a centroid phase vector over time, respectively. The phase of neuron j is obtained by means of

where t_{j, m} corresponds to the time of the m−th spike of neuron j (t_{j, m} < t < t_{j, m+1}) (Rosenblum et al., 1996, 1997). We consider that the spike occurs for V_j > V_thres. R(t) is equal to 0 for fully desynchronized and 1 for fully synchronized patterns, respectively.

We calculate the time-average order parameter R¯ (Batista et al., 2017) given by

where t_fin − t_ini is the time window. We consider t_fin = 200s and t_ini = 180s.

We monitor the instantaneous synaptic conductances arriving at each neuron i through

The instantaneous synaptic input changes over time due to the excitatory and inhibitory inputs received by neuron i. The average instantaneous synaptic conductances is given by

The m−th inter-spike interval ISIim is defined as the difference between two consecutive spikes of neuron i,

where tim is the time of the m−th spike of neuron i.

Using the mean value of ISI_i, ISI¯i, and its standard deviation, σ_ISIi, we calculate the coefficient of variation (CV)

The average of CV (CV¯) is then obtained through

Finally, we use CV¯ to identify spike (when CV¯<0.5) and burst fire patterns (when CV¯≥0.5) (Borges et al., 2017; Protachevicz et al., 2018).

The instantaneous firing-rate in intervals of t_step = 1ms is given by

where t_i is the firing time of neuron i in the time interval (t ≤ t_i ≤ t + 1) ms. This measure allows to obtain the instantaneous population activity in the network. The mean firing-rate can then be calculated by means of

where ISI¯ is the average ISI obtained over all N neurons in the network, that is ISI¯=1N∑i=1NISI¯i.

The balance between excitation and inhibition generates an asynchronous activity in the network (Lundqvist et al., 2010; Ostojic, 2014). However, for the unbalanced network we observe synchronized spikes and bursts. Figures 1A–C show the time-average order parameter (R¯), the mean coefficient of variation (CV¯) and the mean firing-rate (F¯), respectively, for the parameter space (g, r), where g is the ratio between inhibitory (g_inh) and excitatory (g_exc) synaptic conductances, and r the relative rheobase current. For g_exc = 0.4nS and g > 6, we observe that R¯<0.5 and that CV¯<0.5, corresponding to desynchronized spikes. In Figure 1D, we see the raster plot and membrane potential for 2 neurons in the network with a desynchronized spike-pattern for g = 5.5 and r = 2 (blue triangle). For g = 4 and r = 1.5 (magenta square), the dynamics exhibits synchronized spikes (Figure 1E), as a result of setting R¯>0.9 and CV¯<0.5. Figure 1F shows synchronized bursts of activity for g = 2.5 and r = 2 (green circle), where R¯>0.9 and CV¯≥0.5. Within this framework, we have verified the existence of transitions from desynchronized spikes to synchronized bursting activities without significant changes in the mean firing-rate.

The appearance of synchronous behavior cannot only be related to the decrease of the inhibitory synaptic strength, but also to a loss of inhibitory neurons. In particular, we show this in Figure 1G which illustrates a network composed of excitatory (red) and inhibitory (blue) neurons, where some inhibitory neurons were removed (dashed circles). In Figure 1H, we see that the synchronous behavior depends on g and the percentage of removed inhibitory neurons. Figure 1I shows the transition from spiking dynamics (CV¯<0.5) to bursting dynamics (CV¯≥0.5), and Figure 1J shows the instantaneous firing-rate F(t). For g = 2.9 and g_exc = 0.4nS (green dashed line), the transition to synchronized bursts occurs when 10% of inhibitory neurons are removed from the network, and as a consequence F(t) reaches the maximum value of 0.2.

Concluding, alterations in the inhibitory synaptic strength or in the number of inhibitory neurons can induce transition to synchronous patterns. Wang et al. (2017) presented results where synchronization transition occurs as a result of small changes in the topology of the network, whereas here, we study transitions caused due to changes in the inhibitory synaptic strength and the emergence of a bistable regime.

Next, we analyse synchronization in the parameter space (g, g_exc). In particular, Figure 2A shows R¯ with values depicted in the color bar. The black region corresponds to desynchronized spike activity, while the remaining colored regions are associated with burst activities. The white region represents the bistable regime, where desynchronized spikes or synchronized bursts are possible depending on the initial conditions. In the bistable regime, decreasing g_exc (backward direction), R¯ is higher than increasing g_exc (forward direction), as shown in Figure 2B for g = 3, r = 2, and g_exc = [0.35, 0.45] nS (green dashed line in Figure 2A). We identify bistability (white region) in the parameter space when the condition R¯backward-R¯forward>0.4 is fulfilled. The raster plot and instantaneous synaptic input for desynchronized spikes (blue circle) and synchronized bursts (red square) are shown in Figures 2C,D, respectively. When the dynamics on the random network is characterized by desynchronized spikes, the instantaneous synaptic inputs exhibit I_syn(t) ≈ 50pA. For synchronized bursts, I_syn(t) ≈ 0 when a large number of neurons in the network are silent (i.e., not firing), and I_syn(t) > 200 pA during synchronous firing activities. In Figure 2E, we compute the probability of occurrence of excessively high synchronicity within the bistable regime. We observe a small synchronization probability value in the bistable region. This result has a biological importance due to the fact that the seizure state is a relatively small probability event compared with the normal state. DaQing et al. (2017) showed that noise can regulate seizure dynamics in partial epilepsy. Figure 2F displays R¯×gexc for Gaussian noise with mean 0 and standard deviation σ_noise equal to 25 pA and 250 pA. We verify that the bistable region decreases when the noise level increases.

In the bistable regime, we investigate the evolution of a trajectory for a finite time interval in the phase space (w_i, V_i) and the time evolution of w_i shown in Figure 3 for i = 1, where the gray regions correspond to dV_i/dt < 0. The boundary between the gray and white regions (black line) is given by dV_i/dt = 0, the V_i-nullcline (Naud et al., 2008). During spiking activity, the trajectory (see Figure 3A) and time evolution of w_i (see Figure 3B) do not cross the V_i-nullcline. For bursting activities (see Figures 3C,D), we observe that w_i lies in the region enclosed by the V_i-nullcline. The emergence of the bistable behavior is related to changes in the V_i-nullcline caused by the variation of I_syn.

Here, following a similar idea as in Antonopoulos (2016), we investigate the effect of the application of SCP on the bistable regime. We apply SCP considering different values of A_I, T_I, and number of removed inhibitory neurons. The SCP is immediately switched off after T_I and the analysis of the effect on the dynamical behavior is started.

Initially, we apply SCP to all neurons in the network with parameter values in the bistable regime with desynchronous behavior (white region in Figure 2A). Figure 4A displays the time (in color scale) that the neurons show a synchronized pattern after the application of SCP. In the black region, we see that SCP does not change the dynamical behavior, namely the neurons remain in a regime of desynchronized behavior. The yellow region depicts the values of T_I and A_I of the SCP that induce a change in the behavior of the neurons from desynchronized spikes to synchronized bursts. Picking up one point close to the border of the black and blue regions (white circle), we see that the instantaneous firing-rate (F(t)) of Equation 11 (see Figure 4B, blue line) exhibits low-amplitude oscillations corresponding to desynchronized spikes. For T_I and A_I values in the yellow region, F(t) (see Figure 4C, red line) exhibits a high-amplitude oscillation after the application of SCP, corresponding to synchronized bursts. For sufficiently large amplitudes, the change in the behavior induced by SCP does not depend on time. Importantly, perturbations with small amplitudes applied for short times is a sufficient condition for the induction of synchronous burst activity in the bistable regime. Therefore, our results suggest that even small excitatory stimuli in a random neural network arriving from other parts might be sufficient for the initiation of excessively high neural synchronization, related to the onset of epileptic seizures. Thus, further work on other neural networks that resemble brain activity might provide more insights on epileptogenesis.

Similarly, we apply SCP when the neurons in the network show synchronized bursts of firing activity in the bistable regime. Here, we aim to suppressing the synchronous behavior by means of applying SCP. We consider SCP with positive and negative amplitudes applied to 10% of the neurons in the random network. Figure 5A shows how long the bursts remain synchronized after SCP is switched off (color bar). We verified that both negative and positive amplitudes exhibit regions where the synchronous behaviors are suppressed, namely there is a transition from synchronized bursts to desynchronized spikes. In addition, for T_I>0.4 s and considering the absolute value of the amplitudes, the transition occurs for positive values with smaller amplitudes than for negative values. In Figure 5B, we show the dependence of the percentage of the perturbed neurons by the stimulus on the time the neurons remain in the bursting synchronous regime. The black region represents parameters for which the dynamics on the network does not remain synchronous, and therefore, synchronization is suppressed. In this figure, T_I = 1s. These results allow us to conclude that desynchronous behavior is achieved for A_I > 15pA and for at least 10% of the perturbed neurons.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.