We build a neuronal network composed of two areas, where each one has a thousand of adaptive exponential integrate-and-fire neurons (N = 1, 000) (Brette and Gerstner, 2005). Each area has a fraction of excitatory (P_exc = 0.8) and inhibitory (P_inh = 0.2) neurons (Noback et al., 2005; di Volo et al., 2019). In each area, the neurons are randomly coupled by means of excitatory and inhibitory connections. The connection is excitatory (inhibitory) when it occurs from an excitatory (inhibitory) neuron. Inside of each area, the probabilities of connections in the same neuronal populations (excitatory or inhibitory) are given by p_ee = 0.05 and p_ii = 0.2, while between different neuronal populations by p_ei = p_ie = 0.05 (di Volo et al., 2019). Between the areas, the probabilities are given by peeA=0.01 (from excitatory to excitatory neurons), peiA=0.05 (from excitatory to inhibitory neurons), piiA=0.10 (from inhibitory to inhibitory neurons), and pieA=0.05 (from inhibitory to excitatory neurons). Figure 1 shows how the probabilities are distributed in a connection matrix, where k and j correspond to the presynaptic and postsynaptic neurons, respectively. In Figure 1, j, k ∈ [1, 1000] correspond to the neurons in Area 1 and j, k ∈ [1001, 2000] to the neurons in Area 2. For the unidirectional configuration, the connections given by k = [1001, 2000] and j = [1, 1000] are not considered. For both unidirectional and bidirectional configuration, we consider only excitatory or inhibitory connections between the areas in each case.

With regard to the coupling intensities, each one is associated with a probability of connection and denoted by g_ee, g_ii, g_ei, g_ie, geeA, geiA, giiA, and gieA. For instance, g_ei and g_ie are related to the connections with the probabilities p_ei and p_ie, respectively.

The cortex is mainly constituted by excitatory pyramidal neurons and inhibitory interneurons (Atencio and Schreiner, 2008). Excitatory neurons have a relatively lower firing rate than inhibitory ones (Wilson et al., 1994; Inawashiro et al., 1999; Baeg et al., 2001). In the mammalian cortex, excitatory neurons show regular spike (RS), while inhibitory neurons exhibit fast spike (FS) activities (Neske et al., 2015; Wang et al., 2016). In addition, while inhibitory neurons exhibit a negligible adaptation, excitatory neurons show an adaptation mechanism in their firings (Foehring et al., 1991; Mancilla et al., 1998; Hensch and Fagiolini, 2004; Destexhe, 2009; Masia et al., 2018; Borges et al., 2020). The adaptive exponential integrate-and-fire (AEIF) model is able to mimic these different firing patterns, including RS and FS (di Volo et al., 2019). In this work, the dynamics of each neuron j (j = 1, …, N) in the network is given by

The membrane potential V_j and adaptation current w_j represent the state of each neuron j. The capacitance membrane is set to C_m = 200 pF, the leak conductance to g_L = 12 nS, the leak reversal potential to E_L = −70 mV, the slope factor to Δ_T = 2 mV, and the spike threshold to V_T = −50 mV. We consider the injection of current I = 270 pA, which is the intensity above the rheobase current. The application of this constant current allows that the neurons change their potentials from resting potentials to spikes. The level of the subthreshold and triggered adaptation are represented by a_j and b_j, respectively. We consider inhibitory neurons of fast spiking activities without adaptation (a_j = 0 and b_j = 0) and excitatory neurons of regular spiking with adaptation mechanisms (a_j = [1.9, 2.1] nS and b_j = 70 pA). Neuronal adaptation corresponds to the capacity of the neuronal membrane in adapting to its excitability according to the past neuronal activity. A sub- and a triggered-threshold adaptation mechanism can be associate with the parameters a_j and b_j, respectively. The adaptation current also depends on the adaptation time constant τ_w = 300 ms. g_j represents the synaptic conductance of each neuron j with an exponential decay associated with the synaptic time constant τ_s = 2.728 ms. The connections from excitatory and inhibitory neurons are related to the excitatory and inhibitory matrix, M⃗exc and M⃗inh, where each matrix element is identified as Mjkexc and Mjkinh, respectively. A matrix element is equal to 1 when there is a connection from k to j or 0 in the absence of a connection. The excitatory and inhibitory elements of the matrix are associated with the red and blue dots in Figure 1. The chemical current input Ijchem arriving on each neuron j is defined by the expression

where the excitatory and inhibitory currents are given by

and

where Ijxy and Ijxy,A are associated with the excitatory (x=e) or inhibitory connectivity (x=i) arriving at the excitatory (y=e) or inhibitory neurons (y=i). The type of synapse depends on the synaptic reversal potential V_REV. We consider the VREVexc=0 mV for excitatory and VREVinh=-80 mV for inhibitory synapses. N_T is the total number of neurons in the network. When the membrane potential of the neuron j is above the threshold V_j > V_thres (Naud et al., 2008), the state variable is updated by the rule

In our simulations, we consider V_r = −58 mV. The value of b_j depends whether the neuron j is excitatory or inhibitory. Each synaptic current is related to the respective conductance g_s. Inside of each area, g_s is equal to g_ee for synapses between excitatory neurons, g_ei for synapses from excitatory to inhibitory neurons, g_ii for synapses between inhibitory neurons, and g_ie for synapses from inhibitory to excitatory neurons. Between different areas we include the superscript “A.” The time delay in the conductance is d_exc = 1.5 ms for excitatory connections and d_inh = 0.8 ms for inhibitory ones (Borges et al., 2020).

Table 1 gives the values of the parameters used in our simulations. We consider g_ee = 0.5 nS, g_ii = 2 nS, and g_ie = 1.5 nS. The areas exhibit synchronous and desynchronous behaviour when uncoupled between them for g_ei = 1 nS and g_ei = 2 nS, respectively. The area 1 is considered synchronised and desynchronised when uncoupled, while the area 2 is always considered desynchronised when uncoupled.

The initial values of V_j are randomly distributed in the interval [−70, −50] mV for all neurons. The initial values of w_j are randomly distributed in the interval [0, 300] pA for excitatory neurons and equal to 0 for inhibitory ones. The initial value of g_j is equal to 0 for all neurons. To solve the delayed differential equations, we consider that the excitatory and inhibitory neurons in the network are not spiking before the beginning of the simulation (t = 0). To integrate the set of ordinary differential equations, we use the 4th Runge-Kutta method with the time-step of integration equal to 10⁻² ms.

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

where R(t) is the amplitude of a centroid phase vector over time. The phase of each neuron j is obtained through

where t_{j, m} corresponds to the time of the m−th spike of the neuron j (t_j,m < t < t_{j, m+1}) (Rosenblum et al., 1997). We consider that the spike occurs when V_j > V_thres. The value of R(t) is equal to 0 for completely desynchronised behaviour and equal to 1 for fully synchronised patterns.

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

where t_fin−t_ini is the time window with t_fin = 100 s and t_ini = 20 s.

The order parameter for each area is given by

where A denotes the area number. The mean value of R_A(t) (R¯A) is computed by Equation (8). The resultant phase angle of each area A is defined as

The real RAx and complex RAy components of the order parameter can be described as

and

respectively. Θ_A(t) evolves in the counter-clockwise direction, since each individual neuron evolves in this direction.

We define a relative phase angle for each area ΘA′(t) (Varela et al., 2001) as

The phase of the area 1 changes over time and its relative value Θ1′(t) is equal to 0. For the area 2, the relative phase angle Θ2′(t) can change over time. We calculate the mean value of the relative phase angle of the area 2 by means of

To compute the predominant rotation direction of the area 2, we consider the first-order derivative of their relative phase angle, which corresponds to the instantaneous velocity of the relative phase angle (rad/s), that is given by

The mean value of the relative velocity of the area 2 is obtained via

The value is close to 0 for non-preponderant direction of rotation and positive (negative) for predominant counter-clockwise (clockwise) direction.

The variability of Θ2′ is given by

Small and high deviations are given by σΘA′≈0 and σΘA′>0.5, respectively. Figure 2 shows a schematic representation of the mean order parameter of the area 1 and area 2 for Figure 2A desynchronised patterns of both areas out-of-phase, Figure 2B in-phase (Θ2′¯=0), Figure 2C anti-phase (Θ2′¯=π), and Figure 2D shift-phase synchronisation (Θ2′¯≈4π/3). In Figure 2A, the clockwise (blue) and counter-clockwise (red) arrows indicate that the relative phase angle of area 2 change over time. The amplitude and direction of the relative phase angle of the area A can be described by R¯AeiΘA′. In our simulations, we observe results in which desynchronised activities can be related to high variability of the relative phase angle of the area 2. When the areas are synchronised, the variability of the relative phase angle can go to 0. The mean value of relative phase angle is efficient for small variability of the relative phase angle.

We analyse a neuronal network separated into two areas (sender-receiver) coupled by means of the excitatory connections. Figure 3 displays a schematic representation of the sender area 1 to receiver area 2 via excitatory connections, that are related to the geiA and geeA conductances.

Firstly, we consider the case in which the neurons in the sender and receiver area are desynchronised. Through the time evolution of Θ2′, we verify that, depending on the conductance values, it can occur a relative counter-clockwise rotation, clockwise rotation, or neither of them in the area 2. We observe that geeA contributes to generate positive Θ∙2′¯, while geiA to negative one. Both areas show order parameters with small values, i.e., the neurons remain desynchronised.

Secondly, we consider that the neurons in the sender area are synchronised while the neurons in the receiver area are initially desynchronised. In Figure 4A, the parameter space geeA×geiA exhibits values of Θ∙2′¯ in [−0.1, 0.1] rad/s, where the Θ∙2′¯ values approximately or less than −0.1 are in a small black diagonal region. For a large set of parameters, Θ∙2′¯ is close to zero. Figure 4B displays σΘ2′ values about 0, except for a small red region where the values are greater than 0.5. We compute the mean relative phase angle of the area 2, as shown in Figure 4C. We verify that geeA leads the area 2 to a value of the phase angle equal to 0, while geiA leads to a shifted phase angle, Θ2′¯≈4π/3. The stabilisation of the phase angle is associated with the synchronisation of the area 2, as shown in Figure 4D. We see a large region in which R¯2 is close to 1, meaning that the neurons in the area 2 are synchronised.

Figure 5 displays the raster plot (top), relative phase angle (middle), and order parameter (bottom) for (a) geeA=2.7 nS and geiA=0.6 nS, (b) geeA=1.5 nS and geiA=3 nS, and (c) geeA=0.3 nS and geiA=5.4 nS, that are indicated in Figure 4 through the circle, square, and triangle symbols, respectively. In Figure 5A, we see that Θ2′ is closed to the sender phase angle and the receiver area has synchronised neurons due to a high geeA conductance. Figure 5B shows that for a combination of geeA and geiA, neurons in the area 2 are not synchronised. In Figure 5C, we observe that due to high geiA conductance, neurons in the area 2 have a shift-phase synchronisation, corresponding to Θ2′¯≈4π/3, as indicated in Figure 2D.

Figure 6 displays a schematic representation of the sender area to the receiver area via the giiA and gieA conductances.

For desynchronised neurons in the sender and receiver areas, increasing giiA, we see that Θ∙2′¯ is positive for small values of gieA, as shown in Figure 7A. On the other hand, gieA contribute to negative values of Θ∙2′¯. Figure 7B exhibits a high variability (σΘ2′>0.5) of Θ2′. In Figures 7C,D, we compute the mean order parameters for the areas 1 and 2, respectively. The neurons are desynchronised in the area 1, while in the area 2, we see a small region in the parameter space where there is synchronous behaviour (R¯2>0.8). The increase of R¯2 is due to the giiA parameter, that induces the inhibition of the inhibitory neurons of the receiver area (circle). For high gieA values (triangle), a large number of excitatory neurons do not fire. A similar result was reported by Zhou et al. (2010), where silent activities of excitatory neurons were observed due to a strong inhibition. Urban-Ciecko et al. (2015) found which inhibition can silence excitatory synapses in the neocortex. Pals et al. (2020) demonstrated that activity-silence maintenance can be related to a working memory process. The silence of neurons has received great attention in the last years (Mochol et al., 2015; Wiegert et al., 2015; Barbosa et al., 2020; Xu et al., 2020).

Figure 8 displays the raster plot (top), time evolution of Θ2′ (middle), and time evolution of the order parameter (bottom) for the values of giiA and gieA indicated by circle, square, and triangle symbols in Figure 7. For giiA=3.6 nS and gieA=0.3 nS (Figure 8A), the neurons in the area 2 synchronise and Θ2′ denotes a relative counter-clockwise rotation. For giiA=2 nS and gieA=1 nS (Figure 8B), the neurons in the area 2 are desynchronised and Θ2′ denotes a relative clockwise rotation. In Figure 8C, there is no synchronous behaviour and we see a large number of excitatory neurons that remain in silence for a long time.

We also consider the case in which the sender area has synchronised neurons. Figure 9A shows that the Θ∙2′¯ values are positive due to the giiA conductance with small values of gieA for small values of gieA. This result is similar to the situation in which the neurons in the sender area are desynchronised. However, due to the synchronised neurons in the area 1, we verify the existence of regions in the parameter space giiA×gieA with small values of the variability, as shown in Figure 9B. σΘ2′<0.5 corresponds to a certain stabilisation of the relative phase angle of the area 2. Figure 9C displays the mean relative phase angle of the area 2. For the lowest σΘ2′, we find a region where Θ2′¯≈π. In Figure 9D, we see that synchronous behaviour in the area 2 for large giiA and small gieA values arise, where there is a circle symbol. Partial anti-phase synchronisation is observed in the region indicated by the square. The triangle denotes the region in which a high inhibition of the excitatory neurons occurs, and as a consequence a great quantity of excitatory neurons in the receiver area are silenced.

Figure 10 exhibits a schematic representation of bidirectional interactions via excitatory connections with geiA and geeA conductances. Without an interaction between the areas (geeA=geiA=0), the neurons exhibit desynchronised activities.

Figure 11A displays the mean order parameter of the neuronal network. The region, where the circle is located, has a larger value of R¯, due to the fact that the neurons in the two areas are synchronised, namely phase synchronisation among neurons. In Figure 11B, we verify that the reduction of the variability can be associated with the synchronised activities between the areas. Figures 11C,D shows the mean order parameters of the areas 1 and 2, respectively. We can see that the regions of the small variability of Θ2′ (circle and triangle symbols) correspond to the synchronised activities. For the region with large variability values (square symbol), the neurons of the areas are desynchronised. In the region where the triangle symbol is located, there is an anti-phase synchronisation between the neurons of the areas 1 and 2.

Figure 12 displays a schematic representation of the bidirectional configuration interacting through inhibitory connections associated with the giiA and gieA conductances. Without interaction between the areas (giiA=gieA=0), the neurons exhibit desynchronised activities.

In the parameter space giiA×gieA, the mean order parameter for the neuronal network (Figure 13A) shows a region in which the neurons in the areas are synchronised, where a circle symbol is located. For the standard deviation of relative phase angle rotation of the area 2, we identify three regions with values about 0, as shown in Figure 13B. The inhibitory connections are responsible for decreasing the relative phase angle variability (circle, square, and triangle symbols). Figures 13C,D exhibits the mean order parameter of the areas 1 and 2, respectively. The regions of small variability of Θ2′ can correspond to the synchronous behaviour in each area or silence of some excitatory neurons.

In Figure 14, we show the raster plot (top), the relative phase angle (middle), and the order parameter (bottom) for (a) giiA=3.6 nS and gieA=0.2 nS, (b) giiA=1.5 nS and gieA=1.0 nS, and (c) giiA=0.2 nS and geiA=2.5 nS, according to the circle, square, and triangle symbols, respectively, denoted in Figure 13. Figure 14A displays the occurrence of phase synchronisation among the neurons between the areas. Depending on the values of giiA and gieA, it is possible to observe partial anti-phase synchronisation and silenced excitatory neurons, as shown in Figures 14B,C.