A prefrontal network model operating near steady and oscillatory states links spike desynchronization and synaptic deficits in schizophrenia
Abstract
Schizophrenia results in part from a failure of prefrontal networks but we lack full understanding of how disruptions at a synaptic level cause failures at the network level. This is a crucial gap in our understanding because it prevents us from discovering how genetic mutations and environmental risks that alter synaptic function cause prefrontal network to fail in schizophrenia. To address that question, we developed a recurrent spiking network model of prefrontal local circuits that can explain the link between NMDAR synaptic and 0lag spike synchrony deficits we recently observed in a pharmacological monkey model of prefrontal network failure in schizophrenia. We analyze how the balance between AMPA and NMDA components of recurrent excitation and GABA inhibition in the network influence oscillatory spike synchrony to inform the biological data. We show that reducing recurrent NMDAR synaptic currents prevents the network from shifting from a steady to oscillatory state in response to extrinsic inputs such as might occur during behavior. These findings strongly parallel dynamic modulation of 0lag spike synchrony we observed between neurons in monkey prefrontal cortex during behavior, as well as the suppression of this 0lag spiking by administration of NMDAR antagonists. As such, our cortical network model provides a plausible mechanism explaining the link between NMDAR synaptic and 0lag spike synchrony deficits observed in a pharmacological monkey model of prefrontal network failure in schizophrenia.
Editor's evaluation
This valuable modeling study proposes a local circuit mechanism based on a network of recurrently connected excitatory and inhibitory neurons for the recently reported effect that NMDA receptor antagonists cause a drastic reduction of prefrontal neural synchronization in preparation for motor responses in a cognitive task. This mechanism is convincingly supported by simulations of spiking networks and a thorough analysis of the parameter dependency of network dynamics using meanfield theory. The work will be of general interest to computational neuroscientists, and especially for those interested in computational psychiatry.
https://doi.org/10.7554/eLife.79352.sa0eLife digest
Schizophrenia is a longterm mental health condition that can cause a person to see, hear or believe things that are not real. Although researchers do not fully understand the causes of schizophrenia, it is known to disrupt synapses, which connect neurons in the brain to form circuits that carry out a specific function when activated. This disruption alters the pattern of activity among the neurons, distorting the way that information is processed and leading to symptoms. Development of schizophrenia is thought to be due to interactions between many factors, including genetic makeup, changes in how the brain matures during development, and environmental stress.
Despite animal studies revealing how neural circuits can fail at the level of individual cells, it remains difficult to predict or understand the complex ways that this damage affects advanced brain functions. Previous research in monkeys showed that mimicking schizophrenia using a drug that blocks a particular type of synapse prevented neurons from coordinating their activity. However, this did not address how synaptic and cellular changes lead to disrupted neural circuits.
To better understand this, Crowe et al. developed a computational model of neural circuits to study how they respond to synapse disruption. To replicate the brain, the model consisted of two types of neurons – those that activate connecting cells in response to received signals and those that suppress them. This model could replicate the complex network behavior that causes brain cells to respond to sensory inputs. Increasing the strength of inputs to the network caused it to switch from a state in which the cells fired independently to one where the cells fired at the same time. As was previously seen in monkeys, blocking a particular type of synapse thought to be involved in schizophrenia prevented the cells from coordinating their signaling.
The findings suggest that schizophreniacausing factors can reduce the ability of neurons to fire at the same instant. Disrupting this process could lead to weaker and fewer synapses forming during brain development or loss of synapses in adults. If that is the case, and scientists can understand how factors combine to trigger this process, the mechanism of coordinated activity failure revealed by the model could help identify treatments that prevent or reverse the synapse disruption seen in schizophrenia.
Introduction
NMDAR synaptic malfunction has been implicated as causal in schizophrenia (Fromer et al., 2014; Schizophrenia Working Group of the Psychiatric Genomics Consortium, 2014; Timms et al., 2013), and loss of NMDAR synaptic function in prefrontal networks is believed to contribute to cognitive deficits as well as clinical symptoms in the disease (GoldmanRakic, 1999; Javitt et al., 2012; Wang et al., 2013). However, we do not have a complete understanding of how NMDAR synaptic mechanisms influence neural dynamics in prefrontal networks, nor how the disruption of NMDAR synaptic mechanisms might cause prefrontal networks to malfunction. To address these questions, we recently investigated how blocking NMDAR altered neural dynamics and effective communication between neurons in prefrontal cortex of monkeys performing a cognitive control task measuring deficits in schizophrenia (Blackman et al., 2013; Jones et al., 2010; Kummerfeld et al., 2020; Zick et al., 2018). We found that reducing NMDAR synaptic communication reduced the frequency of synchronous (‘0lag’) spiking between neurons, as well as effective communication between neurons on timescales consistent with monosynaptic interactions between them (Kummerfeld et al., 2020; Zick et al., 2022; Zick et al., 2018). Whereas these studies suggested that NMDAR synaptic function and spike timing in prefrontal networks were linked, they did not elucidate the circuit mechanisms responsible.
In the current study, we developed a spiking neural network model to understand mechanisms that might mediate the link between NMDAR synaptic malfunction and neural dynamics (reduced 0lag synchronous spiking) we observed in biological data (Kummerfeld et al., 2020; Zick et al., 2022; Zick et al., 2018). The network is comprised of leaky integrateandfire neurons embedded in a sparsely connected recurrent network employing realistic NMDAR, GABAR, and AMPAR mediated synaptic currents. We use network stability and mean field analyses to investigate how the balance between NMDA and AMPA components of recurrent excitatory and GABA inhibitory currents influence regimes of network dynamics and spiking synchrony.
For cortical neurons synchrony can occur naturally due to the local recurrent network connectivity, even when external afferent inputs are entirely uncorrelated. Theoretical studies have shown that such synchrony can arise in randomly connected recurrent networks operating in asynchronous irregular (Amit, 1989; Amit and Brunel, 1997; Brunel, 2000; Renart et al., 2010; van Vreeswijk and Sompolinsky, 1996; Vicente et al., 2008) and synchronous irregular regimes (Brunel, 2000; Brunel and Hakim, 1999; Brunel and Wang, 2003; Ledoux and Brunel, 2011). In both regimes, individual neurons fire spikes highly irregularly at low rates, a typical situation in a cortex. The major distinction is that in an asynchronous regime population spike rate is steady in time, whereas in a synchronous regime it becomes oscillatory.
We show that simulated prefrontal networks operating near the boundary between steady (asynchronous irregular) and oscillatory (synchronous irregular) regimes in the synaptic parameter space can explain several key experimental observations. First, such networks achieve biologically realistic stochastic spike trains and firing rates of excitatory and inhibitory neurons in prefrontal cortex. Second, increased extrinsic inputs, such as those that might occur during behavior, shift these networks from a steady to an oscillatory regime that causes the emergence of 0lag spiking between neurons as they stochastically entrain to oscillatory population activity. Third, and perhaps most importantly, we show that reducing recurrent NMDAR synaptic currents prevents these networks from transitioning into oscillatory activity in response to extrinsic inputs, thereby preventing the emergence of 0lag spike synchrony. Although prior modeling studies have addressed the relationship between NMDAR function and oscillatory activity in prefrontal networks (Brunel and Wang, 2003; Compte et al., 2000; Kirli et al., 2014; Wang, 1999), none account for this range of experimental observations. The current results allow us to establish strong parallels between simulated and biological data, including the emergence of 0lag synchronous spiking via recurrent synaptic interactions between neurons during behavior, the association between synchronous spiking and oscillatory population activity, as well as their joint dependence on NMDAR synaptic mechanisms, both in our current simulation and in the neural data (Zick et al., 2018).
Results
Summary of experimental results
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 dotpattern 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: drugnaive 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 drugnaive (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 mslong) 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 drugnaive 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 drugnaive 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.
Network model and theoretical framework
To understand the phenomenon of druginduced 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 integrateandfire 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 $\nu}_{\mathrm{X}$.
Recurrent synaptic currents of excitatory connections are twocomponent, 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 $g}_{\mathrm{X},\alpha$, $g}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A},\alpha$, $g}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A},\alpha$, $g}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A},\alpha$ corresponding to the external AMPA, recurrent AMPA, NMDA, and GABA currents ($\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: $I}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}/{I}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}$, $I}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}/{I}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}$, and $I}_{\mathrm{X},\mathrm{E}}/{I}_{\theta ,\mathrm{E}$, where ${I}_{R}$ is the mean current of the $R$receptor mediated synapse ($R=\mathrm{X},\phantom{\rule{thinmathspace}{0ex}}\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A},\phantom{\rule{thinmathspace}{0ex}}\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A},\phantom{\rule{thinmathspace}{0ex}}\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 $\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 $\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 $\left({I}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}/{I}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}},\phantom{\rule{thinmathspace}{0ex}}{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 $\left({I}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}/{I}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}},\phantom{\rule{thinmathspace}{0ex}}{I}_{\mathrm{X},\mathrm{E}}/{I}_{\theta ,\mathrm{E}}\right)$ parameter space. The asynchronous stationary state corresponds to the region where $\lambda <0$, whereas the synchronous oscillation state is realized in the region where $\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.
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).
Integration of DPX task context and drug condition into the model
To study spike synchrony in asynchronous and synchronous networks in the context of the DPX task performed in drugnaive 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 taskspecific 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.
Selection of primary networks
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 $\left({I}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}/{I}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}},\phantom{\rule{thinmathspace}{0ex}}{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 $\left({I}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}/{I}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}},\phantom{\rule{thinmathspace}{0ex}}{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 $\left({I}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}^{\ast}/{I}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}^{\ast}=0.4,\phantom{\rule{thinmathspace}{0ex}}{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 ($\lambda <0$). The point marked by a red asterisk in Figure 2A located at $\left({I}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}^{\ast}/{I}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}^{\ast}=0.2,\phantom{\rule{thinmathspace}{0ex}}{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}^{*}$ ($\alpha =\mathrm{E},\mathrm{I}$) by numerically solving the mean field equations.
Correlation of spiking activity and synchrony in the asynchronous and synchronous states
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).
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 drugnaive and drug conditions studied in Zick et al., 2018 and with the purpose of elucidating the mechanism of druginduced 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 (0lag 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).
Circuit mechanisms of spike synchronization modulation
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 networkwide synchronization of spiking activity, we carried out linear stability analysis in the $\left({v}_{\mathrm{X}}/{v}_{\mathrm{X}}^{\ast},\phantom{\rule{thinmathspace}{0ex}}{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 $\left({v}_{\mathrm{X}}/{v}_{\mathrm{X}}^{\ast},\phantom{\rule{thinmathspace}{0ex}}{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 ($\lambda <0$) and synchronous oscillatory ($\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).
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 networkwide synchronization of spiking activity in response to these modulations.
Explaining the effects of blocking of NMDAR observed in primate PFC by the prefrontal circuit model
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 drugnaive (black) and drug (magenta) conditions. It can be now appreciated that in drugnaive condition, population activity during the preresponse 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).
The presence of strong spike synchrony (0ms lag) together with the correlation peaks at ±18 ms lags in the preresponse 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 preresponse 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 drugdependent desynchronization can be illustrated by considering the system’s behavior in the $\left({v}_{\mathrm{X}}/{v}_{\mathrm{X}}^{\ast},\phantom{\rule{thinmathspace}{0ex}}{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 $\left({v}_{\mathrm{X}}/{v}_{\mathrm{X}}^{\ast}=1,\phantom{\rule{thinmathspace}{0ex}}{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 drugnaive (${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 (${v}_{\mathrm{X}}/{v}_{\mathrm{X}}^{\ast}=0.97$) to the network state corresponding to the preresponse period (${v}_{\mathrm{X}}/{v}_{\mathrm{X}}^{\ast}=1.03$).
In drugnaive 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 preresponse 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 drugnaive 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 trapezoidlike 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 drugnaive 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.
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 drugnaive and drug conditions.
Role of the balance between NMDAR mediated recurrent excitation and GABA inhibition
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 $\left({I}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}/{I}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}},\phantom{\rule{thinmathspace}{0ex}}{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 $\left({v}_{\mathrm{X}}/{v}_{\mathrm{X}}^{\ast},\phantom{\rule{thinmathspace}{0ex}}{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 $\left({v}_{\mathrm{X}}/{v}_{\mathrm{X}}^{\ast},\phantom{\rule{thinmathspace}{0ex}}{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 (${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 (${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 drugnaive 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).
Dependence of oscillatory instability growth rate on synaptic parameters
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 $\u2206{g}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}/{g}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}^{*}$, $\u2206{g}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}/{g}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}^{*}$, $\u2206{g}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}/{g}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}^{*}$, and $\u2206{v}_{\mathrm{X}}/{v}_{\mathrm{X}}^{*}$ of the synaptic parameters around the critical point $\left\{{g}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A},\left\{\mathrm{E},\mathrm{I}\right\}}^{\ast},\phantom{\rule{thinmathspace}{0ex}}{g}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A},\left\{\mathrm{E},\mathrm{I}\right\}}^{\ast},\phantom{\rule{thinmathspace}{0ex}}{g}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A},\left\{\mathrm{E},\mathrm{I}\right\}}^{\ast},\phantom{\rule{thinmathspace}{0ex}}{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 ($R=\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A},\phantom{\rule{thinmathspace}{0ex}}\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A},\phantom{\rule{thinmathspace}{0ex}}\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}$) are again scaled with the same factors and, thus, their relative changes are equal: $\u2206{g}_{R,\mathrm{E}}/{g}_{R,\mathrm{E}}^{*}=\u2206{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 currentfrequency response function at the critical state, and $\u2206{\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 ($\alpha =\mathrm{E},\mathrm{I}$ for excitatory and inhibitory neurons, respectively). The deviations of the synaptic conductances $\u2206{g}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}$, $\u2206{g}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}$, $\u2206{g}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}$, and external rate $\u2206{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 $\u2206{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 currentfrequency 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 integrateandfire neuron model, $\varphi $ is a monotonically increasing nonlinear 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 $\u2206{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. $\u2206\varphi $), but also in the change of the slope of the response function $\u2206{\varphi}_{{I}_{\mathrm{s}\mathrm{y}\mathrm{n}}}^{`}$. The latter can be calculated by linearizing the selfconsistent mean field equations (see Materials and methods). As a result, $\u2206{\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 $\u2206{\varphi}_{{I}_{\mathrm{s}\mathrm{y}\mathrm{n}}}^{`}$ is proportional to the change $\u2206{I}_{\mathrm{s}\mathrm{y}\mathrm{n}}$ :
where ${I}_{0}$ is a positive constant. Hence, from Equation 2 it follows that $\u2206{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 $\u2206{v}_{\mathrm{X}}$. However, $\lambda $ depends on $\u2206{v}_{\mathrm{X}}$ indirectly via the terms involving the change in the slope $\u2206{\varphi}_{{I}_{\mathrm{s}\mathrm{y}\mathrm{n}}}^{`}$ due to the change in the average total synaptic current $\u2206{I}_{\mathrm{s}\mathrm{y}\mathrm{n}}$ (Equations 3, 4). In fact, $\u2206{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 currentfrequency 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 $\u2206{g}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}$, $\u2206{g}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}$, and $\u2206{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 $\u2206{g}_{R}$ ($R=\mathrm{A}\mathrm{M}\mathrm{A}\mathrm{P}\mathrm{A},\phantom{\rule{thinmathspace}{0ex}}\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A},\phantom{\rule{thinmathspace}{0ex}}\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}$) is determined only by the corresponding factor ${\mathrm{\Lambda}}_{R}$ via the term ${\mathrm{\Lambda}}_{R}\u2206{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 $\u2206{\varphi}_{{I}_{\mathrm{s}\mathrm{y}\mathrm{n},\mathrm{E}}}^{`}$ and $\u2206{\varphi}_{{I}_{\mathrm{s}\mathrm{y}\mathrm{n},\mathrm{I}}}^{`}$, which depend on $\u2206{g}_{R}$ (Equation 2). For example, the strength of direct contribution to $\lambda $ due to the change in the GABAR conductance $\u2206{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}}\u2206{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 $\u2206{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}}\u2206{\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}}\u2206{\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}}\u2206{\varphi}_{{I}_{\mathrm{s}\mathrm{y}\mathrm{n},\mathrm{I}}}^{`}/{\varphi}_{{I}_{\mathrm{s}\mathrm{y}\mathrm{n},\mathrm{I}}}^{`}$ in Equation 1 because $\u2206{\varphi}_{{I}_{\mathrm{s}\mathrm{y}\mathrm{n},\mathrm{E}}}^{`}$ and $\u2206{\varphi}_{{I}_{\mathrm{s}\mathrm{y}\mathrm{n},\mathrm{I}}}^{`}$ themselves depend on $\u2206{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 currentfrequency 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).
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 $\u2206{g}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}$ and $\u2206{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 $\u2206{\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 $\u2206{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 (drugnaive 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 $\u2206{v}_{\mathrm{X}}$ and NMDAR conductance $\u2206{g}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}$ both affect synchrony via indirect mechanism by changing the excitatory drive $I}_{\mathrm{s}\mathrm{y}\mathrm{n}$ and, therefore, shifting the operating point of the neuron’s response function. In the drugnaive condition, increase in external rate ($\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 ($\mathrm{\Delta}{g}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}<0$), the initial excitatory drive is reduced compared to the drugnaive condition, and now the same increase in external rate $\u2206{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.
Discussion
To better understand how synaptic mechanisms influence neural synchrony in recurrent local circuits in monkey prefrontal cortex, we developed a theoretical framework employing a sparsely connected recurrent network model accounting for AMPAR, NMDAR, and GABAR mediated synaptic currents. This allowed us to examine how varying combinations of synaptic transmission in the recurrent network influenced spike timing at the level of pairs of neurons and oscillatory dynamics at the level of neural populations. Our motivation to pursue this question derives from recent neurophysiological experiments investigating the impact of pharmacological NMDAR blockade on spike timing dynamics in monkey prefrontal cortex (Kummerfeld et al., 2020; Zick et al., 2022; Zick et al., 2018). These studies were initiated to investigate how risk factors associated with schizophrenia alter neural dynamics in prefrontal cortex. Those studies found that pharmacological and genetic factors associated with schizophrenia convergently reduce 0lag synchronous spiking between pairs of prefrontal neurons in monkeys and mice (Zick et al., 2022). The spiking network model we develop in the present study provides a circuit mechanism capable of explaining the biological data. The principal features of this circuit mechanism are as follows: (i) synaptic conductance parameters of the underlying circuit are such that it is in an asynchronous state near a critical boundary in the (NMDAR conductance – external input) parameter plane separating asynchronous and synchronous network states, (ii) small increases in extrinsic inputs push the circuit past this critical boundary into the region of a synchronous state, causing emergence of gamma oscillations in population activity, (iii) 0lag synchronous spiking between neurons emerges as they stochastically entrain to the gamma population rhythm, (iv) blocking NMDAR currents prevents the circuit from switching to a synchronous regime in response to external inputs, (v) thereby precluding emergence of 0lag synchronous spiking in neurons.
This circuit mechanism offers a reasonable explanation accounting for the tasklocked increase in 0lag spike synchrony that occurs in monkey prefrontal cortex just before the motor response in the cognitive control task (Zick et al., 2018): the increase in synchrony could reflect increased synaptic input to prefrontal networks at around this time, potentially from mediodorsal nucleus of thalamus (DeNicola et al., 2020). It also explains why pharmacological blockade of NMDAR attenuates 0lag spike synchrony before the motor response: the deficit in NMDAR mediated synaptic currents prevents prefrontal networks from switching to a synchronous regime in response to external inputs.
In the circuit model, the balance between the AMPA component of recurrent excitation and GABA inhibition controls the network frequency at the onset of oscillation, consistent with results in Brunel and Wang, 2003. This frequency is virtually independent of the balance between the tonic component of recurrent excitation mediated by the NMDAR and GABA inhibition. However, the balance between the NMDA and GABA currents determines the strength of modulation of the external synaptic input needed for switching between the asynchronous stationary and synchronous oscillatory states in the absence and presence of NMDAR antagonist.
Firing rate and synaptic mechanisms jointly influence synchronous spiking
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 $\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 currentoutput frequency response function. The direct dependence manifests the essential influence of the AMPAR synaptic conductance on the strength of an excitatoryinhibitory 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 currentfrequency 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.
Relation to prior studies of NMDAR function and oscillatory dynamics
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 excitatoryinhibitory 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.
Relation to prior studies of NMDAR function and working memory
In monkeys performing a memoryguided saccade task, prefrontal neurons exhibit persistent activity that is associated with the maintenance of information in working memory (Chafee and GoldmanRakic, 1998; Funahashi et al., 1989; GoldmanRakic, 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 GABARmediated 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 (GoldmanRakic, 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.
Spike timing disruptions and rewiring of prefrontal local circuits via STDP
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 spiketiming 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 (LitwinKumar 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.
Potential Ushaped relation between NMDAR function and spike synchrony
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, 0lag 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 Ushaped 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 0lag spiking between prefrontal neurons (Zick et al., 2022; Zick et al., 2018). An inverted Ushaped 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 Ushaped relation to NMDAR function, insofar as our prior neural recording studies did not test a Ushaped 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 0lag 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).
Materials and methods
Experimental data
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.
Experimental task
Request a detailed protocolMale rhesus macaque monkeys (8–10 kg) were trained to perform the dotpattern expectancy (DPX) task. This task is closely related to the AXCPT (continuous performance task) except that dot patterns replace letters as stimuli. During each trial of the DPX tasks, monkeys maintained gaze fixated on a central target as a cue stimulus (1,000ms), followed by a delay period (1,000ms), and a probe stimulus (500ms) were presented. Monkeys were rewarded for moving a joystick to the left if the cueprobe sequence had been AX (69% of trials), or to the right if any other cueprobe sequence had been presented (AY, BX, BY, collectively 31% of trials). Since the correct response to the Xprobe depended on the preceding cue (A or B), the task required both working memory and cognitive control. Both The DPX and AXCPT measure specific cognitive control impairments in schizophrenia (Barch et al., 2003; Jones et al., 2010).
Neurophysiological recording
Request a detailed protocolIn our previous study (Zick et al., 2018), we recorded neural activity from the region of the principal sulcus (centered on Brodmann’s areas 46) in the dorsolateral prefrontal cortex of two macaques performing the DPX task. We found that 0lag synchrony while present in both monkeys was much stronger in one than the other animal. For comparison to spiking dynamics in the present neural network simulation, we used neurophysiological recording data from the monkey that exhibited the strongest 0lag spike correlation during task performance (Zick et al., 2018). For neurophysiological recording, we used a computercontrolled electrode drive (System Eckhorn, Thomas Recording, GmbH) advancing 16, closely spaced, independently movable glass coated platinum/tungsten microelectrodes into the prefrontal cortex. Electrodes were spaced 400 µm apart, and interelectrode distances in the array spanned 400–1,400 µm. Moving the electrodes in depth and the position of the array within recording chambers over days made it possible to isolate the spiking activity of different neural ensembles, each containing 15–30 individually isolated, simultaneously recorded neurons. The database included in the present study consisted of 47 neural ensembles containing a total of 893 prefrontal neurons. Spike correlation was evaluated within ensembles of simultaneously recorded ensembles using spike trains recorded during DPX task performance (Zick et al., 2018).
NMDAR antagonist regimen
Request a detailed protocolWe examined the effect of systemic administration of an NMDAR antagonist (phencyclidine, 0.25–0.30 mg/kg IM) on spike timing dynamics in prefrontal local circuits. Neural activity was recorded in a Naive condition (before first exposure to drug), and a Drug condition (following systemic drug administration) (Zick et al., 2018).
Spike correlation and synchrony
Request a detailed protocolTo 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 $\u2206T$ around a given instant of time $t$. The window size $\u2206T$, thus, defines the temporal resolution of time resolved correlation. The interval $\u2206T$ is subdivided into small time bins of width $\u2206t$. Activity of neuron $i$ in a given trial at a time bin $t}^{\prime$ is represented by a binary variable ${\xi}_{i}({t}^{\prime})$ that can take on two values: 1 if in the time bin $t}^{\prime$ one or more spikes are present, and 0 if there are no spikes. Correspondingly, timelagged joint spike activity of neurons $i$ and $j$ is described by the product ${\xi}_{i}({t}^{\prime})\times {\xi}_{j}({t}^{\prime}+\tau )$: it is 1 if neuron $i$ fired a spike in the time bin $t}^{\prime$ and neuron $j$ fired a spike in the time bin ${t}^{\prime}+\tau$; otherwise, it is 0. The duration of the bin $\u2206t$, thus, defines the spike coincidence window. We assume that spike firing statistics of neurons do not change during the interval $\u2206T$, so that low order moments of the binary variables, such as the mean spike frequencies $\nu}_{i}=\overline{{\xi}_{i}({t}^{\prime})$ and $\nu}_{j}=\overline{{\xi}_{j}({t}^{\prime})$ and the mean joint spike frequency $\rho}_{ij}(\tau )=\overline{{\xi}_{i}({t}^{\prime})\times {\xi}_{j}({t}^{\prime}+\tau )$, can be reliably estimated by averaging over $\u2206T/\u2206t$ time bins (bars $\overline{.}$ above the expressions denote time averaging operation). To avoid a contribution to correlation from possible crosstrial nonstationarity (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 ${\rho}_{ij}(\tau )$ normalized by the expected joint spike frequency ${\nu}_{i}\times {\nu}_{j}$ if activity of the neurons were independent: ${\rho}_{ij}(\tau )/({\nu}_{i}\times {\nu}_{j})$. We then average this ratio over the trials to obtain timelagged correlation of spiking activity as ${c}_{ij}(\tau )=\u27e8{\rho}_{ij}(\tau )/({\nu}_{i}\times {\nu}_{j})\u27e9$, where angular brackets $\u2329\bullet \u232a$ denote trial averaging operation. Finally, ${c}_{ij}(\tau )$ is averaged over the population of simultaneously recorded pairs resulting in the population average spike correlation $C(\tau )$. Spike synchrony is defined as 0lag correlation.
To accurately estimate spike synchrony and timelagged correlation in PFC circuits, it is necessary to keep the value of time bin $\u2206t$, 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 $\u2206t$ small (and, thus, spike synchrony resolution sufficiently high), one needs to increase the duration of time window $\u2206T$ and/or the number of trials $K$ However, $\u2206T$ 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 $({v}_{i}\mathrm{\Delta}t)({v}_{j}\mathrm{\Delta}t)(\mathrm{\Delta}T/\mathrm{\Delta}t)K$ observed in $\u2206T/\u2206t$ 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, $\overline{v}}_{ij}=\sqrt{{v}_{i}{v}_{j}$, of the firing rates of neuron pairs: $\overline{v}}_{ij}>1/\sqrt{K\mathrm{\Delta}T\mathrm{\Delta}t$. The typical values for the time window and spike coincidence window are $\u2206T~100$ ms and $\u2206t~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.
Network model
Request a detailed protocolThe 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 $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 ${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): $\beta =0.062\phantom{\rule{thinmathspace}{0ex}}\mathrm{m}{\mathrm{V}}^{1}$ , $\gamma =3.57$ mM, $[\mathrm{M}{\mathrm{g}}^{2+}]=1$ mM. The gating variable ${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 $\{{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 GABARmediated currents (Gupta et al., 2000). The time integral of $s(t)$ in response to a presynaptic spike equals ${\tau}_{\mathrm{*}}$ and, thus, is independent of the temporal shape of $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 $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.
Network simulations
Request a detailed protocolIn 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 RungeKutta method with interpolation of spike firing times between integration time steps $\u2206t$ (Hansel et al., 1998). In most simulations $\u2206t=0.1$ ms.
Mean field approximation
Request a detailed protocolTo 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 $\u2329V\u232a$:
where $\kappa =1+\left[{\mathrm{M}\mathrm{g}}^{2+}\right]/\gamma \mathrm{e}\mathrm{x}\mathrm{p}(\beta \u2329V\u232a)$. After these simplifications, average components of synaptic currents for excitatory ($\alpha =\text{E}$) and inhibitory ($\alpha =\text{I}$) populations can be written as
where $\u2329{V}_{\alpha}\u232a$ 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 ${V}_{\mathrm{E}}(t)$ and ${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 ${\eta}_{\alpha}(t)$ is the time course of these fluctuations:
In the absence of spiking and fluctuations, the average membrane potential would equal ${\mu}_{\alpha}+{V}_{\text{L}}$ (Equation 21). The average membrane potential $\u2329{V}_{\alpha}\u232a$ of spiking neuron in the presence of synaptic noise can be calculated from the distribution of potentials obtained in Brunel and Hakim, 1999 and is given by (Renart et al., 2003)
The total synaptic noise ${\sigma}_{\alpha}^{2}$ characterizing fluctuations in the input that result from random arrival of spikes is approximated as the sum of the fluctuations in the external and recurrent inputs (Fourcaud and Brunel, 2002):
where
${\eta}_{\alpha}(t)$ is a Gaussian process with zero mean, $\u27e8{\eta}_{\alpha}(t)\u27e9=0$, and an exponentially decaying correlation function, $\u27e8{\eta}_{\alpha}(t){\eta}_{\alpha}({t}^{\prime})\u27e9\propto \mathrm{e}\mathrm{x}\mathrm{p}\left(\leftt{t}^{\prime}\right/{\tau}_{\mathrm{s}\mathrm{y}\mathrm{n},\alpha}\right)$, which is due to synaptic filtering with effective time constant ${\tau}_{\mathrm{s}\mathrm{y}\mathrm{n},\alpha}$ (Fourcaud and Brunel, 2002):
where ${{\tau}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}=\tau}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A},\mathrm{l}}+{\tau}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A},\mathrm{r}}+{\tau}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A},\mathrm{d}}$, ${{\tau}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}=\tau}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A},\mathrm{l}}+{\tau}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A},\mathrm{r}}+{\tau}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A},\mathrm{d}}$, ${{\tau}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}=\tau}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A},\mathrm{l}}+{\tau}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A},\mathrm{r}}+{\tau}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A},\mathrm{d}}$ are effective synaptic time constants for AMPAR, NMDAR, and GABARmediated currents, respectively. In addition, because of sparse connectivity, the correlation of the fluctuations in the synaptic inputs of excitatory and inhibitory populations is neglected: $\u27e8{\eta}_{\mathrm{E}}(t){\eta}_{\mathrm{I}}({t}^{\prime})\u27e9=0$. The firing rate ${\nu}_{\alpha}$ of a neuron, whose potential is governed by Equation 21, is given by a currentfrequency relationship ${\varphi}_{\alpha}\left({\mu}_{\alpha},{\sigma}_{\alpha}\right)$ that is a function of the mean and fluctuating part of synaptic input (Brunel and Sergi, 1998; Fourcaud and Brunel, 2002):
where
Since ${\mu}_{\alpha}$ and ${\sigma}_{\alpha}$ themselves depend on the population firing rates ${\nu}_{\text{E}}$ and ${\nu}_{\text{I}}$, the two coupled frequencycurrent equations
provide a selfconsistent description of the network in stationary states, that is regimes of network dynamics when the population average quantities such as firing rates and synaptic inputs are constant in time. In the framework of our model, synaptic conductances ${g}_{\text{X},\alpha}$, ${g}_{\text{AMPA},\alpha}$, ${g}_{\text{NMDA},\alpha}$, ${g}_{\text{GABA},\alpha}$ ($\alpha =\text{E},\text{I}$) and the external spike rate ${\nu}_{\text{X}}$ are system parameters controlling the regime of network dynamics; they enter to the mean field analysis through expressions for ${\mu}_{\alpha}$, and ${\sigma}_{\alpha}$,. If these parameters are given, one can solve the selfconsistent equations to obtain predicted by the mean field approximation population firing rates ${v}_{\mathrm{E}}^{0}$ and ${v}_{\mathrm{I}}^{0}$ in a stationary state of the network. Conversely, once external ${\nu}_{\text{X}}$ and population spike rates ${v}_{\mathrm{E}}^{0}$ and ${v}_{\mathrm{I}}^{0}$ are specified, the selfconsistent equations could be solved to find the values of synaptic conductance parameters ${g}_{\text{X},\alpha}$, ${g}_{\text{AMPA},\alpha}$, ${g}_{\text{NMDA},\alpha}$, ${g}_{\text{GABA},\alpha}$ ($\alpha =\text{E},\text{I}$) that correspond to these spike rates. However, because there are eight unknown parameters and only two equations, to find a unique solution one would need six additional equations imposing constraints on conductance parameters.
Model parametrization
Request a detailed protocolWe derive three of these equations by implementing a commonly used constraint (e.g. Brunel and Wang, 2003; Compte et al., 2000) that equalizes the ratio of synaptic conductance parameters for component currents in excitatory and inhibitory neurons. Since each component current is proportional to its respective synaptic conductance, this constraint implies that the balance between different components of average synaptic currents ${I}_{\text{X},\alpha}^{0}$, ${I}_{\text{AMPA},\alpha}^{0}$, ${I}_{\text{NMDA},\alpha}^{0}$, ${I}_{\text{GABA},\alpha}^{0}$ for excitatory ($\alpha =\text{E}$) and inhibitory ($\alpha =\text{I}$) populations is the same, thus providing the following three equations:
As a result, whenever the ratio of synaptic conductances and/or component currents is involved, the index $\alpha $ designating the type of the neuron can be dropped.
Two additional equations are obtained by fixing the balance between inhibition and twocomponent recurrent excitation at certain values:
The last constraint is provided in terms of the relative magnitude of average external current of excitatory neurons, ${I}_{\mathrm{X},\mathrm{E}}^{0}$:
where ${I}_{\theta ,\mathrm{E}}^{0}$ is the current that is needed for an excitatory neuron to reach firing threshold $\theta $ in absence of recurrent feedback. This approach allowed to parametrize network dynamics in terms of three parameters expressed as ratios of absolute values of average synaptic currents, ${I}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}/{I}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}$, ${I}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}/{I}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}$, and $I}_{\mathrm{X},\mathrm{E}}/{I}_{\theta ,\mathrm{E}$, characterizing the balance between components of recurrent excitation and inhibition, and the balance between external input and firing threshold. For a given external spike rate ${\nu}_{\mathrm{X}}$ and fixed values of these three parameters, we are now able to solve the selfconsistent equations for the eight synaptic conductances that provide the prescribed population firing rates ${v}_{\mathrm{E}}^{0}$ and ${v}_{\mathrm{I}}^{0}$ in a stationary state of the network.
We are interested in the asynchronous stationary state in which neurons fire spikes irregularly and at low rates, like neurons in prefrontal cortex. When mean synaptic inputs ${\mu}_{\alpha}$ are well below threshold $\theta $, firing is driven by the synaptic fluctuations ${\sigma}_{\alpha}$ around the mean input, therefore, resulting in irregular spike trains and low rates (Renart et al., 2003). Given that the number of synaptic connections received by individual neurons is large and network connectivity is sparse, solutions of selfconsistent equations providing the subthreshold regime for ${\mu}_{\alpha}$ and, thus, low rate asynchronous network dynamics, arise when inhibition strongly dominates recurrent excitation and the mean external inputs are around or above threshold $\theta $ (Brunel, 2000; Renart et al., 2003; van Vreeswijk and Sompolinsky, 1996). Thus, for the network to be in asynchronous irregular state the three system parameters characterizing the balance between recurrent excitation and inhibition, and the relative strength of external inputs should be within certain bounds: ${I}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}/{I}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}+{I}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}/{I}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}<1$, and ${I}_{\mathrm{X},\mathrm{E}}/{I}_{\theta ,\mathrm{E}}\gtrsim 1$.
Linear stability analysis
Request a detailed protocolWe perform a linear stability analysis of the asynchronous state (Abbott and van Vreeswijk, 1993; Brunel and Hakim, 1999) on the basis of an analytical consideration in Brunel and Wang, 2003. To understand if the network develops instability caused by fluctuations in population firing rates, we consider small deviations from the stationary population rates ${v}_{\mathrm{E}}^{0}$ and ${v}_{\mathrm{I}}^{0}$. In order to analyze the resulting network behavior, the mean field approach and selfconsistent equations providing population mean firing rates ${v}_{\mathrm{E}}^{0}$ and ${v}_{\mathrm{I}}^{0}$ are extended to describe the dynamics of population rates ${\nu}_{\mathrm{E}}(t)$ and ${\nu}_{\mathrm{I}}(t)$.
In the framework of mean field approximation, each component of synaptic current is determined by the product of effective synaptic strength $J$ and average gating variable $S$ (Equations 17–20 for the steady state consideration). The dynamics of $S$ is governed by the same type of equations as for the gating variable $s$ of an individual synapse in a given postsynaptic neuron (Equation 14), except that the instantaneous rate of spikes $\sum _{k}\delta \left(t{t}_{k}{\tau}_{\mathrm{l}}\right)$ arriving from the presynaptic cell is replaced by the instantaneous average rate of spikes, ${C}_{{\alpha}_{R}}{\nu}_{{\alpha}_{R}}(t{\tau}_{\mathrm{l}})$, arriving from all presynaptic cells making the same type of synapse in the postsynaptic neuron:
where $R$ designates the type of the synapse ($R=\mathrm{X},\phantom{\rule{thinmathspace}{0ex}}\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A},\phantom{\rule{thinmathspace}{0ex}}\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A},\phantom{\rule{thinmathspace}{0ex}}\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}$), and ${\alpha}_{R}$ designates the presynaptic population establishing these synapses ($\alpha}_{R}=\mathrm{X},\mathrm{E$ for glutamatergic and ${\alpha}_{R}=\mathrm{I}$ for GABAergic synapse). Since external firing rate ${\nu}_{\text{X}}$ is stationary, the gating variable for external current is constant in time: ${S}_{\mathrm{X}}={C}_{\mathrm{X}}{\nu}_{\text{X}}$ . For recurrent currents, the temporal course of ${S}_{R}$ is dependent on the instantaneous presynaptic population activity ${\nu}_{{\alpha}_{R}}(t)$. Consequently, the total synaptic input current ${I}_{\mathrm{s}\mathrm{y}\mathrm{n}}(t)$, given as a sum of contributions from external and recurrent components
depends on the population firing rates ${\nu}_{\mathrm{E}}(t)$ and ${\nu}_{\mathrm{I}}(t)$. The output firing rate of population neurons, in turn, is determined by the input current and can be modeled in terms of an inputoutput response function $F$.
In general, the inputoutput relationship $v(t)=F\left({I}_{\mathrm{s}\mathrm{y}\mathrm{n}}(t)\right)$ depends on the spectral characteristics of the input current, resulting in frequency dependent phase shifts and/or amplitude modulations between the oscillatory components of ${I}_{\mathrm{s}\mathrm{y}\mathrm{n}}$ and $v$. However, it has been shown (Brunel et al., 2001; Fourcaud and Brunel, 2002) that the output rate in the leaky integrate and fire neuron model follows instantaneously the temporal variations in its synaptic input current given that synaptic noise is sufficiently strong and synaptic time constant is comparable with membrane time constant. That is, in these conditions, the response does not exhibit a phase shift, and its amplitude is independent of the frequency of oscillatory components of the input current. As a result, even if the input current is varying in time, the inputoutput function $F$ can be approximated by the currentfrequency response function $\varphi $, given by Equation 34, describing the output due to the steady input current.
In the framework of mean field approximation, the output rates ${\varphi}_{\text{E}}\left({I}_{\mathrm{s}\mathrm{y}\mathrm{n},\mathrm{E}}(t)\right)$ and ${\varphi}_{\text{I}}\left({I}_{\mathrm{s}\mathrm{y}\mathrm{n},\mathrm{I}}(t)\right)$ for excitatory and inhibitory populations must be the same as the instantaneous presynaptic population rates ${\nu}_{\mathrm{E}}(t)$ and ${\nu}_{\mathrm{I}}(t)$ because both presynaptic and output rates are of the same populations. This requirement results in two selfconsistent equations:
Since the amplitudes of firing rate deviations from the rates in asynchronous steady state are small, $\varphi \left({I}_{\mathrm{s}\mathrm{y}\mathrm{n}}(t)\right)$ can be linearized about the input current ${I}_{\mathrm{s}\mathrm{y}\mathrm{n}}^{0}$ in asynchronous state as:
With this approximation, the selfconsistent equations for excitatory and inhibitory populations become
where ${A}_{\alpha}=\frac{{I}_{\mathrm{s}\mathrm{y}\mathrm{n},\alpha}^{0}}{{v}_{\alpha}^{0}}\frac{d{\varphi}_{\alpha}\left({I}_{\mathrm{s}\mathrm{y}\mathrm{n},\alpha}^{0}\right)}{d{I}_{\mathrm{s}\mathrm{y}\mathrm{n},\alpha}}$ is the dimensionless slope of the currentfrequency response function at the current value in asynchronous state, expressed as the ratio between the relative changes in the firing rate and the input current (Brunel and Wang, 2003).
The selfconsistent equations Equation 45 together with Equation 41 for the gating variables and Equation 42 for the total synaptic current describe approximate firing rate dynamics of excitatory and inhibitory populations. To determine if the network develops oscillatory instability caused by small fluctuations in population firing rates, we seek solutions for the rates ${\nu}_{\mathrm{E}}(t)$ and ${\nu}_{\mathrm{I}}(t)$ in which initially small (with relative amplitudes $\left{\epsilon}_{\mathrm{E}}\right\ll 1$ and $\left{\epsilon}_{\mathrm{I}}\right\ll 1$) oscillatory perturbations that can change exponentially with time are added to the stationary rates ${v}_{\mathrm{E}}^{0}$ and ${v}_{\mathrm{I}}^{0}$ such that: ${\nu}_{\alpha}(t)={v}_{\alpha}^{0}\left(1+\left{\epsilon}_{\alpha}\right\mathrm{exp}(\lambda t)\mathrm{c}\mathrm{o}\mathrm{s}(\omega t+{\phi}_{\alpha})\right)$ or, equivalently, in complex form
where $\lambda $ is the rate of perturbation growth, $\omega $ is the oscillation frequency, and ${\epsilon}_{\alpha}$ is complex accounting for a possible shift in oscillation phase ${\phi}_{\alpha}$ between the two populations. We can now replace the firing rates in Equation 41 with these expressions to solve the two equations and determine the synaptic variables ${S}_{R}(t)$ for recurrent currents mediated by $R=\text{AMPA, NMDA, GABA}$ receptors:
where
and
The components of synaptic currents and the total currents ${I}_{\mathrm{s}\mathrm{y}\mathrm{n},\mathrm{E}}(t)$ and ${I}_{\mathrm{s}\mathrm{y}\mathrm{n},\mathrm{I}}(t)$ can now be calculated and inserted into the linearized selfconsistent Equation 45 for population firing rates. Taking into account that the balance between the components of synaptic currents in excitatory and inhibitory populations is equal, we arrive at the following set of two equations
and the relationship between the relative amplitudes:
where
Solving Equation 50, we obtain the rate of perturbation growth $\lambda $ and the oscillation frequency $\omega $. Because both ${A}_{\mathrm{E}}$ and ${A}_{\mathrm{I}}$ are real, the linear relationship between the amplitudes ${\epsilon}_{\mathrm{E}}$ and ${\epsilon}_{\mathrm{I}}$ given by Equation 51 means that there is no phase lag between firing rates of excitatory and inhibitory populations.
Analytical consideration of the dependence of oscillation growth rate on network parameters
To further elucidate how specifically 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}}$ affect synchrony, we linearize the mean field equations Equation 37 and equations Equation 50 for the stability analysis around the point $\{{g}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A},\left\{\mathrm{E},\mathrm{I}\right\}}^{\ast},\text{}{g}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A},\left\{\mathrm{E},\mathrm{I}\right\}}^{\ast},\text{}{g}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A},\left\{\mathrm{E},\mathrm{I}\right\}}^{\ast},\text{}{v}_{\mathrm{X}}^{\ast}\text{}\}$ corresponding to the critical state network where $\lambda =0$. We then derive an analytical approximation for the oscillation growth rate $\lambda $ describing its dependence on the synaptic conductances and external rate in the vicinity of this point.
Linearization of mean field equations
Request a detailed protocolApproximate analytic description of the changes in the population firing rates $\u2206{v}_{\mathrm{E}}$ and $\u2206{v}_{\mathrm{I}}$ due to small changes in the synaptic conductances and external rate can be obtained by linearizing the currentfrequency response function $\varphi $, providing population firing rates ${v}_{\mathrm{E}}$ and ${v}_{\mathrm{I}}$ as a function of synaptic conductances and external rate. We note that the function $\varphi $ (Equations 34–36) explicitly depends on the mean effective synaptic input μ, synaptic noise $\sigma $, membrane time constant $\tau $, and synaptic time constant ${\tau}_{\mathrm{s}\mathrm{y}\mathrm{n}}$, which in turn depend on the synaptic conductances and external rate (Equations 22–33). Thus, changes in the firing rates $\u2206{v}_{\mathrm{E}}$ and $\u2206{v}_{\mathrm{I}}$ in response to small changes in the synaptic conductances $\u2206{g}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}$, $\u2206{g}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}$, $\u2206{g}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}$, and external rate $\u2206{v}_{\mathrm{X}}$ can be approximated as:
The dominant contribution to $\u2206{v}_{\alpha}$ is due to the change in synaptic input, $\mathrm{\Delta}{\mu}_{\alpha}$. Contributions from the remaining terms are relatively small, with the largest contribution being due to the change in the effective membrane time constant, $\mathrm{\Delta}{\tau}_{\alpha}$. Therefore, the expression for $\u2206{v}_{\alpha}$ can be simplified by retaining only the terms involving $\mathrm{\Delta}{\mu}_{\alpha}$ and $\mathrm{\Delta}{\tau}_{\alpha}$:
$\mathrm{\Delta}{\mu}_{\alpha}$ and $\mathrm{\Delta}{\tau}_{\alpha}$ are expressed through the relative changes in synaptic conductances $\u2206{g}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}/{g}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}^{*}$, $\u2206{g}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}/{g}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}^{*}$, $\u2206{g}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}/{g}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}^{*}$, external rate $\u2206{v}_{\mathrm{X}}/{v}_{\mathrm{X}}^{*}$, and the changes in population rates $\u2206{v}_{\mathrm{E}}$ and $\u2206{v}_{\mathrm{I}}$:
where $R=\left\{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A},\phantom{\rule{thinmathspace}{0ex}}\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A},\phantom{\rule{thinmathspace}{0ex}}\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}\right\}$, and
Inserting expressions for $\mathrm{\Delta}{\mu}_{\alpha}$ and $\mathrm{\Delta}{\tau}_{\alpha}$ into Equation 56, we obtain a closed system of linear equations for the changes in the firing rates of excitatory and inhibitory populations in response to small changes in the synaptic conductances and external rates. In matrix form these equations can be written as
where
Here, the elements of matrices $\mathbf{a}$ and $\mathbf{b}$ are constants defined by the point in the network parameter space around which the mean field equations are linearized. Components of the vector $\u2206\mathbf{v}$ are the changes in the firing rates of excitatory and inhibitory populations due to the changes in the synaptic conductances and external rate given by the components of vector $\u2206\mathbf{p}$. Taking into account that ${\varphi}_{\tau ,\alpha}^{`}{b}_{R,\alpha}^{\tau}\ll {\varphi}_{\mu ,\alpha}^{`}{b}_{R,\alpha}^{\mu}$ and that ${\mu}_{\alpha}+{V}_{\text{L}}\approx \u2329{V}_{\alpha}\u232a$, we neglect the ${\varphi}_{\tau ,\alpha}^{`}{b}_{R,\alpha}^{\tau}$ terms in $\mathbf{b}$ and replace ${\mu}_{\alpha}+{V}_{\text{L}}$ with $\u2329{V}_{\alpha}\u232a$. With these approximations $\mathbf{b}$ simplifies to:
Equation Equation 71 can now be rewritten as
where $\mathbf{I}$ is the identity matrix. Solving this equation for $\u2206\mathbf{v}$ we obtain
or in component form
where $\mathbf{W}={\left[\begin{array}{cc}{W}_{\mathrm{E}}& {W}_{\mathrm{I}}\end{array}\right]}^{T}$ is given by
In summary, equations Equations 77, 78 describe changes in the excitatory, $\u2206{v}_{\mathrm{E}}$, and inhibitory, $\u2206{v}_{\mathrm{I}}$, population firing rates due to the small relative changes in the synaptic conductances $\u2206{g}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}/{g}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}^{*}$, $\u2206{g}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}/{g}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}^{*}$, $\u2206{g}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}/{g}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}^{*}$, and external rate $\u2206{v}_{\mathrm{X}}/{v}_{\mathrm{X}}^{*}$.
Linearization of equations for oscillatory instability analysis
Request a detailed protocolChanges in synaptic parameters result not only in the changes of population firing rates, but also affect the stability of population dynamics. To understand the precise role played by the synaptic conductances and external input in the destabilization of the steady dynamics and emergence of network oscillation near the boundary between asynchronous and synchronous states, we derive an approximate analytic description of the change in the rate of oscillatory instability growth $\u2206\lambda $ and the change in the oscillation frequency $\u2206\omega $ caused by small changes in the synaptic conductances and external rate. For this purpose, we linearize equations Equation 50 for $\lambda $ and $\omega $ around the point $\{{g}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A},\left\{\mathrm{E},\mathrm{I}\right\}}^{\ast},\phantom{\rule{thinmathspace}{0ex}}{g}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A},\left\{\mathrm{E},\mathrm{I}\right\}}^{\ast},\phantom{\rule{thinmathspace}{0ex}}{g}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A},\left\{\mathrm{E},\mathrm{I}\right\}}^{\ast},\phantom{\rule{thinmathspace}{0ex}}{v}_{\mathrm{X}}^{\ast}\}$ corresponding to the critical state network that is on the boundary between steady and oscillatory states where $\lambda =0$. We do this by taking the differentials with respect to the synaptic variables ${\mathrm{\Phi}}_{R}$ and ${X}_{R}$, ($R=\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A},\phantom{\rule{thinmathspace}{0ex}}\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A},\phantom{\rule{thinmathspace}{0ex}}\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}$) that, in turn, depend on $\lambda $ and $\omega $:
The parameter ${X}_{R}$ (see Equations 52–54) characterizes the relative attenuation in the strength of the underlying synapse due to the $R$current dynamics. In addition to the dependency on $\lambda $ and $\omega $ through ${Q}_{R}$ (Equation 48), ${X}_{R}$ depends directly on its corresponding synaptic conductance ${g}_{R}$ and indirectly on all the synaptic conductances and external rate through its dependency on the slope ${\varphi}_{{I}_{\mathrm{s}\mathrm{y}\mathrm{n}},{\alpha}_{R}}^{`}$ of the currentfrequency response function. The change $\mathrm{\Delta}{X}_{R}$ due to small variations in the synaptic conductances and external rate is given by
where ${\varphi}_{{I}_{\mathrm{s}\mathrm{y}\mathrm{n}},{\alpha}_{R}}^{`}$, ${X}_{R}$, ${Q}_{R}$, and ${\mathrm{\Phi}}_{R}$ are constants whose values are defined by the point $\{{g}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A},\left\{\mathrm{E},\mathrm{I}\right\}}^{\ast},\phantom{\rule{thinmathspace}{0ex}}{g}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A},\left\{\mathrm{E},\mathrm{I}\right\}}^{\ast},\text{}{g}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A},\left\{\mathrm{E},\mathrm{I}\right\}}^{\ast},\phantom{\rule{thinmathspace}{0ex}}{v}_{\mathrm{X}}^{\ast}\}$ in the synaptic parameter space around which the stability analysis equations are linearized.
The relative change ${\mathrm{\Delta}Q}_{R}/{Q}_{R}$ can be obtained from Equation 48:
and the change in ${\mathrm{\Phi}}_{R}$ from Equation 49:
where
and $\omega $ is the oscillation frequency at the critical state. Inserting expressions for ${\mathrm{\Delta}Q}_{R}/{Q}_{R}$ into equations Equation 81 for $\mathrm{\Delta}{X}_{R}$ and, subsequently, expressions for ${\mathrm{\Delta}\mathrm{\Phi}}_{R}$ and $\mathrm{\Delta}{X}_{R}$ into equations Equation 80, we obtain a system of two linear equations for $\mathrm{\Delta}\lambda $ and $\mathrm{\Delta}\mathrm{\omega}$:
where
and $R=\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A},\phantom{\rule{thinmathspace}{0ex}}\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A},\phantom{\rule{thinmathspace}{0ex}}\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}$. Solving the system of equations Equation 86 for $\mathrm{\Delta}\lambda $ and $\mathrm{\Delta}\mathrm{\omega}$ we obtain:
Here, ${\mathrm{\Lambda}}_{R}$ and ${\mathrm{\Omega}}_{R}$ are constants defined by the parameters of the critical state network around which the stability analysis equations are linearized:
or, equivalently,
where ${T}_{0}=\sqrt{{T}_{+}^{2}+{T}_{}^{2}}$ and ${\mathrm{\Phi}}_{0}=\mathrm{a}\mathrm{t}\mathrm{a}\mathrm{n}\left({T}_{+}/{T}_{}\right)$.
Note that while $\u2206\lambda $ and $\u2206\omega $ given by equations Equations 93 and 94 depend directly on the changes in the synaptic conductances, they also depend indirectly on these parameters and the change in external rate through the terms involving ${\mathrm{\Delta}\varphi}_{{I}_{\mathrm{s}\mathrm{y}\mathrm{n},\mathrm{E}}}^{`}$ and ${\mathrm{\Delta}\varphi}_{{I}_{\mathrm{s}\mathrm{y}\mathrm{n},\mathrm{I}}}^{`}$ characterizing changes in the slopes of the currentfrequency response functions of excitatory and inhibitory neurons. To calculate these changes due to the changes in the synaptic conductances and external rate, we note that
Taking into account the linear relationship $\mu ~{I}_{\mathrm{s}\mathrm{y}\mathrm{n}}/{g}_{\mathrm{m}}S$ between the effective synaptic input μ and total synaptic current ${I}_{\mathrm{s}\mathrm{y}\mathrm{n}}$, we arrive at
As in the case of the change in the currentfrequency response function ${\mathrm{\Delta}\varphi}_{\mu ,\alpha}$, the dominant contribution to the change in the slope of the response function ${\mathrm{\Delta}\varphi}_{\mu ,\alpha}^{`}$ is coming from the change in the synaptic input $\u2206{\mu}_{\alpha}$, while the change in the effective membrane time constant $\u2206{\tau}_{\alpha}$, similarly, is the next largest contribution. Therefore, ${\mathrm{\Delta}\varphi}_{\mu ,\alpha}^{`}$ can be approximated as
Note also that using Equation 22 one can express the relative change $\mathrm{\Delta}{S}_{\alpha}/{S}_{\alpha}$ through the change $\u2206{\tau}_{\alpha}$ as $\mathrm{\Delta}{S}_{\alpha}/{S}_{\alpha}=\mathrm{\Delta}{\tau}_{\alpha}/{\tau}_{\alpha}$. Inserting expressions for ${\mathrm{\Delta}\varphi}_{\mu ,\alpha}^{`}$ and $\mathrm{\Delta}{S}_{\alpha}/{S}_{\alpha}$ into Equation 100 we obtain
Equations for $\u2206{\mu}_{\alpha}$ and $\u2206{\tau}_{\alpha}$ in terms of the changes in the synaptic conductances, external rate, and the resulting changes in the population firing rates $\u2206{v}_{\mathrm{E}}$ and $\u2206{v}_{\mathrm{I}}$ have been already derived and are given by Equations 57 and 58. We replace $\u2206{v}_{\mathrm{E}}$ and $\u2206{v}_{\mathrm{I}}$ in these equations with the solution obtained from the linearization of the mean field equations given, respectively, by Equations 77, 78. Next, by inserting the resulting $\u2206{\mu}_{\alpha}$ and $\u2206{\tau}_{\alpha}$ into Equation 102, we obtain expressions describing the relative changes in the slopes of the response functions for excitatory and inhibitory neurons due to the small changes in the synaptic conductances and external rate. In matrix form, these expressions can be written as
where
The elements of matrices $\stackrel{~}{\mathbf{a}}$ and $\stackrel{~}{\mathbf{b}}$ are constants defined by the parameters of the critical state network. Noting that $\left(\frac{{\varphi}_{\tau \mu ,\alpha}^{``}}{{\varphi}_{\mu ,\alpha}^{`}}+\frac{1}{{\tau}_{\alpha}}\right){b}_{R,\alpha}^{\tau}\ll \frac{{\varphi}_{\mu \mu ,\alpha}^{``}}{{\varphi}_{\mu ,\alpha}^{`}}{b}_{R,\alpha}^{\mu}$, we neglect the $\left(\frac{{\varphi}_{\tau \mu ,\alpha}^{``}}{{\varphi}_{\mu ,\alpha}^{`}}+\frac{1}{{\tau}_{\alpha}}\right){b}_{R,\alpha}^{\tau}$ terms in $\stackrel{~}{\mathbf{b}}$ and, as in the calculation of the change in the currentfrequency response function, replace ${\mu}_{\alpha}+{V}_{\text{L}}$ with $\u2329{V}_{\alpha}\u232a$. With these approximations $\stackrel{~}{\mathbf{b}}$ simplifies to:
Equation Equation 103 can now be written as
or in component form
where $\mathbf{U}={\left[\begin{array}{cc}{U}_{\mathrm{E}}& {U}_{\mathrm{I}}\end{array}\right]}^{T}$ is given by
Equations Equations 108 and 109 describing the relative changes in the slopes of the response functions for excitatory and inhibitory neurons can now be combined with Equations 93 and 94 to account for both direct and indirect dependence of the change in the oscillation growth rate $\mathrm{\Delta}\lambda $ and change in the oscillation frequency $\mathrm{\Delta}\omega $ on the small relative changes in the synaptic conductances and external rate.
Numerical solutions
Request a detailed protocolSelfconsistent mean field equations for the eight conductance parameters, and linear stability equations for the perturbation growth rate $\lambda $ and the oscillation frequency $\omega $ were both solved numerically using custom codes written in MATLAB (The MathWorks) with the aid of fsolve function.
Appendix 1
Accuracy of analytical approximation for the instability growth rate
Appendix 1—figure 1shows the comparison of the exact value of $\lambda$ obtained by numerically solving the mean field and stability analysis equations, with the analytical approximation given by Equations 1, 2. The comparison is performed separately for each of the three synaptic conductances and external rate by varying one of the underlying parameters while the remaining three keeping constant at their critical values. It is seen that in all cases, the analytical expression for $\lambda$ provides a very good approximation of the exact relationship.
Appendix 2
Contributions from $\mathbf{\Lambda}}_{\mathbf{A}\mathbf{M}\mathbf{P}\mathbf{A}$, $\mathbf{\Lambda}}_{\mathbf{N}\mathbf{M}\mathbf{D}\mathbf{A}$, and $\mathbf{\Lambda}}_{\mathbf{G}\mathbf{A}\mathbf{B}\mathbf{A}$ terms to the instability growth rate
The reasons of why ${\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}}$ can be understood by considering analytical expressions of 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}}$. These quantities have a form ${\mathrm{\Lambda}}_{R}={X}_{R}\mathrm{c}\mathrm{o}\mathrm{s}({\mathrm{\Phi}}_{R}+{\mathrm{\Phi}}_{0})/{T}_{0}$, where ${X}_{R}$ characterizes the relative attenuation in the strength of the corresponding synapse due to the $R$current dynamics, ${\mathrm{\Phi}}_{R}$ is the phase lag introduced by $R$current synaptic filtering, and ${T}_{0}$ and ${\mathrm{\Phi}}_{0}$ are constants (see Materials and methods). For illustration purposes, it is convenient to use a geometric interpretation of ${\mathrm{\Lambda}}_{R}$ in which it is associated with a vector emanating from the origin of a cartesian plane, with the length ${X}_{R}/{T}_{0}$, and the angle ${\mathrm{\Phi}}_{R}+{\mathrm{\Phi}}_{0}$ between the vector and the $x$axis. Then, the projection of this vector on the $x$axis (i.e. its $x$component) is ${\mathrm{\Lambda}}_{R}$.
Appendix 2—figure 1 shows vectors $\mathit{L}}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}$, $\mathit{L}}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}$, and $\mathit{L}}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}$ corresponding to the values of $\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 the crttical network. It is seen that the length of the vector $\mathit{L}}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}$ is nearly zero and, therefore, its projection $\mathrm{\Lambda}}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}$ is nearly zero as well. The length of $\mathit{L}}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}$ being very short means that the relative attenuation in the strength of NMDA synapse ${X}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}\ll 1$. This is because the characteristic synaptic time constant ${\tau}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}\sim 100$ ms of NMDA synapse is much longer than the period of the $\gamma$band oscillation $1/{f}_{\mathrm{n}\mathrm{t}\mathrm{w}\mathrm{r}\mathrm{k}}\sim 15$ ms emerging in the critical network at the boundary between the steady and oscillatory states. This results in a strong attenuation of the oscillatory component of NMDAR mediated current that is quantified by the $X}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}$ numeric value (Equations 48, 53).
The vector ${\mathit{L}}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}$, in contrast, has a length comparable with the length of ${\mathit{L}}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}$. However, because its angle ${\mathrm{\Phi}}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}+{\mathrm{\Phi}}_{0}$ is close to 90°, the ${\mathit{L}}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}$ projection ${\mathrm{\Lambda}}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}$ is much smaller than the projection ${\mathrm{\Lambda}}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}$ of the vector ${\mathit{L}}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}$. This occurs because the phase shift ${\mathrm{\Phi}}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}$ introduced by GABAcurrent synaptic filtering is lagging the phase shift ${\mathrm{\Phi}}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}$ introduced by AMPAcurrent filtering. This difference in the phase delays, in turn, is due to the fact that the GABA synapse (${\tau}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}~5$ ms) is slower than the AMPA synapse (${\tau}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}~2$ ms). Thus, while both ${\mathrm{\Lambda}}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}$ and ${\mathrm{\Lambda}}_{\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}$ are much smaller than ${\mathrm{\Lambda}}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}$, the primary reasons are different.
Appendix 3
Direct and indirect pathways for synaptic conductances and external rate affecting synchrony
Equation 6 for the oscillatory instability growth rate $\lambda $ illuminates specific pathways by which modulations of synaptic conductances and external rate affect synchrony. In general, when fast excitatory feedback becomes sufficiently strong and is followed by slow and powerful inhibition the asynchronous state becomes unstable and synchronous oscillation emerges (Brunel and Wang, 2003; Compte et al., 2000; Tsodyks et al., 1997; Wang, 1999). The intuitive explanation is as follows. An increase of activity in the excitatory population caused by a fluctuation rapidly surges due to the fast excitatory feedback. The buildup of excess excitation, in turn, activates the inhibitory population. The growth of excitation continues until slow rising inhibition eventually suppresses excitatory population activity. As a result, activity of inhibitory population decays as well, and the next oscillation cycle starts. Equation 5 explicitly shows that the strong excitatory feedback necessary to induce oscillation can emerge via two pathways. One is manifested in the indirect dependence of $\lambda $ on the synaptic parameters via their effect on the operating point of the currentfrequency response function. Changes in the synaptic conductances and external rate producing an increase in the average total synaptic current ($\mathrm{\Delta}{I}_{\mathrm{s}\mathrm{y}\mathrm{n}}>0$) shift the operating point of excitatory population toward a steeper slope of the currentfrequency response function ($\mathrm{\Delta}{\varphi}_{{I}_{\mathrm{s}\mathrm{y}\mathrm{n}}}^{\prime}>0$). As a result, the amplitude of the neuron’s response to dynamically varying input increases, giving rise to stronger excitatory feedback. If AMPAR conductance is fixed, $\mathrm{\Delta}{I}_{\mathrm{s}\mathrm{y}\mathrm{n}}>0$ results in $\lambda >0$ (Equation 5). Therefore, changes in the synaptic parameters that increase the average total synaptic currents, in general, have a destabilizing effect on the asynchronous activity and favor synchronous oscillations.
It should be noted that the change in the slope with a shift in the operating point is critical for the indirect mechanism functionality. If the response function were linear and, thus, the slope is constant, independent of ${I}_{\mathrm{s}\mathrm{y}\mathrm{n}}$ (${U}_{\mathrm{E}}=0$, Equation 5), a change in the operating point $\u2206{I}_{\mathrm{s}\mathrm{y}\mathrm{n}}$ will still affect the neuron’s average firing frequency. However, the neuron’s response amplitude to dynamically varying input will be independent of the operating point, i.e., the average input current ${I}_{\mathrm{s}\mathrm{y}\mathrm{n}}$ . As a result, shifts in the operating point will not change the strength of excitatory feedback. In the case of integrateandfire neurons, the input currentoutput frequency relationship is a monotonically increasing nonlinear function (${U}_{\mathrm{E}}>0$). Consequently, when the operating point shifts due to increasing average input current ($\mathrm{\Delta}{I}_{\mathrm{s}\mathrm{y}\mathrm{n}}>0$), the amplitude of the neuron’s response to the same variation in the input current increases since the slope of the response function becomes steeper, giving rise to stronger excitatory feedback.
The second pathway is manifested in the direct dependence of $\lambda $ on the AMPAR conductance. An increase in the conductance results in stronger excitatory to excitatory connections and excitatory to inhibitory connections that, in turn, amplify fast AMPAR mediated excitatory feedback and slow GABAR mediated inhibition. If other parameters remain unchanged, $\mathrm{\Delta}{g}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}>0$ also results in $\mathrm{\Delta}{I}_{\mathrm{s}\mathrm{y}\mathrm{n}}>0$ (Equation 4). This means that because of the AMPAR conductance increase, both direct and indirect pathways tend to make $\lambda >0$ (Equation 5). Thus, increased fast excitation, due to the stronger AMPAR mediated current, and enhanced response to varying input, due to the shift in the operating point, destabilize the asynchronous dynamic and push the network toward synchronous oscillations.
Appendix 4
The mechanism by which NMDAR influences network synchronization
In this section we provide a more formal consideration of the mechanism for the transition between the steady and oscillatory states, and the lack thereof when the NMDAR is blocked. For convenience, we rewrite here Equation 7 for the instability growth rate $\lambda $:
The expression in the brackets is proportional to the change $\u2206{I}_{\mathrm{s}\mathrm{y}\mathrm{n}}$ in the average total synaptic current (Equation 4). Without loss of generality, we assume that the drugnaive condition corresponds to the network for which $\u2206{g}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}=0$ and $\mathrm{\Delta}{v}_{\mathrm{X}}=\delta <0$. Thus, initially $\mathrm{\Delta}{I}_{\mathrm{s}\mathrm{y}\mathrm{n}}<0$ and the response function operating point is below the point of the critical network. As a result, $\lambda <0$ and the network is in the asynchronous steady state. When external input increases by $2\delta $, so that $\mathrm{\Delta}{v}_{\mathrm{X}}=+\delta >0$, the total synaptic current increases as well, so that $\u2206{I}_{\mathrm{s}\mathrm{y}\mathrm{n}}$ becomes positive and the operating point shifts above the point of the critical network. This results in $\lambda >0$ and, therefore, the network transitions to the synchronous oscillatory state. This transition is due to the indirect dependence of $\lambda $ on the external rate change $\u2206{v}_{\mathrm{X}}$ via its effect on the operating point of the currentfrequency response function. In the drug condition, the NMDAR conductance is reduced and, if it becomes sufficiently small so that $\frac{\mathrm{\Delta}{g}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}}{{g}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}^{\ast}}<\frac{\delta}{{v}_{\mathrm{X}}^{\ast}}\frac{{I}_{\mathrm{X}}}{{I}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}}$, the same increase in external rate will not result in positive $\lambda $. As a result, no oscillatory instability develops, and the network remains in the asynchronous steady state. This is because the decrease in the NMDAR conductance, which like $\u2206{v}_{\mathrm{X}}$ influences $\lambda $ indirectly, reduces the average total synaptic current ${I}_{\mathrm{s}\mathrm{y}\mathrm{n}}$ and brings down the operating point to such locus that it cannot be shifted above the operating point of the critical network by the same $2\delta $ increase in the external rate.
When the oscillatory instability growth rate $\lambda $ on the lefthand side of Equation 111 is fixed at zero, it becomes the equation for the critical hyperplane in the synaptic parameter space separating the asynchronous and synchronous states. In the context of our model, we are interested in the network states within the $\left({v}_{\mathrm{X}}/{v}_{\mathrm{X}}^{\ast},\phantom{\rule{thinmathspace}{0ex}}{g}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}/{g}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}^{\ast}\right)$ subspace and, thus, Equation 111 becomes the equation for the line
This critical line is the analytical approximation of the exact critical line shown in the state diagram depicted in Figure 6C that was obtained by solving the mean field and stability analysis equations numerically. Appendix 4—figure 1 shows the comparison between the exact critical line (dotted line) and its approximation given by Equation 112 (solid line). It is seen that these lines virtually overlap and, thus, the accuracy of approximation is very good.
Appendix 5
Special case of the zero net effect of the variation of synaptic parameters on the average total synaptic current
Here we consider an interesting case (see a specific example below) when variations in the synaptic conductances and external rate are such that they have no net effect on the average total synaptic current because the sum of their contributions cancels out. In this case, $\u2206{I}_{\mathrm{s}\mathrm{y}\mathrm{n}}=0$ and there will be no shifts in the operating points and, therefore, no change in the firing rates and slopes of the response functions (Equation 3). As a result, despite the synaptic parameter changes there will not be any indirect contribution to $\lambda $ (Equation 1). However, even though the background firing rates remain unchanged, the oscillatory instability growth rate $\lambda $ will be affected via its direct dependence on the changes in the synaptic conductances.
Such case was studied in Compte et al., 2000. Specifically, Compte and colleagues observed in simulations that the concurrent increase in the AMPAR and decrease in the NMDAR conductances, so that the total strength of excitation remains fixed, leads to the destabilization of the asynchronous steady state and emergence of synchronous oscillations. Assuming that the network is close to the boundary between the asynchronous and synchronous states, this phenomenon can be explained in terms of Equation 6 for the oscillatory instability growth rate $\lambda $. In this case, it simplifies to
The expression in the brackets is proportional to the change $\u2206{I}_{\mathrm{s}\mathrm{y}\mathrm{n}}$ in the average total synaptic current (Equation 4). It can be seen that decreasing the NMDAR conductance by $\mathrm{\Delta}{g}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}^{0}<0$ and simultaneously increasing the AMPAR conductance in the amount of $\mathrm{\Delta}{g}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}^{0}=\mathrm{\Delta}{g}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}^{0}\frac{{I}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}}{{I}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}}\frac{{g}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}^{\ast}}{{g}_{\mathrm{N}\mathrm{M}\mathrm{D}\mathrm{A}}^{\ast}}>0$ results in zero for the expression in the brackets. This means that the contributions from the changes in the NMDAR and AMPAR mediated currents to the average total synaptic current ${I}_{\mathrm{s}\mathrm{y}\mathrm{n}}$ cancel each other and $\u2206{I}_{\mathrm{s}\mathrm{y}\mathrm{n}}=0$. Thus, in this case, although the AMPAR and NMDAR conductances change, the operating point of the response function does not. However, because $\lambda $ also depends directly on the AMPAR conductance (the first term in the r.h.s. of Equation 113) and $\mathrm{\Delta}{g}_{\mathrm{A}\mathrm{M}\mathrm{P}\mathrm{A}}^{0}>0$, the instability growth rate $\lambda $ becomes positive and synchronous oscillations emerge due to the stronger excitatoryinhibitory feedback loop.
Data availability
The current manuscript is a computational study. Neural data analyzed in this paper are available at the Mendeley Data repository. Custom MATLAB codes used for simulations, mean field and stability analysis are available at the GitHub: https://github.com/amirikian/SCZSynapticCircuitFailureModel (copy archived at Amirikian, 2024).

Mendeley DataNeural data accompanying Crowe et al. 'A prefrontal network operating near steady and oscillatory states links spike desynchronization and synaptic deficits in schizophrenia'.https://doi.org/10.17632/bt6j9gyz5t.1
References

Asynchronous states in networks of pulsecoupled oscillatorsPhysical Review. E, Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics 48:1483–1490.https://doi.org/10.1103/physreve.48.1483

BookCorticonics: Neural Circuits of the Cerebral CortexCambridge: Cambridge University Press.

SoftwareSCZsynapticcircuitfailuremodel, version swh:1:rev:bd7a267d8852fa7cea26358775ae5908fbce17d0Software Heritage.

BookModeling Brain Function: The World of Attractor Neural NetworksCambridge University Press.https://doi.org/10.1017/CBO9780511623257

Pairwise analysis can account for network structures arising from spiketiming dependent plasticityPLOS Computational Biology 9:e1002906.https://doi.org/10.1371/journal.pcbi.1002906

Contextprocessing deficits in schizophrenia: diagnostic specificity, 4week course, and relationships to clinical symptomsJournal of Abnormal Psychology 112:132–143.

Controlling the false discovery rate: a practical and powerful approach to multiple testingJournal of the Royal Statistical Society 57:289–300.https://doi.org/10.1111/j.25176161.1995.tb02031.x

BookCortex: Statistics and Geometry of Neuronal ConnectivityBerlin, Heidelberg: Springer.https://doi.org/10.1007/9783662037331

Firing frequency of leaky intergrateandfire neurons with synaptic current dynamicsJournal of Theoretical Biology 195:87–95.https://doi.org/10.1006/jtbi.1998.0782

Dynamics of sparsely connected networks of excitatory and inhibitory spiking neuronsJournal of Computational Neuroscience 8:183–208.https://doi.org/10.1023/a:1008925309027

Effects of synaptic noise and filtering on the frequency response of spiking neuronsPhysical Review Letters 86:2186–2189.https://doi.org/10.1103/PhysRevLett.86.2186

Effects of neuromodulation in a cortical network model of object working memory dominated by recurrent inhibitionJournal of Computational Neuroscience 11:63–85.https://doi.org/10.1023/a:1011204814320

What determines the frequency of fast network oscillations with irregular neural discharges? ISynaptic Dynamics and ExcitationInhibition Balance. J Neurophysiol 90:415–430.https://doi.org/10.1152/jn.01095.2002

BookTheoretical Neuroscience Computational and Mathematical Modeling of Neural SystemsCambridge, MA: MIT Press.

Dynamics of the firing probability of noisy integrateandfire neuronsNeural Computation 14:2057–2110.https://doi.org/10.1162/089976602320264015

Mnemonic coding of visual space in the monkey’s dorsolateral prefrontal cortexJournal of Neurophysiology 61:331–349.https://doi.org/10.1152/jn.1989.61.2.331

Decreased dendritic spine density on prefrontal cortical pyramidal neurons in schizophreniaArchives of General Psychiatry 57:65–73.https://doi.org/10.1001/archpsyc.57.1.65

On numerical simulations of integrateandfire neural networksNeural Computation 10:467–483.https://doi.org/10.1162/089976698300017845

Spiketiming dynamics of neuronal groupsCerebral Cortex 14:933–944.https://doi.org/10.1093/cercor/bhh053

Voltage dependence of NMDAactivated macroscopic conductances predicted by singlechannel kineticsThe Journal of Neuroscience 10:3178–3182.https://doi.org/10.1523/JNEUROSCI.100903178.1990

The dot pattern expectancy task: Reliability and replication of deficits in schizophreniaPsychological Assessment 22:131–141.https://doi.org/10.1037/a0017828

Computational study of NMDA conductance and cortical oscillations in schizophreniaFrontiers in Computational Neuroscience 8:133.https://doi.org/10.3389/fncom.2014.00133

BookBiophysics of Computation: Information Processing in Single NeuronsOxford University Press.

Cognitive control errors in nonhuman primates resembling those in schizophrenia reflect opposing effects of NMDA receptor blockade on causal interactions between cells and circuits in prefrontal and parietal corticesBiological Psychiatry. Cognitive Neuroscience and Neuroimaging 5:705–714.https://doi.org/10.1016/j.bpsc.2020.02.013

Dynamics of networks of excitatory and inhibitory neurons in response to timedependent inputsFrontiers in Computational Neuroscience 5:25.https://doi.org/10.3389/fncom.2011.00025

Formation and maintenance of neuronal assemblies through synaptic plasticityNature Communications 5:5319.https://doi.org/10.1038/ncomms6319

A dynamical systems hypothesis of schizophreniaPLOS Computational Biology 3:e228.https://doi.org/10.1371/journal.pcbi.0030228

Selective loss of smaller spines in schizophreniaThe American Journal of Psychiatry 174:586–594.https://doi.org/10.1176/appi.ajp.2017.16070814

Spiketimingdependent plasticity in balanced random networksNeural Computation 19:1437–1467.https://doi.org/10.1162/neco.2007.19.6.1437

BookMeanfield theory of recurrent cortical networks: working memory circuits with irregularly Spiking neurons inIn: Feng J, editors. Computational Neuroscience: A Comprehensive Approach. Boca Raton, FL: Chapman and Hall/ CRC Press. pp. 431–490.

Paradoxical effects of external modulation of inhibitory interneuronsThe Journal of Neuroscience 17:4382–4388.https://doi.org/10.1523/JNEUROSCI.171104382.1997

InvertedU dopamine D1 receptor actions on prefrontal neurons engaged in working memoryNature Neuroscience 10:376–384.https://doi.org/10.1038/nn1846

Synaptic basis of cortical persistent activity: the importance of NMDA receptors to working memoryThe Journal of Neuroscience 19:9587–9603.https://doi.org/10.1523/JNEUROSCI.192109587.1999
Decision letter

Albert CompteReviewing Editor; IDIBAPS, Spain

John R HuguenardSenior Editor; Stanford University School of Medicine, United States

Albert CompteReviewer; IDIBAPS, Spain

Yashar AhmadianReviewer; University of Cambridge, United Kingdom
Our editorial process produces two outputs: (i) public reviews designed to be posted alongside the preprint for the benefit of readers; (ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.
Decision letter after peer review:
Thank you for submitting your article "A prefrontal network model operating near steady and oscillatory states links spike desynchronization and synaptic deficits in schizophrenia" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, including Albert Compte as the Reviewing Editor and Reviewer #1, and the evaluation has been overseen by John Huguenard as the Senior Editor. The following individual involved in the review of your submission has agreed to reveal their identity: Yashar Ahmadian (Reviewer #2).
The manuscript is currently not suitable for publication in eLife, but we are willing to consider a revised version that would fully address all concerns of the reviewers. The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.
Essential revisions:
The reviewers consider that the manuscript has potential but it requires extensive revisions and new simulations to support its claims. Because these are substantial revisions, a revised manuscript will be critically reviewed in depth again and sent back out to the original reviewers.
The interest in the computational results depends critically on the ability of the model to provide a convincing explanation of the experimental data (Zick et al. 2018, 2021). For this, it is essential to:
1) consider not only the dynamics of synchronization (Figure 1) but also the dynamics of firing rate, and not only the period immediately preceding motor action but also the dynamics of the effects following motor action and around the stimulus cue. The mechanisms of the computational model make specific predictions as to how the dynamics of firing rate and oscillatory synchronization will be associated, and this should be specifically validated against the experimental data.
2) consider how the model accommodates a Ushape relationship between NMDAR modulation and network synchrony. Existing evidence shows that 0lag synchrony in prefrontal networks is affected by manipulations that reduce NMDAR function (Zick et al. 2018) and by manipulations that enhance NMDAR function (Zick et al. 2021). The computational model presented in this manuscript does not show this Ushape behavior and the discussion does not mention this.
In addition, the theoretical understanding of the model should be enhanced by:
3) expanding the theoretical analysis to provide more intuition on the way NMDA works the way it does, and clarification of its distinct role in synchronization (in contrast to AMPA).
4) investigating how synchrony in the model reacts to transient external inputs and how the model can produce weakly oscillatory synchrony.
Finally, all reviewers were concerned about the emphasis on the association with schizophrenia, so we recommend:
5) moving speculative links to schizophrenia to the Discussion section. Also, please broaden the scope of this discussion with available literature on how synchrony and STDP interact in computational models, and with alternative explanations based on changes in NMDARdependent synaptic plasticity irrespective of synchronization.
Reviewer #1 (Recommendations for the authors):
Here are some specific recommendations for improvement of the manuscript:
1) The last two paragraphs in the introduction are largely discussional in spirit. The paragraph before last does not make much sense here, at a point where the reader is not yet familiar with the details of the model proposed. The very last paragraph is also mostly speculative in how the results in this manuscript would be important for schizophrenia and would appear more natural in the Discussion section.
2) The data presented in Figure 1 could be presented more fully in order to provide more constraints for the model so as to enhance its biological plausibility. For instance, it would be most interesting to have a parallel plot of the firing rate in these same neurons to see how neuronal activity evolves through the task. Neurons show a response to the cue stimulus (Zick et al. 2018), which may be similar to the response at the time of the response, and that is something that the model would have to explain: how does the control network react with strong synchronization when it receives inputs prior to the response but instead does not show any synchronization when it receives inputs at the time of the cue? Also, it would be good to show panel b extended for a time window after motor response time t=0 to compare the time course of firing rate and synchrony increases. This could also constrain further the network model.
3) In Figure 1 it would be good to have a dashed horizontal line at y=0 to separate correlated (y>0) from anticorrelated (y<0) activity more clearly. I would also recommend changing the xaxis label ("time lag") because this is now used in crosscorrelation functions in Figure 4 and they are intrinsically different. The synchrony plot should also include error bars and proper inference about differences between the two conditions at different time points. For clarity, it would be nice to indicate with a shaded area or horizontal bar the period of cue stimulus presentation, too.
4) Briefly mention in line 167 what spiking neuron model you are using
5) Mention of "GABA receptors" in line 173 suggests that the model includes both GABA_A and GABA_B receptors.
6) Verb missing in "and mediated" in line 174.
7) In Figure 3 the color code red/blue means two different things: excitatory vs. inhibitory neurons, and steady vs. critical networks. This is confusing and should be clarified with different color codes for each thing.
Reviewer #2 (Recommendations for the authors):
1) I assume elucidation of the precise mechanism by which external inputs or NMDA affect synchrony can be achieved by analytically inspecting the direct and indirect dependence of the real part of the relevant complex eigenvalue on the synaptic conductances and external inputs, either directly or indirectly via the dependence of linearized neural gains in the background operating point (which itself depends on the external input and the synaptic parameters).
2) The authors do refer to the Compte et al. 2000 study that found that when the total strength of excitatory connectivity is fixed, tilting the relative strength of NMDAR and AMPAR towards the former disfavors synchronous oscillations (with relatively high frequency in or near the γ band). Previous studies (e.g. Wang and Brunel 2003, Tsodyks et al. 1997, and Compte et al. 2000) have pointed out the role of the negative feedback loop between excitatory and inhibitory populations mediated by the fast receptors AMPA and GABA in γband oscillations. Again I think the meanfield theory and the stability analysis based on it could be used to elucidate the differences and similarities in the pathways by which NMDAR and AMPAR affect the transition to synchrony. For example, the role of NMDAR could be via the explicit dependence of the growth rate of oscillations (the eigenvalue real part) on g_NMDA, or (most likely) via the effect of NMDARbased excitation on background rates (which indirectly affect the eigenvalue).
3) Regarding the link to "disconnection" via STDP. Various studies have found that STDP in (initially randomly connected) networks in the asynchronous state can maintain random connectivity (Morrison et al. 2007 – DOI: 10.1162/neco.2007.19.6.1437) or even lead to the emergence of structured connectivity (Izhikevich et al. 2004, Babadi and Abbott 2013, and LitwinKumar and Doiron 2014  DOIs: 10.1093/cercor/bhh053, 10.1371/journal.pcbi.1002906, 10.1038/ncomms6319). In lieu of evidence (either previous studies or new simulations/modeling by the authors) for this claim, I suggest weakening the claim to saying that the NMDARrelated reduction in synchrony, as a specific disruption in the patterns of precise spike timing, can lead to abnormal changes in patterns of synaptic connectivity via STDP. I also suggest citing the abovementioned (or similar) studies on previous studies of STDP in asynchronous (or synchronous) states, or any other studies that support specific claims by the authors in this regard.
Other suggestions:
– Figure 6 top: here the crosscorrelation plot based on monkey data (as in the model) exhibits nonzero lag "γ" peaks. But these seem to be absent in the data shown in the Neuron paper (e.g. Figure 4D therein). Can the authors comment on this seeming discrepancy? Is this because data only from one of the monkeys (showing stronger effects) was included in the current paper? Can the authors also comment on the flatness of the CCH in the "initial probe period" in the monkey data in contrast to the small 0lag peak in the model CCH?
– Lines 104105: I found the meaning of *balance between* in this sentence vague/unclear. What does the balance between NMDA inputs and oscillatory activity mean?
– I suggest the authors also mention (e.g. around lines 476477 where they discuss their assumption 1, and as further support for this assumption) that in their empirical data (Neuron paper) they do find an increase in PFC neural firing rates in the preresponse period, which will presumably also arise in their model from the increase in external inputs.
– Line 781: "themes" → "times".
– Line 838: I would explain where this equation comes from (I understand it accounts for the effect of reset on mean membrane potential), and also what v_α denotes.
Reviewer #3 (Recommendations for the authors):
Here are issues that I believe the authors should address before the paper can be published:
1. In Figure 1, the authors only show the synchrony before the motor response. I believe it would be instructive to show also what happens after motor response, to visualize how temporally precise is this increase. It would also be useful to show the time course of average firing rates. The explanation of the authors relies on an increase in external input to generate an increase in synchrony – this increase in external input should be reflected in an increase in the average firing rate as well.
2. Figure 6 shows a strong difference between the network model and data – side peaks in spike correlations are much more prominent in the model than in the data. Thus, synchronization in the data is much less oscillatory than in the model. Is there a way in the model to get such weakly oscillatory synchrony? If not, can the authors speculate as to what additional mechanisms would lead to a reduction of γ power, preserving zero lag synchrony?
3. In the data, the increase in synchrony is transient. In the simulations, the author considers the case of stationary external inputs. I believe the authors should also investigate the case of transient increases of external inputs leading to a more realistic time course of synchrony. Such a scenario could also help in solving
the issue mentioned in 2.
[Editors' note: further revisions were suggested prior to acceptance, as described below.]
Thank you for resubmitting your work entitled "A prefrontal network model operating near steady and oscillatory states links spike desynchronization and synaptic deficits in schizophrenia" for further consideration by eLife. Your revised article has been evaluated by John Huguenard (Senior Editor) and a Reviewing Editor.
The manuscript has been improved but there are some remaining issues that need to be addressed, as outlined below:
1. The dynamics of average firing rates in Figure 1 should mark clearly the various relevant periods of the task on the time axis: probe interval, motor response, etc.
2. The new theoretical analyses in this revision have greatly extended the length of the main text, but the additional insight provided is not clearly explained and in its current form it can lead to confusion. It is suggested that most of this material is moved to a supplementary document and a more concise and clearly delineated section in the main text provides the fundamental message of how NMDARs influence network synchronization.
Reviewer #3 (Recommendations for the authors):
The paper has been considerably improved and I would like to congratulate the authors for their hard work. There are still a couple of issues though that I think should be addressed to improve the paper further.
1. In response to one of the main concerns of the reviewers, the authors have expanded Figure 1 to include the dynamics of average firing rates in the data and the effect of NMDAR blockade on this dynamic. This is a welcome addition, but I think this figure could be improved to better understand the time course of the task. Currently, there is no discussion at all about when the motor response is happening and so the left and right panels seem totally disconnected. One has to go back to the original paper by Zick et al. to understand when the motor response happens. I would suggest extending the left panels to include the whole trial (including the `probe' interval), adding where the mean response time is similar to Figure 4 in the Zick et al. 2018 paper. Or at the very least indicate when is the mean response time (time 0 in right panels) with respect to the time axis of the left panels.
2. In response to another concern, the authors have considerably expanded their theoretical analysis with a detailed description of an approximation for the instability growth rate. However, while the approximation itself represents an interesting new derivation, I did not find that it brings much insight, and in fact in some ways introduces additional confusion, for the following reasons. (i) The main new result of this new analysis is the new Equation (1). In this equation, the LambdaRs play critical roles and are the subject of much of the later discussion, however, their meaning is not really explained in the main text and one has to go to the Methods to understand what they really are. (ii) The new Figure (10) indicates that the growth rate is primarily determined by the AMPA term and that the GABA term plays little role. This is quite confusing as oscillations in such networks relies on GABAR mediated negative feedback. How can the authors explain this fact? (iii) After reading the new part, I ended up not having much more insight about mechanisms. My guess at this point is that NMDAR influences the synchronisation properties of the network by providing an additional excitatory drive, thereby playing a similar role as the external excitatory input. I have the feeling that is what the authors are trying to say when they say that the NMDAR term acts through changes in phiprime (the slope of the transfer function), which is exactly how the external inputs act, but it feel like it could be explained in a much clearer way.
https://doi.org/10.7554/eLife.79352.sa1Author response
Essential revisions:
The reviewers consider that the manuscript has potential but it requires extensive revisions and new simulations to support its claims. Because these are substantial revisions, a revised manuscript will be critically reviewed in depth again and sent back out to the original reviewers.
The interest in the computational results depends critically on the ability of the model to provide a convincing explanation of the experimental data (Zick et al. 2018, 2021). For this, it is essential to:
1) consider not only the dynamics of synchronization (Figure 1) but also the dynamics of firing rate, and not only the period immediately preceding motor action but also the dynamics of the effects following motor action and around the stimulus cue. The mechanisms of the computational model make specific predictions as to how the dynamics of firing rate and oscillatory synchronization will be associated, and this should be specifically validated against the experimental data.
We have extended our analysis of both the simulated and recorded neural data to include modulation of both synchrony and firing rate during the DPX task (in the neural data) and in response to transient inputs to the PFC network (in the simulated data), both in the baseline case, and following reduction of NMDAR synaptic transmission. We added panels C and D to Figure 1 to illustrate modulation of firing rate in the neural data and the influence of NMDAR blockade on neural firing rates. This new analysis shows (1) firing rate in monkey prefrontal cortex ramped gradually up during the delay period considerably before the motor response, (2) spike synchrony increased more abruptly and dramatically immediately before the motor response, and (3) blocking NMDAR attenuated both the increase in spike synchrony and the increase in firing rate (new Figure 1C, D). In addition, we have extended the time of analysis of both synchrony and firing rate to encompass the response as well as the time period following the response. To facilitate the comparison of the dynamics of firing rate and synchrony between experimental data and the computational model we modified the location of the critical network in the parameter space by changing the balance (ratio) between NMDAR and GABAR mediated currents from 0.2 to 0.15. We, correspondingly, repeated all previous simulations for the new value of the underlying parameter and updated the manuscript figures accordingly. These new figures are very similar to the previous ones. Next, we extended our consideration to transient external inputs to the spiking network to emulate timevarying inputs to prefrontal networks that were likely to drive tasklocked changes in firing rates (Figure 1C, D), as suggested by Reviewers (see Point 4, below), and provide the results of the simulations in the new Figure 7. Figure 7 demonstrates that modulation of firing rate and spike synchrony in the network simulation exhibit several features that were evident in the recorded neural data. Namely, increases in spike rate were early and gradual following transient increase in network inputs, whereas increases in spike synchrony were delayed and abrupt in relation to transient increase in network inputs in the network simulation (Figure 7B, C; black), as we saw in the neural data (Figure 1C, D; black). In addition, the increase in both firing rate and spike synchrony were attenuated by reducing NMDAR synaptic currents in the model (Figure 7B, C; magenta) as they were in the neural data (Figure 1C, D; magenta). Thus, there was considerable agreement between the simulated and recorded neural data, although the neural data additionally exhibited several unique features, including additional dynamics after the response (Figure 1C, D) that likely reflected complex modulation of synaptic inputs during this time. The new simulation results characterizing the network response to transient inputs both in terms of spike synchrony and spike rate are provided in the Results section titled ‘Explaining the effects of blocking NMDAR observed in primate PFC by the prefrontal circuit model’, starting at Line 367.
2) consider how the model accommodates a Ushape relationship between NMDAR modulation and network synchrony. Existing evidence shows that 0lag synchrony in prefrontal networks is affected by manipulations that reduce NMDAR function (Zick et al. 2018) and by manipulations that enhance NMDAR function (Zick et al. 2021). The computational model presented in this manuscript does not show this Ushape behavior and the discussion does not mention this.
We are grateful to the Reviewers for raising this possibility, which we consider in a new section we have added to the revised Discussion entitled ‘Potential Ushaped relation between NMDAR function and spike synchrony’ (starting at Line 777). There we note the possibility that either too little or too much NMDAR synaptic activity could impair spike synchrony, in a manner analogous to the Ushaped relationship between D1 receptor activation and prefrontal persistent activity as reported by Vijayraghavan, Arnsten and colleagues (2007), for example. However, existing neurophysiological data do not establish yet that the relationship between spike synchrony and NMDAR function has an inverted Ushape. Specifically, we did not test the doseresponse relationship between spike synchrony and NMDAR modifying drugs in the primate drug model, nor establish that NMDAR upregulation was responsible for reduced spike synchrony in the mouse genetic model (Deleting Dgcr8 would alter expression of many downstream proteins). The absence of a Ushaped relationship between NDMAR function and spike synchrony in the model from this perspective does not seem to represent a disagreement with available biological data.
In addition, the theoretical understanding of the model should be enhanced by:
3) expanding the theoretical analysis to provide more intuition on the way NMDA works the way it does, and clarification of its distinct role in synchronization (in contrast to AMPA).
We have now extended the theoretical consideration and obtained an analytic approximation for the oscillatory instability growth rate λ describing the dependence on the AMAPR, NMDAR, GABAR synaptic conductances and external rate. Based on this consideration, we have now provided a substantially more detailed theoretical account of the precise mechanism implemented in our model for the transition between the steady and oscillatory states and the lack thereof when the NMDAR conductance is blocked. Also, the closed form analytic expression for λ allowed us to reveal 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. The results of this new theoretical analysis are presented in the new section “Dependence of oscillatory instability growth rate on synaptic parameters” in the Results and include four new figures (Figures 912). The mathematical details are given in the new section “Analytical consideration of the dependence of oscillation growth rate on network parameters” in the Methods.
4) investigating how synchrony in the model reacts to transient external inputs and how the model can produce weakly oscillatory synchrony.
In the revised manuscript we extended our consideration to transient external inputs. To this end, we carried out a new series of simulations in which stationary external input received by network neurons was replaced with a timevarying input that had a trapezoidlike temporal profile. Results of these new simulations are summarized in the new Figure 7 mentioned above. Also, as it was suggested by Reviewer #3, we found that the oscillatory synchrony became slightly weaker, which could be appreciated from relatively smaller offzero lag peaks in the population average spike correlation function as compared with the case in which external input was stationary. More generally, as one can appreciate from Figure 4 panels A2 and C2 and Figure 5B the strength of oscillatory synchrony depends on how far from the oscillatory instability boundary the synchronous network state is: the closer it is, the weaker oscillatory synchrony is, as it can be appreciated by decreasing offzero lag peaks while the 0lag peak remains present.
Finally, all reviewers were concerned about the emphasis on the association with schizophrenia, so we recommend:
5) moving speculative links to schizophrenia to the Discussion section. Also, please broaden the scope of this discussion with available literature on how synchrony and STDP interact in computational models, and with alternative explanations based on changes in NMDARdependent synaptic plasticity irrespective of synchronization.
These are good points and we thank the reviewers for raising them. We have made the suggested changes by limiting the link between the current study and schizophrenia in the Introduction to the motivation for the original neurophysiological experiments (as this link dictated the pharmacological and genetic manipulations we employed in the animal models). We have also added a new section to the Discussion with the heading ‘Spike timing disruptions and rewiring of prefrontal local circuits via STDP’ where we discuss the complexity of the interaction between STDP, synchrony, and connectivity in prior modeling studies. Namely, it is difficult to predict whether loss of synchronous spiking would cause disconnection via STDP without additional data. We acknowledge this constraint on our original hypothesis that asynchrony would cause disconnection considering these prior theoretical studies in this new section.
Reviewer #1 (Recommendations for the authors):
Here are some specific recommendations for improvement of the manuscript:
1) The last two paragraphs in the introduction are largely discussional in spirit. The paragraph before last does not make much sense here, at a point where the reader is not yet familiar with the details of the model proposed. The very last paragraph is also mostly speculative in how the results in this manuscript would be important for schizophrenia and would appear more natural in the Discussion section.
These are good points and we thank the reviewers for raising them. We have made the suggested changes by limiting the link between the current study and schizophrenia in the Introduction to the motivation for the original neurophysiological experiments (as this link dictated the pharmacological and genetic manipulations we employed in the animal models). Specifically, we have removed speculative links between the model results we present and synaptic disconnection in schizophrenia. In addition, we have removed the second to last paragraph in the original submission describing the relationship to prior modeling studies. Detailed comparison to prior work is now restricted to the Discussion as recommended.
2) The data presented in Figure 1 could be presented more fully in order to provide more constraints for the model so as to enhance its biological plausibility. For instance, it would be most interesting to have a parallel plot of the firing rate in these same neurons to see how neuronal activity evolves through the task. Neurons show a response to the cue stimulus (Zick et al. 2018), which may be similar to the response at the time of the response, and that is something that the model would have to explain: how does the control network react with strong synchronization when it receives inputs prior to the response but instead does not show any synchronization when it receives inputs at the time of the cue? Also, it would be good to show panel b extended for a time window after motor response time t=0 to compare the time course of firing rate and synchrony increases. This could also constrain further the network model.
We completely agree that the neural data used in Figure 1 should be presented more fully.
1) In the revised Figure 1, we have now added parallel plots displaying the time course of population average firing rate observed in neural data under the two conditions. We have also added a new figure, Figure 7, showing how population spike rate and synchrony obtained from spike trains of simulated networks in response to transient inputs evolve in time and compare the results of simulations with neural data presented in the revised Figure 1.
2) “how does the control network react with strong synchronization when it receives inputs prior to the response but instead does not show any synchronization when it receives inputs at the time of the cue?” This is an important point, that we address now in a new Figure 6—figure supplement 2. First, it can be now appreciated from Figure 1C and D that the population firing rate response to the cue and initial probe periods is weaker than the response at the time of the motor response. The likely reason is that during the cue/probe periods, as it is commonly assumed in modeling studies, only a small group of neurons (typically <10% in network models) receive the external/selective inputs. Second, to understand whether our model is consistent with the experimental observation that during the cue and initial probe periods no increase in the spike synchrony occurs, we carried out a series of simulations in which the fraction of neurons that received the increase in external spike rate varied from 0 to 1. The plot in Figure 6—figure supplement 2 shows spike synchrony obtained from the spike trains of simulated networks as a function of the fraction f. In this approach, the cue and probe presentations correspond to the case when relatively small fraction of neurons (f<0.2) receive the increase in external input, in contrast to the case of the preresponse period, when in the framework of our model all neurons receive increased input (f=1). It is seen that for small f the spike synchrony remains week, which is in agreement with the experimental observations.
3) We have now extended the panel B in Figure 1, as well as the new panel D displaying the time evolution of population firing rate, beyond the motor response time t=0 allowing the comparison of the time courses of synchrony and firing rate observed in neural data. In the revised Results section “Explaining the effects of blocking of NMDAR observed in primate PFC by the prefrontal circuit model” we also carry out a similar comparison for data obtained from simulated networks presented in new Figure 7 mentioned above. In addition, we compare the simulated data with neural data presented in the revised Figure 1.
3) In Figure 1 it would be good to have a dashed horizontal line at y=0 to separate correlated (y>0) from anticorrelated (y<0) activity more clearly. I would also recommend changing the xaxis label ("time lag") because this is now used in crosscorrelation functions in Figure 4 and they are intrinsically different. The synchrony plot should also include error bars and proper inference about differences between the two conditions at different time points. For clarity, it would be nice to indicate with a shaded area or horizontal bar the period of cue stimulus presentation, too.
We thank the reviewer for the valuable suggestions.
1) In the revised Figure 1 we have now added a dashed horizontal line at y=0.
2) Thank you for noticing the mislabeled axis: we have now changed the “time lag” label to “time”.
3) We have now included the standard errors depicted as shaded bands in the synchrony plots.
4) We have now marked the time points in Figure 1 when the two conditions are statistically different by asterisks and provide the details in the figure legend.
5) We have now indicated by a yellow shaded area the period of cue stimulus presentation.
4) Briefly mention in line 167 what spiking neuron model you are using
We have added a statement saying that we are using leaky integrateandfire neuron model (line 124).
5) Mention of "GABA receptors" in line 173 suggests that the model includes both GABA_A and GABA_B receptors.
We have appended the subscript A to GABA to remove the uncertainty (line 130).
6) Verb missing in "and mediated" in line 174.
We have inserted the missing verb (line 131).
7) In Figure 3 the color code red/blue means two different things: excitatory vs. inhibitory neurons, and steady vs. critical networks. This is confusing and should be clarified with different color codes for each thing.
We agree that the choice of the colors is confusing. We have now changed the color code for the firing rates of excitatory/inhibitory neurons in Figure 3 from red/blue to black/green, respectively.
Reviewer #2 (Recommendations for the authors):
1) I assume elucidation of the precise mechanism by which external inputs or NMDA affect synchrony can be achieved by analytically inspecting the direct and indirect dependence of the real part of the relevant complex eigenvalue on the synaptic conductances and external inputs, either directly or indirectly via the dependence of linearized neural gains in the background operating point (which itself depends on the external input and the synaptic parameters).
Thank you for this valuable suggestion. We have now extended the theoretical consideration and obtained an analytic approximation for the oscillatory instability growth rate λ describing the dependence on the AMAPR, NMDAR, GABAR synaptic conductances and external rate. We have accomplished this by linearizing the mean field and linear stability analysis equations near the boundary between the asynchronous and synchronous states where λ=0 . We have shown that λ essentially directly depends only on the AMPAR conductance. However, λ 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 currentoutput frequency response function. The results of this new theoretical analysis are presented in the new section “Dependence of oscillatory instability growth rate on synaptic parameters” in the Results and include three new figures. The mathematical details are given in the new section “Analytical consideration of the dependence of oscillation growth rate on network parameters” in the Methods.
In the new Results section, we have now provided a more detailed theoretical account of the precise mechanism implemented in our model for the transition between the steady and oscillatory states and the lack thereof when the NMDAR conductance is blocked. We do this in terms of the equation for the instability growth rate $\lambda $ describing the dependence on the NMDAR conductance and external rate in closed form. We also show that the critical line separating the asynchronous and synchronous states in the synaptic parameter space predicted by the analytical approximation is in excellent agreement with the exact line obtained by solving the meanfiled and linear stability analysis equations numerically. We believe that these and other related results described in the new Results section considerably improved the manuscript, and we are very grateful to the reviewer for giving us the opportunity to provide a deeper analytical insight into the involved synaptic mechanisms.
2) The authors do refer to the Compte et al. 2000 study that found that when the total strength of excitatory connectivity is fixed, tilting the relative strength of NMDAR and AMPAR towards the former disfavors synchronous oscillations (with relatively high frequency in or near the γ band). Previous studies (e.g. Wang and Brunel 2003, Tsodyks et al. 1997, and Compte et al. 2000) have pointed out the role of the negative feedback loop between excitatory and inhibitory populations mediated by the fast receptors AMPA and GABA in γband oscillations. Again I think the meanfield theory and the stability analysis based on it could be used to elucidate the differences and similarities in the pathways by which NMDAR and AMPAR affect the transition to synchrony. For example, the role of NMDAR could be via the explicit dependence of the growth rate of oscillations (the eigenvalue real part) on g_NMDA, or (most likely) via the effect of NMDARbased excitation on background rates (which indirectly affect the eigenvalue).
As implied by the reviewer, the analytic expression for the oscillation growth rate λ allowed us to reveal the differences and similarities in how AMPAR and NMDAR, both of which mediate recurrent excitation, influence the stability of the asynchronous state and the transition to synchronous oscillations. Both AMPAR and NMDAR conductances affect λ indirectly by influencing the amplitude of the neuron’s response to varying input current. However, AMPAR conductance also affects λ directly by influencing the strength of fast excitatory feedback. We have now added the corresponding text in the Discussion section.
In the new Results section mentioned in point (1), we explain that the direct dependence manifests the essential influence of the AMPAR synaptic conductance on the strength of excitatoryinhibitory feedback loops via fast excitatory to excitatory and excitatory to inhibitory recurrent connections. The indirect dependence manifests the influence of the synaptic parameters, including NMDAR and external rate, on the location of the operating point on the neuron’s currentfrequency 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 the feedback loop.
Also, in the new Results section, we have provided a theoretical account of the simulation results observed in Compte et al., 2000, when concurrent increase in the AMPAR and decrease in the NMDAR conductances leads to synchronous oscillations. As in the case of our model, the explanation is given in terms of the equation for the instability growth rate λ that describes the direct (AMPAR) and indirect (both AMPAR and NMDAR) synaptic conductance contributions to the emergence of oscillations. In addition, we have also revised the Discussion section where we address the apparent contradiction between Compte et al., 2000, simulations and Zick et al., 2018, experiments, emphasizing the critical role played by the direct and indirect mechanisms by which NMDAR and AMPAR affect the transition from the asynchronous to the synchronous dynamic.
3) Regarding the link to "disconnection" via STDP. Various studies have found that STDP in (initially randomly connected) networks in the asynchronous state can maintain random connectivity (Morrison et al. 2007 – DOI: 10.1162/neco.2007.19.6.1437) or even lead to the emergence of structured connectivity (Izhikevich et al. 2004, Babadi and Abbott 2013, and LitwinKumar and Doiron 2014  DOIs: 10.1093/cercor/bhh053, 10.1371/journal.pcbi.1002906, 10.1038/ncomms6319). In lieu of evidence (either previous studies or new simulations/modeling by the authors) for this claim, I suggest weakening the claim to saying that the NMDARrelated reduction in synchrony, as a specific disruption in the patterns of precise spike timing, can lead to abnormal changes in patterns of synaptic connectivity via STDP. I also suggest citing the abovementioned (or similar) studies on previous studies of STDP in asynchronous (or synchronous) states, or any other studies that support specific claims by the authors in this regard.
We thank the reviewer for raising this important point. It motivated us to reexamine our initial hypothesis that asynchrony at the level of spiking neurons would necessarily be expected to disconnect recurrent networks considering prior theoretical studies examining STDP. We consider constraints on our hypothesis imposed by the prior theoretical studies in a new section we have added to the Discussion (‘Spike timing disruptions and rewiring of prefrontal local circuits via STDP’, see new text starting on line 739). There we acknowledge that the relationship between spike timing dynamics and synaptic connectivity in networks incorporating STDP is hard to predict because changes in spike synchrony and synaptic connectivity influence one another in complex ways after STDP is implemented and spiking patterns are allowed to sculpt synaptic connectivity. We thank the reviewer for pointing out this caveat and we now acknowledge it as a constraint on our prior hypothesis in the Discussion. However, to our understanding, prior network stimulations incorporating STDP have obtained divergent results as to how synchrony and connectivity influence one another that appear to depend strongly on the initial assumptions of the models. These crucial factors appear to include delays in signal conduction imposed by axons, synapses, and dendrites, and how network simulations model these delays. For example, delays in signal conduction imposed by the length of axons in recurrent connections between pyramidal neurons can be on the order of tens of ms, e.g., very large in comparison to the STDP window. These delays reflecting anatomical details about cortical network architecture therefore have profound effects on synaptic plasticity via STDP mechanisms (as demonstrated in the Izhikevich et al., 2004 paper). Indeed, how these delays are modeled appears to contribute to divergent results between theoretical studies (comparing the Izhikevich et al.,2004 and Morrison et al., 2007 papers). Synchronous spiking is not required to maintain structure in some STDP networks (Izhikevich 2004) yet correlated synchronous spiking can induce structure in other STDP networks (LitwinKumar and Doiron, 2014), according to the specific assumptions adopted. We are grateful to the reviewer for posing these questions, they seem essential. We hope these points and the data we provide will motivate further integrated computational and experimental studies of cortical local circuit dynamics.
Other suggestions:
– Figure 6 top: here the crosscorrelation plot based on monkey data (as in the model) exhibits nonzero lag "γ" peaks. But these seem to be absent in the data shown in the Neuron paper (e.g. Figure 4D therein). Can the authors comment on this seeming discrepancy? Is this because data only from one of the monkeys (showing stronger effects) was included in the current paper? Can the authors also comment on the flatness of the CCH in the "initial probe period" in the monkey data in contrast to the small 0lag peak in the model CCH?
1) The absence of nonzero “γ” lag peaks in Figure 4, Neuron paper, is likely due to the fact that the crosscorrelation depicted in Figure 4D is the average of the crosscorrelation function over the trial time. Since the strongest spike synchrony is restricted to the end of the trial around the time of the motor response, flanking peaks in the CCH may be washed when averaging over the trial, with only the stronger, central peak remaining. In addition, restriction of the current data to the animal with the strongest spike synchrony may have contributed to this difference.
2) Flatness of the CCH in the "initial probe period" in the monkey data in contrast to the small 0lag peak in the model CCH. This question is now addressed in Figure 6—figure supplement 1. We note that in a large (network size $N\gg 1$) and sparsely connected network (connection probability $p\ll 1$), when the population average firing rate is stationary and the network is in the asynchronous state, spike synchrony vanishes in the limit $N\to \mathrm{\infty},\phantom{\rule{thinmathspace}{0ex}}p\to 0$, while the average number of connections of individual neurons $C=pN$ remains finite. The small 0lag peak seen for the asynchronous network in panel Figure 6 A2 is due to the effect of the finite size of the network (Brunel and Hakim, 1999) in the simulations. The plot in Figure 6—figure supplement 1 shows spike synchrony obtained from the spike trains of the asynchronous steady (red circles) and synchronous oscillatory networks (blue circles) as a function of the simulated network size (from left to right, $N=$ 20,000, 10,000, 5,000, and 3,333). The red and blue lines represent linear fits to the corresponding simulated data. It can be seen that, as the sizes of the networks increase, the synchrony of the asynchronous network extrapolates to zero in the limit $N\to \mathrm{\infty}$, whereas for the synchronous network it remains finite.
– Lines 104105: I found the meaning of *balance between* in this sentence vague/unclear. What does the balance between NMDA inputs and oscillatory activity mean?
We agree this was a vague statement and we have removed it from the Introduction in the revision.
– I suggest the authors also mention (e.g. around lines 476477 where they discuss their assumption 1, and as further support for this assumption) that in their empirical data (Neuron paper) they do find an increase in PFC neural firing rates in the preresponse period, which will presumably also arise in their model from the increase in external inputs.
Thank you for the suggestion. We added a corresponding statement and referred to a new panel D in Figure 1 showing the time course of population firing rate in empirical data (line 344).
– Line 781: "themes" → "times".
Corrected (line 936).
– Line 838: I would explain where this equation comes from (I understand it accounts for the effect of reset on mean membrane potential), and also what v_α denotes.
Thank you for alerting us to this oversight. We now explain that this equation is for the average membrane potential $\u27e8{V}_{\alpha}\u27e9$. It was calculated in Renart et al., 2003, using the distribution of potentials obtained in Brunel and Hakim, 1999, and accounts for neuron spiking (i.e., the effect of reset and refractory period) and synaptic noise (line 1000).
Reviewer #3 (Recommendations for the authors):
Here are issues that I believe the authors should address before the paper can be published:
1. In Figure 1, the authors only show the synchrony before the motor response. I believe it would be instructive to show also what happens after motor response, to visualize how temporally precise is this increase. It would also be useful to show the time course of average firing rates. The explanation of the authors relies on an increase in external input to generate an increase in synchrony – this increase in external input should be reflected in an increase in the average firing rate as well.
Thank you for these valuable suggestions. We have now added new panels C, and D to Figure 1 that show the time courses of average firing rates. We also extended panel B as well as the new panel D beyond the time of the motor response. It can be now appreciated from Figure 1B that the peak of the synchrony modulation matches the onset of the motor response. Comparing Figures 1B and 1D one can also appreciate that, consistent with the model’s assumption about the increase in external input, the empirical firing rates gradually increase before the motor response. We have added a corresponding statement (see new text at Lines 110115) to support our assumption about the increase in external input by referring to the increase in the empirical firing rates shown in the new panel D, Figure 1. Additionally, in a new Figure 7 we have now provided a similar comparison between the time courses of synchrony, population spike rate, and external input for the network model based on spike trains obtained in a series of new simulations in which the network received transient external inputs (see point 3 below).
2. Figure 6 shows a strong difference between the network model and data – side peaks in spike correlations are much more prominent in the model than in the data. Thus, synchronization in the data is much less oscillatory than in the model. Is there a way in the model to get such weakly oscillatory synchrony? If not, can the authors speculate as to what additional mechanisms would lead to a reduction of γ power, preserving zero lag synchrony?
In general, as one can appreciate from Figure 4 panels A2 and C2 and Figure 5B the strength of oscillatory synchrony depends on how far from the oscillatory instability boundary the synchronous network state is: the closer it is, the weaker oscillatory synchrony is, as it can be appreciated by decreasing offzero lag peaks while the 0lag peak remains present. This suggests that inputs to prefrontal networks in the neural data push the network just past the state boundary into the domain of oscillatory activity, whereas we consider larger excursions through the synaptic parameter space in our simulations to better characterize the full range of effects.
3. In the data, the increase in synchrony is transient. In the simulations, the author considers the case of stationary external inputs. I believe the authors should also investigate the case of transient increases of external inputs leading to a more realistic time course of synchrony. Such a scenario could also help in solving the issue mentioned in 2.
In the revised manuscript we extended our consideration to transient external inputs. To this end, we carried out a new series of simulations in which stationary external input received by network neurons was replaced with a timevarying input that had a trapezoidlike temporal profile. Results of these new simulations are presented in the Results section “Explaining the effects of blocking of NMDAR observed in primate PFC by the prefrontal circuit model,” and are summarized in the new Figure 7 mentioned above providing the time courses for synchrony and firing rate obtained from spike trains of simulated networks in response to transient external inputs. In the revised Results section, we now compare the dynamics of firing rate and synchrony and association between them observed in simulations with those observed in experiments. Figure 7 demonstrates that modulation of firing rate and spike synchrony in response to a transient input to the network exhibit several features that were evident in the recorded neural data. Namely, increases in spike rate were early and gradual following transient increase in network inputs, whereas increases in spike synchrony were delayed and abrupt in relation to transient increase in network inputs in the network simulation (Figure 7B, C; black), as we saw in the neural data (Figure 1C, D; black). We acknowledge, however, that while our relatively simple model qualitatively is consistent with dynamical features of the firing rate and synchrony observed in primate PFC, there are some quantitative discrepancies in firing rates.
[Editors’ note: what follows is the authors’ response to the second round of review.]
The manuscript has been improved but there are some remaining issues that need to be addressed, as outlined below:
1. The dynamics of average firing rates in Figure 1 should mark clearly the various relevant periods of the task on the time axis: probe interval, motor response, etc.
We expanded the time axis in the left panels of Figure 1 to include the whole trial period. The probe presentation time is now indicated by green shaded area whereas the mean and standard deviation of the motor response time by color coded horizontal bars.
2. The new theoretical analyses in this revision have greatly extended the length of the main text, but the additional insight provided is not clearly explained and in its current form it can lead to confusion. It is suggested that most of this material is moved to a supplementary document and a more concise and clearly delineated section in the main text provides the fundamental message of how NMDARs influence network synchronization.
We have now moved large parts of this material, including three figures, from Results to newly added Appendices 15. We have also revised the remaining part in the Results section to improve the clarity of the text.
Reviewer #3 (Recommendations for the authors):
The paper has been considerably improved and I would like to congratulate the authors for their hard work. There are still a couple of issues though that I think should be addressed to improve the paper further.
Thank you for the positive comments and valuable suggestions.
1. In response to one of the main concerns of the reviewers, the authors have expanded Figure 1 to include the dynamics of average firing rates in the data and the effect of NMDAR blockade on this dynamic. This is a welcome addition, but I think this figure could be improved to better understand the time course of the task. Currently, there is no discussion at all about when the motor response is happening and so the left and right panels seem totally disconnected. One has to go back to the original paper by Zick et al. to understand when the motor response happens. I would suggest extending the left panels to include the whole trial (including the `probe' interval), adding where the mean response time is similar to Figure 4 in the Zick et al. 2018 paper. Or at the very least indicate when is the mean response time (time 0 in right panels) with respect to the time axis of the left panels.
As suggested, we expanded the time axis in the left panels to include the whole trial period. The probe presentation time is now indicated by green shaded area whereas the mean and standard deviation of the motor response time by color coded horizontal bars. We also added accompanying explanation in the figure legend (lines 14911493).
2. In response to another concern, the authors have considerably expanded their theoretical analysis with a detailed description of an approximation for the instability growth rate. However, while the approximation itself represents an interesting new derivation, I did not find that it brings much insight, and in fact in some ways introduces additional confusion, for the following reasons. (i) The main new result of this new analysis is the new Equation (1). In this equation, the LambdaRs play critical roles and are the subject of much of the later discussion, however, their meaning is not really explained in the main text and one has to go to the Methods to understand what they really are. (ii) The new Figure (10) indicates that the growth rate is primarily determined by the AMPA term and that the GABA term plays little role. This is quite confusing as oscillations in such networks relies on GABAR mediated negative feedback. How can the authors explain this fact? (iii) After reading the new part, I ended up not having much more insight about mechanisms. My guess at this point is that NMDAR influences the synchronisation properties of the network by providing an additional excitatory drive, thereby playing a similar role as the external excitatory input. I have the feeling that is what the authors are trying to say when they say that the NMDAR term acts through changes in phiprime (the slope of the transfer function), which is exactly how the external inputs act, but it feel like it could be explained in a much clearer way.
(i) We thank the reviewer for noting that the meaning of the factors $\mathrm{\Lambda}}_{\text{AMPA}$, $\mathrm{\Lambda}}_{\text{NMDA}$, and $\mathrm{\Lambda}}_{\text{GABA}$ have not been clearly explained in the main text. We believe that this fact also caused the confusion expressed in point (ii), which we address below. We have now revised the corresponding section of the manuscript (lines 465477). The factors $\mathrm{\Lambda}}_{\text{AMPA}$, $\mathrm{\Lambda}}_{\text{NMDA}$, and $\mathrm{\Lambda}}_{\text{GABA}$ govern the strength of the direct and indirect contributions of the changes in the synaptic conductances $\mathrm{\Delta}{g}_{\text{AMPA}}$, $\mathrm{\Delta}{g}_{\text{NMDA}}$, and $\mathrm{\Delta}{g}_{\text{GABA}}$ to the oscillatory instability. The direct contribution of the conductance change $\mathrm{\Delta}{g}_{R}$ is due to the explicit dependence of $\lambda$ on the $\mathrm{\Delta}{g}_{R}$ in Equation 1. The indirect contribution is due to the implicit dependence of $\lambda$ on the $\mathrm{\Delta}{g}_{R}$ in Equation 1 via the changes in the slopes of the currentfrequency response functions $\mathrm{\Delta}{\varphi}_{{I}_{syn,E}}^{\prime}$ and $\mathrm{\Delta}{\varphi}_{{I}_{syn,I}}^{\prime}$ of excitatory and inhibitory neurons, respectively. These changes in the slopes occur because the change in the conductance $\mathrm{\Delta}{g}_{R}$ produces changes in the average total synaptic currents of excitatory and inhibitory neurons. These, in turn, change the operating points of the response functions and, hence, the slopes of these functions at the operating points (lines 436443). We explain that the strength of the direct contribution of the change in the synaptic conductance $\mathrm{\Delta}{g}_{R}$ to the instability growth rate $\lambda$ is determined solely by the corresponding factor $\mathrm{\Lambda}}_{R$ via the term $\mathrm{\Lambda}}_{R}\frac{\mathrm{\Delta}{g}_{R}}{{g}_{R}^{\ast}$ in Equation 1. However, the strength of the indirect contribution of the $\mathrm{\Delta}{g}_{R}$ is determined by all three, $\mathrm{\Lambda}}_{\text{AMPA}$, $\mathrm{\Lambda}}_{\text{NMDA}$, and $\mathrm{\Lambda}}_{\text{GABA}$, factors via the terms involving changes in the slopes $\mathrm{\Delta}{\varphi}_{{I}_{syn,E}}^{\prime}$ and $\mathrm{\Delta}{\varphi}_{{I}_{syn,I}}^{\prime}$. In addition, we now provide a specific example for the change in the GABAR conductance $\mathrm{\Delta}{g}_{\text{GABA}}$ (lines 471477) and explain the resulting contributions to the instability growth rate $\lambda$ determined by each of the $\mathrm{\Lambda}}_{\text{AMPA}$, $\mathrm{\Lambda}}_{\text{NMDA}$, and $\mathrm{\Lambda}}_{\text{GABA}$ factors.
(ii) Figure 9 (Figure 10 in previous version) demonstrates that contributions to the oscillatory instability growth rate $\lambda$ from the terms involving $\mathrm{\Lambda}}_{\text{NMDA}$ and $\mathrm{\Lambda}}_{\text{GABA}$ in Equation 1 are small compared to the contribution involving the term $\mathrm{\Lambda}}_{\text{AMPA}$. Figure 9 also demonstrates that this does not mean that changes in the NMDAR and GABAR conductances do not affect the oscillatory instability. Indeed, one can see from Figure 9 C and D that changes in the instability growth rate $\lambda$ (black line) in response to the changes in the NMDAR (panel C) and GABAR (panel D) conductances are comparable to the changes in $\lambda$ in response to the change in the AMPAR conductance (panel B). This is because changes in the NMDAR and GABAR conductances, while do not affect directly, influence the oscillatory instability growth rate $\lambda$ indirectly by changing the average total synaptic current (excitatory drive) and, therefore, the operating point of the neuron’s currentfrequency response function. We provide this explanation in a new paragraph (lines 485490).
(iii) The reviewer’s understanding of our explanation is correct: NMDAR influences the synchronisation properties of the network by providing an additional excitatory drive, thereby playing a similar role as the external excitatory input. To improve the overall readability and clarity of this section we have now moved some parts of it, including three figures, from the main text to newly added Appendices 15. We have also revised the main text to provide a more concise and clear explanation of how the NMDAR conductance and external input modulations affect synchrony (lines 506524).
https://doi.org/10.7554/eLife.79352.sa2Article and author information
Author details
Funding
National Institute of Mental Health (R01MH107491)
 Matthew V Chafee
National Institute of Mental Health (P50MH119569)
 Matthew V Chafee
National Institute of Mental Health (R25 MH101076)
 Rachael K Blackman
National Institute of General Medical Sciences (T32 GM008244)
 Rachael K Blackman
National Institute of General Medical Sciences (T32 HD007151)
 Rachael K Blackman
University of Minnesota Foundation
 Bagrat Amirikian
U.S. Department of Veterans Affairs
 Bagrat Amirikian
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We would like to provide our special thanks to Apostolos P Georgopoulos for his intellectual contributions to the study of the temporal correlations in networks and the temporal dynamics of neural ensembles that provided part of the motivation for this work, as well as for his constant support and encouragement. We thank Dean Evans for lab and project management as well as his assistance with surgeries, animal care, and neural recordings; Dale Boeff for his assistance with neurophysiological recording system design and construction, as well as computer programming for signal processing and data analysis; Sofia Sakellaridi for her assistance with neural recordings; Olivia Newman for preliminary analysis. This work was performed while RKB was employed at the University of Minnesota. The opinions expressed in this article are the author's own and do not reflect the views of the Food and Drug Administration, the Department of Health and Human Services, or the United States Government. Support for this work was provided by the Wilfred Wetzel Graduate Fellowship (to RKB), the National Institute of General Medical Sciences (T32 GM008244 and T32 HD007151 to RKB), and the National Institute of Mental Health (R25 MH101076 to RKB, and R01MH107491 and P50MH119569 to MVC). Partial funding for this study was provided by the University of Minnesota Foundation and the U.S. Department of Veterans Affairs. The sponsors had no role in the current study design, analysis or interpretation, or in the writing of this paper. The contents do not represent the views of the U.S. Department of Veterans Affairs, the National Institutes of Health, the Department of Health and Human Services, or the United States Government.
Ethics
All animal care and experimental procedures conformed to National Institutes of Health (NIH) guidelines and were approved by the Institutional Animal Care and Use Committees of the Minneapolis Veterans Administration Medical Center (VA Animal Component of Research Protocol Number 170601).
Senior Editor
 John R Huguenard, Stanford University School of Medicine, United States
Reviewing Editor
 Albert Compte, IDIBAPS, Spain
Reviewers
 Albert Compte, IDIBAPS, Spain
 Yashar Ahmadian, University of Cambridge, United Kingdom
Version history
 Received: April 8, 2022
 Preprint posted: June 13, 2022 (view preprint)
 Accepted: November 24, 2023
 Version of Record published: February 6, 2024 (version 1)
Copyright
© 2024, Crowe et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 188
 Page views

 20
 Downloads

 0
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Neuroscience
Basal forebrain cholinergic neurons modulate how organisms process and respond to environmental stimuli through impacts on arousal, attention, and memory. It is unknown, however, whether basal forebrain cholinergic neurons are directly involved in conditioned behavior, independent of secondary roles in the processing of external stimuli. Using fluorescent imaging, we found that cholinergic neurons are active during behavioral responding for a reward – even prior to reward delivery and in the absence of discrete stimuli. Photostimulation of basal forebrain cholinergic neurons, or their terminals in the basolateral amygdala (BLA), selectively promoted conditioned responding (licking), but not unconditioned behavior nor innate motor outputs. In vivo electrophysiological recordings during cholinergic photostimulation revealed rewardcontingencydependent suppression of BLA neural activity, but not prefrontal cortex. Finally, ex vivo experiments demonstrated that photostimulation of cholinergic terminals suppressed BLA projection neuron activity via monosynaptic muscarinic receptor signaling, while also facilitating firing in BLA GABAergic interneurons. Taken together, we show that the neural and behavioral effects of basal forebrain cholinergic activation are modulated by reward contingency in a targetspecific manner.

 Neuroscience
Orbitofrontal cortex (OFC) is classically linked to inhibitory control, emotion regulation, and reward processing. Recent perspectives propose that the OFC also generates predictions about perceptual events, actions, and their outcomes. We tested the role of the OFC in detecting violations of prediction at two levels of abstraction (i.e., hierarchical predictive processing) by studying the eventrelated potentials (ERPs) of patients with focal OFC lesions (n = 12) and healthy controls (n = 14) while they detected deviant sequences of tones in a local–global paradigm. The structural regularities of the tones were controlled at two hierarchical levels by rules defined at a local (i.e., between tones within sequences) and at a global (i.e., between sequences) level. In OFC patients, ERPs elicited by standard tones were unaffected at both local and global levels compared to controls. However, patients showed an attenuated mismatch negativity (MMN) and P3a to local prediction violation, as well as a diminished MMN followed by a delayed P3a to the combined local and global level prediction violation. The subsequent P3b component to conditions involving violations of prediction at the level of global rules was preserved in the OFC group. Comparable effects were absent in patients with lesions restricted to the lateral PFC, which lends a degree of anatomical specificity to the altered predictive processing resulting from OFC lesion. Overall, the altered magnitudes and time courses of MMN/P3a responses after lesions to the OFC indicate that the neural correlates of detection of auditory regularity violation are impacted at two hierarchical levels of rule abstraction.