In this section, we summarize main experimental findings reported previously by our group (Zick et al., 2018). In that study, spike trains of ensembles of single neurons were recorded simultaneously from PFC of monkeys while they performed the dot-pattern expectancy (DPX) task, a task that measures specific deficits in cognitive control in schizophrenia (Jones et al., 2010). In the DPX task, the correct response (left or right joystick movement) to a probe stimulus depends on a preceding cue followed by a delay period (Materials and methods).

In the present study, we focus on PFC population spike dynamics recorded in the DPX task under two conditions: drug-naive and drug. The drug naive data were collected before monkeys were administered drug, phencyclidine, which is an NMDA receptor antagonist. Figure 1 shows the population average pairwise correlation between spike trains of neurons recorded in drug-naive (black) and drug (magenta) conditions. The strength of spike correlation was quantified by the ratio between the observed frequency of synchronous spikes (1ms resolution) and the frequency expected if the spike trains were uncorrelated (we subtracted 1 from this ratio so that correlation value is zero for uncorrelated, positive for correlated, and negative for anticorrelated spike activity, Materials and methods). The frequency of spike synchrony was determined from activity observed during a short (100 ms-long) window that was slid across time of task performance. Figure 1A shows that spike synchrony obtained from trials aligned to the cue onset (time 0) remained relatively weak and unchanged during the cue and delay periods, until the probe onset, in both drug-naive and drug conditions. The corresponding population average spike rates during these periods are shown in Figure 1C. Because the instant of response after probe presentation varied from trial to trial, to appreciate the time course of synchrony and spike rate after the delay period immediately preceding the response, in Figure 1B and D we aligned trials to response time (time 0). It is seen that synchrony started to increase sharply about 200ms before the motor response in the drug-naive condition and reached its peak at the time of the response (Figure 1B, black). The spike rate also started to increase before the response but more gradually and starting earlier before the response (Figure 1D, black). Both spike synchrony (Figure 1B) and spike rate (Figure 1D) exhibited secondary peaks occurring approximately 150–250ms after the response. In the drug condition, however, the increase in spike synchrony at the time of the response was markedly weakened (Figure 1B, magenta). The increase in spike rate was also reduced, although less dramatically (Figure 1D, magenta). We term this effect as NMDAR blockage induced desynchronization of spiking activity.

Figure 1

Download asset Open asset

To understand the phenomenon of drug-induced desynchronization of spiking activity and the role played by various components of synaptic currents, we considered a spiking network model representing a local circuit of monkey PFC. Details of the model and the theoretical framework are given in Materials and methods. Here, we only highlight their main aspects.

The network comprises excitatory and inhibitory neurons representing populations of pyramidal cells and interneurons, respectively. All neurons are modeled as leaky integrate-and-fire units (see, e.g., Dayan and Abbott, 2001). Synaptic connections are random and sparse, but the number of connections received by individual neurons is large. In addition to the recurrent local connections, each neuron also receives external connections from excitatory neurons outside of the network that fire spikes with rate ${\displaystyle {\nu }_{\mathrm{X}}}$.

Recurrent synaptic currents of excitatory connections are two-component, mediated by AMPA and NMDA receptors, whereas currents of inhibitory connections are mediated by GABA_A receptors (GABA thereafter). External currents represent the noisy inputs due to the background synaptic activity and are mediated by AMPA receptors. Thus, the model entails eight maximal synaptic conductance parameters ${\displaystyle g_{\mathrm{X},\alpha }}$, ${\displaystyle g_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A},\alpha }}$, ${\displaystyle g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A},\alpha }}$, ${\displaystyle g_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A},\alpha }}$ corresponding to the external AMPA, recurrent AMPA, NMDA, and GABA currents (${\displaystyle \alpha =\mathrm{E},\mathrm{I}}$ for excitatory and inhibitory neurons, respectively).

To produce a desired regime of network dynamics (asynchronous or synchronous) with a given firing rate of excitatory and inhibitory neurons ${\nu }_{\text{E}}$ and ${\nu }_{\text{I}}$, respectively, the values of the conductance parameters should be properly adjusted. For this purpose, we used mean field analysis. In this framework, population mean firing rates $v_{\mathrm{E}}^0$ and $v_{\mathrm{I}}^0$ in the asynchronous stationary state of the network can be effectively parametrized by three parameters expressed as ratios of component synaptic currents: ${\displaystyle I_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}/I_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}}$, ${\displaystyle I_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}/I_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}}$, and ${\displaystyle I_{\mathrm{X},\mathrm{E}}/I_{\theta ,\mathrm{E}}}$, where $I_R$ is the mean current of the $R$-receptor mediated synapse (${\displaystyle R=\mathrm{X},\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A},\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A},\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}$), and $I_{\theta ,\mathrm{E}}$ is the current that is needed for an excitatory neuron to reach firing threshold $\theta$ in absence of recurrent feedback. These parameters characterize the balance between recurrent excitation and inhibition, and the balance between external input and firing threshold. Once they are specified, for a given external spike rate ${\nu }_{\mathrm{X}}$ one can solve the mean field equations to obtain the underlying eight synaptic conductances providing the desired population mean firing rates $v_{\mathrm{E}}^0$ and $v_{\mathrm{I}}^0$ in asynchronous state of the network.

While the mean field analysis allows us to determine synaptic conductances that achieve desired firing rates of neurons, whether these rates remain stable over time is another issue. To address it, we conduct a linear stability analysis of the asynchronous state to understand if the network develops oscillatory instability caused by small fluctuations in population firing rates. This analysis entails two parameters, $\lambda$ and $\omega$, describing the rate of instability growth and the oscillation frequency. The asynchronous state is stable when ${\displaystyle \lambda <0}$; in this case small perturbations of firing rates cause exponentially damped oscillation of network activity. The case $\lambda =0$ corresponds to the onset of instability of the asynchronous state and the emergence of sustained sinusoidal oscillations of population average firing rates with frequency $\omega$; in the oscillatory regime spike trains remain sparse and irregular but at each oscillation cycle a random subset of network neurons fire synchronously giving rise to the synchronous irregular state. Lastly, when ${\displaystyle \lambda >0}$, small fluctuations in the stationary rates develop oscillatory instability with the amplitude of oscillations growing exponentially in time; however, higher order terms neglected in linear analysis can eventually saturate the instability growth (Brunel and Hakim, 1999), resulting in a stable oscillation with a finite amplitude.

To examine the boundary between the regions of asynchronous and synchronous states, we fix the balance of tonic NMDA current relative to GABA current, $I_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}/I_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}$, and vary the remaining two parameters: the balance between recurrent excitation and inhibition, $I_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}/I_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}$, and the balance between external excitation and firing threshold, $I_{\mathrm{X},\mathrm{E}}/I_{\theta ,\mathrm{E}}$. For a given point in this ${\displaystyle \left(I_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}/I_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}},I_{\mathrm{X},\mathrm{E}}/I_{\theta ,\mathrm{E}}\right)}$ parameter plane we solve the mean field equations to find the underlying set of eight synaptic conductances that provide the prescribed rates $v_{\mathrm{E}}^0$ and $v_{\mathrm{I}}^0$ given external spike rate ${\nu }_{\mathrm{X}}$, and then carry out linear stability analysis to find out if these rates are stable. Figure 2A shows a state diagram of the system for which external spike rate is set to ${\nu }_{\mathrm{X}}=5$ Hz, the rates of excitatory and inhibitory populations are set to $v_{\mathrm{E}}^0=5$ Hz, $v_{\mathrm{I}}^0=20$ Hz, and the NMDA current balance is fixed at $I_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}/I_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}=0.15$. The diagram shows solutions for $\lambda$ obtained from the linear stability analysis in the ${\displaystyle \left(I_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}/I_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}},I_{\mathrm{X},\mathrm{E}}/I_{\theta ,\mathrm{E}}\right)}$ parameter space. The asynchronous stationary state corresponds to the region where ${\displaystyle \lambda <0}$, whereas the synchronous oscillation state is realized in the region where ${\displaystyle \lambda >0}$. The asynchronous and synchronous states are separated by a “critical” or instability line on which $\lambda =0$ (shown in white color in Figure 2A). This boundary is the locus where the stationary network dynamics becomes unstable, and the sinusoidal oscillation of network activity develops. The oscillation frequency on the critical line, $f_{\mathrm{n}\mathrm{t}\mathrm{w}\mathrm{r}\mathrm{k}}=\omega /2\pi$, as a function of the balance between the recurrent AMPA and GABA currents, $I_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}/I_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}$, is shown in Figure 2B.

Figure 2 with 1 supplement see all

Download asset Open asset

The characteristic features of the state diagram qualitatively remain unchanged when the balance between the NMDA and GABA currents is varied (Figure 2—figure supplement 1A). Furthermore, the network frequency at the onset of oscillation, $f_{\mathrm{n}\mathrm{t}\mathrm{w}\mathrm{r}\mathrm{k}}$, essentially is independent of the $I_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}/I_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}$ balance (Figure 2—figure supplement 1B).

To study spike synchrony in asynchronous and synchronous networks in the context of the DPX task performed in drug-naive and drug conditions (Zick et al., 2018), we make two assumptions regarding neural and synaptic activity: (1) the increase in spike synchronization observed before the monkey’s response in Zick et al., 2018 is due to task-specific external afferent signals received by PFC neurons after probe presentation; (2) administration of NMDAR antagonist results in blocking NMDAR mediated synaptic currents. In the framework of our model, we implemented these assumptions as follows: task specific external signals were accounted for by an increase in the external spike rate from its background level ${\nu }_{\mathrm{X}}$, whereas the effect of drug administration was modeled by setting NMDAR conductances $g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A},\mathrm{E}}$ and $g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A},\mathrm{I}}$ to zero.

Next, to investigate how spike synchrony in asynchronous and synchronous networks depends on the modulations of $v_{\mathrm{X}}$ and $g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A},\alpha }$, for each network regime we proceed with the following three steps. First, we choose proper values for conductances, so that the underlying network operates in a desired regime providing the prescribed population firing rates $v_{\mathrm{E}}^{*}$ and $v_{\mathrm{I}}^{*}$ for a given external spike rate $v_{\mathrm{X}}^{*}$. We shall designate this network as the primary network relating to the underlying regime and distinguish the corresponding values of all its parameters by the asterisk (*). Second, we carry out a series of network simulations, in which external spike rate $v_{\mathrm{X}}$ and NMDAR conductance $g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A},\alpha }$ are varied relative to their standard values $v_{\mathrm{X}}^{*}$ and $g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A},\alpha }^{*}$, respectively. Lastly, for each simulated network, we compute population average pairwise correlation between spike trains of neurons and analyze how this correlation depends on the external spike rate and NMDAR conductance.

To perform a comparison between the primary networks, we need to choose appropriate values for their parameters. We begin with the parameters that are common to both networks. First, we set the excitatory and inhibitory population mean firing rates to $v_{\mathrm{E}}^{*}=5$ Hz and $v_{\mathrm{I}}^{*}=20$ Hz, respectively, which are on the order of magnitude of spontaneous rates observed for PFC neurons. Second, since external inputs represent activity of excitatory neurons outside the PFC circuit model, we choose the background external rate $v_{\mathrm{X}}^{*}$ to be the same as the excitatory population rate $v_{\mathrm{E}}^{*}$ inside the model and, thus, set $v_{\mathrm{X}}^{*}=5$ Hz. Lastly, for both networks, we fix the balance between NMDA and GABA currents at $I_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}^{*}/I_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}^{*}=0.15$. Note that the state diagram in the ${\displaystyle \left(I_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}/I_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}},I_{\mathrm{X},\mathrm{E}}/I_{\theta ,\mathrm{E}}\right)}$ space shown in Figure 2A was obtained exactly for these values of the above listed parameters. We use this state diagram for selecting the primary networks and determining the remaining parameters that are network specific.

In this regard, we note that each point in the ${\displaystyle \left(I_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}/I_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}},I_{\mathrm{X},\mathrm{E}}/I_{\theta ,\mathrm{E}}\right)}$ plane corresponds to a network with a specific set of synaptic conductances. For synchronous regime, we look for a network on the critical line ($\lambda =0$, white line in Figure 2A), at the onset of oscillatory instability with a frequency in the $\gamma$-band (a frequency band associated with the LFPs recorded from prefrontal areas [Bastos et al., 2018; Lundqvist et al., 2016]). For instance, the point marked by a blue asterisk in Figure 2A located at ${\displaystyle \left(I_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}^{\ast }/I_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}^{\ast }=0.4,I_{\mathrm{X},\mathrm{E}}^{\ast }/I_{\theta ,\mathrm{E}}^{\ast }=1.09\right)}$ corresponds to such a network with oscillation frequency $f_{\mathrm{n}\mathrm{t}\mathrm{w}\mathrm{r}\mathrm{k}}^{*}~58$ Hz (Figure 2B). In the following, we refer to this network as the critical state primary network.

Correspondingly, for the asynchronous regime, we need to select a network that is far from the critical line and deep in the region of stable network dynamics (${\displaystyle \lambda <0}$). The point marked by a red asterisk in Figure 2A located at ${\displaystyle \left(I_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}^{\ast }/I_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}^{\ast }=0.2,I_{\mathrm{X},\mathrm{E}}^{\ast }/I_{\theta ,\mathrm{E}}^{\ast }=1.09\right)}$ is an example of such a network. We shall refer to this network as the steady state primary network. For each primary network, we obtain the underlying set of eight synaptic conductance parameters $g_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A},\alpha }^{*},g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A},\alpha }^{*},g_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A},\alpha }^{*},g_{\mathrm{x}\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A},\alpha }^{*}$ (${\displaystyle \alpha =\mathrm{E},\mathrm{I}}$) by numerically solving the mean field equations.

To investigate characteristic features of spiking dynamics in asynchronous and synchronous regimes, we carried out direct simulations of the primary networks. Both networks comprise $N=\mathrm{5,000}$ neurons, of which $N_{\text{E}}=\mathrm{4,000}$ are excitatory and $N_{\text{I}}=\mathrm{1,000}$ inhibitory. Neurons are connected randomly with a probability $p=0.2$. Figure 3 illustrate the behavior of simulated networks with synaptic conductance parameters corresponding to the steady and critical primary networks indicated by the red and blue asterisks, respectively, in the state diagram presented in Figure 2A. The dynamic behavior is shown at the level of individual cell activity (spike rasters, top of panels in Figure 3), as well as whole population activity (bottom of panels in Figure 3).

Figure 3 with 1 supplement see all

Download asset Open asset

In simulations shown in Figure 3 panels A1 and A2 external spike rate ${\nu }_{\mathrm{X}}$ was fixed at the level of $v_{\mathrm{X}}^{*}=5$ Hz chosen for the primary networks. It is seen that excitatory and inhibitory neurons exhibit highly irregular firing with average rates, ${\nu }_{\text{E}}$ and ${\nu }_{\text{I}}$, about 5.2 Hz and 20 Hz in the steady state primary network (Figure 3A1) and 5.5 Hz and 21 Hz in the critical state primary network (Figure 3A2). These observed in simulations rates ${\nu }_{\text{E}}$ and ${\nu }_{\text{I}}$ are in good agreement with the prescribed rates $v_{\mathrm{E}}^{*}=5$ Hz and $v_{\mathrm{I}}^{*}=20$ Hz that were used to derive the synaptic conductance parameters of the simulated networks. Moreover, Figure 3A1 demonstrates that population activity of the steady state primary network is rather stationary in time, whereas activity of the critical primary network shown in Figure 3A2 exhibits signs of developing of oscillatory instability (compare Figure 3—figure supplement 1A1 vs Figure 3—figure supplement 1A2). Thus, spiking dynamics observed in the simulated steady state primary network displays basic characteristics of the asynchronous regime—irregular firing of individual neurons and stationary population activity. Correspondingly, the behavior of the simulated critical state primary network exhibits similarity with the boundary regime on which the asynchronous stationary state destabilizes and oscillatory behavior of the population activity emerges.

Panels B1 and B2 in Figure 3 demonstrate results of simulations in which external spike rate ${\nu }_{\mathrm{X}}$ was increased by 5% relative to the rate $v_{\mathrm{X}}^{*}$ used in simulations illustrated in Figure 3 panels A1 and A2. For the steady state primary network (Figure 3B1), the firing rates of excitatory and inhibitory neurons increase with the external drive. However, the regime of network dynamics qualitatively does not change and remains asynchronous (compare Figure 3—figure supplement 1A1 vs Figure 3—figure supplement 1B1). In contrast, stronger external inputs received by the critical state primary network synchronize population activity (Figure 3B2). It is seen that while individual neurons continue to fire irregularly, population activity now clearly exhibits oscillatory behavior, indicating that the network is in synchronous irregular regime in which the average firing frequency of neurons is low, about 20 Hz, compared to the frequency of network oscillation, which is about 50 Hz (see Figure 3—figure supplement 1B2). This frequency is close to the theoretically predicted network frequency of 58 Hz near the onset of oscillation.

Thus, direct simulations confirm that analytically derived network parameters for both steady and critical primary networks provide the anticipated regimes of network dynamics.

To facilitate the comparison of characteristic features exhibited by a simulated network with experimentally measurable quantities, we compute temporal correlation of spiking activity that quantifies average pairwise correlation between spike trains of excitatory neurons. In the context of the DPX task performed in drug-naive and drug conditions studied in Zick et al., 2018 and with the purpose of elucidating the mechanism of drug-induced desynchronization of spiking activity, we investigated how temporal correlations depend on the strength of external drive and the NMDAR mediated synaptic current. To this end, we varied external input rate $v_{\mathrm{X}}$ and the NMDAR conductance parameters $g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A},\mathrm{E}}$ and $g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A},\mathrm{I}}$ relative to their respective standard values $v_{\mathrm{X}}^{*}$, and $g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A},\mathrm{E}}^{*}$ and $g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A},\mathrm{I}}^{*}$ , while keeping all other system parameters fixed, and performed simulations of the ensuing networks. Conductances for excitatory and inhibitory neurons were scaled with the same factor and, therefore, their relative values $g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A},\mathrm{E}}/g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A},\mathrm{E}}^{*}$ and $g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A},\mathrm{I}}/g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A},\mathrm{I}}^{*}$ are the same; in the following we drop the $\text{E},\text{I}$ designation.

Figure 4 displays correlation of spiking activity (panels A1, A2, C1, C2) and synchrony (0-lag correlation, panels B1, B2, D1, D2) obtained from spike trains of simulated steady (panels A1, B1, C1, D1) and critical (panels A2, B2, C2, D2) networks for a range of $v_{\mathrm{X}}/v_{\mathrm{X}}^{*}$ (panels A1, A2, B1, B2) and $g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}/g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}^{*}$ (panels C1, C2, D1, D2) values. It is seen that in the steady state primary network correlations are weak and insensitive to the modulations of external input rate or NMDAR conductance (Figure 4 panels A1, B1, C1, D1). In contrast, in the critical state primary network temporal correlations show sharp dependence on these parameters (Figure 4 panels A2, C2), and with decreasing external drive or decreasing NMDAR conductance profoundly attenuating spike synchrony (Figure 4 panels B2, D2).

Figure 4

Download asset Open asset

What are the network mechanisms of external drive and NMDA conductance dependent spike synchronization? Why in the network close to the boundary between asynchronous and synchronous regimes, are spike correlations strongly affected by the modulations of external inputs and recurrent NMDA currents, but in the network far from this boundary and deep in the region of the asynchronous regime, correlations are essentially independent of these modulations? How does the interplay between synchronous and asynchronous regimes at their boundary lead to spike synchronization when external input rate ${\nu }_{\mathrm{X}}$ increases, and to desynchronization when the NMDA conductance $g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}$ decreases?

To answer these questions and to illuminate the role of asynchronous and synchronous regimes in the shaping of network-wide synchronization of spiking activity, we carried out linear stability analysis in the ${\displaystyle \left(v_{\mathrm{X}}/v_{\mathrm{X}}^{\ast },g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}/g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}^{\ast }\right)}$ parameter plane while keeping the remaining parameters fixed. For both steady and critical state primary networks, stability is investigated in the vicinity of the standard values of the external input spike rate and NMDAR conductances corresponding to the respective networks.

Figure 5 illustrates state diagrams in the ${\displaystyle \left(v_{\mathrm{X}}/v_{\mathrm{X}}^{\ast },g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}/g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}^{\ast }\right)}$ plane in the neighborhood of the steady (Figure 5A) and critical (Figure 5B) state primary networks. As in Figure 2A, the critical line ($\lambda =0$) separating the asynchronous stationary (${\displaystyle \lambda <0}$) and synchronous oscillatory (${\displaystyle \lambda >0}$) states is shown in white color. Asterisks correspond to the loci of the steady (Figure 5A) and critical (Figure 5B) state primary networks in these parameter planes. It is seen that the modulations of ${\nu }_{\mathrm{X}}$ and $g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}$ in the steady state primary network (Figure 5A) do not change the network state; these modulations have no impact on the spike correlation and the strength of synchrony (Figure 4B1 and D1 and Figure 5A insets).

Figure 5

Download asset Open asset

In contrast, the modulations of ${\nu }_{\mathrm{X}}$ and $g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}$ in the critical state primary network (Figure 5B) induce transitions between the network states. Specifically, as the external input spike rate ${\nu }_{\mathrm{X}}$ increases (horizontal yellow arrow in Figure 5B) the system crosses the boundary between asynchronous and synchronous regimes and the network state changes from stationary to oscillatory; this transition is accompanied by a sharp increase in spike synchrony (Figure 4B2 and Figure 5B bottom inset). The decrease of NMDAR conductance $g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}$ (vertical magenta arrow in Figure 5B) causes the system to cross the boundary again, and the network state changes now from oscillatory to stationary; this transition is accompanied by a sharp decrease in spike synchrony (Figure 4D2 and Figure 5B right inset).

Thus, this analysis reveals that networks that are close to the boundary between asynchronous and synchronous regimes, in contrast to asynchronous networks that are far from this boundary, have a rich dynamic behavior. The dynamic states of these networks could be easily switched around by modulations in the external drive and the strength of recurrent excitation by NMDAR mediated currents. Switching between the network states, in turn, results in sharp changes in the degree of network-wide synchronization of spiking activity in response to these modulations.

As illustrated in Figure 1B, spiking activity observed in monkey PFC in the DPX task (Zick et al., 2018) remains practically desynchronized after probe presentation for about 200ms but it begins to increase sharply about 200ms before the motor response. To get a deeper insight into the properties of spike timing dynamics, we show in Figure 6 temporal correlations of spiking activity during the 200ms period following probe presentation (Figure 6A1) and during the 200ms period preceding the motor response (Figure 6B1) in drug-naive (black) and drug (magenta) conditions. It can be now appreciated that in drug-naive condition, population activity during the pre-response period develops characteristics of synchronized oscillation behavior, as signaled by the appearance of time lagged peaks of correlation (blue arrows, Figure 6B1, black). However, administration of a drug blocking NMDAR desynchronizes neuronal activity during this period (Figure 6B1, magenta).

Figure 6 with 2 supplements see all

Download asset Open asset

The presence of strong spike synchrony (0ms lag) together with the correlation peaks at ±18 ms lags in the pre-response period (Figure 6B1), and the absence of these characteristics in the initial probe period (Figure 6A1) suggest that after probe presentation but before motor response network dynamics switches from the asynchronous stationary state to the synchronous oscillation state with a $\gamma$-frequency around 55 Hz. Desynchronization of neuronal activity produced by drug administration implies that NMDAR blockage prevents PFC circuits operating in the asynchronous regime from switching to synchronous dynamics.

These experimental findings could be readily explained by a prefrontal network model that operates on the boundary between asynchronous and synchronous regimes. We start by recalling that in the framework of our approach the pre-response afferent signals, which we assume are received by PFC neurons before the monkey’s response, are modeled as an increase in the external spike rate from its background level ${\nu }_{\mathrm{X}}$. This assumption is supported by the increase in the population spike rate preceding the monkey’s response observed in neural data shown in Figure 1D. Secondly, the effect of drug administration is modeled by setting NMDAR conductances $g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A},\mathrm{E}}$ and $g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A},\mathrm{I}}$ to zero. The capacity of the prefrontal network model to provide a circuit mechanism for the emergence of synchrony in spiking activity and drug-dependent desynchronization can be illustrated by considering the system’s behavior in the ${\displaystyle \left(v_{\mathrm{X}}/v_{\mathrm{X}}^{\ast },g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}/g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}^{\ast }\right)}$ state plane around the point ${\displaystyle \left(v_{\mathrm{X}}/v_{\mathrm{X}}^{\ast }=1,g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}/g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}^{\ast }=1\right)}$ corresponding to the critical state primary network (Figure 6C). In this space, the effects of probe presentation on the spiking dynamics of the prefrontal circuit model under drug-naive (${\displaystyle g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}/g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}^{\ast }=1.25}$) and drug ($g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}/g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}^{*}=0$) conditions are represented, respectively, by black and magenta horizontal arrows (Figure 6C). The arrows are pointing from the state of the network corresponding to the initial probe period (${\displaystyle v_{\mathrm{X}}/v_{\mathrm{X}}^{\ast }=0.97}$) to the network state corresponding to the pre-response period (${\displaystyle v_{\mathrm{X}}/v_{\mathrm{X}}^{\ast }=1.03}$).

In drug-naive condition, increase in the external spike rate $v_{\mathrm{X}}$ switches the circuit model from asynchronous to synchronous regime (Figure 6C, black arrow crosses the boundary between the regimes). The oscillation frequency is about 50 Hz, which is manifested in the temporal correlations of spiking activity as a sharp increase in synchrony and appearance of peaks at ±20 ms lags (Figure 6A2 vs B2, black line). This is very similar to what is observed in monkey PFC during the initial probe and pre-response periods in the DPX task (Fig. 6A1 vs A2 and B1 vs B2, black line). In the drug condition, setting NMDAR conductance to zero prevents the circuit model from switching to the synchronous regime in response to an increase in the external spike rate $v_{\mathrm{X}}$ (Figure 6C, magenta arrow does not cross the boundary between the regimes). This, in turn, considerably reduces the degree of spike synchrony compared to drug-naive condition (Figure 6B2, magenta vs black line), similar to the desynchronizing effect of NMDAR antagonist administration on spiking activity in monkey PFC (Figure 6 B1, magenta vs black line).

In the consideration above, we investigated the network spiking dynamics in the asynchronous and synchronous states during stationary external input at a decreased ($v_{\mathrm{X}}/v_{\mathrm{X}}^{*}=0.97$) and increased ($v_{\mathrm{X}}/v_{\mathrm{X}}^{*}=1.03$) external rate. To simulate a more biologically realistic scenario, we also examined the network behavior in response to transient external input. In this analysis, external input rate had a trapezoid-like temporal profile (Figure 7A). First, external rate was fixed at a lower level ($v_{\mathrm{X}}/v_{\mathrm{X}}^{*}=0.97$) setting the network in the asynchronous state. Then, throughout 100ms period the rate was linearly increased to a higher level ($v_{\mathrm{X}}/v_{\mathrm{X}}^{*}=1.05$) and kept constant for 400ms, pushing the network across the boundary to the synchronous state. Finally, the rate was decreased to the initial level during the next 100ms to switch the network back to the asynchronous state. Figure 7 shows time evolution of population spike rate (Figure 7B) and synchrony (Figure 7C) in response to such transient external input (Figure 7A) for $g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}/g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}^{*}=1.25$ (black) and $g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}/g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}^{*}=0$ (magenta) corresponding to drug-naive and drug conditions, respectively. These simulated temporal profiles can be compared with the temporal profiles shown in Figure 1 for population average spike rate (Figure 1D) and synchrony (Figure 1B) obtained from experimental data. Parallels between the simulated and recorded neural data are evident. Spike synchrony and spike rate peak at about the same time both in simulated (Figure 7B and C) and recorded (Figure 1B and D) neural activity. Further, the increase in spike rate is early and gradual in comparison to the increase in spike synchrony which is delayed and abrupt both in simulated (Figure 7B and C) and recorded (Figure 1B and D) neural activity. While our relatively simple model qualitatively is consistent with dynamical features of the firing rate and synchrony observed in primate PFC, there are, however, some quantitative discrepancies in firing rates. In addition, recorded neural activity exhibits complex dynamics following the response (Figure 1B and D), that are not evident in the simulation (Figure 7B and D). This presumably reflects temporal modulation of synaptic inputs to the recorded neurons in the biological data that are more complex than the ramp transient we implemented in the simulation.

Figure 7

Download asset Open asset

In summary, the analyses of simulations with stationary and transient external inputs suggest that when the prefrontal network model operates close to the boundary between asynchronous stationary and synchronous oscillatory regimes it has a considerable capacity to capture experimentally observed aspects of spike synchrony in both drug-naive and drug conditions.

So far, in most of our analyses, we did not vary the balance between the tonic component of recurrent excitation mediated by NMDA and GABA inhibition, keeping it fixed at $I_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}^{*}/I_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}^{*}=0.15$. We have only shown that the network frequency at the onset of oscillation essentially is independent of the $I_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}/I_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}$ balance (Figure 2—figure supplement 1B), and that the characteristic features of the ${\displaystyle \left(I_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}/I_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}},I_{\mathrm{X},\mathrm{E}}/I_{\theta ,\mathrm{E}}\right)}$ state diagram qualitatively remain unchanged when this balance is varied (Figure 2—figure supplement 1A). Could, however, the $I_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}^{*}/I_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}^{*}$ balance be crucial for the prefrontal circuit model capacity to provide the underlying mechanism for external input and NMDA conductance dependent spike synchronization? To investigate this issue, we analyzed how characteristic features of the ${\displaystyle \left(v_{\mathrm{X}}/v_{\mathrm{X}}^{\ast },g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}/g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}^{\ast }\right)}$ state diagram shown in Figure 6C depend on the $I_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}^{*}/I_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}^{*}$ balance.

Figure 8 shows state diagrams in the ${\displaystyle \left(v_{\mathrm{X}}/v_{\mathrm{X}}^{\ast },g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}/g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}^{\ast }\right)}$ plane obtained for several $I_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}^{*}/I_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}^{*}$ balance values. It is seen that the orientation of the critical line in the state space depends on the $I_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}^{*}/I_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}^{*}$ balance. When the balance is shifted toward stronger inhibition (${\displaystyle I_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}^{\ast }/I_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}^{\ast }<0.15}$, Figure 8A), the critical line becomes too steep: in the drug condition, blocking NMDA current may not necessarily lead to spike desynchronization because the external spike modulation could trigger the network to switch to the synchronous regime (magenta arrow in Figure 8A). On the other hand, when the balance is shifted toward stronger tonic excitation (${\displaystyle I_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}^{\ast }/I_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}^{\ast }>0.15}$, Figure 8C), the critical line becomes too flat: in the drug-naive condition the external spike modulation may not be able to produce strong enough synchrony because the system would be too close to the critical line, and not shift deep enough into the region of the synchronous regime (black arrow in Figure 8C).

Figure 8

Download asset Open asset

Further insights into how synaptic conductances and external rate affect synchrony can be achieved by obtaining an analytic expression describing the dependence of the rate of oscillatory instability growth $\lambda$ on these parameters near the boundary between the asynchronous and synchronous states. Such expression can be derived by linearizing the stability analysis equations in the limit of small relative changes $∆g_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}/g_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}^{*}$, $∆g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}/g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}^{*}$, $∆g_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}/g_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}^{*}$, and $∆v_{\mathrm{X}}/v_{\mathrm{X}}^{*}$ of the synaptic parameters around the critical point ${\displaystyle \left\{g_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A},\left\{\mathrm{E},\mathrm{I}\right\}}^{\ast },g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A},\left\{\mathrm{E},\mathrm{I}\right\}}^{\ast },g_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A},\left\{\mathrm{E},\mathrm{I}\right\}}^{\ast },v_{\mathrm{X}}^{\ast }\right\}}$ corresponding to the onset of oscillatory instability where $\lambda =0$ (conductances $g_{R,\mathrm{E}}$ and $g_{R,\mathrm{I}}$ of excitatory and inhibitory neurons (${\displaystyle R=\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A},\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A},\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}$) are again scaled with the same factors and, thus, their relative changes are equal: $∆g_{R,\mathrm{E}}/g_{R,\mathrm{E}}^{*}=∆g_{R,\mathrm{I}}/g_{R,\mathrm{I}}^{*}$). The calculation is detailed in the Materials and methods section. The result is that $\lambda$ in the vicinity of the critical point on the boundary between the steady and oscillatory states can be approximated by

where ${\mathrm{\Lambda }}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}$, ${\mathrm{\Lambda }}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}$, and ${\mathrm{\Lambda }}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}$ are quantities defined by the parameters of the critical state network around which the equations are linearized, ${\varphi }_{I_{\mathrm{s}\mathrm{y}\mathrm{n},\alpha }}^{`}$ is the slope of the neuron’s current-frequency response function at the critical state, and $∆{\varphi }_{I_{\mathrm{s}\mathrm{y}\mathrm{n},\alpha }}^{`}$ is the change in the slope of the response function due to the deviations of the synaptic parameters from their critical values (${\displaystyle \alpha =\mathrm{E},\mathrm{I}}$ for excitatory and inhibitory neurons, respectively). The deviations of the synaptic conductances $∆g_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}$, $∆g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}$, $∆g_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}$, and external rate $∆v_{\mathrm{X}}$ give rise to the changes in the corresponding average recurrent $I_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}$, $I_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}$, $I_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}$ and external $I_{\mathrm{X}}$ synaptic currents. This produces the change $∆I_{\mathrm{s}\mathrm{y}\mathrm{n}}$ in the average total current $I_{\mathrm{s}\mathrm{y}\mathrm{n}}=I_{\mathrm{X}}+I_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}+I_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}-I_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}$ and shifts the operating point of the current-frequency response function $v=\varphi \left(I_{\mathrm{s}\mathrm{y}\mathrm{n}}\right)$ that describes the relationship between the average total input current $I_{\mathrm{s}\mathrm{y}\mathrm{n}}$ and the output firing frequency of the neuron $v$. For the leaky integrate-and-fire neuron model, $\varphi$ is a monotonically increasing non-linear function (see, e.g., Renart et al., 2003). Thus, the shift of the operating point of the neuron’s response function $\varphi$ due to the change $∆I_{\mathrm{s}\mathrm{y}\mathrm{n}}$ in the total average synaptic current results not only in the change of the firing rate (i.e. $∆\varphi$), but also in the change of the slope of the response function $∆{\varphi }_{I_{\mathrm{s}\mathrm{y}\mathrm{n}}}^{`}$. The latter can be calculated by linearizing the self-consistent mean field equations (see Materials and methods). As a result, $∆{\varphi }_{I_{\mathrm{s}\mathrm{y}\mathrm{n}}}^{`}$ is approximated as

where $U_{\alpha }$ is a positive constant defined by the parameters of the critical state network around which the mean field equations are linearized. The analytical expression for $\lambda$ given by Equation 1, 2 provides a very good approximation of the exact relationship (see Appendix 1).

Within the linear approximation, the change $∆{\varphi }_{I_{\mathrm{s}\mathrm{y}\mathrm{n}}}^{`}$ is proportional to the change $∆I_{\mathrm{s}\mathrm{y}\mathrm{n}}$ :

where $I_0$ is a positive constant. Hence, from Equation 2 it follows that $∆I_{\mathrm{s}\mathrm{y}\mathrm{n}}$ is proportional to the expression in the brackets:

From Equation 1 it follows that the rate of oscillatory instability growth $\lambda$ directly depends on the changes in the synaptic conductances but does not explicitly depend on the external rate variation $∆v_{\mathrm{X}}$. However, $\lambda$ depends on $∆v_{\mathrm{X}}$ indirectly via the terms involving the change in the slope $∆{\varphi }_{I_{\mathrm{s}\mathrm{y}\mathrm{n}}}^{`}$ due to the change in the average total synaptic current $∆I_{\mathrm{s}\mathrm{y}\mathrm{n}}$ (Equations 3, 4). In fact, $∆I_{\mathrm{s}\mathrm{y}\mathrm{n}}$ is affected by the variations of the synaptic conductances as well. Thus, the rate of instability growth $\lambda$ not only directly depends on the synaptic conductances, but also indirectly via the effect of the recurrent excitatory and inhibitory currents mediated by them on the average total synaptic current and, therefore, the operating point of the current-frequency response function.

The factors ${\mathrm{\Lambda }}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}$, ${\mathrm{\Lambda }}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}$, and ${\mathrm{\Lambda }}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}$ govern the strength of the direct and indirect contributions of the changes in the synaptic conductances $∆g_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}$, $∆g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}$, and $∆g_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}$ to the oscillatory instability. By inspecting Equation 1, one can see that the strength of the direct contribution of the change $∆g_R$ (${\displaystyle R=\mathrm{A}\mathrm{M}\mathrm{A}\mathrm{P}\mathrm{A},\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A},\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}$) is determined only by the corresponding factor ${\mathrm{\Lambda }}_R$ via the term ${\mathrm{\Lambda }}_R∆g_R/g_R^{*}$. However, the strength of its indirect contribution is determined by all three factors, ${\mathrm{\Lambda }}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}$, ${\mathrm{\Lambda }}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}$, and ${\mathrm{\Lambda }}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}$, through the changes in the slopes $∆{\varphi }_{I_{\mathrm{s}\mathrm{y}\mathrm{n},\mathrm{E}}}^{`}$ and $∆{\varphi }_{I_{\mathrm{s}\mathrm{y}\mathrm{n},\mathrm{I}}}^{`}$, which depend on $∆g_R$ (Equation 2). For example, the strength of direct contribution to $\lambda$ due to the change in the GABAR conductance $∆g_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}$ is determined only by ${\mathrm{\Lambda }}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}$ via the term ${\mathrm{\Lambda }}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}∆g_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}/g_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}^{*}$ in Equation 1. However, the strength of indirect contribution from $∆g_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}$ is determined by all three factors ${\mathrm{\Lambda }}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}$, ${\mathrm{\Lambda }}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}$, and ${\mathrm{\Lambda }}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}$ via the terms ${\mathrm{\Lambda }}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}∆{\varphi }_{I_{\mathrm{s}\mathrm{y}\mathrm{n},\mathrm{E}}}^{`}/{\varphi }_{I_{\mathrm{s}\mathrm{y}\mathrm{n},\mathrm{E}}}^{`}$, ${\mathrm{\Lambda }}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}∆{\varphi }_{I_{\mathrm{s}\mathrm{y}\mathrm{n},\mathrm{E}}}^{`}/{\varphi }_{I_{\mathrm{s}\mathrm{y}\mathrm{n},\mathrm{E}}}^{`}$, and ${\mathrm{\Lambda }}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}∆{\varphi }_{I_{\mathrm{s}\mathrm{y}\mathrm{n},\mathrm{I}}}^{`}/{\varphi }_{I_{\mathrm{s}\mathrm{y}\mathrm{n},\mathrm{I}}}^{`}$ in Equation 1 because $∆{\varphi }_{I_{\mathrm{s}\mathrm{y}\mathrm{n},\mathrm{E}}}^{`}$ and $∆{\varphi }_{I_{\mathrm{s}\mathrm{y}\mathrm{n},\mathrm{I}}}^{`}$ themselves depend on $∆g_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}$ (Equation 2). As noted above, this indirect contribution is due to the change in the average total synaptic current and, therefore, the change in the operating point of the current-frequency response function.

Figure 9 illustrates the contributions of individual terms involving ${\mathrm{\Lambda }}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}$, ${\mathrm{\Lambda }}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}$, and ${\mathrm{\Lambda }}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}$ in Equation 1 to the oscillatory instability growth rate $\lambda$. The panels display separately four cases in which one of the synaptic parameters is varied while the remaining three are kept constant at their critical values. It is seen that in all four cases the dominant contribution to $\lambda$ is coming from the term involving ${\mathrm{\Lambda }}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}$. The contribution related to ${\mathrm{\Lambda }}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}$ is nearly zero, whereas the contribution from ${\mathrm{\Lambda }}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}$ term is much smaller than the one from ${\mathrm{\Lambda }}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}$. While both ${\mathrm{\Lambda }}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}$ and ${\mathrm{\Lambda }}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}\ll {\mathrm{\Lambda }}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}$, the primary reasons are different (see Appendix 2).

Figure 9

Download asset Open asset

It should be noted that even though ${\mathrm{\Lambda }}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}$ and ${\mathrm{\Lambda }}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}$ are negligibly small, this does not mean that changes in the NMDAR and GABAR conductances do not affect oscillatory instability (black lines, panels C and D, Figure 9). The fact that ${\mathrm{\Lambda }}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}$ and ${\mathrm{\Lambda }}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}$ are small only means that $∆g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}$ and $∆g_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}$ do not affect the oscillatory instability directly. However, the changes in the NMDAR and GABAR conductances still affect the instability growth rate $\lambda$ indirectly via the term involving the product of ${\mathrm{\Lambda }}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}$ and $∆{\varphi }_{I_{\mathrm{s}\mathrm{y}\mathrm{n},\mathrm{E}}}^{`}/{\varphi }_{I_{\mathrm{s}\mathrm{y}\mathrm{n},\mathrm{E}}}^{`}$ in Equation 1, as mentioned above (and summarized below).

Since ${\mathrm{\Lambda }}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}},{\mathrm{\Lambda }}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}\ll {\mathrm{\Lambda }}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}$, we can neglect the terms involving ${\mathrm{\Lambda }}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}$ and ${\mathrm{\Lambda }}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}$ in Equation 1 for the oscillatory instability growth rate $\lambda$. With this approximation, the equation for $\lambda$ simplifies to

Inserting the expression for $∆I_{\mathrm{s}\mathrm{y}\mathrm{n}}$ from Equation 4, we obtain

Thus, the instability growth rate $\lambda$, in essence, directly depends only on the AMPAR conductance via the first term in Equation 6. The term in the brackets describes the dependence on the NMDAR mediated excitation, GABAR mediated inhibition, and external rate $v_{\mathrm{X}}$ that affect $\lambda$ only indirectly through their effect on the operating point of the response function. In addition, $\lambda$ also depends indirectly on the AMPAR conductance. For the critical state network $U_{\mathrm{E}}=2.5$ and $I_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}/I_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}=0.4$. Therefore, half of the contribution to $\lambda$ is due to the indirect and the second half due to the direct dependence on the AMPAR conductance. A more detailed consideration of the direct and indirect pathways by which modulations of synaptic conductances and external rate affect synchrony is given in Appendix 3.

Since in our network model we vary only the NMDAR conductance and external rate, Equation 6 for $\lambda$ simplifies to

The expression in the brackets is proportional to the change in the average total synaptic current $I_{\mathrm{s}\mathrm{y}\mathrm{n}}$ (Equation 4). The transition to synchrony in the model simulations is achieved by increasing external input (drug-naive condition in Zick et al., 2018), whereas reducing the NMDAR conductance prevents the network from such transition (drug condition in Zick et al., 2018). These simulation results and the mechanism implemented in our model for the transition between the steady and oscillatory states, and the lack thereof when the NMDAR conductance is blocked can be explained in terms of Equation 7 for the instability growth rate $\lambda$. As explained above, changes in external rate $∆v_{\mathrm{X}}$ and NMDAR conductance $∆g_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}$ both affect synchrony via indirect mechanism by changing the excitatory drive ${\displaystyle I_{\mathrm{s}\mathrm{y}\mathrm{n}}}$ and, therefore, shifting the operating point of the neuron’s response function. In the drug-naive condition, increase in external rate (${\displaystyle \mathrm{\Delta }{v}_{\mathrm{X}}>0}$) increases the excitatory drive. As a result, $\lambda$ becomes positive (see Equation 7) and the network switches to the synchronous regime. However, in the drug condition, when NMDAR is blocked (${\displaystyle \mathrm{\Delta }{g}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}<0}$), the initial excitatory drive is reduced compared to the drug-naive condition, and now the same increase in external rate $∆v_{\mathrm{X}}$ becomes insufficient to offset the reduced excitatory drive caused by the NMDAR blockage. As a result, $\lambda$ stays negative and the network remains in asynchronous regime.

A more formal consideration of the mechanism implemented in our model for the transition between the steady and oscillatory states as well as an analytical approximation for the critical line separating these two states are given in Appendix 4. In Appendix 5, we provide theoretical explanations in terms of the equation for the oscillatory instability growth rate $\lambda$ for some other simulation results obtained earlier.

To gain further insights into how specifically synaptic conductances and external rate affect emergence of synchronous oscillations, we obtained an analytic approximation for the oscillatory instability growth rate $\lambda$ describing the dependence on these parameters near the boundary between the asynchronous and synchronous states where ${\displaystyle \lambda =0}$. We showed that $\lambda$, in essence, directly depends only on the AMPAR synaptic conductance; it is virtually independent of the NMDAR conductance due to the slow synaptic decay time constant, while the dependence on the GABAR conductance is much weaker compared to AMPAR because of nearly 90° effective phase lag introduced by synaptic filtering. However, $\lambda$ depends on the NMDAR, GABAR as well as AMPAR conductances and external rate indirectly via their effect on the operating point of the neuron’s input current-output frequency response function. The direct dependence manifests the essential influence of the AMPAR synaptic conductance on the strength of an excitatory-inhibitory feedback loop via fast excitatory to excitatory and excitatory to inhibitory recurrent connections. The indirect dependence manifests the influence of the synaptic conductances and external rate on the location of the operating point on the current-frequency response curve and, therefore, the slope of the response function. The steepness of the slope, in turn, determines the amplitude of the neuron’s response to dynamically varying input current and, therefore, affects the strength of excitatory feedback.

The analytic expression for the oscillation growth rate $\lambda$ also reveals the differences and similarities in how AMPAR and NMDAR, both of which mediate recurrent excitation, influence the stability of asynchronous state and transition to synchronous oscillations. Both AMPAR and NMDAR conductances affect $\lambda$ indirectly by influencing the amplitude of the neuron’s response to varying input current. However, because AMPA currents are much faster than NMDA currents, unlike NMDAR, AMPAR conductance also affects $\lambda$ directly by influencing the strength of fast excitatory feedback.

Previous work (Wang, 1999) suggested that NMDAR mediated recurrent currents have a stabilizing effect on the network activity. Compte and colleagues (Compte et al., 2000) carried out spiking network simulations with different relative contributions of the NMDAR and AMPAR mediated currents to the recurrent excitation and showed that with less NMDA but more AMPA currents, the asynchronous steady state becomes unstable and neurons begin to synchronize, leading to network oscillations in the gamma band. At first glance, these simulation results seem to contradict the experimental findings in Zick et al., 2018. Indeed, in the neural recording experiments blocking NMDAR caused desynchronization of neurons, whereas in the simulations (Compte et al., 2000) the reduction of NMDAR currents provoked strong synchronization. Our model and theoretical analysis allows to explain this apparent paradox. In general, the asynchronous state becomes unstable and oscillation emerges when an excitatory feedback from the fast AMPA currents becomes sufficiently strong and is followed by a strong inhibitory feedback from the slower GABA currents (Brunel and Wang, 2003; Compte et al., 2000; Tsodyks et al., 1997; Wang, 1999). As explained above, the excitatory feedback can be enhanced via different mechanisms involving direct and indirect influence of synaptic parameters on the instability growth rate. In Compte et al., 2000, the concurrent increasing AMPAR and decreasing NMDAR conductances nullifies the indirect effect because contributions from the changes in the NMDAR and AMPAR mediated currents to the average total synaptic current, in essence, cancel each other. As a result, the operating point of the response function, defined by the average total current, does not change. However, due to the direct effect of the AMPAR on the instability growth rate $\lambda$, increasing AMPAR conductance enhances the excitatory-inhibitory feedback loop leading to the destabilization of the asynchronous activity and emergence of synchronous oscillations. In our model, by contrast, there is no direct effect on the instability growth rate because the AMPAR conductance is kept fixed, and the enhancement of recurrent excitatory feedback is entirely due to the indirect mechanism. It is achieved through external rate increase at a certain strength of the NMDAR conductance resulting in the neuron’s operating point shift toward a steeper slope above the point of the critical network. This induces network oscillation and synchronization of neurons as observed in monkey PFC when NMDAR is not blocked (Zick et al., 2018). However, when the NMDAR conductance is set to zero, the average total synaptic current is reduced, and the operating point moves down to such locus that it cannot be shifted above the point of the critical network by the same increase in external rate. As a result, external rate increase no longer provides a strong enough excitatory feedback, the network remains in asynchronous state, and no increase in synchrony occurs, consistent with observations in Zick et al., 2018 when NMDAR is blocked.

In monkeys performing a memory-guided saccade task, prefrontal neurons exhibit persistent activity that is associated with the maintenance of information in working memory (Chafee and Goldman-Rakic, 1998; Funahashi et al., 1989; Goldman-Rakic, 1995). Prior theoretical studies have investigated circuit and synaptic mechanisms that can generate persistent activity in recurrent prefrontal networks, specifically addressing how reducing NMDAR function destabilizes attractor states (patterns of stable neural activity) in these networks during a delay period (when the memory of the stimulus must be retained) leading to working memory deficits (Calvin and Redish, 2021; Compte et al., 2000; Funahashi et al., 1989; Goldman‐Rakic, 1987; Loh et al., 2007; Murray et al., 2014). In one seminal study by Compte et al., 2000, the authors investigated the robustness of working memory storage against external synaptic noise and distraction stimuli in attractor networks. They showed that a concomitant increase of NMDAR- and GABAR-mediated currents leads to an increase of persistent activity and to a decrease of spontaneous activity, thereby enhancing the resistance of the network to distractors (Brunel and Wang, 2001; Compte et al., 2000). In another prominent work, Murray et al., 2014, employing an attractor network model, investigated the neural and behavioral effects of synaptic disinhibition induced by the malfunction of NMDAR mediated synapses targeting inhibitory neurons. They demonstrated that disinhibition resulted in a broadening of stimulus selective persistent activity at the neural level, with a concomitant loss of precision, increase in variability over time, and increase in distractibility of stored information at the behavioral level. Although these modeling studies provide important mechanistic insight into prefrontal network dynamics underlying working memory, and potentially, working memory deficits in schizophrenia (Goldman-Rakic, 1999), they do not address the topic of the current study, which is how slow NMDAR recurrent excitation and external input received by the network jointly influence spike timing dynamics at the neuron level and oscillatory dynamics at the population level in the presence of fast AMPA excitation and GABA inhibition. Thus, no prior modeling study captures the relationship between NMDAR synaptic mechanisms, spike timing, and network oscillations that we have observed in neural recordings (Kummerfeld et al., 2020; Zick et al., 2018), and for which we provide a theoretical explanation in the current report.

We previously hypothesized that reduced synchrony at the level of spiking neurons (Zick et al., 2022; Zick et al., 2018) could disconnect prefrontal local circuits via spike-timing dependent synaptic plasticity (STDP; Dan and Poo, 2004; Feldman, 2012), contributing to the reduction in dendritic spine density that has been observed in postmortem analysis of prefrontal cortex in schizophrenia (Glantz and Lewis, 2000; MacDonald et al., 2017). However, the interaction between neural synchrony and synaptic connectivity in networks incorporating STDP is hard to predict, as changes in connectivity patterns and neural dynamics are mutually dependent and interact in complex ways as connectivity and synchrony influence each other over time. Perhaps for this reason, prior theoretical studies incorporating STDP into spiking networks have obtained divergent results with respect to how STDP changes the pattern of synaptic connections between neurons in networks, and whether synchronous inputs to the neurons are required for STDP to influence the pattern of synaptic connections. For example, STDP operating on random spiking in neurons can either lead to the formation of structured stable connections between neurons in the absence of synchronous inputs (forming neural ‘groups’) (Izhikevich et al., 2004), or not (Morrison et al., 2007), depending on the assumptions incorporated into the models. Similarly, correlated external input to recurrent networks incorporating STDP can either fail to produce structured synaptic connections between neurons (Morrison et al., 2007), or it can lead to the formation of such structured connections (Litwin-Kumar and Doiron, 2014) depending on the specifics of the simulations. Key parameters that could influence the diversity of outcomes among studies include whether (Izhikevich et al., 2004) or not (Morrison et al., 2007) axonal conduction delays and the geometry of recurrent connections are incorporated into the models (since circuit architecture and associated signal conduction delays powerfully influences when action potentials arrive at pre- and postsynaptic elements), as well as the specific form of the STDP rule employed (Babadi and Abbott, 2013; Bono and Clopath, 2017; Izhikevich et al., 2004; Morrison et al., 2007). Based on these results, it seems reasonable that distortions of spike timing dynamics in prefrontal networks may alter the pattern of neural connections via STDP in schizophrenia. However, the diversity of results obtained from theoretical studies of STDP outlined above make it difficult to conclude that the reduction in synchronous spiking we observed would lead to synaptic disconnection via STDP, imposing important constraints on our prior hypothesis (Zick et al., 2022; Zick et al., 2018), although this remains a possibility. Network simulations that accurately incorporate as many of these biological variables as possible may be useful in predicting how spike timing changes that may emerge downstream of schizophrenia risk factors would be likely to influence synaptic connectivity in the human cortex. In addition, as noted, genetic linkage studies have implicated altered NMDAR function in schizophrenia (Fromer et al., 2014; Schizophrenia Working Group of the Psychiatric Genomics Consortium, 2014). Since NMDAR play a central role in the molecular mechanisms that implement STDP in the brain, disruption of NMDAR synaptic transmission in schizophrenia may alter STDP directly, independently of the impact of disrupted NMDAR function on neural spiking dynamics in the disease state.

We had previously reported that blocking NMDAR in monkeys (Zick et al., 2018) and deleting a schizophrenia risk gene (Dgcr8) in mice (Zick et al., 2022), both reduced the frequency of synchronous, 0-lag spiking between prefrontal neurons. Dgcr8 encodes a protein involved in the synthesis of miRNA, which in turn bind to mRNA and suppress their translation into proteins, including mRNA coding for NMDAR subunits (Corbel et al., 2015). Deleting Dgcr8 would therefore be expected to reduce miRNA synthesis and increase translation of mRNA coding for NMDAR subunits. Given these considerations, the convergent spike desynchronization we observed in monkey drug and mouse genetic models could be explained by an inverted U-shaped relationship wherein either too little NMDAR function (as produced by NMDAR blockade in monkeys) or too much NMDAR function (as predicted to result from deletion of Dgcr8 in mice) decreases the frequency of 0-lag spiking between prefrontal neurons (Zick et al., 2022; Zick et al., 2018). An inverted U-shaped relationship has been reported between the level of D1 dopamine receptor stimulation and the strength of persistent neural activity in prefrontal neurons during working memory tasks wherein small doses of an agonist amplify persistent activity, and larger doses degrade it (Vijayraghavan et al., 2007). However, additional experimental data are needed to establish that spike synchrony exhibits a similar inverted U-shaped relation to NMDAR function, insofar as our prior neural recording studies did not test a U-shaped relationship directly (Zick et al., 2022; Zick et al., 2018). These studies did not for example contrast the effect of low versus high doses of an NMDAR agonist (such as NMDA) on spike synchrony in the monkey model, nor relate reduction in spike synchrony specifically to the upregulation of NMDAR subunit expression in the mouse model (rather than the many other proteins regulated by miRNA that are dependent on Dgcr8).

Results we present in the current study establish a theoretical basis and circuit mechanism explaining how reduction of NMDAR synaptic function implicated in schizophrenia could lead to the desynchronization of neural activity in prefrontal recurrent circuits. We provide evidence that spiking networks situated close to a boundary in the synaptic parameter space separating asynchronous and synchronous activity states can explain a variety of biological observations. These include the emergence of 0-lag synchronous spiking between individual prefrontal neurons when external inputs to the network push it across this state boundary, and failure of synchronous spiking to emerge between prefrontal neurons when NMDAR synaptic currents are reduced, as we have observed in neural recordings in primate prefrontal cortex (Kummerfeld et al., 2020; Zick et al., 2022; Zick et al., 2018).

For the present theoretical study, we used experimental data obtained in our previous work (Zick et al., 2018). Here, we provide brief descriptions of the experimental task, NMDAR antagonist regimen, and neurophysiological recording methodology employed in that work; details have been reported in Blackman et al., 2016; Zick et al., 2018.

To estimate correlation between spiking activity of simultaneously recorded neuron pairs as a function of time, we used a similar approach described in Zick et al., 2018. Correlation is evaluated from spiking activity observed during a time window $∆T$ around a given instant of time $t$. The window size $∆T$, thus, defines the temporal resolution of time resolved correlation. The interval $∆T$ is subdivided into small time bins of width $∆t$. Activity of neuron $i$ in a given trial at a time bin ${\displaystyle t^{\prime }}$ is represented by a binary variable ${\displaystyle {\xi }_i(t^{\prime })}$ that can take on two values: 1 if in the time bin ${\displaystyle t^{\prime }}$ one or more spikes are present, and 0 if there are no spikes. Correspondingly, time-lagged joint spike activity of neurons $i$ and $j$ is described by the product ${\displaystyle {\xi }_i(t^{\prime })\times {\xi }_j(t^{\prime }+\tau )}$: it is 1 if neuron $i$ fired a spike in the time bin ${\displaystyle t^{\prime }}$ and neuron $j$ fired a spike in the time bin ${\displaystyle t^{\prime }+\tau }$; otherwise, it is 0. The duration of the bin $∆t$, thus, defines the spike coincidence window. We assume that spike firing statistics of neurons do not change during the interval $∆T$, so that low order moments of the binary variables, such as the mean spike frequencies ${\displaystyle {\nu }_i=\overline{{\xi }_i(t^{\prime })}}$ and ${\displaystyle {\nu }_j=\overline{{\xi }_j(t^{\prime })}}$ and the mean joint spike frequency ${\displaystyle {\rho }_{ij}(\tau )=\overline{{\xi }_i(t^{\prime })\times {\xi }_j(t^{\prime }+\tau )}}$, can be reliably estimated by averaging over $∆T/∆t$ time bins (bars ${\displaystyle \overline{.}}$ above the expressions denote time averaging operation). To avoid a contribution to correlation from possible cross-trial non-stationarity (slow covariation) of neural activity, for each neuron pair correlation is estimated from single trials and then averaged over all trials. Spiking correlation between neurons $i$ and $j$ in a single trial is characterized by the observed frequency of joint spikes ${\displaystyle {\rho }_{ij}(\tau )}$ normalized by the expected joint spike frequency ${\nu }_i\times {\nu }_j$ if activity of the neurons were independent: ${\displaystyle {\rho }_{ij}(\tau )/({\nu }_i\times {\nu }_j)}$. We then average this ratio over the trials to obtain time-lagged correlation of spiking activity as ${\displaystyle c_{ij}(\tau )=⟨{\rho }_{ij}(\tau )/({\nu }_i\times {\nu }_j)⟩}$, where angular brackets $〈\bullet 〉$ denote trial averaging operation. Finally, ${\displaystyle c_{ij}(\tau )}$ is averaged over the population of simultaneously recorded pairs resulting in the population average spike correlation ${\displaystyle C(\tau )}$. Spike synchrony is defined as 0-lag correlation.

To accurately estimate spike synchrony and time-lagged correlation in PFC circuits, it is necessary to keep the value of time bin $∆t$, controlling the spike coincidence window, sufficiently small, within 1–2ms (no more than one spike occurred in a bin). On the other hand, the firing rates of PFC neurons are relatively low, on the order of 10 Hz. Therefore, to increase the number of counts of joint spike events and improve the estimate of spike synchrony while keeping $∆t$ small (and, thus, spike synchrony resolution sufficiently high), one needs to increase the duration of time window $∆T$ and/or the number of trials $K$ However, $∆T$ should be kept sufficiently short so that during this interval spiking activity remains nearly stationary, whereas $K$ cannot be made arbitrarily large because it is limited by practical considerations.

These experimental restrictions, as a result, impose constraints on the firing rates of the neurons in the pair. To derive a meaningful criterion for selecting ‘good’ neuron pairs, we note that for a reliable estimation of the mean joint spike firing frequency, which is a second order statistic, one needs quadratically more experimental samples than for a reliable estimation of the mean spike frequency, which is a first order statistic. We also note that the expected joint spike frequency if neurons in the pair were independent is simply given as the product of their mean spike frequencies. It is this quantity that is used as a reference (normalization) for the quantification of spike correlation strength. Therefore, to reliably estimate the joint spike firing frequency from available samples of a given pair, one should be confident that at least when assuming that neurons fire independently, a sufficiently accurate estimation of the expected joint spike frequency from these samples is possible. This, in turn, means that, given the neuron firing rates $v_i$ and $v_j$, the average total number of counts of joint spikes ${\displaystyle (v_i\mathrm{\Delta }t)(v_j\mathrm{\Delta }t)(\mathrm{\Delta }T/\mathrm{\Delta }t)K}$ observed in $∆T/∆t$ bins in $K$ trials predicted under the assumption of independence and calculated from experimental samples should be ‘detectable’, that is, it should be at least greater than 1. This condition results in a constraint for the geometric mean, ${\displaystyle {\overline{v}}_{ij}=\sqrt{v_i{v}_j}}$, of the firing rates of neuron pairs: ${\displaystyle {\overline{v}}_{ij}>1/\sqrt{K\mathrm{\Delta }T\mathrm{\Delta }t}}$. The typical values for the time window and spike coincidence window are $∆T~100$ ms and $∆t~1$ ms. Given that the number of correct trials in the DPX task were on the order of $K~200$, this means that the geometric mean firing rate of neuron pairs, for which a reliable estimation of synchrony can be achieved, should be at least 7 Hz.

The network consists of $N$ leaky integrate and fire neurons (see, e.g., Dayan and Abbott, 2001), of which $N_{\text{E}}=0.8N$ are excitatory and $N_{\text{I}}=0.2N$ are inhibitory (Abeles, 1991; Braitenberg and Schüz, 1998). Neurons are connected randomly with a probability $p$, so that, on average, each neuron receives $C_{\text{E}}=p{N}_{\text{E}}$ connections from excitatory and $C_{\text{I}}=p{N}_{\text{I}}$ from inhibitory neurons. In the framework of mean field consideration, the network is large ($N\gg 1$) and connections are sparse ($p\ll 1$) but the average number of connections received by individual neurons, $C$, is large ($C=pN\gg 1$). In most simulations, networks consisted of $N=5\bullet {10}^3$ neurons that were randomly connected with the probability $p=0.2$ and, therefore each neuron, on average, received $C={10}^3$ connections. In addition, each neuron also receives $C_{\text{X}}$ external connections from excitatory neurons outside of the network that fire spikes independently according to a Poisson process with rate ${\nu }_{\text{X}}$.

The dynamics of the membrane potential ${\displaystyle V(t)}$ of a neuron below the spike firing potential threshold $\theta$ obeys the standard leaky integrate and fire equation:

where $C_{\mathrm{m}}$ is the cell membrane capacitance, $g_{\mathrm{m}}$ is the membrane leak conductance, $V_{\mathrm{L}}$ is the resting potential, and ${\displaystyle i_{\mathrm{s}\mathrm{y}\mathrm{n}}(t)}$ is the total synaptic current. When the membrane potential reaches the threshold $\theta$, the neuron fires a spike, the potential is reset to $V_{\mathrm{r}\mathrm{s}\mathrm{t}}$, and the neuron becomes insensitive to its input for the duration of a refractory period ${\tau }_{\mathrm{r}\mathrm{p}}$. Both excitatory and inhibitory neurons have $\theta =-50$ mV, $V_{\mathrm{L}}=-70$ mV, and $V_{\mathrm{r}\mathrm{s}\mathrm{t}}=-55$ mV. For excitatory neurons $C_{\mathrm{m}}=0.5$ nF, $g_{\mathrm{m}}=25$ nS, ${\tau }_{\mathrm{r}\mathrm{p}}=2$ ms, and for inhibitory neurons $C_{\mathrm{m}}=0.2$ nF, $g_{\mathrm{m}}=20$ nS, ${\tau }_{\mathrm{r}\mathrm{p}}=1$ ms (see, e.g., Koch, 2004).

The total synaptic input for each neuron is a linear sum of four components:

where $i_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}$ and $i_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}$ correspond to recurrent excitatory currents mediated by AMPA and NMDA receptors, respectively, $i_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}$ corresponds to inhibitory currents mediated by GABA receptors, and $i_{\mathrm{X}}$ corresponds to external currents mediated by AMPA receptors. The purpose of external currents is twofold: (i) to represent the noisy inputs due to the background synaptic activity and (ii) to convey neural signals from outside of the network.

The description of component synaptic currents of a postsynaptic neuron follows Wang, 1999:

where synaptic reversal potentials $V_{\mathrm{E}}=0$ mV and $V_{\mathrm{I}}=-70$ mV. NMDAR mediated currents have voltage dependence controlled by the extracellular magnesium concentration (Jahr and Stevens, 1990): ${\displaystyle \beta =0.062\mathrm{m}{\mathrm{V}}^{-1}}$ , $\gamma =3.57$ mM, ${\displaystyle [\mathrm{M}{\mathrm{g}}^{2+}]=1}$ mM. The gating variable ${\displaystyle s_{R,j}(t)}$, describes the temporal course of postsynaptic currents received from the presynaptic neuron $j$ mediated by the receptor $R$, where $R=$ X, AMPA, NMDA, GABA. For a spike train generated by a presynaptic neuron with emission times ${\displaystyle \{t_k\}}$, the temporal dynamics of the gating variable obeys the equations

where ${\tau }_{\mathrm{l}}$, ${\tau }_{\mathrm{r}}$ and ${\tau }_{\mathrm{d}}$ are, respectively, latency, rising, and decay time constants. Their values are ${\tau }_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A},\mathrm{l}}=1$ mS, ${\tau }_{\mathrm{A}\mathrm{M}\mathrm{A}\mathrm{P},\mathrm{r}}=0.2$ ms, ${\tau }_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A},\mathrm{d}}=2$ ms for AMPAR mediated currents (Zhou and Hablitz, 1998), ${\tau }_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A},\mathrm{l}}=1$ mS, ${\tau }_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A},\mathrm{r}}=2$ ms, ${\tau }_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A},\mathrm{d}}=100$ ms for NMDAR mediated currents (Hestrin et al., 1990), and ${\tau }_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A},\mathrm{l}}=1$ mS, ${\tau }_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A},\mathrm{r}}=0.5$ ms, ${\tau }_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A},\mathrm{d}}=5$ ms for GABAR-mediated currents (Gupta et al., 2000). The time integral of ${\displaystyle s(t)}$ in response to a presynaptic spike equals ${\tau }_{\mathrm{*}}$ and, thus, is independent of the temporal shape of ${\displaystyle s(t)}$, which is determined by the rising and decay time constants that are specific to each receptor type. Because the charge flowing to the cell is determined by the product of the time integral of ${\displaystyle s(t)}$ and the maximal conductance, we set ${\tau }_{\mathrm{*}}$ to be the same for all types of receptors, so that the charge entry mediated by each type of receptor is parametrized, in essence, solely by the corresponding maximal conductance parameter.

In all direct network simulations, the numerical integration of the coupled differential equations describing the dynamics of membrane potentials and synaptic variables of all cells and synapses were carried out using a custom MATALAB (The MathWorks) code implementing a second order Runge-Kutta method with interpolation of spike firing times between integration time steps $∆t$ (Hansel et al., 1998). In most simulations $∆t=0.1$ ms.

To derive maximal synaptic conductance parameters $g_{\text{X},\alpha }$, $g_{\text{AMPA},\alpha }$, $g_{\text{NMDA},\alpha }$, $g_{\text{GABA},\alpha }$ ($\alpha =\text{E},\text{I}$) providing prescribed neural firing rates ${\nu }_{\text{E}}$ and ${\nu }_{\text{I}}$, we used mean field analysis (Amit and Brunel, 1997; Brunel, 2000; van Vreeswijk and Sompolinsky, 1996) extended to networks of neurons with realistic, conductance based synapses (Brunel and Wang, 2001; Renart et al., 2003). For simplicity, we disregard the heterogeneity of synaptic connectivity and assume that each neuron receives $C_{\text{E}}$ excitatory and $C_{\text{I}}$ inhibitory connections. In the mean field approximation synaptic inputs are described in terms of their average and their fluctuations arising from both external and recurrent inputs. To this end, the sums of gating variables in Equations 10–13 are replaced by their respective population averages ${\tau }_{*}{S}_R^0$, where $R$ designates the type of the synapse, and

The voltage dependence of NMDAR conductance is linearized around the mean value of the potential $〈V〉$:

where $\kappa =1+\left[{\mathrm{M}\mathrm{g}}^{2+}\right]/\gamma \mathrm{e}\mathrm{x}\mathrm{p}(-\beta 〈V〉)$. After these simplifications, average components of synaptic currents for excitatory ($\alpha =\text{E}$) and inhibitory ($\alpha =\text{I}$) populations can be written as

where $〈V_{\alpha }〉$ is the average membrane potential, and $J_{R,\alpha }$ is the effective strength of the $R$-receptor mediated synapse, expressed as the total charge entering the postsynaptic neuron due to a single presynaptic spike. In this framework, the system of equations describing the dynamics of membrane potentials for each of $N_{\text{E}}$ excitatory and $N_{\text{I}}$ inhibitory neurons is reduced to equations describing the dynamics of membrane potentials ${\displaystyle V_{\mathrm{E}}(t)}$ and ${\displaystyle V_{\mathrm{I}}(t)}$ of just two neurons representing, respectively, excitatory, $\text{E}$, and inhibitory, $\text{I}$, populations (Brunel and Wang, 2001; Renart et al., 2003):

where $V_{\text{L}}$ is the resting potential, ${\tau }_{\alpha }$ is the effective membrane time constant, ${\mu }_{\alpha }$ is the effective mean synaptic input, ${\sigma }_{\alpha }$ is the magnitude of the fluctuations in the synaptic input, and ${\displaystyle {\eta }_{\alpha }(t)}$ is the time course of these fluctuations: