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.

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

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

  • 349
    views
  • 35
    downloads
  • 2
    citations

Views, downloads and citations are aggregated across all versions of this paper published by eLife.

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. Immunology and Inflammation
    2. Neuroscience
    Yuedan Wang, Ying Li ... Xuan Xiao
    Research Article

    Acute retinal ischemia and ischemia-reperfusion injury are the primary causes of retinal neural cell death and vision loss in retinal artery occlusion (RAO). The absence of an accurate mouse model for simulating the retinal ischemic process has hindered progress in developing neuroprotective agents for RAO. We developed a unilateral pterygopalatine ophthalmic artery occlusion (UPOAO) mouse model using silicone wire embolization combined with carotid artery ligation. The survival of retinal ganglion cells and visual function were evaluated to determine the duration of ischemia. Immunofluorescence staining, optical coherence tomography, and haematoxylin and eosin staining were utilized to assess changes in major neural cell classes and retinal structure degeneration at two reperfusion durations. Transcriptomics was employed to investigate alterations in the pathological process of UPOAO following ischemia and reperfusion, highlighting transcriptomic differences between UPOAO and other retinal ischemia-reperfusion models. The UPOAO model successfully replicated the acute interruption of retinal blood supply observed in RAO. 60 min of Ischemia led to significant loss of major retinal neural cells and visual function impairment. Notable thinning of the inner retinal layer, especially the ganglion cell layer, was evident post-UPOAO. Temporal transcriptome analysis revealed various pathophysiological processes related to immune cell migration, oxidative stress, and immune inflammation during the non-reperfusion and reperfusion periods. A pronounced increase in microglia within the retina and peripheral leukocytes accessing the retina was observed during reperfusion periods. Comparison of differentially expressed genes (DEGs) between the UPOAO and high intraocular pressure models revealed specific enrichments in lipid and steroid metabolism-related genes in the UPOAO model. The UPOAO model emerges as a novel tool for screening pathogenic genes and promoting further therapeutic research in RAO.

    1. Computational and Systems Biology
    2. Neuroscience
    Megan E Farquhar, Qianqian Yang, Viktor Vegh
    Research Article

    Diffusional kurtosis imaging (DKI) is a methodology for measuring the extent of non-Gaussian diffusion in biological tissue, which has shown great promise in clinical diagnosis, treatment planning, and monitoring of many neurological diseases and disorders. However, robust, fast, and accurate estimation of kurtosis from clinically feasible data acquisitions remains a challenge. In this study, we first outline a new accurate approach of estimating mean kurtosis via the sub-diffusion mathematical framework. Crucially, this extension of the conventional DKI overcomes the limitation on the maximum b-value of the latter. Kurtosis and diffusivity can now be simply computed as functions of the sub-diffusion model parameters. Second, we propose a new fast and robust fitting procedure to estimate the sub-diffusion model parameters using two diffusion times without increasing acquisition time as for the conventional DKI. Third, our sub-diffusion-based kurtosis mapping method is evaluated using both simulations and the Connectome 1.0 human brain data. Exquisite tissue contrast is achieved even when the diffusion encoded data is collected in only minutes. In summary, our findings suggest robust, fast, and accurate estimation of mean kurtosis can be realised within a clinically feasible diffusion-weighted magnetic resonance imaging data acquisition time.