A prefrontal network model operating near steady and oscillatory states links spike desynchronization and synaptic deficits in schizophrenia

  1. David A Crowe
  2. Andrew Willow
  3. Rachael K Blackman
  4. Adele L DeNicola
  5. Matthew V Chafee  Is a corresponding author
  6. Bagrat Amirikian  Is a corresponding author
  1. Department of Biology, Augsburg University, United States
  2. Department of Neuroscience, University of Minnesota, United States
  3. Medical Scientist Training Program (MD/PhD), University of Minnesota, United States
  4. Brain Sciences Center, VA Medical Center, United States
  5. Center for Cognitive Sciences, University of Minnesota, United States

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 0-lag 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 0-lag spike synchrony we observed between neurons in monkey prefrontal cortex during behavior, as well as the suppression of this 0-lag spiking by administration of NMDAR antagonists. As such, our cortical network model provides a plausible mechanism explaining the link between NMDAR synaptic and 0-lag 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 mean-field 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.sa0

eLife digest

Schizophrenia is a long-term 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 schizophrenia-causing 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 (Goldman-Rakic, 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 (‘0-lag’) 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 0-lag 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 integrate-and-fire 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 0-lag 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 0-lag 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 0-lag 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 dot-pattern expectancy (DPX) task, a task that measures specific deficits in cognitive control in schizophrenia (Jones et al., 2010). In the DPX task, the correct response (left or right joystick movement) to a probe stimulus depends on a preceding cue followed by a delay period (Materials and methods).

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

Population average spike rate and synchrony between spike trains of neuron pairs recorded during the DPX task as a function of time.

Plots show time evolution of spike synchrony (A, B) and spike rate (C, D) estimated with 100 ms temporal resolution for drug-naive (black) and drug (magenta) conditions. Spike synchrony was measured with 1 ms resolution, and only neuron pairs for which a reasonably reliable estimation of synchrony could be achieved contributed to the plots (see Materials and methods). (A, C): Trials are aligned to the cue onset (t=0 ms); in all trials, the cue was presented until t=1,000 ms (yellow shaded area), followed by a 1,000 ms delay period, after which the probe was presented at t=2,000 ms for 500 ms (green shaded area). Color-coded horizontal error-bars indicate the mean and standard deviation of the motor response time for the corresponding drug condition. The numbers of contributing pairs for drug-naive and drug conditions are 524 and 195 (A), and the number of neurons are, correspondingly, 514 and 343 (C). (B, D): Trials are aligned to the time of motor response (t=0 ms) to show the temporal modulation of synchrony and spike rate during the last 600 ms immediately preceding the response. Color-coded horizontal error-bars indicate the mean and standard deviation of the probe presentation time for the corresponding drug condition. The numbers of contributing pairs for drug-naive and drug conditions are 661 and 223 (B), and the number of neurons are, correspondingly, 538 and 343 (D). Shaded grey and magenta bands show the standard errors for spike synchrony (A, B) and rate (C, D). Green asterisks show the instances of times when the drug-naive and drug conditions are statistically different (false discovery rate 0.05 [Benjamini and Hochberg, 1995] using two-sample t-test p-values).

Network model and theoretical framework

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

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

Recurrent synaptic currents of excitatory connections are two-component, mediated by AMPA and NMDA receptors, whereas currents of inhibitory connections are mediated by GABAA 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 gX,α, gAMPA,α, gNMDA,α, gGABA,α corresponding to the external AMPA, recurrent AMPA, NMDA, and GABA currents (α=E,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 νE and ν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 vE0 and vI0 in the asynchronous stationary state of the network can be effectively parametrized by three parameters expressed as ratios of component synaptic currents: IAMPA/IGABA, INMDA/IGABA, and IX,E/Iθ,E, where IR is the mean current of the R-receptor mediated synapse (R=X,AMPA,NMDA,GABA), and Iθ,E is the current that is needed for an excitatory neuron to reach firing threshold θ 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 νX one can solve the mean field equations to obtain the underlying eight synaptic conductances providing the desired population mean firing rates vE0 and vI0 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, λ and ω, describing the rate of instability growth and the oscillation frequency. The asynchronous state is stable when λ<0; in this case small perturbations of firing rates cause exponentially damped oscillation of network activity. The case λ=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 ω; 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 λ>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, INMDA/IGABA, and vary the remaining two parameters: the balance between recurrent excitation and inhibition, IAMPA/IGABA, and the balance between external excitation and firing threshold, IX,E/Iθ,E. For a given point in this (IAMPA/IGABA,IX,E/Iθ,E) parameter plane we solve the mean field equations to find the underlying set of eight synaptic conductances that provide the prescribed rates vE0 and vI0 given external spike rate ν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 νX=5 Hz, the rates of excitatory and inhibitory populations are set to vE0=5 Hz, vI0=20 Hz, and the NMDA current balance is fixed at INMDA/IGABA=0.15. The diagram shows solutions for λ obtained from the linear stability analysis in the (IAMPA/IGABA,IX,E/Iθ,E) parameter space. The asynchronous stationary state corresponds to the region where λ<0, whereas the synchronous oscillation state is realized in the region where λ>0. The asynchronous and synchronous states are separated by a “critical” or instability line on which λ=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, fntwrk=ω/2π, as a function of the balance between the recurrent AMPA and GABA currents, IAMPA/IGABA, is shown in Figure 2B.

Figure 2 with 1 supplement see all
Characteristics of the system predicted by the linear stability analysis.

Parameters are as follows: prescribed firing rates of excitatory and inhibitory populations are 5 Hz and 20 Hz, respectively; external input spike rate is 5 Hz; and the balance between NMDA and GABA currents is fixed at 0.15. (A:) State diagram in the (IAMPA/IGABA,IX,E/Iθ,E) parameter plane showing color coded value of the rate of instability growth λ: in the region of the parameter space where λ<0 the asynchronous state is stable, whereas the region where λ>0 corresponds to the synchronous oscillation state. The two regimes are separated by a critical line on which λ=0. This boundary, shown by a white line, is the locus where the stationary network dynamic becomes unstable, and oscillatory population activity develops. Each point in this parameter plane corresponds to a network with a specific set of eight synaptic conductances provided by the mean field approximation. Red and blue asterisks are the points in the state diagram corresponding to the steady and critical primary networks, respectively (see Selection of Primary Networks in Results). (B:) Network oscillation frequency that develops on the critical line as a function of the balance between AMPA component of recurrent excitation and inhibition.

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, fntwrk, essentially is independent of the INMDA/IGABA 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 drug-naive and drug conditions (Zick et al., 2018), we make two assumptions regarding neural and synaptic activity: (1) the increase in spike synchronization observed before the monkey’s response in Zick et al., 2018 is due to task-specific external afferent signals received by PFC neurons after probe presentation; (2) administration of NMDAR antagonist results in blocking NMDAR mediated synaptic currents. In the framework of our model, we implemented these assumptions as follows: task specific external signals were accounted for by an increase in the external spike rate from its background level νX, whereas the effect of drug administration was modeled by setting NMDAR conductances gNMDA,E and gNMDA,I to zero.

Next, to investigate how spike synchrony in asynchronous and synchronous networks depends on the modulations of vX and gNMDA,α, 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 vE* and vI* for a given external spike rate vX*. 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 vX and NMDAR conductance gNMDA,α are varied relative to their standard values vX* and gNMDA,α*, 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 vE*=5 Hz and vI*=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 vX* to be the same as the excitatory population rate vE* inside the model and, thus, set vX*=5 Hz. Lastly, for both networks, we fix the balance between NMDA and GABA currents at INMDA*/IGABA*=0.15. Note that the state diagram in the (IAMPA/IGABA,IX,E/Iθ,E) 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 (IAMPA/IGABA,IX,E/Iθ,E) plane corresponds to a network with a specific set of synaptic conductances. For synchronous regime, we look for a network on the critical line (λ=0, white line in Figure 2A), at the onset of oscillatory instability with a frequency in the γ-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 (IAMPA/IGABA=0.4,IX,E/Iθ,E=1.09) corresponds to such a network with oscillation frequency fntwrk*~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 (λ<0). The point marked by a red asterisk in Figure 2A located at (IAMPA/IGABA=0.2,IX,E/Iθ,E=1.09) 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 gGABA,α*,gNMDA,α*,gAMPA,α*,gxAMPA,α* (α=E,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=5,000 neurons, of which NE=4,000 are excitatory and NI=1,000 inhibitory. Neurons are connected randomly with a probability p=0.2. Figure 3 illustrate the behavior of simulated networks with synaptic conductance parameters corresponding to the steady and critical primary networks indicated by the red and blue asterisks, respectively, in the state diagram presented in Figure 2A. The dynamic behavior is shown at the level of individual cell activity (spike rasters, top of panels in Figure 3), as well as whole population activity (bottom of panels in Figure 3).

Figure 3 with 1 supplement see all
Simulations of networks composed of 4,000 excitatory and 1,000 inhibitory neurons connected randomly with probability 0.2.

Conductance parameters are solutions of mean field equations for the steady state primary network (A1, B1) and the critical state primary network (A2, B2) corresponding to the red and blue asterisks, respectively, in the (IAMPA/IGABA,IX,E/Iθ,E) state plane shown in Figure 2A and inset. (A1, B1), (A2, B2): Top, spike rasters (sorted by rate) of 200 excitatory (black) and 50 inhibitory (green) neurons. Bottom, time-varying activity (1ms resolution) of excitatory (black) and inhibitory (green) populations. (A1, A2): External input spike rate νX=5 Hz. Excitatory and inhibitory neurons display average firing rates of, respectively, 5.3 Hz and 20 Hz (A1), and 6.3 Hz and 22 Hz (A2). (B1, B2): In these simulations νX was increased by 5%. Excitatory and inhibitory neurons display average firing rates of, respectively, 7.5 Hz and 25 Hz (B1), and 12 Hz and 34 Hz (B2).

In simulations shown in Figure 3 panels A1 and A2 external spike rate νX was fixed at the level of vX*=5 Hz chosen for the primary networks. It is seen that excitatory and inhibitory neurons exhibit highly irregular firing with average rates, νE and ν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 νE and νI are in good agreement with the prescribed rates vE*=5 Hz and vI*=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 νX was increased by 5% relative to the rate vX* used in simulations illustrated in Figure 3 panels A1 and A2. For the steady state primary network (Figure 3B1), the firing rates of excitatory and inhibitory neurons increase with the external drive. However, the regime of network dynamics qualitatively does not change and remains asynchronous (compare Figure 3—figure supplement 1A1 vs Figure 3—figure supplement 1B1). In contrast, stronger external inputs received by the critical state primary network synchronize population activity (Figure 3B2). It is seen that while individual neurons continue to fire irregularly, population activity now clearly exhibits oscillatory behavior, indicating that the network is in synchronous irregular regime in which the average firing frequency of neurons is low, about 20 Hz, compared to the frequency of network oscillation, which is about 50 Hz (see Figure 3—figure supplement 1B2). This frequency is close to the theoretically predicted network frequency of 58 Hz near the onset of oscillation.

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

To facilitate the comparison of characteristic features exhibited by a simulated network with experimentally measurable quantities, we compute temporal correlation of spiking activity that quantifies average pairwise correlation between spike trains of excitatory neurons. In the context of the DPX task performed in drug-naive and drug conditions studied in Zick et al., 2018 and with the purpose of elucidating the mechanism of drug-induced desynchronization of spiking activity, we investigated how temporal correlations depend on the strength of external drive and the NMDAR mediated synaptic current. To this end, we varied external input rate vX and the NMDAR conductance parameters gNMDA,E and gNMDA,I relative to their respective standard values vX*, and gNMDA,E* and gNMDA,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 gNMDA,E/gNMDA,E* and gNMDA,I/gNMDA,I* are the same; in the following we drop the E,I designation.

Figure 4 displays correlation of spiking activity (panels A1, A2, C1, C2) and synchrony (0-lag correlation, panels B1, B2, D1, D2) obtained from spike trains of simulated steady (panels A1, B1, C1, D1) and critical (panels A2, B2, C2, D2) networks for a range of vX/vX* (panels A1, A2, B1, B2) and gNMDA/gNMDA* (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).

Spiking activity correlation and synchrony computed from spike trains of simulated networks.

Conductance parameters are solutions of mean field equations for the steady state primary network (A1, B1, C1, D1) and the critical state primary network (A2, B2, C2, D2) corresponding to the red and blue asterisks, respectively, in the (IAMPA/IGABA,IX,E/Iθ,E) state plane shown in Figure 2A and inset. For the steady state network, correlation and synchrony are weak and insensitive to the modulation of external input spike rate νX (A1, B1) and NMDAR conductance gNMDA (C1, D1). In contrast, for the critical state network spike correlation depends strongly on the external spike rate (A2) and NMDAR conductance (C2) and the degree of spike synchrony could be modulated from relatively weak to strong (B2, D2). Results shown in (C1, D1, C2, D2) are obtained from simulations in which νX is increased by 5%. The magnitudes of modulation of νX and gNMDA are normalized by their standard values vX* and gNMDA*, respectively. The numbers next to color-coded lines for spike correlation plots show the normalized magnitudes of external input spike rates, vX/vX*, (A1, A2) and NMDAR conductance, gNMDA/gNMDA*, (C1, C2).

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 νX increases, and to desynchronization when the NMDA conductance gNMDA decreases?

To answer these questions and to illuminate the role of asynchronous and synchronous regimes in the shaping of network-wide synchronization of spiking activity, we carried out linear stability analysis in the (vX/vX,gNMDA/gNMDA) 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 (vX/vX,gNMDA/gNMDA) plane in the neighborhood of the steady (Figure 5A) and critical (Figure 5B) state primary networks. As in Figure 2A, the critical line (λ=0) separating the asynchronous stationary (λ<0) and synchronous oscillatory (λ>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 νX and gNMDA 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).

Network state diagrams in th e (vX/vX,gNMDA/gNMDA) plane.

The critical line (λ=0, white line) separates the parameter plane into regions of asynchronous stationary (λ<0) and synchronous oscillatory (λ>0) regimes. In the state diagram for the steady state network (A) the critical line is beyond the area covered by the diagram. Asterisks correspond to the steady (A) and critical (B) state primary networks in these planes. Color-coded arrows show the range of modulation of νX (yellow) and gNMDA (magenta) corresponding to the range of modulation of these parameters for which temporal correlations of spiking activity and synchrony are shown in Figure 4. The insets show how spike synchrony changes along the corresponding arrows in the state diagrams. These insets display the same plots for spike synchrony that are shown in panels B1 and B2 (bottom insets in A and B, correspondingly) and D1 and D2 (right insets in A and B, correspondingly) in Figure 4.

In contrast, the modulations of νX and gNMDA in the critical state primary network (Figure 5B) induce transitions between the network states. Specifically, as the external input spike rate ν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 gNMDA (vertical magenta arrow in Figure 5B) causes the system to cross the boundary again, and the network state changes now from oscillatory to stationary; this transition is accompanied by a sharp decrease in spike synchrony (Figure 4D2 and Figure 5B right inset).

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

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 drug-naive (black) and drug (magenta) conditions. It can be now appreciated that in drug-naive condition, population activity during the pre-response period develops characteristics of synchronized oscillation behavior, as signaled by the appearance of time lagged peaks of correlation (blue arrows, Figure 6B1, black). However, administration of a drug blocking NMDAR desynchronizes neuronal activity during this period (Figure 6B1, magenta).

Figure 6 with 2 supplements see all
Comparison of the effects of blocking of NMDAR in primate PFC and in the prefrontal circuit model.

(A1, B1) Plots show population average temporal correlations between spiking activity of neuron pairs recorded from PFC during the 200 ms period immediately following probe presentation (A1) and the 200 ms period immediately preceding the motor response (B1) in the DPX task (Zick et al., 2018). In the drug-naive condition (black line), population activity during the pre-response period develops characteristics of synchronous oscillation with a frequency of ∼55 Hz (peaks at time lags ±18 ms, blue arrows, B1). Administration of a drug blocking NMDAR (magenta line) desynchronizes neuronal activity during the pre-response period (B1). (A2, B2, C) Temporal correlations (A2, B2) computed from spike trains of simulated networks corresponding to four conditions shown in the (νX/vX,gNMDA/gNMDA) state plane (C) by bold dots and arrow heads: initial probe ( vX/vX*=0.97, A2) and pre-response ( vX/vX*=1.03, B2) periods for drug-naive ( gNMDA/gNMDA=1.25, black line) and drug ( gNMDA/gNMDA*=0, magenta line) conditions. The critical line (λ=0, white line in panel C) separates the parameter plane into regions of asynchronous stationary (λ<0) and synchronous oscillation (λ>0) regimes. The locus of the blue asterisk corresponds to the critical state primary network in this plane.

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

These experimental findings could be readily explained by a prefrontal network model that operates on the boundary between asynchronous and synchronous regimes. We start by recalling that in the framework of our approach the pre-response afferent signals, which we assume are received by PFC neurons before the monkey’s response, are modeled as an increase in the external spike rate from its background level ν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 gNMDA,E and gNMDA,I to zero. The capacity of the prefrontal network model to provide a circuit mechanism for the emergence of synchrony in spiking activity and drug-dependent desynchronization can be illustrated by considering the system’s behavior in the (vX/vX,gNMDA/gNMDA) state plane around the point (vX/vX=1,gNMDA/gNMDA=1) corresponding to the critical state primary network (Figure 6C). In this space, the effects of probe presentation on the spiking dynamics of the prefrontal circuit model under drug-naive (gNMDA/gNMDA=1.25) and drug (gNMDA/gNMDA*=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 (vX/vX=0.97) to the network state corresponding to the pre-response period (vX/vX=1.03).

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

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

Network spiking dynamics in response to transient input.

(A) Temporal profile of external rate. Initially, external rate is fixed at a lower level (vX/vX*=0.97) and the network is in the asynchronous state. At time t=0 the rate begins to increase and in the next 100 ms it crosses the boundary between the asynchronous and synchronous states reaching a higher level (vX/vX*=1.05). The rate is kept constant for the next 400 ms and, afterwards, it decreases within 100 ms and returns to its initial level corresponding to the asynchronous state. (B, C) Average population spike rate (B) and synchrony (C) obtained from spike trains of 100 network simulations that received the transient external input shown in A for drug-naive (gNMDA/gNMDA=1.25, black line) and drug conditions (gNMDA/gNMDA*=0, magenta line). Shaded grey and magenta bands show the standard errors for spike rate (B) and synchrony (C). (D) Population average temporal correlations between spiking activity of neuron pairs obtained in simulations during the 200 ms period shown in (C) by yellow shaded area.

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

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 INMDA*/IGABA*=0.15. We have only shown that the network frequency at the onset of oscillation essentially is independent of the INMDA/IGABA balance (Figure 2—figure supplement 1B), and that the characteristic features of the (IAMPA/IGABA,IX,E/Iθ,E) state diagram qualitatively remain unchanged when this balance is varied (Figure 2—figure supplement 1A). Could, however, the INMDA*/IGABA* 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 (vX/vX,gNMDA/gNMDA) state diagram shown in Figure 6C depend on the INMDA*/IGABA* balance.

Figure 8 shows state diagrams in the (vX/vX,gNMDA/gNMDA) plane obtained for several INMDA*/IGABA* balance values. It is seen that the orientation of the critical line in the state space depends on the INMDA*/IGABA* balance. When the balance is shifted toward stronger inhibition (INMDA/IGABA<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 (INMDA/IGABA>0.15, Figure 8C), the critical line becomes too flat: in the drug-naive condition the external spike modulation may not be able to produce strong enough synchrony because the system would be too close to the critical line, and not shift deep enough into the region of the synchronous regime (black arrow in Figure 8C).

State diagrams in the (vX/vX,gNMDA/gNMDA) plane obtained for several values of the balance between the NMDA and GABA currents.

Notations are the same as in Figure 6C. (A) INMDA*/IGABA*=0.05; (B) INMDA*/IGABA*=0.15; (C) INMDA*/IGABA*=0.25. Note that the critical line orientation depends on the INMDA*/IGABA* balance.

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 λ 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 gAMPA/gAMPA*, gNMDA/gNMDA*, gGABA/gGABA*, and vX/vX* of the synaptic parameters around the critical point {gAMPA,{E,I},gNMDA,{E,I},gGABA,{E,I},vX} corresponding to the onset of oscillatory instability where λ=0 (conductances gR,E and gR,I of excitatory and inhibitory neurons (R=AMPA,NMDA,GABA) are again scaled with the same factors and, thus, their relative changes are equal: gR,E/gR,E*=gR,I/gR,I*). The calculation is detailed in the Materials and methods section. The result is that λ in the vicinity of the critical point on the boundary between the steady and oscillatory states can be approximated by

(1) λ=ΛAMPA(ΔgAMPAgAMPA+ΔϕIsyn,EϕIsyn,E)+ΛNMDA(ΔgNMDAgNMDA+ΔϕIsyn,EϕIsyn,E)ΛGABA(ΔgGABAgGABA+ΔϕIsyn,IϕIsyn,I),

where ΛAMPA, ΛNMDA, and ΛGABA are quantities defined by the parameters of the critical state network around which the equations are linearized, ϕIsyn,α` is the slope of the neuron’s current-frequency response function at the critical state, and ϕIsyn,α` is the change in the slope of the response function due to the deviations of the synaptic parameters from their critical values (α=E,I for excitatory and inhibitory neurons, respectively). The deviations of the synaptic conductances gAMPA, gNMDA, gGABA, and external rate vX give rise to the changes in the corresponding average recurrent IAMPA, INMDA, IGABA and external IX synaptic currents. This produces the change Isyn in the average total current Isyn=IX+IAMPA+INMDA-IGABA and shifts the operating point of the current-frequency response function v=ϕIsyn that describes the relationship between the average total input current Isyn and the output firing frequency of the neuron v. For the leaky integrate-and-fire neuron model, ϕ is a monotonically increasing non-linear function (see, e.g., Renart et al., 2003). Thus, the shift of the operating point of the neuron’s response function ϕ due to the change Isyn in the total average synaptic current results not only in the change of the firing rate (i.e. ϕ), but also in the change of the slope of the response function ϕIsyn`. The latter can be calculated by linearizing the self-consistent mean field equations (see Materials and methods). As a result, ϕIsyn` is approximated as

(2) ΔϕIsyn,αϕIsyn,α=Uα(IXIGABAΔvXvX+IAMPAIGABAΔgAMPAgAMPA+INMDAIGABAΔgNMDAgNMDAΔgGABAgGABA),α=E,I,

where Uα 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 λ given by Equation 1, 2 provides a very good approximation of the exact relationship (see Appendix 1).

Within the linear approximation, the change ϕIsyn` is proportional to the change Isyn :

(3) ΔϕIsyn,α`ϕIsyn,α`=UαIsynI0,

where I0 is a positive constant. Hence, from Equation 2 it follows that Isyn is proportional to the expression in the brackets:

(4) Isyn=I0IXIGABAvXvX*+IAMPAIGABAgAMPAgAMPA*+INMDAIGABAgNMDAgNMDA*-gGABAgGABA*.

From Equation 1 it follows that the rate of oscillatory instability growth λ directly depends on the changes in the synaptic conductances but does not explicitly depend on the external rate variation vX. However, λ depends on vX indirectly via the terms involving the change in the slope ϕIsyn` due to the change in the average total synaptic current Isyn (Equations 3, 4). In fact, Isyn is affected by the variations of the synaptic conductances as well. Thus, the rate of instability growth λ not only directly depends on the synaptic conductances, but also indirectly via the effect of the recurrent excitatory and inhibitory currents mediated by them on the average total synaptic current and, therefore, the operating point of the current-frequency response function.

The factors ΛAMPA, ΛNMDA, and ΛGABA govern the strength of the direct and indirect contributions of the changes in the synaptic conductances gAMPA, gNMDA, and gGABA to the oscillatory instability. By inspecting Equation 1, one can see that the strength of the direct contribution of the change gR (R=AMAPA,NMDA,GABA) is determined only by the corresponding factor ΛR via the term ΛRgR/gR*. However, the strength of its indirect contribution is determined by all three factors, ΛAMPA, ΛNMDA, and ΛGABA, through the changes in the slopes ϕIsyn,E` and ϕIsyn,I`, which depend on gR (Equation 2). For example, the strength of direct contribution to λ due to the change in the GABAR conductance gGABA is determined only by ΛGABA via the term ΛGABAgGABA/gGABA* in Equation 1. However, the strength of indirect contribution from gGABA is determined by all three factors ΛAMPA, ΛNMDA, and ΛGABA via the terms ΛAMPAϕIsyn,E`/ϕIsyn,E`, ΛNMDAϕIsyn,E`/ϕIsyn,E`, and ΛGABAϕIsyn,I`/ϕIsyn,I` in Equation 1 because ϕIsyn,E` and ϕIsyn,I` themselves depend on gGABA (Equation 2). As noted above, this indirect contribution is due to the change in the average total synaptic current and, therefore, the change in the operating point of the current-frequency response function.

Figure 9 illustrates the contributions of individual terms involving ΛAMPA, ΛNMDA, and ΛGABA in Equation 1 to the oscillatory instability growth rate λ. 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 λ is coming from the term involving ΛAMPA. The contribution related to ΛNMDA is nearly zero, whereas the contribution from ΛGABA term is much smaller than the one from ΛAMPA. While both ΛNMDA and ΛGABAΛAMPA, the primary reasons are different (see Appendix 2).

Contributions from various terms in the analytical approximation of the oscillatory instability growth rate λ.

The plots show separately the rate λ and its individual terms ΛAMPA, ΛNMDA, and ΛGABA (Equation 1) as functions of the relative deviations from the critical value of external rate (A), AMPAR conductance (B), NMDAR conductance (C), and GABAR conductance (D). The comparison is performed by varying the underlying parameter while keeping the other parameters at their critical values. Black lines correspond to the rate λ, whereas red, green, and blue lines correspond to the terms involving ΛAMPA, ΛNMDA, and ΛGABA, respectively. In each plot, the values corresponding to red, blue, and green lines add up to the values of black lines. Note that red lines run very close to black lines, and blue and green lines are nearly horizontal. This indicates that the term ΛAMPA alone approximates the dependence of λ on the synaptic parameters rather accurately and that the contributions from the remaining terms ΛNMDA and ΛGABA are rather small.

It should be noted that even though ΛNMDA and ΛGABA 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 ΛNMDA and ΛGABA are small only means that gNMDA and gGABA do not affect the oscillatory instability directly. However, the changes in the NMDAR and GABAR conductances still affect the instability growth rate λ indirectly via the term involving the product of ΛAMPA and ϕIsyn,E`/ϕIsyn,E` in Equation 1, as mentioned above (and summarized below).

Since ΛNMDA,ΛGABAΛAMPA, we can neglect the terms involving ΛNMDA and ΛGABA in Equation 1 for the oscillatory instability growth rate λ. With this approximation, the equation for λ simplifies to

(5) λΛAMPA=ΔgAMPAgAMPA*+ΔϕIsyn,E`ϕIsyn,E`=ΔgAMPAgAMPA*+UEIsynI0.

Inserting the expression for Isyn from Equation 4, we obtain

(6) λΛAMPA=ΔgAMPAgAMPA*+UEIXIGABAvXvX*+IAMPAIGABAgAMPAgAMPA*+INMDAIGABAgNMDAgNMDA*-gGABAgGABA*,

Thus, the instability growth rate λ, 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 vX that affect λ only indirectly through their effect on the operating point of the response function. In addition, λ also depends indirectly on the AMPAR conductance. For the critical state network UE=2.5 and IAMPA/IGABA=0.4. Therefore, half of the contribution to λ 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 λ simplifies to

(7) λΛAMPA=UEIXIGABAvXvX*+INMDAIGABAgNMDAgNMDA*.

The expression in the brackets is proportional to the change in the average total synaptic current Isyn (Equation 4). The transition to synchrony in the model simulations is achieved by increasing external input (drug-naive condition in Zick et al., 2018), whereas reducing the NMDAR conductance prevents the network from such transition (drug condition in Zick et al., 2018). These simulation results and the mechanism implemented in our model for the transition between the steady and oscillatory states, and the lack thereof when the NMDAR conductance is blocked can be explained in terms of Equation 7 for the instability growth rate λ. As explained above, changes in external rate vX and NMDAR conductance gNMDA both affect synchrony via indirect mechanism by changing the excitatory drive Isyn and, therefore, shifting the operating point of the neuron’s response function. In the drug-naive condition, increase in external rate (ΔvX>0) increases the excitatory drive. As a result, λ becomes positive (see Equation 7) and the network switches to the synchronous regime. However, in the drug condition, when NMDAR is blocked (ΔgNMDA<0), the initial excitatory drive is reduced compared to the drug-naive condition, and now the same increase in external rate vX becomes insufficient to offset the reduced excitatory drive caused by the NMDAR blockage. As a result, λ 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 λ 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 0-lag 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) 0-lag 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 0-lag synchronous spiking in neurons.

This circuit mechanism offers a reasonable explanation accounting for the task-locked increase in 0-lag 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 0-lag 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 λ describing the dependence on these parameters near the boundary between the asynchronous and synchronous states where λ=0. We showed that λ, 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, λ depends on the NMDAR, GABAR as well as AMPAR conductances and external rate indirectly via their effect on the operating point of the neuron’s input current-output frequency response function. The direct dependence manifests the essential influence of the AMPAR synaptic conductance on the strength of an excitatory-inhibitory feedback loop via fast excitatory to excitatory and excitatory to inhibitory recurrent connections. The indirect dependence manifests the influence of the synaptic conductances and external rate on the location of the operating point on the current-frequency response curve and, therefore, the slope of the response function. The steepness of the slope, in turn, determines the amplitude of the neuron’s response to dynamically varying input current and, therefore, affects the strength of excitatory feedback.

The analytic expression for the oscillation growth rate λ 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 λ 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 λ 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 λ, increasing AMPAR conductance enhances the excitatory-inhibitory feedback loop leading to the destabilization of the asynchronous activity and emergence of synchronous oscillations. In our model, by contrast, there is no direct effect on the instability growth rate because the AMPAR conductance is kept fixed, and the enhancement of recurrent excitatory feedback is entirely due to the indirect mechanism. It is achieved through external rate increase at a certain strength of the NMDAR conductance resulting in the neuron’s operating point shift toward a steeper slope above the point of the critical network. This induces network oscillation and synchronization of neurons as observed in monkey PFC when NMDAR is not blocked (Zick et al., 2018). However, when the NMDAR conductance is set to zero, the average total synaptic current is reduced, and the operating point moves down to such locus that it cannot be shifted above the point of the critical network by the same increase in external rate. As a result, external rate increase no longer provides a strong enough excitatory feedback, the network remains in asynchronous state, and no increase in synchrony occurs, consistent with observations in Zick et al., 2018 when NMDAR is blocked.

Relation to prior studies of NMDAR function and working memory

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

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

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

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

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 protocol

Male rhesus macaque monkeys (8–10 kg) were trained to perform the dot-pattern expectancy (DPX) task. This task is closely related to the AX-CPT (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 cue-probe sequence had been AX (69% of trials), or to the right if any other cue-probe sequence had been presented (AY, BX, BY, collectively 31% of trials). Since the correct response to the X-probe depended on the preceding cue (A or B), the task required both working memory and cognitive control. Both The DPX and AX-CPT measure specific cognitive control impairments in schizophrenia (Barch et al., 2003; Jones et al., 2010).

Neurophysiological recording

Request a detailed protocol

In 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 0-lag 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 0-lag spike correlation during task performance (Zick et al., 2018). For neurophysiological recording, we used a computer-controlled 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 protocol

We 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 protocol

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

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

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

Network model

Request a detailed protocol

The network consists of N leaky integrate and fire neurons (see, e.g., Dayan and Abbott, 2001), of which NE=0.8N are excitatory and NI=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 CE=pNE connections from excitatory and CI=pNI from inhibitory neurons. In the framework of mean field consideration, the network is large (N1) and connections are sparse (p1) but the average number of connections received by individual neurons, C, is large (C=pN1). In most simulations, networks consisted of N=5103 neurons that were randomly connected with the probability p=0.2 and, therefore each neuron, on average, received C=103 connections. In addition, each neuron also receives CX external connections from excitatory neurons outside of the network that fire spikes independently according to a Poisson process with rate νX.

The dynamics of the membrane potential V(t) of a neuron below the spike firing potential threshold θ obeys the standard leaky integrate and fire equation:

(8) CmdV(t)dt=gm(V(t)VL)isyn(t),

where Cm is the cell membrane capacitance, gm is the membrane leak conductance, VL is the resting potential, and isyn(t) is the total synaptic current. When the membrane potential reaches the threshold θ, the neuron fires a spike, the potential is reset to Vrst, and the neuron becomes insensitive to its input for the duration of a refractory period τrp. Both excitatory and inhibitory neurons have θ=-50 mV, VL=-70 mV, and Vrst=-55 mV. For excitatory neurons Cm=0.5 nF, gm=25 nS, τrp=2 ms, and for inhibitory neurons Cm=0.2 nF, gm=20 nS, τrp=1 ms (see, e.g., Koch, 2004).

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

(9) isyn(t)=iAMPA(t)+iNMDA(t)+iGABA(t)+iX(t),

where iAMPA and iNMDA correspond to recurrent excitatory currents mediated by AMPA and NMDA receptors, respectively, iGABA corresponds to inhibitory currents mediated by GABA receptors, and iX 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:

(10) iAMPA(t)=gAMPA(V(t)VE)jsAMPA,j(t)
(11) iNMDA(t)=gNMDA(V(t)VE)1+[Mg2+]/γexp(βV(t))jsNMDA,j(t)
(12) iGABA(t)=gGABA(V(t)VI)jsGABA,j(t)
(13) iX(t)=gX(V(t)VE)jsX,j(t),

where synaptic reversal potentials VE=0 mV and VI=-70 mV. NMDAR mediated currents have voltage dependence controlled by the extracellular magnesium concentration (Jahr and Stevens, 1990): β=0.062mV1 , γ=3.57 mM, [Mg2+]=1 mM. The gating variable sR,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 {tk}, the temporal dynamics of the gating variable obeys the equations

(14) {τrdx(t)dt=x(t)+τkδ(ttkτl)τdds(t)dt=s(t)+x(t),

where τl, τr and τd are, respectively, latency, rising, and decay time constants. Their values are τAMPA,l=1 mS, τAMAP,r=0.2 ms, τAMPA,d=2 ms for AMPAR mediated currents (Zhou and Hablitz, 1998), τNMDA,l=1 mS, τNMDA,r=2 ms, τNMDA,d=100 ms for NMDAR mediated currents (Hestrin et al., 1990), and τGABA,l=1 mS, τGABA,r=0.5 ms, τGABA,d=5 ms for GABAR-mediated currents (Gupta et al., 2000). The time integral of s(t) in response to a presynaptic spike equals τ* 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 τ* 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 protocol

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

Mean field approximation

Request a detailed protocol

To derive maximal synaptic conductance parameters gX,α, gAMPA,α, gNMDA,α, gGABA,α (α=E,I) providing prescribed neural firing rates νE and ν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 CE excitatory and CI 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 τ*SR0, where R designates the type of the synapse, and

(15) SX0=CXνX,SAMPA0=CEνE,SNMDA0=CEνE,SGABA0=CIνI.

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

(16) (V(t)VE)1+[Mg2+]/γexp(βV(t))V(t)VEκ+β(V(t)V)(VVE)(κ1)κ2,

where κ=1+Mg2+/γexp(-βV). After these simplifications, average components of synaptic currents for excitatory (α=E) and inhibitory (α=I) populations can be written as

(17) IX,α0=gX,αVα-VEτ*SX0=JX,αSX0
(18) IAMPA,α0=gAMPA,αVα-VEτ*SAMPA0=JAMPA,αSAMPA0
(19) INMDA,α0=gNMDA,α/κVα-VEτ*SNMDA0=JNMDA,αSNMDA0
(20) IGABA,α0=gGABA,αVα-VIτ*SGABA0=JGABA,αSGABA0,

where Vα is the average membrane potential, and JR,α 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 NE excitatory and NI inhibitory neurons is reduced to equations describing the dynamics of membrane potentials VE(t) and VI(t) of just two neurons representing, respectively, excitatory, E, and inhibitory, I, populations (Brunel and Wang, 2001; Renart et al., 2003):

(21) ταdVα(t)dt=(Vα(t)VL)+μα+σαταηα(t),α=E,I,

where VL is the resting potential, τα is the effective membrane time constant, μα is the effective mean synaptic input, σα is the magnitude of the fluctuations in the synaptic input, and ηα(t) is the time course of these fluctuations:

(22) τα=Cm,αgm,αSα
(23) Sα=1+TX,ανX+TAMPA,ανE+TNMDA1,α+TNMDA2,ανE+TGABA,ανI
(24) TX,α=gX,αCXτ*gm,α
(25) TAMPA,α=gAMPA,αCEτ*gm,α
(26) TNMDA1,α=gNMDA,αCEτgm,ακ
(27) TNMDA2,α=βgNMDA,αCEτ(VαVE)(κ1)gm,ακ2
(28) TGABA,α=gGABA,αCIτ*gm,α
(29) μα=(TX,ανX+TAMPA,ανE+TNMDA1,ανE)(VEVL)Sα+TNMDA2,ανE(VαVL)+TGABA,ανI(VIVL)Sα.

In the absence of spiking and fluctuations, the average membrane potential would equal μα+VL (Equation 21). The average membrane potential Vα 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)

(30) Vα=μα+VL-θ-Vrstvατα-μα+VL-Vrstvατrp,α.

The total synaptic noise σα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):

(31) σα2=σX,α2+σAMPA,α2+σNMDA,α2+σGABA,α2,

where

(32) σR,α2=JR,α2SR0ταCm,α2,R=X,AMPA,NMDA,GABA.

ηα(t) is a Gaussian process with zero mean, ηα(t)=0, and an exponentially decaying correlation function, ηα(t)ηα(t)exp(|tt|/τsyn,α), which is due to synaptic filtering with effective time constant τsyn,α (Fourcaud and Brunel, 2002):

(33) τsyn,α=σα2σX,α2τAMPA+σAMPA,α2τAMPA+σNMDA,α2τNMDA+σGABA,α2τGABA,

where τAMPA=τAMPA,l+τAMPA,r+τAMPA,d, τNMDA=τNMDA,l+τNMDA,r+τNMDA,d, τGABA=τGABA,l+τGABA,r+τGABA,d are effective synaptic time constants for AMPAR, NMDAR, and GABAR-mediated currents, respectively. In addition, because of sparse connectivity, the correlation of the fluctuations in the synaptic inputs of excitatory and inhibitory populations is neglected: ηE(t)ηI(t)=0. The firing rate να of a neuron, whose potential is governed by Equation 21, is given by a current-frequency relationship ϕαμα,σα that is a function of the mean and fluctuating part of synaptic input (Brunel and Sergi, 1998; Fourcaud and Brunel, 2002):

(34) ϕα(μα,σα)=( τrp,α+ταa(μα,σα)b(μα,σα)dxπexp(x2)(1+erf(x)))1,

where

(35) aμα,σα=Vrst-VL-μασα
(36) b(μα,σα)=θVLμασα(1+0.5τsyn,ατα)+1.03τsyn,ατα0.5τsyn,ατα.

Since μα and σα themselves depend on the population firing rates νE and νI, the two coupled frequency-current equations

(37) {νE=ϕE(μE(νE,νI),σE(νE,νI))νI=ϕI(μI(νE,νI),σI(νE,νI))

provide a self-consistent 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 gX,α, gAMPA,α, gNMDA,α, gGABA,α (α=E,I) and the external spike rate νX are system parameters controlling the regime of network dynamics; they enter to the mean field analysis through expressions for μα, and σα,. If these parameters are given, one can solve the self-consistent equations to obtain predicted by the mean field approximation population firing rates vE0 and vI0 in a stationary state of the network. Conversely, once external νX and population spike rates vE0 and vI0 are specified, the self-consistent equations could be solved to find the values of synaptic conductance parameters gX,α, gAMPA,α, gNMDA,α, gGABA,α (α=E,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 protocol

We 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 IX,α0, IAMPA,α0, INMDA,α0, IGABA,α0 for excitatory (α=E) and inhibitory (α=I) populations is the same, thus providing the following three equations:

(38) INMDA,E0IGABA,E0=INMDA,I0IGABA,I0,IAMPA,E0IGABA,E0=IAMPA,I0IGABA,I0,IX,E0IGABA,E0=IX,I0IGABA,I0.

As a result, whenever the ratio of synaptic conductances and/or component currents is involved, the index α designating the type of the neuron can be dropped.

Two additional equations are obtained by fixing the balance between inhibition and two-component recurrent excitation at certain values:

(39) INMDA0IGABA0=q1,IAMPA0IGABA0=q2

The last constraint is provided in terms of the relative magnitude of average external current of excitatory neurons, IX,E0:

(40) IX,E0Iθ,E0=q3,

where Iθ,E0 is the current that is needed for an excitatory neuron to reach firing threshold θ 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, IAMPA/IGABA, INMDA/IGABA, and IX,E/Iθ,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 νX and fixed values of these three parameters, we are now able to solve the self-consistent equations for the eight synaptic conductances that provide the prescribed population firing rates vE0 and vI0 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 μα are well below threshold θ, firing is driven by the synaptic fluctuations σα 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 self-consistent equations providing the subthreshold regime for μα and, thus, low rate asynchronous network dynamics, arise when inhibition strongly dominates recurrent excitation and the mean external inputs are around or above threshold θ (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: IAMPA/IGABA+INMDA/IGABA<1, and IX,E/Iθ,E1.

Linear stability analysis

Request a detailed protocol

We 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 vE0 and vI0. In order to analyze the resulting network behavior, the mean field approach and self-consistent equations providing population mean firing rates vE0 and vI0 are extended to describe the dynamics of population rates νE(t) and ν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 kδt-tk-τl arriving from the presynaptic cell is replaced by the instantaneous average rate of spikes, CαRναR(t-τl), arriving from all presynaptic cells making the same type of synapse in the postsynaptic neuron:

(41) {τrdx(t)dt=x(t)+CαRναR(tτl)τddSR(t)dt=SR(t)+x(t),

where R designates the type of the synapse (R=X,AMPA,NMDA,GABA), and αR designates the presynaptic population establishing these synapses (αR=X,E for glutamatergic and αR=I for GABAergic synapse). Since external firing rate νX is stationary, the gating variable for external current is constant in time: SX=CXνX . For recurrent currents, the temporal course of SR is dependent on the instantaneous presynaptic population activity ναR(t). Consequently, the total synaptic input current Isyn(t), given as a sum of contributions from external and recurrent components

(42) Isyn(t)=JXSX+JAMPASAMPA(t)+JNMDASNMDA(t)+JGABASGABA(t),

depends on the population firing rates νE(t) and ν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 input-output response function F.

In general, the input-output relationship v(t)=F(Isyn(t)) depends on the spectral characteristics of the input current, resulting in frequency dependent phase shifts and/or amplitude modulations between the oscillatory components of Isyn 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 input-output function F can be approximated by the current-frequency response function ϕ, given by Equation 34, describing the output due to the steady input current.

In the framework of mean field approximation, the output rates ϕE(Isyn,E(t)) and ϕI(Isyn,I(t)) for excitatory and inhibitory populations must be the same as the instantaneous presynaptic population rates νE(t) and νI(t) because both presynaptic and output rates are of the same populations. This requirement results in two self-consistent equations:

(43) {νE(t)=ϕE(Isyn,E(νE(t),νI(t)))νI(t)=ϕI(Isyn,I(νE(t),νI(t))).

Since the amplitudes of firing rate deviations from the rates in asynchronous steady state are small, ϕ(Isyn(t)) can be linearized about the input current Isyn0 in asynchronous state as:

(44) ϕ(Isyn(t))ϕ(Isyn0)+dϕ(Isyn0)dIsyn(Isyn(t)Isyn0).

With this approximation, the self-consistent equations for excitatory and inhibitory populations become

(45) {νE(t)=vE0(1+AEIsyn,E(t)Isyn,E0Isyn,E0)νI(t)=vI0(1+AIIsyn,I(t)Isyn,I0Isyn,I0),

where Aα=Isyn,α0vα0dϕαIsyn,α0dIsyn,α is the dimensionless slope of the current-frequency 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 self-consistent 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 νE(t) and νI(t) in which initially small (with relative amplitudes εE1 and εI1) oscillatory perturbations that can change exponentially with time are added to the stationary rates vE0 and vI0 such that: να(t)=vα0(1+|εα|exp(λt)cos(ωt+φα)) or, equivalently, in complex form

(46) να(t)=vα0(1+εαexp(λt+iωt)),α=E,I,

where λ is the rate of perturbation growth, ω is the oscillation frequency, and εα is complex accounting for a possible shift in oscillation phase φα 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 SR(t) for recurrent currents mediated by R=AMPA, NMDA, GABA receptors:

(47) SR(t)=SR0[1+εαRQR(λ,ω)exp(λt+iωtiΦR(λ,ω))],

where

(48) QR(λ,ω)=exp(λτR,l)((1+λτR,r)2+ω2τR,r2)((1+λτR,d)2+ω2τR,d2)

and

(49) ΦR(λ,ω)=ωτR,l+atan(ωτR,r1+λτR,r)+atan(ωτR,d1+λτR,d).

The components of synaptic currents and the total currents Isyn,E(t) and Isyn,I(t) can now be calculated and inserted into the linearized self-consistent 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

(50) {XAMPA(λ,ω)cos(ΦAMPA(λ,ω))+XNMDA(λ,ω)cos(ΦNMDA(λ,ω))XGABA(λ,ω)cos(ΦGABA(λ,ω))=1XAMPA(λ,ω)sin(ΦAMPA(λ,ω))+XNMDA(λ,ω)sin(ΦNMDA(λ,ω))XGABA(λ,ω)sin(ΦGABA(λ,ω))=0

and the relationship between the relative amplitudes:

(51) εEAI=εIAE,

where

(52) XAMPA(λ,ω)=AEIAMPAIsynQAMPA(λ,ω)
(53) XNMDA(λ,ω)=AEINMDAIsynQNMDA(λ,ω)
(54) XGABA(λ,ω)=AIIGABAIsynQGABA(λ,ω).

Solving Equation 50, we obtain the rate of perturbation growth λ and the oscillation frequency ω. Because both AE and AI are real, the linear relationship between the amplitudes εE and ε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 gAMPA, gNMDA, gGABA, and external rate vX affect synchrony, we linearize the mean field equations Equation 37 and equations Equation 50 for the stability analysis around the point {gAMPA,{E,I}, gNMDA,{E,I}, gGABA,{E,I}, vX } corresponding to the critical state network where λ=0. We then derive an analytical approximation for the oscillation growth rate λ describing its dependence on the synaptic conductances and external rate in the vicinity of this point.

Linearization of mean field equations

Request a detailed protocol

Approximate analytic description of the changes in the population firing rates vE and vI due to small changes in the synaptic conductances and external rate can be obtained by linearizing the current-frequency response function ϕ, providing population firing rates vE and vI as a function of synaptic conductances and external rate. We note that the function ϕ (Equations 34–36) explicitly depends on the mean effective synaptic input μ, synaptic noise σ, membrane time constant τ, and synaptic time constant τsyn, which in turn depend on the synaptic conductances and external rate (Equations 22–33). Thus, changes in the firing rates vE and vI in response to small changes in the synaptic conductances gAMPA, gNMDA, gGABA, and external rate vX can be approximated as:

(55) Δvα=dϕαdμαΔμα+dϕαdσαΔσα+dϕαdταΔτα+dϕαdτsyn,αΔτsyn,α,α=E,I.

The dominant contribution to vα is due to the change in synaptic input, Δμα. Contributions from the remaining terms are relatively small, with the largest contribution being due to the change in the effective membrane time constant, Δτα. Therefore, the expression for vα can be simplified by retaining only the terms involving Δμα and Δτα:

(56) ΔvαdϕαdμαΔμα+dϕαdταΔτα,α=E,I.

Δμα and Δτα are expressed through the relative changes in synaptic conductances gAMPA/gAMPA*, gNMDA/gNMDA*, gGABA/gGABA*, external rate vX/vX*, and the changes in population rates vE and vI:

(57) Δμα=aαEμνE+aαIμvI+bX,αμvXvX*+RbR,αμgRgR*
(58) Δτα=aαEτνE+aαIτvI+bX,ατvXvX*+RbR,ατgRgR*,

where R={AMPA,NMDA,GABA}, and

(59) aαEμ=(TAMPA,α+TNMDA1,α)(VEVL)+TNMDA2,α(V0VL)Sαμα(TAMPA,α+TNMDA1,α+TNMDA2,α)Sα
(60) aαIμ=TGABA,α(VI-μα-VL)Sα
(61) bAMPA,αμ=TAMPA,αVE-μα-VLSανE
(62) bNMDA,αμ=TNMDA1,αVE-μα-VL+TNMDA2,α(V0-μα-VL)SανE
(63) bGABA,αμ=TGABA,αVI-μα-VLSανI
(64) bX,αμ=TX,αVE-μα-VLSανX
(65) aαEτ=-TAMPA,α+TNMDA1,α+TNMDA2,αSατα
(66) aαIτ=-TGABA,αSατα
(67) bAMPA,ατ=-TAMPA,αSατανE
(68) bNMDA,ατ=-TNMDA1,α+TNMDA2,αSατανE
(69) bGABA,ατ=-TGABA,αSατανI
(70) bX,ατ=-TX,αSατανX.

Inserting expressions for Δμα and Δτα 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

(71) v=av+bp,

where

(72) a=[ϕμ,EaEEμ+ϕτ,EaEEτϕμ,EaEIμ+ϕτ,EaEIτϕμ,IaIEμ+ϕτ,IaIEτϕμ,IaIIμ+ϕτ,IaIIτ],Δv=[ΔνEΔvI]
(73) bT=[ϕμ,EbX,Eμ+ϕτ,EbX,Eτϕμ,IbX,Iμ+ϕτ,IbX,Iτϕμ,EbAMPA,Eμ+ϕτ,EbAMPA,Eτϕμ,IbAMPA,Iμ+ϕτ,IbAMPA,Iτϕμ,EbNMDA,Eμ+ϕτ,EbNMDA,Eτϕμ,IbNMDA,Iμ+ϕτ,IbNMDA,Iτϕμ,EbGABA,Eμ+ϕτ,EbGABA,Eτϕμ,IbGABA,Iμ+ϕτ,IbGABA,Iτ],Δp=[ΔvX/vXΔgAMPA/gAMPAΔgNMDA/gNMDAΔgGABA/gGABA].

Here, the elements of matrices a and b are constants defined by the point in the network parameter space around which the mean field equations are linearized. Components of the vector 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 p. Taking into account that ϕτ,α`bR,ατϕμ,α`bR,αμ and that μα+VLVα, we neglect the ϕτ,α`bR,ατ terms in b and replace μα+VL with Vα. With these approximations b simplifies to:

(74) bb0[IX/IGABAIAMPA/IGABAINMDA/IGABA1]T,b0=[ϕμ,EIGABA,E/gm,ESEϕμ,IIGABA,I/gm,ISI].

Equation Equation 71 can now be rewritten as

(75) a-Iv+b0IX/IGABAIAMPA/IGABAINMDA/IGABA-1Tp=0,

where I is the identity matrix. Solving this equation for v we obtain

(76) v=WIX/IGABAIAMPA/IGABAINMDA/IGABA-1Tp,

or in component form

(77) vE=WEIXIGABAvXvX*+IAMPAIGABAgAMPAgAMPA*+INMDAIGABAgNMDAgNMDA*-gGABAgGABA*
(78) vI=WIIXIGABAvXvX*+IAMPAIGABAgAMPAgAMPA*+INMDAIGABAgNMDAgNMDA*-gGABAgGABA*,

where W=WEWIT is given by

(79) W=-a-I-1b0.

In summary, equations Equations 77, 78 describe changes in the excitatory, vE, and inhibitory, vI, population firing rates due to the small relative changes in the synaptic conductances gAMPA/gAMPA*, gNMDA/gNMDA*, gGABA/gGABA*, and external rate vX/vX*.

Linearization of equations for oscillatory instability analysis

Request a detailed protocol

Changes 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 λ and the change in the oscillation frequency ω caused by small changes in the synaptic conductances and external rate. For this purpose, we linearize equations Equation 50 for λ and ω around the point {gAMPA,{E,I},gNMDA,{E,I},gGABA,{E,I},vX} corresponding to the critical state network that is on the boundary between steady and oscillatory states where λ=0. We do this by taking the differentials with respect to the synaptic variables ΦR and XR, (R=AMPA,NMDA,GABA) that, in turn, depend on λ and ω:

(80) {ΔXAMPAcos(ΦAMPA)XAMPAsin(ΦAMPA)ΔΦAMPA+ΔXNMDAcos(ΦNMDA)XNMDAsin(ΦNMDA)ΔΦNMDAΔXGABAcos(ΦGABA)+XGABAsin(ΦGABA)ΔΦGABA=0ΔXAMPAsin(ΦAMPA)+XAMPAcos(ΦAMPA)ΔΦAMPA+ΔXNMDAsin(ΦNMDA)+XNMDAcos(ΦNMDA)ΔΦNMDAΔXGABAsin(ΦGABA)XGABAcos(ΦGABA)ΔΦGABA=0,

The parameter XR (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 λ and ω through QR (Equation 48), XR depends directly on its corresponding synaptic conductance gR and indirectly on all the synaptic conductances and external rate through its dependency on the slope ϕIsyn,αR` of the current-frequency response function. The change ΔXR due to small variations in the synaptic conductances and external rate is given by

(81) ΔXR=XRΔgRgR*+ΔϕIsyn,αR`ϕIsyn,αR`+ΔQRQR,

where ϕIsyn,αR`, XR, QR, and ΦR are constants whose values are defined by the point {gAMPA,{E,I},gNMDA,{E,I}, gGABA,{E,I},vX} in the synaptic parameter space around which the stability analysis equations are linearized.

The relative change ΔQR/QR can be obtained from Equation 48:

(82) ΔQRQR=-τR1Δλ-τR2Δω,

and the change in ΦR from Equation 49:

(83) ΔΦR=τR1Δω-τR2Δλ,

where

(84) τR(1)=τR,l+τR,r1+(ωτR,r)2+τR,d1+(ωτR,d)2
(85) τR2=ωτR,r21+(ωτR,r)2+τR,d21+(ωτR,d)2,

and ω is the oscillation frequency at the critical state. Inserting expressions for ΔQR/QR into equations Equation 81 for ΔXR and, subsequently, expressions for ΔΦR and ΔXR into equations Equation 80, we obtain a system of two linear equations for Δλ and Δω:

(86) {T+Δω+TΔλ=ΔξAMPA+ΔξNMDAΔξGABATΔωT+Δλ=ΔζAMPA+ΔζNMDAΔζGABA,

where

(87) T+=XAMPAτAMPA++XNMDAτNMDA+-XGABAτGABA+
(88) T-=XAMPAτAMPA-+XNMDAτNMDA--XGABAτGABA-
(89) τR+=τR(1)sin(ΦR)+τR(2)cos(ΦR)
(90) τR=τR(1)cos(ΦR)τR(2)sin(ΦR)
(91) ΔξR=XRcos(ΦR)ΔgRgR*+ΔϕIsyn,αR`ϕIsyn,αR`
(92) ΔζR=-XRsin(ΦR)ΔgRgR*+ΔϕIsyn,αR`ϕIsyn,αR`

and R=AMPA,NMDA,GABA. Solving the system of equations Equation 86 for Δλ and Δω we obtain:

(93) Δλ=ΛAMPA(ΔgAMPAgAMPA+ΔϕIsyn,EϕIsyn,E)+ΛNMDA(ΔgNMDAgNMDA+ΔϕIsyn,EϕIsyn,E)ΛGABA(ΔgGABAgGABA+ΔϕIsyn,IϕIsyn,I)
(94) Δω=ΩAMPA(ΔgAMPAgAMPA+ΔϕIsyn,EϕIsyn,E)+ΩNMDA(ΔgNMDAgNMDA+ΔϕIsyn,EϕIsyn,E)ΩGABA(ΔgGABAgGABA+ΔϕIsyn,IϕIsyn,I).

Here, ΛR and ΩR are constants defined by the parameters of the critical state network around which the stability analysis equations are linearized:

(95) ΛR=XRT+2+T-2T+sinΦR+T-cosΦR
(96) ΩR=XRT+2+T-2T+cosΦR-T-sinΦR,

or, equivalently,

(97) ΛR=XRT0cos(ΦR+Φ0)
(98) ΩR=XRT0sin(ΦR+Φ0),

where T0=T+2+T-2 and Φ0=-atan(T+/T-).

Note that while λ and ω 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 ΔϕIsyn,E` and ΔϕIsyn,I` characterizing changes in the slopes of the current-frequency 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

(99) ΔϕIsynϕIsyn=Δ(dϕdμdμdIsyn)dϕdμdμdIsyn=dμdIsynΔ(dϕdμ)+dϕdμΔ(dμdIsyn)dϕdμdμdIsyn=Δϕμϕμ+Δ(dμdIsyn)dμdIsyn.

Taking into account the linear relationship μ~-Isyn/gmS between the effective synaptic input μ and total synaptic current Isyn, we arrive at

(100) ΔϕIsyn,αϕIsyn,α=Δϕμ,αϕμ,αΔSαSα,α=E,I.

As in the case of the change in the current-frequency response function Δϕμ,α, the dominant contribution to the change in the slope of the response function Δϕμ,α` is coming from the change in the synaptic input μα, while the change in the effective membrane time constant τα, similarly, is the next largest contribution. Therefore, Δϕμ,α` can be approximated as

(101) Δϕμ,α`d2ϕαdμα2μα+d2ϕαdταdματα.

Note also that using Equation 22 one can express the relative change ΔSα/Sα through the change τα as ΔSα/Sα=-Δτα/τα. Inserting expressions for Δϕμ,α` and ΔSα/Sα into Equation 100 we obtain

(102) ΔϕIsyn,α`ϕIsyn,α`=ϕμμ,α``ϕμ,α`μα+ϕτμ,α``ϕμ,α`+1τατα.

Equations for μα and τα in terms of the changes in the synaptic conductances, external rate, and the resulting changes in the population firing rates vE and vI have been already derived and are given by Equations 57 and 58. We replace vE and vI 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 μα and τα 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

(103) [ΔϕIsyn,EϕIsyn,EΔϕIsyn,IϕIsyn,I]=aW[IX/IGABAIAMPA/IGABAINMDA/IGABA1]TΔp+bΔp,

where

(104) a~=ϕμμ,E``ϕμ,E`aEEμ+ϕτμ,E``ϕμ,E`+1τEaEEτϕμμ,E``ϕμ,E`aEIμ+ϕτμ,E``ϕμ,E`+1τEaEIτϕμμ,I``ϕμ,I`aIEμ+ϕτμ,I``ϕμ,I`+1τIaIEτϕμμ,I``ϕμ,I`aIIμ+ϕτμ,I``ϕμ,I`+1τIaIIτ
(105) b~T=ϕμμ,E``ϕμ,E`bX,Eμ+ϕτμ,E``ϕμ,E`+1τEbX,Eτϕμμ,I``ϕμ,I`bX,Iμ+ϕτμ,I``ϕμ,I`+1τIbX,Iτϕμμ,E``ϕμ,E`bAMPA,Eμ+ϕτμ,E``ϕμ,E`+1τEbAMPA,Eτϕμμ,I``ϕμ,I`bAMPA,Iμ+ϕτμ,I``ϕμ,I`+1τIbAMPA,Iτϕμμ,E``ϕμ,E`bNMDA,Eμ+ϕτμ,E``ϕμ,E`+1τEbNMDA,Eτϕμμ,I``ϕμ,I`bNMDA,Iμ+ϕτμ,I``ϕμ,I`+1τIbNMDA,Iτϕμμ,E``ϕμ,E`bGABA,Eμ+ϕτμ,E``ϕμ,E`+1τEbGABA,Eτϕμμ,I``ϕμ,I`bGABA,Iμ+ϕτμ,I``ϕμ,I`+1τIbGABA,Iτ.

The elements of matrices a~ and b~ are constants defined by the parameters of the critical state network. Noting that ϕτμ,α``ϕμ,α`+1ταbR,ατϕμμ,α``ϕμ,α`bR,αμ, we neglect the ϕτμ,α``ϕμ,α`+1ταbR,ατ terms in b~ and, as in the calculation of the change in the current-frequency response function, replace μα+VL with Vα. With these approximations b~ simplifies to:

(106) bb0[IX/IGABAIAMPA/IGABAINMDA/IGABA1]T,b0=[ϕμμ,Eϕμ,EIGABA,Egm,ESEϕμμ,Iϕμ,IIGABA,Igm,ISI].

Equation Equation 103 can now be written as

(107) ΔϕIsyn,E`ϕIsyn,E`ΔϕIsyn,I`ϕIsyn,I`=UIX/IGABAIAMPA/IGABAINMDA/IGABA-1Tp,

or in component form

(108) ΔϕIsyn,E`ϕIsyn,E`=UEIXIGABAvXvX*+IAMPAIGABAgAMPAgAMPA*+INMDAIGABAgNMDAgNMDA*-gGABAgGABA*
(109) ΔϕIsyn,I`ϕIsyn,I`=UIIXIGABAvXvX*+IAMPAIGABAgAMPAgAMPA*+INMDAIGABAgNMDAgNMDA*-gGABAgGABA*,

where U=UEUIT is given by

(110) U=a~W+b0~.

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 Δλ and change in the oscillation frequency Δω on the small relative changes in the synaptic conductances and external rate.

Numerical solutions

Request a detailed protocol

Self-consistent mean field equations for the eight conductance parameters, and linear stability equations for the perturbation growth rate λ and the oscillation frequency ω 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 λ 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 λ provides a very good approximation of the exact relationship.

Appendix 1—figure 1
Comparison between exact solution and analytical approximation for the oscillatory instability growth rate λ.

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 λ provides a very good approximation of the exact relationship.

The plots show the rate λ as a function of the relative deviation from the critical value of external rate (A), AMPAR conductance (B), NMDAR conductance (C), and GABAR conductance (D). The comparison is performed by varying the underlying parameter while keeping the other parameters at their critical values. The solid lines show analytical approximations for λ given by Equations 1, 2, and the dotted lines are the exact solutions obtained from the mean field and stability analysis equations.

Appendix 2

Contributions from ΛAMPA, ΛNMDA, and ΛGABA terms to the instability growth rate

The reasons of why ΛNMDA,ΛGABAΛAMPA can be understood by considering analytical expressions of the factors ΛAMPA, ΛNMDA, and ΛGABA. These quantities have a form ΛR=XRcos(ΦR+Φ0)/T0, where XR characterizes the relative attenuation in the strength of the corresponding synapse due to the R-current dynamics, ΦR is the phase lag introduced by R-current synaptic filtering, and T0 and Φ0 are constants (see Materials and methods). For illustration purposes, it is convenient to use a geometric interpretation of ΛR in which it is associated with a vector emanating from the origin of a cartesian plane, with the length XR/T0, and the angle ΦR+Φ0 between the vector and the x-axis. Then, the projection of this vector on the x-axis (i.e. its x-component) is ΛR.

Appendix 2—figure 1 shows vectors LAMPA, LNMDA, and LGABA corresponding to the values of ΛAMPA, ΛNMDA, and ΛGABA in the crttical network. It is seen that the length of the vector LNMDA is nearly zero and, therefore, its projection ΛNMDA is nearly zero as well. The length of LNMDA being very short means that the relative attenuation in the strength of NMDA synapse XNMDA1. This is because the characteristic synaptic time constant τNMDA100 ms of NMDA synapse is much longer than the period of the γ-band oscillation 1/fntwrk15 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 XNMDA numeric value (Equations 48, 53).

The vector LGABA, in contrast, has a length comparable with the length of LAMPA. However, because its angle ΦGABA+Φ0 is close to 90°, the LGABA projection ΛGABA is much smaller than the projection ΛAMPA of the vector LAMPA. This occurs because the phase shift ΦGABA introduced by GABA-current synaptic filtering is lagging the phase shift ΦAMPA introduced by AMPA-current filtering. This difference in the phase delays, in turn, is due to the fact that the GABA synapse (τGABA~5 ms) is slower than the AMPA synapse (τAMPA~2 ms). Thus, while both ΛNMDA and ΛGABA are much smaller than ΛAMPA, the primary reasons are different.

Appendix 2—figure 1
Geometric interpretation of the factors ΛAMPA, ΛNMDA, and ΛGABA.

For illustration purposes, the factors ΛAMPA, ΛNMDA, and ΛGABA (Equation 1) are associated with corresponding vectors LAMPA, LNMDA, and LGABA (see text for details). The projections of these vectors on the x-axis are equal to their respective factors. The vectors LAMPA and LGABA corresponding to the values of ΛAMPA and ΛGABA in the critical network are shown in red and blue, respectively. The green arrowhead at the origin relates to the vector LNMDA corresponding to ΛNMDA. Due to the small size of the vector, it cannot be depicted in scale on this plot.

Appendix 3

Direct and indirect pathways for synaptic conductances and external rate affecting synchrony

Equation 6 for the oscillatory instability growth rate λ 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 λ on the synaptic parameters via their effect on the operating point of the current-frequency response function. Changes in the synaptic conductances and external rate producing an increase in the average total synaptic current (ΔIsyn>0) shift the operating point of excitatory population toward a steeper slope of the current-frequency response function (ΔϕIsyn>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, ΔIsyn>0 results in λ>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 Isyn (UE=0, Equation 5), a change in the operating point Isyn 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 Isyn . As a result, shifts in the operating point will not change the strength of excitatory feedback. In the case of integrate-and-fire neurons, the input current-output frequency relationship is a monotonically increasing non-linear function (UE>0). Consequently, when the operating point shifts due to increasing average input current (ΔIsyn>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 λ 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, ΔgAMPA>0 also results in ΔIsyn>0 (Equation 4). This means that because of the AMPAR conductance increase, both direct and indirect pathways tend to make λ>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 λ:

(111) λΛAMPA=UEIXIGABAvXvX*+INMDAIGABAgNMDAgNMDA*.

The expression in the brackets is proportional to the change Isyn in the average total synaptic current (Equation 4). Without loss of generality, we assume that the drug-naive condition corresponds to the network for which gNMDA=0 and ΔvX=δ<0. Thus, initially ΔIsyn<0 and the response function operating point is below the point of the critical network. As a result, λ<0 and the network is in the asynchronous steady state. When external input increases by 2δ, so that ΔvX=+δ>0, the total synaptic current increases as well, so that Isyn becomes positive and the operating point shifts above the point of the critical network. This results in λ>0 and, therefore, the network transitions to the synchronous oscillatory state. This transition is due to the indirect dependence of λ on the external rate change vX via its effect on the operating point of the current-frequency response function. In the drug condition, the NMDAR conductance is reduced and, if it becomes sufficiently small so that ΔgNMDAgNMDA<δvXIXINMDA, the same increase in external rate will not result in positive λ. 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 vX influences λ indirectly, reduces the average total synaptic current Isyn 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δ increase in the external rate.

When the oscillatory instability growth rate λ on the left-hand 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 (vX/vX,gNMDA/gNMDA) subspace and, thus, Equation 111 becomes the equation for the line

(112) 0=IXvXvX*+INMDAgNMDAgNMDA*.

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 4—figure 1
Comparison between exact and analytical approximation of the critical line separating asynchronous and synchronous states.

The solid line shows the analytical approximation of the critical line given by Equation 112. The dotted line corresponds to the exact solution obtained from the mean field and stability analysis equations. On the scale of this plot, the two curves virtually overlap.

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, Isyn=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 λ (Equation 1). However, even though the background firing rates remain unchanged, the oscillatory instability growth rate λ 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 λ. In this case, it simplifies to

(113) λΛAMPA=ΔgAMPAgAMPA*+UEIAMPAIGABAgAMPAgAMPA*+INMDAIGABAgNMDAgNMDA*.

The expression in the brackets is proportional to the change Isyn in the average total synaptic current (Equation 4). It can be seen that decreasing the NMDAR conductance by ΔgNMDA0<0 and simultaneously increasing the AMPAR conductance in the amount of ΔgAMPA0=ΔgNMDA0INMDAIAMPAgAMPAgNMDA>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 Isyn cancel each other and Isyn=0. Thus, in this case, although the AMPAR and NMDAR conductances change, the operating point of the response function does not. However, because λ also depends directly on the AMPAR conductance (the first term in the r.h.s. of Equation 113) and ΔgAMPA0>0, the instability growth rate λ becomes positive and synchronous oscillations emerge due to the stronger excitatory-inhibitory 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/SCZ-Synaptic-Circuit-Failure-Model (copy archived at Amirikian, 2024).

The following data sets were generated
    1. Chafee MV
    (2024) Mendeley Data
    Neural 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

    1. Abbott LF
    2. van Vreeswijk C
    (1993) Asynchronous states in networks of pulse-coupled oscillators
    Physical Review. E, Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics 48:1483–1490.
    https://doi.org/10.1103/physreve.48.1483
  1. Book
    1. Abeles M
    (1991)
    Corticonics: Neural Circuits of the Cerebral Cortex
    Cambridge: Cambridge University Press.
    1. Barch DM
    2. Carter CS
    3. MacDonald AW
    4. Braver TS
    5. Cohen JD
    (2003)
    Context-processing deficits in schizophrenia: diagnostic specificity, 4-week course, and relationships to clinical symptoms
    Journal of Abnormal Psychology 112:132–143.
  2. Book
    1. Dayan P
    2. Abbott LF
    (2001)
    Theoretical Neuroscience Computational and Mathematical Modeling of Neural Systems
    Cambridge, MA: MIT Press.
  3. Book
    1. Koch C
    (2004)
    Biophysics of Computation: Information Processing in Single Neurons
    Oxford University Press.
  4. Book
    1. Renart A
    2. Brunel N
    3. Wang XJ
    (2003)
    Mean-field theory of recurrent cortical networks: working memory circuits with irregularly Spiking neurons in
    In: Feng J, editors. Computational Neuroscience: A Comprehensive Approach. Boca Raton, FL: Chapman and Hall/ CRC Press. pp. 431–490.

Decision letter

  1. Albert Compte
    Reviewing Editor; IDIBAPS, Spain
  2. John R Huguenard
    Senior Editor; Stanford University School of Medicine, United States
  3. Albert Compte
    Reviewer; IDIBAPS, Spain
  4. Yashar Ahmadian
    Reviewer; University of Cambridge, United Kingdom

Our editorial process produces two outputs: (i) public reviews designed to be posted alongside the preprint for the benefit of readers; (ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Decision letter after peer review:

Thank you for submitting your article "A prefrontal network model operating near steady and oscillatory states links spike desynchronization and synaptic deficits in schizophrenia" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, including Albert Compte as the Reviewing Editor and Reviewer #1, and the evaluation has been overseen by John Huguenard as the Senior Editor. The following individual involved in the review of your submission has agreed to reveal their identity: Yashar Ahmadian (Reviewer #2).

The manuscript is currently not suitable for publication in eLife, but we are willing to consider a revised version that would fully address all concerns of the reviewers. The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Essential revisions:

The reviewers consider that the manuscript has potential but it requires extensive revisions and new simulations to support its claims. Because these are substantial revisions, a revised manuscript will be critically reviewed in depth again and sent back out to the original reviewers.

The interest in the computational results depends critically on the ability of the model to provide a convincing explanation of the experimental data (Zick et al. 2018, 2021). For this, it is essential to:

1) consider not only the dynamics of synchronization (Figure 1) but also the dynamics of firing rate, and not only the period immediately preceding motor action but also the dynamics of the effects following motor action and around the stimulus cue. The mechanisms of the computational model make specific predictions as to how the dynamics of firing rate and oscillatory synchronization will be associated, and this should be specifically validated against the experimental data.

2) consider how the model accommodates a U-shape relationship between NMDAR modulation and network synchrony. Existing evidence shows that 0-lag synchrony in prefrontal networks is affected by manipulations that reduce NMDAR function (Zick et al. 2018) and by manipulations that enhance NMDAR function (Zick et al. 2021). The computational model presented in this manuscript does not show this U-shape behavior and the discussion does not mention this.

In addition, the theoretical understanding of the model should be enhanced by:

3) expanding the theoretical analysis to provide more intuition on the way NMDA works the way it does, and clarification of its distinct role in synchronization (in contrast to AMPA).

4) investigating how synchrony in the model reacts to transient external inputs and how the model can produce weakly oscillatory synchrony.

Finally, all reviewers were concerned about the emphasis on the association with schizophrenia, so we recommend:

5) moving speculative links to schizophrenia to the Discussion section. Also, please broaden the scope of this discussion with available literature on how synchrony and STDP interact in computational models, and with alternative explanations based on changes in NMDAR-dependent synaptic plasticity irrespective of synchronization.

Reviewer #1 (Recommendations for the authors):

Here are some specific recommendations for improvement of the manuscript:

1) The last two paragraphs in the introduction are largely discussional in spirit. The paragraph before last does not make much sense here, at a point where the reader is not yet familiar with the details of the model proposed. The very last paragraph is also mostly speculative in how the results in this manuscript would be important for schizophrenia and would appear more natural in the Discussion section.

2) The data presented in Figure 1 could be presented more fully in order to provide more constraints for the model so as to enhance its biological plausibility. For instance, it would be most interesting to have a parallel plot of the firing rate in these same neurons to see how neuronal activity evolves through the task. Neurons show a response to the cue stimulus (Zick et al. 2018), which may be similar to the response at the time of the response, and that is something that the model would have to explain: how does the control network react with strong synchronization when it receives inputs prior to the response but instead does not show any synchronization when it receives inputs at the time of the cue? Also, it would be good to show panel b extended for a time window after motor response time t=0 to compare the time course of firing rate and synchrony increases. This could also constrain further the network model.

3) In Figure 1 it would be good to have a dashed horizontal line at y=0 to separate correlated (y>0) from anti-correlated (y<0) activity more clearly. I would also recommend changing the x-axis label ("time lag") because this is now used in cross-correlation functions in Figure 4 and they are intrinsically different. The synchrony plot should also include error bars and proper inference about differences between the two conditions at different time points. For clarity, it would be nice to indicate with a shaded area or horizontal bar the period of cue stimulus presentation, too.

4) Briefly mention in line 167 what spiking neuron model you are using

5) Mention of "GABA receptors" in line 173 suggests that the model includes both GABA_A and GABA_B receptors.

6) Verb missing in "and mediated" in line 174.

7) In Figure 3 the color code red/blue means two different things: excitatory vs. inhibitory neurons, and steady vs. critical networks. This is confusing and should be clarified with different color codes for each thing.

Reviewer #2 (Recommendations for the authors):

1) I assume elucidation of the precise mechanism by which external inputs or NMDA affect synchrony can be achieved by analytically inspecting the direct and indirect dependence of the real part of the relevant complex eigenvalue on the synaptic conductances and external inputs, either directly or indirectly via the dependence of linearized neural gains in the background operating point (which itself depends on the external input and the synaptic parameters).

2) The authors do refer to the Compte et al. 2000 study that found that when the total strength of excitatory connectivity is fixed, tilting the relative strength of NMDAR and AMPAR towards the former disfavors synchronous oscillations (with relatively high frequency in or near the γ band). Previous studies (e.g. Wang and Brunel 2003, Tsodyks et al. 1997, and Compte et al. 2000) have pointed out the role of the negative feedback loop between excitatory and inhibitory populations mediated by the fast receptors AMPA and GABA in γ-band oscillations. Again I think the mean-field theory and the stability analysis based on it could be used to elucidate the differences and similarities in the pathways by which NMDAR and AMPAR affect the transition to synchrony. For example, the role of NMDAR could be via the explicit dependence of the growth rate of oscillations (the eigenvalue real part) on g_NMDA, or (most likely) via the effect of NMDAR-based excitation on background rates (which indirectly affect the eigenvalue).

3) Regarding the link to "disconnection" via STDP. Various studies have found that STDP in (initially randomly connected) networks in the asynchronous state can maintain random connectivity (Morrison et al. 2007 – DOI: 10.1162/neco.2007.19.6.1437) or even lead to the emergence of structured connectivity (Izhikevich et al. 2004, Babadi and Abbott 2013, and Litwin-Kumar and Doiron 2014 -- DOIs: 10.1093/cercor/bhh053, 10.1371/journal.pcbi.1002906, 10.1038/ncomms6319). In lieu of evidence (either previous studies or new simulations/modeling by the authors) for this claim, I suggest weakening the claim to saying that the NMDAR-related reduction in synchrony, as a specific disruption in the patterns of precise spike timing, can lead to abnormal changes in patterns of synaptic connectivity via STDP. I also suggest citing the above-mentioned (or similar) studies on previous studies of STDP in asynchronous (or synchronous) states, or any other studies that support specific claims by the authors in this regard.

Other suggestions:

– Figure 6 top: here the cross-correlation plot based on monkey data (as in the model) exhibits nonzero lag "γ" peaks. But these seem to be absent in the data shown in the Neuron paper (e.g. Figure 4D therein). Can the authors comment on this seeming discrepancy? Is this because data only from one of the monkeys (showing stronger effects) was included in the current paper? Can the authors also comment on the flatness of the CCH in the "initial probe period" in the monkey data in contrast to the small 0-lag peak in the model CCH?

– Lines 104-105: I found the meaning of *balance between* in this sentence vague/unclear. What does the balance between NMDA inputs and oscillatory activity mean?

– I suggest the authors also mention (e.g. around lines 476-477 where they discuss their assumption 1, and as further support for this assumption) that in their empirical data (Neuron paper) they do find an increase in PFC neural firing rates in the pre-response period, which will presumably also arise in their model from the increase in external inputs.

– Line 781: "themes" → "times".

– Line 838: I would explain where this equation comes from (I understand it accounts for the effect of reset on mean membrane potential), and also what v_α denotes.

Reviewer #3 (Recommendations for the authors):

Here are issues that I believe the authors should address before the paper can be published:

1. In Figure 1, the authors only show the synchrony before the motor response. I believe it would be instructive to show also what happens after motor response, to visualize how temporally precise is this increase. It would also be useful to show the time course of average firing rates. The explanation of the authors relies on an increase in external input to generate an increase in synchrony – this increase in external input should be reflected in an increase in the average firing rate as well.

2. Figure 6 shows a strong difference between the network model and data – side peaks in spike correlations are much more prominent in the model than in the data. Thus, synchronization in the data is much less oscillatory than in the model. Is there a way in the model to get such weakly oscillatory synchrony? If not, can the authors speculate as to what additional mechanisms would lead to a reduction of γ power, preserving zero lag synchrony?

3. In the data, the increase in synchrony is transient. In the simulations, the author considers the case of stationary external inputs. I believe the authors should also investigate the case of transient increases of external inputs leading to a more realistic time course of synchrony. Such a scenario could also help in solving

the issue mentioned in 2.

[Editors' note: further revisions were suggested prior to acceptance, as described below.]

Thank you for resubmitting your work entitled "A prefrontal network model operating near steady and oscillatory states links spike desynchronization and synaptic deficits in schizophrenia" for further consideration by eLife. Your revised article has been evaluated by John Huguenard (Senior Editor) and a Reviewing Editor.

The manuscript has been improved but there are some remaining issues that need to be addressed, as outlined below:

1. The dynamics of average firing rates in Figure 1 should mark clearly the various relevant periods of the task on the time axis: probe interval, motor response, etc.

2. The new theoretical analyses in this revision have greatly extended the length of the main text, but the additional insight provided is not clearly explained and in its current form it can lead to confusion. It is suggested that most of this material is moved to a supplementary document and a more concise and clearly delineated section in the main text provides the fundamental message of how NMDARs influence network synchronization.

Reviewer #3 (Recommendations for the authors):

The paper has been considerably improved and I would like to congratulate the authors for their hard work. There are still a couple of issues though that I think should be addressed to improve the paper further.

1. In response to one of the main concerns of the reviewers, the authors have expanded Figure 1 to include the dynamics of average firing rates in the data and the effect of NMDAR blockade on this dynamic. This is a welcome addition, but I think this figure could be improved to better understand the time course of the task. Currently, there is no discussion at all about when the motor response is happening and so the left and right panels seem totally disconnected. One has to go back to the original paper by Zick et al. to understand when the motor response happens. I would suggest extending the left panels to include the whole trial (including the `probe' interval), adding where the mean response time is similar to Figure 4 in the Zick et al. 2018 paper. Or at the very least indicate when is the mean response time (time 0 in right panels) with respect to the time axis of the left panels.

2. In response to another concern, the authors have considerably expanded their theoretical analysis with a detailed description of an approximation for the instability growth rate. However, while the approximation itself represents an interesting new derivation, I did not find that it brings much insight, and in fact in some ways introduces additional confusion, for the following reasons. (i) The main new result of this new analysis is the new Equation (1). In this equation, the LambdaRs play critical roles and are the subject of much of the later discussion, however, their meaning is not really explained in the main text and one has to go to the Methods to understand what they really are. (ii) The new Figure (10) indicates that the growth rate is primarily determined by the AMPA term and that the GABA term plays little role. This is quite confusing as oscillations in such networks relies on GABAR mediated negative feedback. How can the authors explain this fact? (iii) After reading the new part, I ended up not having much more insight about mechanisms. My guess at this point is that NMDAR influences the synchronisation properties of the network by providing an additional excitatory drive, thereby playing a similar role as the external excitatory input. I have the feeling that is what the authors are trying to say when they say that the NMDAR term acts through changes in phiprime (the slope of the transfer function), which is exactly how the external inputs act, but it feel like it could be explained in a much clearer way.

https://doi.org/10.7554/eLife.79352.sa1

Author response

Essential revisions:

The reviewers consider that the manuscript has potential but it requires extensive revisions and new simulations to support its claims. Because these are substantial revisions, a revised manuscript will be critically reviewed in depth again and sent back out to the original reviewers.

The interest in the computational results depends critically on the ability of the model to provide a convincing explanation of the experimental data (Zick et al. 2018, 2021). For this, it is essential to:

1) consider not only the dynamics of synchronization (Figure 1) but also the dynamics of firing rate, and not only the period immediately preceding motor action but also the dynamics of the effects following motor action and around the stimulus cue. The mechanisms of the computational model make specific predictions as to how the dynamics of firing rate and oscillatory synchronization will be associated, and this should be specifically validated against the experimental data.

We have extended our analysis of both the simulated and recorded neural data to include modulation of both synchrony and firing rate during the DPX task (in the neural data) and in response to transient inputs to the PFC network (in the simulated data), both in the baseline case, and following reduction of NMDAR synaptic transmission. We added panels C and D to Figure 1 to illustrate modulation of firing rate in the neural data and the influence of NMDAR blockade on neural firing rates. This new analysis shows (1) firing rate in monkey prefrontal cortex ramped gradually up during the delay period considerably before the motor response, (2) spike synchrony increased more abruptly and dramatically immediately before the motor response, and (3) blocking NMDAR attenuated both the increase in spike synchrony and the increase in firing rate (new Figure 1C, D). In addition, we have extended the time of analysis of both synchrony and firing rate to encompass the response as well as the time period following the response. To facilitate the comparison of the dynamics of firing rate and synchrony between experimental data and the computational model we modified the location of the critical network in the parameter space by changing the balance (ratio) between NMDAR and GABAR mediated currents from 0.2 to 0.15. We, correspondingly, repeated all previous simulations for the new value of the underlying parameter and updated the manuscript figures accordingly. These new figures are very similar to the previous ones. Next, we extended our consideration to transient external inputs to the spiking network to emulate time-varying inputs to prefrontal networks that were likely to drive task-locked changes in firing rates (Figure 1C, D), as suggested by Reviewers (see Point 4, below), and provide the results of the simulations in the new Figure 7. Figure 7 demonstrates that modulation of firing rate and spike synchrony in the network simulation exhibit several features that were evident in the recorded neural data. Namely, increases in spike rate were early and gradual following transient increase in network inputs, whereas increases in spike synchrony were delayed and abrupt in relation to transient increase in network inputs in the network simulation (Figure 7B, C; black), as we saw in the neural data (Figure 1C, D; black). In addition, the increase in both firing rate and spike synchrony were attenuated by reducing NMDAR synaptic currents in the model (Figure 7B, C; magenta) as they were in the neural data (Figure 1C, D; magenta). Thus, there was considerable agreement between the simulated and recorded neural data, although the neural data additionally exhibited several unique features, including additional dynamics after the response (Figure 1C, D) that likely reflected complex modulation of synaptic inputs during this time. The new simulation results characterizing the network response to transient inputs both in terms of spike synchrony and spike rate are provided in the Results section titled ‘Explaining the effects of blocking NMDAR observed in primate PFC by the prefrontal circuit model’, starting at Line 367.

2) consider how the model accommodates a U-shape relationship between NMDAR modulation and network synchrony. Existing evidence shows that 0-lag synchrony in prefrontal networks is affected by manipulations that reduce NMDAR function (Zick et al. 2018) and by manipulations that enhance NMDAR function (Zick et al. 2021). The computational model presented in this manuscript does not show this U-shape behavior and the discussion does not mention this.

We are grateful to the Reviewers for raising this possibility, which we consider in a new section we have added to the revised Discussion entitled ‘Potential U-shaped relation between NMDAR function and spike synchrony’ (starting at Line 777). There we note the possibility that either too little or too much NMDAR synaptic activity could impair spike synchrony, in a manner analogous to the U-shaped relationship between D1 receptor activation and prefrontal persistent activity as reported by Vijayraghavan, Arnsten and colleagues (2007), for example. However, existing neurophysiological data do not establish yet that the relationship between spike synchrony and NMDAR function has an inverted U-shape. Specifically, we did not test the dose-response relationship between spike synchrony and NMDAR modifying drugs in the primate drug model, nor establish that NMDAR upregulation was responsible for reduced spike synchrony in the mouse genetic model (Deleting Dgcr8 would alter expression of many downstream proteins). The absence of a U-shaped relationship between NDMAR function and spike synchrony in the model from this perspective does not seem to represent a disagreement with available biological data.

In addition, the theoretical understanding of the model should be enhanced by:

3) expanding the theoretical analysis to provide more intuition on the way NMDA works the way it does, and clarification of its distinct role in synchronization (in contrast to AMPA).

We have now extended the theoretical consideration and obtained an analytic approximation for the oscillatory instability growth rate λ describing the dependence on the AMAPR, NMDAR, GABAR synaptic conductances and external rate. Based on this consideration, we have now provided a substantially more detailed theoretical account of the precise mechanism implemented in our model for the transition between the steady and oscillatory states and the lack thereof when the NMDAR conductance is blocked. Also, the closed form analytic expression for λ allowed us to reveal the differences and similarities in how AMPAR and NMDAR, both of which mediate recurrent excitation, influence the stability of asynchronous state and transition to synchronous oscillations. The results of this new theoretical analysis are presented in the new section “Dependence of oscillatory instability growth rate on synaptic parameters” in the Results and include four new figures (Figures 9-12). The mathematical details are given in the new section “Analytical consideration of the dependence of oscillation growth rate on network parameters” in the Methods.

4) investigating how synchrony in the model reacts to transient external inputs and how the model can produce weakly oscillatory synchrony.

In the revised manuscript we extended our consideration to transient external inputs. To this end, we carried out a new series of simulations in which stationary external input received by network neurons was replaced with a time-varying input that had a trapezoid-like temporal profile. Results of these new simulations are summarized in the new Figure 7 mentioned above. Also, as it was suggested by Reviewer #3, we found that the oscillatory synchrony became slightly weaker, which could be appreciated from relatively smaller off-zero lag peaks in the population average spike correlation function as compared with the case in which external input was stationary. More generally, as one can appreciate from Figure 4 panels A2 and C2 and Figure 5B the strength of oscillatory synchrony depends on how far from the oscillatory instability boundary the synchronous network state is: the closer it is, the weaker oscillatory synchrony is, as it can be appreciated by decreasing off-zero lag peaks while the 0-lag peak remains present.

Finally, all reviewers were concerned about the emphasis on the association with schizophrenia, so we recommend:

5) moving speculative links to schizophrenia to the Discussion section. Also, please broaden the scope of this discussion with available literature on how synchrony and STDP interact in computational models, and with alternative explanations based on changes in NMDAR-dependent synaptic plasticity irrespective of synchronization.

These are good points and we thank the reviewers for raising them. We have made the suggested changes by limiting the link between the current study and schizophrenia in the Introduction to the motivation for the original neurophysiological experiments (as this link dictated the pharmacological and genetic manipulations we employed in the animal models). We have also added a new section to the Discussion with the heading ‘Spike timing disruptions and rewiring of prefrontal local circuits via STDP’ where we discuss the complexity of the interaction between STDP, synchrony, and connectivity in prior modeling studies. Namely, it is difficult to predict whether loss of synchronous spiking would cause disconnection via STDP without additional data. We acknowledge this constraint on our original hypothesis that asynchrony would cause disconnection considering these prior theoretical studies in this new section.

Reviewer #1 (Recommendations for the authors):

Here are some specific recommendations for improvement of the manuscript:

1) The last two paragraphs in the introduction are largely discussional in spirit. The paragraph before last does not make much sense here, at a point where the reader is not yet familiar with the details of the model proposed. The very last paragraph is also mostly speculative in how the results in this manuscript would be important for schizophrenia and would appear more natural in the Discussion section.

These are good points and we thank the reviewers for raising them. We have made the suggested changes by limiting the link between the current study and schizophrenia in the Introduction to the motivation for the original neurophysiological experiments (as this link dictated the pharmacological and genetic manipulations we employed in the animal models). Specifically, we have removed speculative links between the model results we present and synaptic disconnection in schizophrenia. In addition, we have removed the second to last paragraph in the original submission describing the relationship to prior modeling studies. Detailed comparison to prior work is now restricted to the Discussion as recommended.

2) The data presented in Figure 1 could be presented more fully in order to provide more constraints for the model so as to enhance its biological plausibility. For instance, it would be most interesting to have a parallel plot of the firing rate in these same neurons to see how neuronal activity evolves through the task. Neurons show a response to the cue stimulus (Zick et al. 2018), which may be similar to the response at the time of the response, and that is something that the model would have to explain: how does the control network react with strong synchronization when it receives inputs prior to the response but instead does not show any synchronization when it receives inputs at the time of the cue? Also, it would be good to show panel b extended for a time window after motor response time t=0 to compare the time course of firing rate and synchrony increases. This could also constrain further the network model.

We completely agree that the neural data used in Figure 1 should be presented more fully.

1) In the revised Figure 1, we have now added parallel plots displaying the time course of population average firing rate observed in neural data under the two conditions. We have also added a new figure, Figure 7, showing how population spike rate and synchrony obtained from spike trains of simulated networks in response to transient inputs evolve in time and compare the results of simulations with neural data presented in the revised Figure 1.

2) “how does the control network react with strong synchronization when it receives inputs prior to the response but instead does not show any synchronization when it receives inputs at the time of the cue?” This is an important point, that we address now in a new Figure 6—figure supplement 2. First, it can be now appreciated from Figure 1C and D that the population firing rate response to the cue and initial probe periods is weaker than the response at the time of the motor response. The likely reason is that during the cue/probe periods, as it is commonly assumed in modeling studies, only a small group of neurons (typically <10% in network models) receive the external/selective inputs. Second, to understand whether our model is consistent with the experimental observation that during the cue and initial probe periods no increase in the spike synchrony occurs, we carried out a series of simulations in which the fraction of neurons that received the increase in external spike rate varied from 0 to 1. The plot in Figure 6—figure supplement 2 shows spike synchrony obtained from the spike trains of simulated networks as a function of the fraction f. In this approach, the cue and probe presentations correspond to the case when relatively small fraction of neurons (f<0.2) receive the increase in external input, in contrast to the case of the pre-response period, when in the framework of our model all neurons receive increased input (f=1). It is seen that for small f the spike synchrony remains week, which is in agreement with the experimental observations.

3) We have now extended the panel B in Figure 1, as well as the new panel D displaying the time evolution of population firing rate, beyond the motor response time t=0 allowing the comparison of the time courses of synchrony and firing rate observed in neural data. In the revised Results section “Explaining the effects of blocking of NMDAR observed in primate PFC by the prefrontal circuit model” we also carry out a similar comparison for data obtained from simulated networks presented in new Figure 7 mentioned above. In addition, we compare the simulated data with neural data presented in the revised Figure 1.

3) In Figure 1 it would be good to have a dashed horizontal line at y=0 to separate correlated (y>0) from anti-correlated (y<0) activity more clearly. I would also recommend changing the x-axis label ("time lag") because this is now used in cross-correlation functions in Figure 4 and they are intrinsically different. The synchrony plot should also include error bars and proper inference about differences between the two conditions at different time points. For clarity, it would be nice to indicate with a shaded area or horizontal bar the period of cue stimulus presentation, too.

We thank the reviewer for the valuable suggestions.

1) In the revised Figure 1 we have now added a dashed horizontal line at y=0.

2) Thank you for noticing the mislabeled axis: we have now changed the “time lag” label to “time”.

3) We have now included the standard errors depicted as shaded bands in the synchrony plots.

4) We have now marked the time points in Figure 1 when the two conditions are statistically different by asterisks and provide the details in the figure legend.

5) We have now indicated by a yellow shaded area the period of cue stimulus presentation.

4) Briefly mention in line 167 what spiking neuron model you are using

We have added a statement saying that we are using leaky integrate-and-fire neuron model (line 124).

5) Mention of "GABA receptors" in line 173 suggests that the model includes both GABA_A and GABA_B receptors.

We have appended the subscript A to GABA to remove the uncertainty (line 130).

6) Verb missing in "and mediated" in line 174.

We have inserted the missing verb (line 131).

7) In Figure 3 the color code red/blue means two different things: excitatory vs. inhibitory neurons, and steady vs. critical networks. This is confusing and should be clarified with different color codes for each thing.

We agree that the choice of the colors is confusing. We have now changed the color code for the firing rates of excitatory/inhibitory neurons in Figure 3 from red/blue to black/green, respectively.

Reviewer #2 (Recommendations for the authors):

1) I assume elucidation of the precise mechanism by which external inputs or NMDA affect synchrony can be achieved by analytically inspecting the direct and indirect dependence of the real part of the relevant complex eigenvalue on the synaptic conductances and external inputs, either directly or indirectly via the dependence of linearized neural gains in the background operating point (which itself depends on the external input and the synaptic parameters).

Thank you for this valuable suggestion. We have now extended the theoretical consideration and obtained an analytic approximation for the oscillatory instability growth rate λ describing the dependence on the AMAPR, NMDAR, GABAR synaptic conductances and external rate. We have accomplished this by linearizing the mean field and linear stability analysis equations near the boundary between the asynchronous and synchronous states where λ=0 . We have shown that λ essentially directly depends only on the AMPAR conductance. However, λ depends on the NMDAR, GABAR as well as AMPAR conductances and external rate indirectly via their effect on the operating point of the neuron’s input current-output frequency response function. The results of this new theoretical analysis are presented in the new section “Dependence of oscillatory instability growth rate on synaptic parameters” in the Results and include three new figures. The mathematical details are given in the new section “Analytical consideration of the dependence of oscillation growth rate on network parameters” in the Methods.

In the new Results section, we have now provided a more detailed theoretical account of the precise mechanism implemented in our model for the transition between the steady and oscillatory states and the lack thereof when the NMDAR conductance is blocked. We do this in terms of the equation for the instability growth rate λ describing the dependence on the NMDAR conductance and external rate in closed form. We also show that the critical line separating the asynchronous and synchronous states in the synaptic parameter space predicted by the analytical approximation is in excellent agreement with the exact line obtained by solving the mean-filed and linear stability analysis equations numerically. We believe that these and other related results described in the new Results section considerably improved the manuscript, and we are very grateful to the reviewer for giving us the opportunity to provide a deeper analytical insight into the involved synaptic mechanisms.

2) The authors do refer to the Compte et al. 2000 study that found that when the total strength of excitatory connectivity is fixed, tilting the relative strength of NMDAR and AMPAR towards the former disfavors synchronous oscillations (with relatively high frequency in or near the γ band). Previous studies (e.g. Wang and Brunel 2003, Tsodyks et al. 1997, and Compte et al. 2000) have pointed out the role of the negative feedback loop between excitatory and inhibitory populations mediated by the fast receptors AMPA and GABA in γ-band oscillations. Again I think the mean-field theory and the stability analysis based on it could be used to elucidate the differences and similarities in the pathways by which NMDAR and AMPAR affect the transition to synchrony. For example, the role of NMDAR could be via the explicit dependence of the growth rate of oscillations (the eigenvalue real part) on g_NMDA, or (most likely) via the effect of NMDAR-based excitation on background rates (which indirectly affect the eigenvalue).

As implied by the reviewer, the analytic expression for the oscillation growth rate λ allowed us to reveal the differences and similarities in how AMPAR and NMDAR, both of which mediate recurrent excitation, influence the stability of the asynchronous state and the transition to synchronous oscillations. Both AMPAR and NMDAR conductances affect λ indirectly by influencing the amplitude of the neuron’s response to varying input current. However, AMPAR conductance also affects λ directly by influencing the strength of fast excitatory feedback. We have now added the corresponding text in the Discussion section.

In the new Results section mentioned in point (1), we explain that the direct dependence manifests the essential influence of the AMPAR synaptic conductance on the strength of excitatory-inhibitory feedback loops via fast excitatory to excitatory and excitatory to inhibitory recurrent connections. The indirect dependence manifests the influence of the synaptic parameters, including NMDAR and external rate, on the location of the operating point on the neuron’s current-frequency response curve and, therefore, the slope of the response function. The steepness of the slope, in turn, determines the amplitude of the neuron’s response to dynamically varying input current and, therefore, affects the strength of the feedback loop.

Also, in the new Results section, we have provided a theoretical account of the simulation results observed in Compte et al., 2000, when concurrent increase in the AMPAR and decrease in the NMDAR conductances leads to synchronous oscillations. As in the case of our model, the explanation is given in terms of the equation for the instability growth rate λ that describes the direct (AMPAR) and indirect (both AMPAR and NMDAR) synaptic conductance contributions to the emergence of oscillations. In addition, we have also revised the Discussion section where we address the apparent contradiction between Compte et al., 2000, simulations and Zick et al., 2018, experiments, emphasizing the critical role played by the direct and indirect mechanisms by which NMDAR and AMPAR affect the transition from the asynchronous to the synchronous dynamic.

3) Regarding the link to "disconnection" via STDP. Various studies have found that STDP in (initially randomly connected) networks in the asynchronous state can maintain random connectivity (Morrison et al. 2007 – DOI: 10.1162/neco.2007.19.6.1437) or even lead to the emergence of structured connectivity (Izhikevich et al. 2004, Babadi and Abbott 2013, and Litwin-Kumar and Doiron 2014 -- DOIs: 10.1093/cercor/bhh053, 10.1371/journal.pcbi.1002906, 10.1038/ncomms6319). In lieu of evidence (either previous studies or new simulations/modeling by the authors) for this claim, I suggest weakening the claim to saying that the NMDAR-related reduction in synchrony, as a specific disruption in the patterns of precise spike timing, can lead to abnormal changes in patterns of synaptic connectivity via STDP. I also suggest citing the above-mentioned (or similar) studies on previous studies of STDP in asynchronous (or synchronous) states, or any other studies that support specific claims by the authors in this regard.

We thank the reviewer for raising this important point. It motivated us to reexamine our initial hypothesis that asynchrony at the level of spiking neurons would necessarily be expected to disconnect recurrent networks considering prior theoretical studies examining STDP. We consider constraints on our hypothesis imposed by the prior theoretical studies in a new section we have added to the Discussion (‘Spike timing disruptions and rewiring of prefrontal local circuits via STDP’, see new text starting on line 739). There we acknowledge that the relationship between spike timing dynamics and synaptic connectivity in networks incorporating STDP is hard to predict because changes in spike synchrony and synaptic connectivity influence one another in complex ways after STDP is implemented and spiking patterns are allowed to sculpt synaptic connectivity. We thank the reviewer for pointing out this caveat and we now acknowledge it as a constraint on our prior hypothesis in the Discussion. However, to our understanding, prior network stimulations incorporating STDP have obtained divergent results as to how synchrony and connectivity influence one another that appear to depend strongly on the initial assumptions of the models. These crucial factors appear to include delays in signal conduction imposed by axons, synapses, and dendrites, and how network simulations model these delays. For example, delays in signal conduction imposed by the length of axons in recurrent connections between pyramidal neurons can be on the order of tens of ms, e.g., very large in comparison to the STDP window. These delays reflecting anatomical details about cortical network architecture therefore have profound effects on synaptic plasticity via STDP mechanisms (as demonstrated in the Izhikevich et al., 2004 paper). Indeed, how these delays are modeled appears to contribute to divergent results between theoretical studies (comparing the Izhikevich et al.,2004 and Morrison et al., 2007 papers). Synchronous spiking is not required to maintain structure in some STDP networks (Izhikevich 2004) yet correlated synchronous spiking can induce structure in other STDP networks (Litwin-Kumar and Doiron, 2014), according to the specific assumptions adopted. We are grateful to the reviewer for posing these questions, they seem essential. We hope these points and the data we provide will motivate further integrated computational and experimental studies of cortical local circuit dynamics.

Other suggestions:

– Figure 6 top: here the cross-correlation plot based on monkey data (as in the model) exhibits nonzero lag "γ" peaks. But these seem to be absent in the data shown in the Neuron paper (e.g. Figure 4D therein). Can the authors comment on this seeming discrepancy? Is this because data only from one of the monkeys (showing stronger effects) was included in the current paper? Can the authors also comment on the flatness of the CCH in the "initial probe period" in the monkey data in contrast to the small 0-lag peak in the model CCH?

1) The absence of nonzero “γ” lag peaks in Figure 4, Neuron paper, is likely due to the fact that the cross-correlation depicted in Figure 4D is the average of the cross-correlation function over the trial time. Since the strongest spike synchrony is restricted to the end of the trial around the time of the motor response, flanking peaks in the CCH may be washed when averaging over the trial, with only the stronger, central peak remaining. In addition, restriction of the current data to the animal with the strongest spike synchrony may have contributed to this difference.

2) Flatness of the CCH in the "initial probe period" in the monkey data in contrast to the small 0-lag peak in the model CCH. This question is now addressed in Figure 6—figure supplement 1. We note that in a large (network size N1) and sparsely connected network (connection probability p1), when the population average firing rate is stationary and the network is in the asynchronous state, spike synchrony vanishes in the limit N,p0, while the average number of connections of individual neurons C=pN remains finite. The small 0-lag peak seen for the asynchronous network in panel Figure 6 A2 is due to the effect of the finite size of the network (Brunel and Hakim, 1999) in the simulations. The plot in Figure 6—figure supplement 1 shows spike synchrony obtained from the spike trains of the asynchronous steady (red circles) and synchronous oscillatory networks (blue circles) as a function of the simulated network size (from left to right, N= 20,000, 10,000, 5,000, and 3,333). The red and blue lines represent linear fits to the corresponding simulated data. It can be seen that, as the sizes of the networks increase, the synchrony of the asynchronous network extrapolates to zero in the limit N, whereas for the synchronous network it remains finite.

– Lines 104-105: I found the meaning of *balance between* in this sentence vague/unclear. What does the balance between NMDA inputs and oscillatory activity mean?

We agree this was a vague statement and we have removed it from the Introduction in the revision.

– I suggest the authors also mention (e.g. around lines 476-477 where they discuss their assumption 1, and as further support for this assumption) that in their empirical data (Neuron paper) they do find an increase in PFC neural firing rates in the pre-response period, which will presumably also arise in their model from the increase in external inputs.

Thank you for the suggestion. We added a corresponding statement and referred to a new panel D in Figure 1 showing the time course of population firing rate in empirical data (line 344).

– Line 781: "themes" → "times".

Corrected (line 936).

– Line 838: I would explain where this equation comes from (I understand it accounts for the effect of reset on mean membrane potential), and also what v_α denotes.

Thank you for alerting us to this oversight. We now explain that this equation is for the average membrane potential Vα. It was calculated in Renart et al., 2003, using the distribution of potentials obtained in Brunel and Hakim, 1999, and accounts for neuron spiking (i.e., the effect of reset and refractory period) and synaptic noise (line 1000).

Reviewer #3 (Recommendations for the authors):

Here are issues that I believe the authors should address before the paper can be published:

1. In Figure 1, the authors only show the synchrony before the motor response. I believe it would be instructive to show also what happens after motor response, to visualize how temporally precise is this increase. It would also be useful to show the time course of average firing rates. The explanation of the authors relies on an increase in external input to generate an increase in synchrony – this increase in external input should be reflected in an increase in the average firing rate as well.

Thank you for these valuable suggestions. We have now added new panels C, and D to Figure 1 that show the time courses of average firing rates. We also extended panel B as well as the new panel D beyond the time of the motor response. It can be now appreciated from Figure 1B that the peak of the synchrony modulation matches the onset of the motor response. Comparing Figures 1B and 1D one can also appreciate that, consistent with the model’s assumption about the increase in external input, the empirical firing rates gradually increase before the motor response. We have added a corresponding statement (see new text at Lines 110-115) to support our assumption about the increase in external input by referring to the increase in the empirical firing rates shown in the new panel D, Figure 1. Additionally, in a new Figure 7 we have now provided a similar comparison between the time courses of synchrony, population spike rate, and external input for the network model based on spike trains obtained in a series of new simulations in which the network received transient external inputs (see point 3 below).

2. Figure 6 shows a strong difference between the network model and data – side peaks in spike correlations are much more prominent in the model than in the data. Thus, synchronization in the data is much less oscillatory than in the model. Is there a way in the model to get such weakly oscillatory synchrony? If not, can the authors speculate as to what additional mechanisms would lead to a reduction of γ power, preserving zero lag synchrony?

In general, as one can appreciate from Figure 4 panels A2 and C2 and Figure 5B the strength of oscillatory synchrony depends on how far from the oscillatory instability boundary the synchronous network state is: the closer it is, the weaker oscillatory synchrony is, as it can be appreciated by decreasing off-zero lag peaks while the 0-lag peak remains present. This suggests that inputs to prefrontal networks in the neural data push the network just past the state boundary into the domain of oscillatory activity, whereas we consider larger excursions through the synaptic parameter space in our simulations to better characterize the full range of effects.

3. In the data, the increase in synchrony is transient. In the simulations, the author considers the case of stationary external inputs. I believe the authors should also investigate the case of transient increases of external inputs leading to a more realistic time course of synchrony. Such a scenario could also help in solving the issue mentioned in 2.

In the revised manuscript we extended our consideration to transient external inputs. To this end, we carried out a new series of simulations in which stationary external input received by network neurons was replaced with a time-varying input that had a trapezoid-like temporal profile. Results of these new simulations are presented in the Results section “Explaining the effects of blocking of NMDAR observed in primate PFC by the prefrontal circuit model,” and are summarized in the new Figure 7 mentioned above providing the time courses for synchrony and firing rate obtained from spike trains of simulated networks in response to transient external inputs. In the revised Results section, we now compare the dynamics of firing rate and synchrony and association between them observed in simulations with those observed in experiments. Figure 7 demonstrates that modulation of firing rate and spike synchrony in response to a transient input to the network exhibit several features that were evident in the recorded neural data. Namely, increases in spike rate were early and gradual following transient increase in network inputs, whereas increases in spike synchrony were delayed and abrupt in relation to transient increase in network inputs in the network simulation (Figure 7B, C; black), as we saw in the neural data (Figure 1C, D; black). We acknowledge, however, that while our relatively simple model qualitatively is consistent with dynamical features of the firing rate and synchrony observed in primate PFC, there are some quantitative discrepancies in firing rates.

[Editors’ note: what follows is the authors’ response to the second round of review.]

The manuscript has been improved but there are some remaining issues that need to be addressed, as outlined below:

1. The dynamics of average firing rates in Figure 1 should mark clearly the various relevant periods of the task on the time axis: probe interval, motor response, etc.

We expanded the time axis in the left panels of Figure 1 to include the whole trial period. The probe presentation time is now indicated by green shaded area whereas the mean and standard deviation of the motor response time by color coded horizontal bars.

2. The new theoretical analyses in this revision have greatly extended the length of the main text, but the additional insight provided is not clearly explained and in its current form it can lead to confusion. It is suggested that most of this material is moved to a supplementary document and a more concise and clearly delineated section in the main text provides the fundamental message of how NMDARs influence network synchronization.

We have now moved large parts of this material, including three figures, from Results to newly added Appendices 1-5. We have also revised the remaining part in the Results section to improve the clarity of the text.

Reviewer #3 (Recommendations for the authors):

The paper has been considerably improved and I would like to congratulate the authors for their hard work. There are still a couple of issues though that I think should be addressed to improve the paper further.

Thank you for the positive comments and valuable suggestions.

1. In response to one of the main concerns of the reviewers, the authors have expanded Figure 1 to include the dynamics of average firing rates in the data and the effect of NMDAR blockade on this dynamic. This is a welcome addition, but I think this figure could be improved to better understand the time course of the task. Currently, there is no discussion at all about when the motor response is happening and so the left and right panels seem totally disconnected. One has to go back to the original paper by Zick et al. to understand when the motor response happens. I would suggest extending the left panels to include the whole trial (including the `probe' interval), adding where the mean response time is similar to Figure 4 in the Zick et al. 2018 paper. Or at the very least indicate when is the mean response time (time 0 in right panels) with respect to the time axis of the left panels.

As suggested, we expanded the time axis in the left panels to include the whole trial period. The probe presentation time is now indicated by green shaded area whereas the mean and standard deviation of the motor response time by color coded horizontal bars. We also added accompanying explanation in the figure legend (lines 1491-1493).

2. In response to another concern, the authors have considerably expanded their theoretical analysis with a detailed description of an approximation for the instability growth rate. However, while the approximation itself represents an interesting new derivation, I did not find that it brings much insight, and in fact in some ways introduces additional confusion, for the following reasons. (i) The main new result of this new analysis is the new Equation (1). In this equation, the LambdaRs play critical roles and are the subject of much of the later discussion, however, their meaning is not really explained in the main text and one has to go to the Methods to understand what they really are. (ii) The new Figure (10) indicates that the growth rate is primarily determined by the AMPA term and that the GABA term plays little role. This is quite confusing as oscillations in such networks relies on GABAR mediated negative feedback. How can the authors explain this fact? (iii) After reading the new part, I ended up not having much more insight about mechanisms. My guess at this point is that NMDAR influences the synchronisation properties of the network by providing an additional excitatory drive, thereby playing a similar role as the external excitatory input. I have the feeling that is what the authors are trying to say when they say that the NMDAR term acts through changes in phiprime (the slope of the transfer function), which is exactly how the external inputs act, but it feel like it could be explained in a much clearer way.

(i) We thank the reviewer for noting that the meaning of the factors ΛAMPA, ΛNMDA, and ΛGABA have not been clearly explained in the main text. We believe that this fact also caused the confusion expressed in point (ii), which we address below. We have now revised the corresponding section of the manuscript (lines 465-477). The factors ΛAMPA, ΛNMDA, and ΛGABA govern the strength of the direct and indirect contributions of the changes in the synaptic conductances ΔgAMPA, ΔgNMDA, and ΔgGABA to the oscillatory instability. The direct contribution of the conductance change ΔgR is due to the explicit dependence of λ on the ΔgR in Equation 1. The indirect contribution is due to the implicit dependence of λ on the ΔgR in Equation 1 via the changes in the slopes of the current-frequency response functions ΔϕIsyn,E and ΔϕIsyn,I of excitatory and inhibitory neurons, respectively. These changes in the slopes occur because the change in the conductance ΔgR produces changes in the average total synaptic currents of excitatory and inhibitory neurons. These, in turn, change the operating points of the response functions and, hence, the slopes of these functions at the operating points (lines 436-443). We explain that the strength of the direct contribution of the change in the synaptic conductance ΔgR to the instability growth rate λ is determined solely by the corresponding factor ΛR via the term ΛRΔgRgR in Equation 1. However, the strength of the indirect contribution of the ΔgR is determined by all three, ΛAMPA, ΛNMDA, and ΛGABA, factors via the terms involving changes in the slopes ΔϕIsyn,E and ΔϕIsyn,I. In addition, we now provide a specific example for the change in the GABAR conductance ΔgGABA (lines 471-477) and explain the resulting contributions to the instability growth rate λ determined by each of the ΛAMPA, ΛNMDA, and ΛGABA factors.

(ii) Figure 9 (Figure 10 in previous version) demonstrates that contributions to the oscillatory instability growth rate λ from the terms involving ΛNMDA and ΛGABA in Equation 1 are small compared to the contribution involving the term ΛAMPA. Figure 9 also demonstrates that this does not mean that changes in the NMDAR and GABAR conductances do not affect the oscillatory instability. Indeed, one can see from Figure 9 C and D that changes in the instability growth rate λ (black line) in response to the changes in the NMDAR (panel C) and GABAR (panel D) conductances are comparable to the changes in λ in response to the change in the AMPAR conductance (panel B). This is because changes in the NMDAR and GABAR conductances, while do not affect directly, influence the oscillatory instability growth rate λ indirectly by changing the average total synaptic current (excitatory drive) and, therefore, the operating point of the neuron’s current-frequency response function. We provide this explanation in a new paragraph (lines 485-490).

(iii) The reviewer’s understanding of our explanation is correct: NMDAR influences the synchronisation properties of the network by providing an additional excitatory drive, thereby playing a similar role as the external excitatory input. To improve the overall readability and clarity of this section we have now moved some parts of it, including three figures, from the main text to newly added Appendices 1-5. We have also revised the main text to provide a more concise and clear explanation of how the NMDAR conductance and external input modulations affect synchrony (lines 506-524).

https://doi.org/10.7554/eLife.79352.sa2

Article and author information

Author details

  1. David A Crowe

    Department of Biology, Augsburg University, Minneapolis, United States
    Contribution
    Software, Formal analysis, Validation, Investigation, Writing – review and editing
    Competing interests
    No competing interests declared
  2. Andrew Willow

    Department of Biology, Augsburg University, Minneapolis, United States
    Contribution
    Software
    Competing interests
    No competing interests declared
  3. Rachael K Blackman

    1. Department of Neuroscience, University of Minnesota, Minneapolis, United States
    2. Medical Scientist Training Program (MD/PhD), University of Minnesota, Minneapolis, United States
    3. Brain Sciences Center, VA Medical Center, Minneapolis, United States
    Contribution
    Validation, Investigation, Writing – review and editing
    Competing interests
    No competing interests declared
  4. Adele L DeNicola

    1. Department of Neuroscience, University of Minnesota, Minneapolis, United States
    2. Brain Sciences Center, VA Medical Center, Minneapolis, United States
    Contribution
    Investigation
    Competing interests
    No competing interests declared
  5. Matthew V Chafee

    1. Department of Neuroscience, University of Minnesota, Minneapolis, United States
    2. Brain Sciences Center, VA Medical Center, Minneapolis, United States
    3. Center for Cognitive Sciences, University of Minnesota, Minneapolis, United States
    Contribution
    Conceptualization, Funding acquisition, Validation, Investigation, Writing – original draft, Writing – review and editing
    For correspondence
    chafe001@umn.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-9289-0239
  6. Bagrat Amirikian

    1. Department of Neuroscience, University of Minnesota, Minneapolis, United States
    2. Brain Sciences Center, VA Medical Center, Minneapolis, United States
    3. Center for Cognitive Sciences, University of Minnesota, Minneapolis, United States
    Contribution
    Conceptualization, Software, Formal analysis, Supervision, Visualization, Methodology, Writing – original draft, Project administration, Writing – review and editing
    For correspondence
    amiri001@umn.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-8080-0902

Funding

National Institute of Mental Health (R01MH107491)

  • Matthew V Chafee

National Institute of Mental Health (P50MH119569)

  • Matthew V Chafee

National Institute of Mental Health (R25 MH101076)

  • Rachael K Blackman

National Institute of General Medical Sciences (T32 GM008244)

  • Rachael K Blackman

National Institute of General Medical Sciences (T32 HD007151)

  • Rachael K Blackman

University of Minnesota Foundation

  • Bagrat Amirikian

U.S. Department of Veterans Affairs

  • Bagrat Amirikian

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We would like to provide our special thanks to Apostolos P Georgopoulos for his intellectual contributions to the study of the temporal correlations in networks and the temporal dynamics of neural ensembles that provided part of the motivation for this work, as well as for his constant support and encouragement. We thank Dean Evans for lab and project management as well as his assistance with surgeries, animal care, and neural recordings; Dale Boeff for his assistance with neurophysiological recording system design and construction, as well as computer programming for signal processing and data analysis; Sofia Sakellaridi for her assistance with neural recordings; Olivia Newman for preliminary analysis. This work was performed while RKB was employed at the University of Minnesota. The opinions expressed in this article are the author's own and do not reflect the views of the Food and Drug Administration, the Department of Health and Human Services, or the United States Government. Support for this work was provided by the Wilfred Wetzel Graduate Fellowship (to RKB), the National Institute of General Medical Sciences (T32 GM008244 and T32 HD007151 to RKB), and the National Institute of Mental Health (R25 MH101076 to RKB, and R01MH107491 and P50MH119569 to MVC). Partial funding for this study was provided by the University of Minnesota Foundation and the U.S. Department of Veterans Affairs. The sponsors had no role in the current study design, analysis or interpretation, or in the writing of this paper. The contents do not represent the views of the U.S. Department of Veterans Affairs, the National Institutes of Health, the Department of Health and Human Services, or the United States Government.

Ethics

All animal care and experimental procedures conformed to National Institutes of Health (NIH) guidelines and were approved by the Institutional Animal Care and Use Committees of the Minneapolis Veterans Administration Medical Center (VA Animal Component of Research Protocol Number 170601).

Senior Editor

  1. John R Huguenard, Stanford University School of Medicine, United States

Reviewing Editor

  1. Albert Compte, IDIBAPS, Spain

Reviewers

  1. Albert Compte, IDIBAPS, Spain
  2. Yashar Ahmadian, University of Cambridge, United Kingdom

Version history

  1. Received: April 8, 2022
  2. Preprint posted: June 13, 2022 (view preprint)
  3. Accepted: November 24, 2023
  4. Version of Record published: February 6, 2024 (version 1)

Copyright

© 2024, Crowe et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 188
    Page views
  • 20
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

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)

  1. David A Crowe
  2. Andrew Willow
  3. Rachael K Blackman
  4. Adele L DeNicola
  5. Matthew V Chafee
  6. Bagrat Amirikian
(2024)
A prefrontal network model operating near steady and oscillatory states links spike desynchronization and synaptic deficits in schizophrenia
eLife 13:e79352.
https://doi.org/10.7554/eLife.79352

Share this article

https://doi.org/10.7554/eLife.79352

Further reading

    1. Neuroscience
    Eyal Y Kimchi, Anthony Burgos-Robles ... Kay M Tye
    Research Article

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

    1. Neuroscience
    Olgerta Asko, Alejandro Omar Blenkmann ... Anne-Kristin Solbakk
    Research Article Updated

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