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
Article 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).
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

 306
 views

 32
 downloads

 1
 citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
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
Maternal choline supplementation (MCS) improves cognition in Alzheimer’s disease (AD) models. However, the effects of MCS on neuronal hyperexcitability in AD are unknown. We investigated the effects of MCS in a wellestablished mouse model of AD with hyperexcitability, the Tg2576 mouse. The most common type of hyperexcitability in Tg2576 mice are generalized EEG spikes (interictal spikes [IIS]). IIS also are common in other mouse models and occur in AD patients. In mouse models, hyperexcitability is also reflected by elevated expression of the transcription factor ∆FosB in the granule cells (GCs) of the dentate gyrus (DG), which are the principal cell type. Therefore, we studied ΔFosB expression in GCs. We also studied the neuronal marker NeuN within hilar neurons of the DG because reduced NeuN protein expression is a sign of oxidative stress or other pathology. This is potentially important because hilar neurons regulate GC excitability. Tg2576 breeding pairs received a diet with a relatively low, intermediate, or high concentration of choline. After weaning, all mice received the intermediate diet. In offspring of mice fed the high choline diet, IIS frequency declined, GC ∆FosB expression was reduced, and hilar NeuN expression was restored. Using the novel object location task, spatial memory improved. In contrast, offspring exposed to the relatively low choline diet had several adverse effects, such as increased mortality. They had the weakest hilar NeuN immunoreactivity and greatest GC ΔFosB protein expression. However, their IIS frequency was low, which was surprising. The results provide new evidence that a diet high in choline in early life can improve outcomes in a mouse model of AD, and relatively low choline can have mixed effects. This is the first study showing that dietary choline can regulate hyperexcitability, hilar neurons, ΔFosB, and spatial memory in an animal model of AD.

 Neuroscience
Brain structural circuitry shapes a richly patterned functional synchronization, supporting for complex cognitive and behavioural abilities. However, how coupling of structural connectome (SC) and functional connectome (FC) develops and its relationships with cognitive functions and transcriptomic architecture remain unclear. We used multimodal magnetic resonance imaging data from 439 participants aged 5.7–21.9 years to predict functional connectivity by incorporating intracortical and extracortical structural connectivity, characterizing SC–FC coupling. Our findings revealed that SC–FC coupling was strongest in the visual and somatomotor networks, consistent with evolutionary expansion, myelin content, and functional principal gradient. As development progressed, SC–FC coupling exhibited heterogeneous alterations dominated by an increase in cortical regions, broadly distributed across the somatomotor, frontoparietal, dorsal attention, and default mode networks. Moreover, we discovered that SC–FC coupling significantly predicted individual variability in general intelligence, mainly influencing frontoparietal and default mode networks. Finally, our results demonstrated that the heterogeneous development of SC–FC coupling is positively associated with genes in oligodendrocyterelated pathways and negatively associated with astrocyterelated genes. This study offers insight into the maturational principles of SC–FC coupling in typical development.