Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC3388156

Formats

Article sections

- Abstract
- 1. Introduction
- 2. Neuromodulation in a spiking network model
- 3. The four-population reduced models
- 4. The two-population reduced model
- 5. A speed-accuracy tradeoff via neuromodulation
- 6. Discussion
- References

Authors

Related links

SIAM J Appl Dyn Syst. Author manuscript; available in PMC 2012 July 3.

Published in final edited form as:

SIAM J Appl Dyn Syst. 2011; 10(1): 148–188.

PMCID: PMC3388156

NIHMSID: NIHMS377931

See other articles in PMC that cite the published article.

Previous models of neuromodulation in cortical circuits have used either physiologically based networks of spiking neurons or simplified gain adjustments in low-dimensional connectionist models. Here we reduce a high-dimensional spiking neuronal network model, first to a four-population mean-field model and then to a two-population model. This provides a realistic implementation of neuromodulation in low-dimensional decision-making models, speeds up simulations by three orders of magnitude, and allows bifurcation and phase-plane analyses of the reduced models that illuminate neuromodulatory mechanisms. As modulation of excitation-inhibition varies, the network can move from unaroused states, through optimal performance to impulsive states, and eventually lose inhibition-driven winner-take-all behavior: all are clear outcomes of the bifurcation structure. We illustrate the value of reduced models by a study of the speed-accuracy tradeoff in decision making. The ability of such models to recreate neuromodulatory dynamics of the spiking network will accelerate the pace of future experiments linking behavioral data to cellular neurophysiology.

Neuromodulation plays an important role in cognitive processes such as decision making [17], but there is currently no direct connection between cellular-level effects of neuromodulators and low-dimensional models of decision dynamics. This paper investigates how neuromodulation at the synaptic level influences performance on a two-alternative, forced-choice (2AFC) decision task. Such tasks are often modeled in cognitive psychology by leaky, competing accumulators [46], each equipped with a decision variable representing the activity of a neural population selective for one of the two stimuli. Choices are signaled when the first of these variables crosses a fixed threshold. In this context, neuromodulation is typically introduced at an abstract level by changing the slope and/or bias point in the functions that characterize how the decision variables respond to stimulus inputs [39, 45, 24, 41, 42].

Such high-level models take the form of stochastic ordinary differential equations (ODEs). These are attractive not only for their analytical and computational tractability but also because, under suitable conditions, their behavior can be approximated by yet simpler one-dimensional Ornstein–Uhlenbeck (OU) and drift-diffusion (DD) processes [25, 8, 18]. As described in [26, 8], DD processes with constant drift rate are optimal in the sense that they render a decision of guaranteed accuracy in the minimum possible time, allowing one to test human and animal behavior against a normative theory, and, specifically, to compute a speed-accuracy tradeoff that optimizes performance [8].

In [19] we began a larger scale computational study of norepinephrine (NE) modulation on decision making. It is known that NE can change cellular excitability and synaptic efficacy, thus altering performance in behavioral tasks [7, 5, 9, 38]. Adapting a spiking neuronal network model from [47], we examined the effects of tonic and phasic NE release when glutamatergic and GABA-ergic synapses are both simultaneously and separately modulated. The biophysical detail in this 2000-neuron model allows comparison with both behavioral and physiological experiments [48, 6], but it contains 9200 ODEs, rendering it computationally expensive and analytically intractable. In order to better understand the underlying mechanisms and to enable faster simulation studies, a low-dimensional reduction that preserves key physiological details is desirable.

Drawing on the studies of [19], here we perform a reduction that includes neuromodulation of synaptic conductances. Spiking networks have previously been reduced to four- and two-population mean-field models for a fixed set of synaptic gains [49]. Here, we show that these models cannot capture extended neuromodulatory effects. In contrast, the models developed in this paper are designed to reproduce behavioral changes over a wide range of synaptic gains. Careful matching of reward rates and analyses of bifurcation diagrams play central roles in our analyses. The extended mean-field version of a working-memory model [13, 15] can capture the effects of neuromodulation and find attractors and bifurcation structures, but it is computationally intensive due to self-consistency calculations. It is therefore difficult to simulate behavior in order to estimate subject reaction times, accuracies, and reward rates under different gain conditions. The present reduced models overcome this problem, speeding up simulations over the spiking model by *O*(10^{3}) for the four-population model and by 3600 for the two-population model, and facilitating bifurcation and nullcline analyses. Moreover, while slices of behavior can be studied and stable states can be determined in the full spiking model, the reduced models allow complete bifurcation studies with computationally determined basins of attraction.

In section 2 we extend the studies of [19] to better map the behavioral effects of different synaptic changes in the original spiking neuronal network, examining task performance metrics and network firing rate behaviors. Section 3 covers the four-population model, describing its derivation, simulated behavioral results, and bifurcation structures. In section 4 we perform the final reduction to a two-population model and compare its behavior to that of the four-population and full spiking models. Bifurcation diagrams for the reduced models illuminate various phenomena of the spiking model, which its high-dimensionality renders opaque. The transition from unarousal, through good performance, to impulsive behavior can be understood from bifurcation diagrams, noise levels, and basins of attraction. The reduced neuromodulation models are then used to implement a speed-accuracy tradeoff, described in section 5, which in principle links behavioral and neurophysiological experiments. The paper closes with a discussion in section 6. Technical details of the reduction process are relegated to appendices.

In this section we extend the computational study of [19], which was motivated by perceptual choice experiments in which subjects are under pressure to respond correctly to as many trials as possible, under a free-response paradigm, in a session of fixed duration. The simplest choices are binary, and such tasks are widely studied. In a visual instantiation, a coherent subset of random dots moves in one of two target directions, which must be correctly signaled to gain a reward. Task difficulty is adjusted by the stimulus strength *E* [0, 1], *E* = 0 indicating complete randomness and *E* = 1 indicating 100% coherence [10, 11, 40].

Direct recordings in monkeys suggest that the lateral intraparietal area (LIP), among others, integrates motion evidence from the middle temporal area (MT) of visual cortex [28, 29, 40, 31], and that its firing rates reach a threshold just prior to response [37]. We apply a similar criterion in the spiking model, thus determining decision times (*DT*) from stimulus onset to threshold crossing, for each trial. The response-to-stimulus interval (*RSI*) imposed by the experimenter is also important in defining performance; it will be studied in section 5. Including nondecision latencies (*NDL*), such as signal processing and motor preparation, the *reward rate* is defined as

$$\mathit{RR}=\frac{\langle \mathit{Acc}\rangle}{\langle \mathit{DT}+\mathit{NDL}+\mathit{RSI}\rangle}:$$

(1)

this is the performance measure used below. In (1) · denote averages over trials and the accuracy *Acc* is the fraction of rewarded trials in the session. If a threshold is crossed before stimulus onset, the trial is defined as an impulsive choice or false alarm, and if no threshold is crossed before 2000 ms, a no-choice trial is declared. No reward is given in either case.

The spiking network of [47, 49, 19] was derived from earlier studies [13] and specifically tuned to represent the MT-LIP system. It contains 2000 leaky integrate-and-fire neurons divided into four groups: two stimulus-selective populations of *N*_{1} = *N*_{2} = 240 pyramidal cells, a nonselective population of *N*_{3} = 1120 pyramidal cells, and a population of *N _{I}* = 400 inhibitory interneurons. The state variables are the (subthreshold) transmembrane voltages

The ODEs describing this leaky integrate-and-fire neuronal model take the form

$${C}_{j}\frac{d{V}_{j}}{\mathit{dt}}=-{g}_{L}({V}_{j}-{V}_{L})+{I}_{\mathit{syn},j}\phantom{\rule{0.1em}{0ex}}(t),$$

(2)

$$\frac{d{S}_{\mathit{type},j}}{\mathit{dt}}=-\frac{{S}_{\mathit{type},j}}{{T}_{\mathit{type}}}+\sum _{l}\delta (t-{t}_{j}^{l}),$$

(3)

where *type* = *AMPA* or *GABA, T _{type}* is the time constant for that synapse type, and the summation of delta functions in (3) is over spikes in the presynaptic neuron. Here

The slower *NMDA* synapses require two ODEs for each *j* [0, 1600] with separate time constants describing the rise and fall of *S _{NMDA,j}*:

$$\frac{d{S}_{\mathit{NMDA},j}\phantom{\rule{0.1em}{0ex}}(t)}{\mathit{dt}}=-\frac{{S}_{\mathit{NMDA},j}\phantom{\rule{0.1em}{0ex}}(t)}{{\tau}_{\mathit{NMDA},\mathit{decay}}}+\alpha {x}_{j}\phantom{\rule{0.1em}{0ex}}(t)\phantom{\rule{0.1em}{0ex}}(1-{S}_{\mathit{NMDA},j}\phantom{\rule{0.1em}{0ex}}(t)),$$

(4)

$$\frac{d{x}_{j}\phantom{\rule{0.1em}{0ex}}(t)}{\mathit{dt}}=-\frac{{x}_{j}\phantom{\rule{0.1em}{0ex}}(t)}{{\tau}_{\mathit{NMDA},\mathit{rise}}}+\sum _{l}\delta (t-{t}_{j}^{l}).$$

(5)

When *V _{j}*(

Spikes from cells external to the local network, with a combined presynaptic firing rate *f _{ext}* = 2.4 kHz due to a 3 Hz background firing rate of approximately 800 external presynaptic connections, as described in [47, 49], and from upstream sensory neurons, are assumed to be filtered by AMPA receptors, and their spike times
${t}_{j}^{l}$, which would normally be taken as independently generated Poisson processes, are approximated by drawing from independent Gaussian distributions with mean and variance

$${\mu}_{S,\mathit{AMPA},\mathit{ext},j}=0.001\phantom{\rule{0.2em}{0ex}}{f}_{\mathit{ext}}\phantom{\rule{0.2em}{0ex}}{\tau}_{\mathit{AMPA}},$$

(6)

$${\sigma}_{S,\mathit{AMPA},\mathit{ext},j}=\sqrt{0.5{\mu}_{S,\mathit{AMPA},\mathit{ext},j}},$$

(7)

derived from the asymptotic values of the conductances. The factor 0.001 in (6) is due to the conversion of *τ* from *ms* to *s*. Thus, the external conductances *S _{AMPA,ext,j}*(

$$d{S}_{\mathit{AMPA},\mathit{ext},j}\phantom{\rule{0.1em}{0ex}}(t)=-({S}_{\mathit{AMPA},\mathit{ext},j}-{\mu}_{S,\mathit{AMPA},\mathit{ext},j})\frac{\mathit{dt}}{{\tau}_{\mathit{AMPA}}}+{\sigma}_{S,\mathit{AMPA},\mathit{ext},j}\sqrt{\frac{2\mathit{dt}}{{\tau}_{\mathit{AMPA}}}}N(0,1),$$

(8)

where *N*(0, 1) is a standard normal distribution. This is the source of variability and noise in the simulation.

The synaptic currents in (2) are given by

$$\begin{array}{ccc}\hfill {I}_{\mathit{syn},k}\phantom{\rule{0.1em}{0ex}}(t)& =& {I}_{\mathit{AMPA},\mathit{ext},k}\phantom{\rule{0.1em}{0ex}}(t)+{I}_{\mathit{AMPA},\mathit{rec},k}\phantom{\rule{0.1em}{0ex}}(t)+{I}_{\mathit{NMDA},\mathit{rec},k}\phantom{\rule{0.1em}{0ex}}(t)+{I}_{\mathit{GABA},\mathit{rec},k}\phantom{\rule{0.1em}{0ex}}(t),\hfill \\ \hfill {I}_{\mathit{AMPA},\mathit{ext},k}\phantom{\rule{0.1em}{0ex}}(t)& =& -{g}_{\mathit{AMPA},\mathit{ext},k}\phantom{\rule{0.1em}{0ex}}({V}_{k}-{V}_{E})\phantom{\rule{0.1em}{0ex}}{S}_{\mathit{AMPA},\mathit{ext},k}\phantom{\rule{0.1em}{0ex}}(t),\hfill \\ \hfill {I}_{\mathit{AMPA},\mathit{rec},k}\phantom{\rule{0.1em}{0ex}}(t)& =& -{g}_{\mathit{AMPA},\mathit{rec},k}\phantom{\rule{0.1em}{0ex}}({V}_{k}-{V}_{E})\sum _{j=1}^{{N}_{E}}{w}_{j,k}\phantom{\rule{0.1em}{0ex}}{S}_{\mathit{AMPA},j}\phantom{\rule{0.1em}{0ex}}(t),\hfill \\ \hfill {I}_{\mathit{NMDA},\mathit{rec},k}\phantom{\rule{0.1em}{0ex}}(t)& =& -\frac{{g}_{\mathit{NMDA},\mathit{rec},k}\phantom{\rule{0.1em}{0ex}}({V}_{k}-{V}_{E})}{(1+[M{g}^{2+}]\phantom{\rule{0.2em}{0ex}}\mathrm{exp}(-0.062{V}_{k}\phantom{\rule{0.1em}{0ex}})/3.57)}\sum _{j=1}^{{N}_{E}}{w}_{j,k}{S}_{\mathit{NMDA},j}\phantom{\rule{0.1em}{0ex}}(t),\hfill \\ \hfill {I}_{\mathit{GABA},\mathit{rec},k}\phantom{\rule{0.1em}{0ex}}(t)& =& -{g}_{\mathit{GABA},\mathit{rec},k}\phantom{\rule{0.1em}{0ex}}({V}_{k}-{V}_{I})\sum _{j=1}^{{N}_{I}}{S}_{\mathit{GABA},j}\phantom{\rule{0.1em}{0ex}}(t).\hfill \end{array}$$

(9)

Here the subscript *rec* (henceforth omitted) indicates recurrent connections within the network, *V _{E}* and

Equations (2)-(9) define a 9200-dimensional dynamical system. Specifically, *k* ranges from 1 to 2000 for voltages *V _{k}* and external inputs

Neuromodulation can be modeled by changing the leak *g _{L}* and synaptic conductances

In [19] we showed that a ratio of GABA-ergic to glutamatergic synaptic modulation of approximately 2:1 is optimal for robustness in that reward rates remain at high levels over the widest range of *γ _{I}*, but only one-dimensional slices of the neuromodulation plane were examined. In Figure 1 we extend those results over the region (

Contour maps of network performance over the glutamatergic (γ_{E}) and GABA-ergic (γ_{I}) gain modulation plane. (*a*) Reward rate RR: ridge of higher rewards has slope ≈ *2* near γ_{E} = γ_{I} = *1.0* that moderates to ≈ **...**

Figure 1(a) shows that no trials are rewarded for *γ _{E}* < 0.65, because low excitation cannot overcome leak conductance, even without inhibition. For

Accuracy (Figure 1(b)) has a similar profile to reward rate, both falling more steeply toward no-choice trials at the upper left than to impulsive ones at lower right (note fast decision times in Figure 1(c)). The optimal reward rate ridge grows wider at higher gains, curving over from a high initial slope to approach a slope of 1. We wish to preserve these key features in our low-dimensional reductions, especially those of Figure 1(a), but, in addition to gross behavioral measures, we also seek to understand the network dynamics that underlie this behavior.

Supplementing Figure 1, Figure 2 delineates regions in the neuromodulation plane in which distinct system behaviors are observed. We describe these moving from top left to bottom right. In the region labeled “below threshold,” inhibition dominates excitation, and, even in the presence of stimulus and noise, a unique stable steady state exists in which the firing rates of both selective populations remain below the decision threshold (here 20 Hz): see the top row of Figure 2(b). Prefiguring the bifurcation analyses in section 3.5, we refer to this as a low-low attractor. As noted above, a certain degree of excitation, primarily from external inputs, is required to overcome leak and allow response. Interneuronal capacitances also differ from those of pyramidal cells, allowing them to spike at lower external currents. Hence as recurrent inhibition *γ _{I}* rises, it takes more excitation

(*a*) Behavior of the spiking model over an expanded neuromodulation plane visualized through a plot of the relative separation between poststimulus firing rates of the two selective populations. In the below threshold region both selective populations **...**

In the peak discriminability region, with stimulus present, the network acquires two additional high-low attractors in which one selective pyramidal population fires above threshold and the other selective and the nonselective populations are suppressed. There are two separate subregions of interest, the first of which includes the high reward rate ridge of Figure 1(a); here both selective populations remain near the low-low attractor with noise present and stimulus off, and both initially ramp up after stimulus onset until one suppresses the other and crosses threshold to approach a high-low attractor: see the second row of Figure 2(b). Moving to the lower right, the basin of attraction of the low-low stimulus-off state shrinks, and the second subregion begins, where more impulsive choices occur due to noise-assisted threshold crossing in the absence of stimulus. Eventually the low-low attractor’s basin shrinks enough that noise rapidly perturbs solutions to one of the high-low attractors: see the third row of Figure 2(b).

Throughout the peak discriminability region firing rates split, endowing one selective population with high activity and the other selective and the nonselective populations with low activity in comparison to the winner, and leaving interneurons with sufficient activity to mediate the pooled inhibition. The winning population’s firing rates remain below 80 Hz, and the nonselective population never wins, even in the impulsive subregion, because its relative recurrent weight factor *ω _{j,k}* = 1 is lower than

The loss of discriminability region marks a transition from splitting to race model behavior [8], in which all three pyramidal populations exhibit increasing firing rates in the presence of stimuli, although one selective population still wins: see the bottom row of Figure 2(b). Moving downward, firing rates of the superthreshold attractors approach 70 Hz at the boundary of the no discriminability, both high, region. Entering this region, pooled inhibition fails to suppress the losing selective population and all pyramidal populations exceed threshold, approaching their maxima at the inverse absolute refractory period 1/*τ _{ref}* = 500 Hz (not shown). Interneuronal firing rates are also high, but pooled inhibition is too weak to suppress any pyramidal population.

In section 3.5 we will investigate the network stability underlying the above behavior by computing bifurcation diagrams for a reduced model.

We now derive a set of four-population models, modifying the mean-field procedure of [49] to incorporate neuromodulation. This yields 11-dimensional systems of stochastic ODEs in which populations become homogeneous units described by average firing rates *ν _{j}*(

Although the synaptic conductances and currents depend on time-varying postsynaptic membrane potentials, we simplify the self-consistency calculations of [13, 35] by employing a fixed average voltage = (*V _{reset}*+

$${J}_{\mathit{type},k}=-{g}_{\mathit{type},k}(\overline{V}-{V}_{\mathit{type}})/1000$$

(10)

(multiplication of *g _{type,k}* in (

Due to all-to-all connectivity, input currents are calculated by summing the presynaptic contributions from each population multiplied by the number of cells in that population:

$${I}_{\mathit{syn},k}\phantom{\rule{0.1em}{0ex}}(t)={I}_{\mathit{NMDA},k}+{I}_{\mathit{AMPA},k}+{I}_{\mathit{GABA},k}+{I}_{\mathit{AMPA},\mathit{ext},k}+{I}_{\mathit{stim},k},$$

(11)

$${I}_{\mathit{AMPA},k}\phantom{\rule{0.1em}{0ex}}(t)=\sum _{j=1}^{3}{N}_{j}\phantom{\rule{0.1em}{0ex}}{J}_{\mathit{AMPA},k}{\omega}_{j,k}\phantom{\rule{0.1em}{0ex}}{S}_{\mathit{AMPA},j}\phantom{\rule{0.1em}{0ex}}(t),$$

(12)

$${I}_{\mathit{NMDA},k}\phantom{\rule{0.1em}{0ex}}(t)=\sum _{j=1}^{3}{N}_{j}\phantom{\rule{0.1em}{0ex}}{J}_{\mathit{NMDA},k}{\omega}_{j,k}\phantom{\rule{0.1em}{0ex}}{S}_{\mathit{NMDA},j}\phantom{\rule{0.1em}{0ex}}(t),$$

(13)

$${I}_{\mathit{GABA},I}\phantom{\rule{0.1em}{0ex}}(t)={N}_{I}\phantom{\rule{0.1em}{0ex}}{J}_{\mathit{GABA},k}\phantom{\rule{0.1em}{0ex}}{S}_{\mathit{GABA},I}\phantom{\rule{0.1em}{0ex}}(t).$$

(14)

Here *ω _{j,k}* is the weight from population

$${I}_{\mathit{AMPA},\mathit{ext},k}={J}_{\mathit{AMPA},\mathit{ext},k}\frac{{T}_{\mathit{AMPA}}}{1000}{N}_{\mathit{ext}}{\nu}_{\mathit{ext}}+{I}_{\mathit{noise},k},$$

(15)

where *N _{ext}* is the number of connections per cell from presynaptic neurons external to the network (here 800),

Following [13, 49], we account for the rise in *S _{NMDA,j}* due to presynaptic spikes by including a term

$${\psi}_{j}=\langle {S}_{\mathit{NMDA},j}\rangle =\frac{0.641{\nu}_{j}{T}_{\mathit{NMDA}}}{1000+0.641{\nu}_{j}{T}_{\mathit{NMDA}}}$$

(16)

[49], which yields *F* (*ψ* (*ν _{j}*)) 0.641

Equipped with these definitions and the currents (11)-(15), we can now state the reduced ODEs. The synaptic variables are determined by

$$\frac{d{S}_{\mathit{NMDA},j}}{\mathit{dt}}=-\frac{{S}_{\mathit{NMDA},j}}{{T}_{\mathit{NMDA}}}+0.641(1-{S}_{\mathit{NMDA},j})\frac{{\nu}_{j}}{1000},$$

(17)

$$\frac{d{S}_{\mathit{AMPA},j}}{\mathit{dt}}=-\frac{{S}_{\mathit{AMPA},j}}{{T}_{\mathit{AMPA}}}+\frac{{\nu}_{j}}{1000},$$

(18)

$$\frac{d{S}_{\mathit{GABA},I}}{\mathit{dt}}=-\frac{{S}_{\mathit{GABA},I}}{{T}_{\mathit{GABA}}}+\frac{{\nu}_{I}}{1000},$$

(19)

and the firing rates obey

$$\frac{d{\nu}_{j}}{\mathit{dt}}=\frac{-({\nu}_{j}-{\varphi}_{j}\phantom{\rule{0.1em}{0ex}}({I}_{\mathit{syn},j}))}{{T}_{\mathit{AMPA}}}.$$

(20)

Here *j* = 1, 2, 3, *I*, and the input-output functions or f-I (firing rate-current) curves *ϕ _{j}*(

We develop a variety of simple analytic expressions for the relationships between input currents and firing rates by following both a semianalytical derivation and an empirical study that better matches the spiking network behavior. For a leaky integrate-and-fire model [14] such as (2), the expected average interspike interval for solutions to travel from reset voltage *V _{reset}* to threshold

$${\varphi}_{\mathit{LIF}}\phantom{\rule{0.1em}{0ex}}({I}_{\mathit{syn}})={\left[{\tau}_{\mathit{ref}}+{\tau}_{m}\sqrt{\pi}{\mathit{\int}}_{\frac{{V}_{\mathit{reset}}-{V}_{\mathit{ss}}}{{\sigma}_{V}}}^{\frac{{V}_{\mathit{thresh}}-{V}_{\mathit{ss}}}{{\sigma}_{V}}}{e}^{{x}^{2}}[1+\mathrm{erf}\phantom{\rule{0.1em}{0ex}}(x)]\phantom{\rule{0.2em}{0ex}}\mathit{dx}\right]}^{-1},$$

(21)

where *V _{ss}* =

At low firing rates *ϕ _{LIF}* (

$${\varphi}_{\mathit{AC}}\phantom{\rule{0.1em}{0ex}}({I}_{\mathit{syn}})\approx \frac{{c}_{p/I}\phantom{\rule{0.1em}{0ex}}({I}_{\mathit{syn}}-{I}_{\mathit{thresh},p/I})}{1-\mathrm{exp}[-{g}_{p/I}\phantom{\rule{0.1em}{0ex}}({c}_{p/I}\phantom{\rule{0.1em}{0ex}}({I}_{\mathit{syn}}-{I}_{\mathit{thresh},p/I}))]}.$$

(22)

This is plotted as the blue curve labeled Abbott–Chance. Here *c _{p/I}* is the asymptotic slope,

$$\begin{array}{l}{\varphi}_{E,a}\phantom{\rule{0.1em}{0ex}}({I}_{\mathit{syn}})\approx {\left[{\tau}_{\mathit{ref}}+\frac{1-\mathrm{exp}[-{g}_{p/I}\phantom{\rule{0.1em}{0ex}}({c}_{p/I}\phantom{\rule{0.1em}{0ex}}({I}_{\mathit{syn}}-{I}_{\mathit{thresh},p/I}))]}{{c}_{p/I}\phantom{\rule{0.1em}{0ex}}({I}_{\mathit{syn}}-{I}_{\mathit{thresh},p/I})}\right]}^{-1}\\ =\frac{{c}_{p/I}\phantom{\rule{0.1em}{0ex}}({I}_{\mathit{syn}}-{I}_{\mathit{thresh},p/I})}{1-\mathrm{exp}[-{g}_{p/I}\phantom{\rule{0.1em}{0ex}}({c}_{p/I}\phantom{\rule{0.1em}{0ex}}({I}_{\mathit{syn}}-{I}_{\mathit{thresh},p/I}))]+{\tau}_{\mathit{ref}}\phantom{\rule{0.2em}{0ex}}{c}_{p/I}\phantom{\rule{0.1em}{0ex}}({I}_{\mathit{syn}}-{I}_{\mathit{thresh},p/I})}.\end{array}$$

(23)

Figures 3(a) and 3(b) show all these f-I curves for both pyramidal cells and interneurons, with the function *ϕ _{E,a}* of (23) in green. With additional modifications to effective synaptic currents

To lower the excessive asymptotic firing rates while retaining currents originally derived from the spiking model, we develop an empirical f-I curve for the pyramidal populations. The full spiking model was simulated under different (*γ _{E}, γ_{I}*) conditions, and the firing rates for each population were averaged over a 500 ms window once the system approached equilibrium. The firing rates of the four populations in the spiking model were then used to compute effective currents to each population using the reduced model currents (11)-(14) calculated above. Comparing these currents to the population firing rates which generated them in the spiking model produces the following function:

$${\varphi}_{E,b}\phantom{\rule{0.1em}{0ex}}({I}_{\mathit{syn},p})={\varphi}_{p,0}+\frac{{c}_{p}\phantom{\rule{0.1em}{0ex}}({I}_{\mathit{syn},p}-{I}_{\mathit{thresh},p})}{1-\mathrm{exp}[-{g}_{p}\phantom{\rule{0.1em}{0ex}}({c}_{p}\phantom{\rule{0.1em}{0ex}}({I}_{\mathit{syn},p}-{I}_{\mathit{thresh},p}))]+{c}_{p}\phantom{\rule{0.1em}{0ex}}({I}_{\mathit{syn},p}-{I}_{\mathit{thresh},p})/{\varphi}_{p,\mathit{max}}},$$

(24)

shown in red in Figure 3(a), compared with empirical results from the spiking model (blue dots: the scatter is due to the fact that firing rates were only averaged over 500 ms windows). While (24) approximates the relationship between the simulated firing rates and the reduced model input currents *I _{syn}* which correspond to the simulated average population firing rates, the function (24) is constrained to have the same initial slope as (23).

For the 2AFC application it is most important to reproduce firing rates up to the decision threshold, i.e., *ϕ*(*I _{syn,p}*) ≤ 20 Hz. Several features of

Contour maps of reward rate over the glutamatergic-GABA-ergic (γ_{E}, γ_{I}) gain modulation plane for the four-population models with f-I curves ϕ_{E,a}(I_{syn,j}) (*a*), and ϕ_{E,b}(I_{syn,p}) and ϕ_{I,b}(I_{syn,I}) (*b*). Peak synaptic **...**

The f-I curve for the interneuronal population must rise sufficiently steeply for inhibition to produce splitting behavior and correctly reproduce the no discriminability region of Figure 2. A piecewise-linear function *ϕ _{I,b}*(

$${\varphi}_{I,b}\phantom{\rule{0.1em}{0ex}}({I}_{\mathit{syn},I})=3+600\phantom{\rule{0.1em}{0ex}}\theta ({I}_{\mathit{syn},I}-{I}_{\mathit{thresh},I}),$$

(25)

in which *θ* is a Heaviside function. The four-population models differ in their use of different f-I curves.

Finally, the two-population model of section 4 is simplified by replacing the trivially linear f-I curve for the nonselective population of [49] by a piecewise-linear function with the same current threshold and initial slope as the selective populations, which will continue to use *ϕ _{E,b}*:

$${\varphi}_{\mathit{TL},3}\phantom{\rule{0.1em}{0ex}}({I}_{\mathit{syn},3})=1+{c}_{p}\phantom{\rule{0.1em}{0ex}}\theta ({I}_{\mathit{syn},3}-{I}_{\mathit{thresh},p}),$$

(26)

but this is not necessary in the four-population models, in which all three pyramidal populations share identical f-I curves.

In typical spiking models each cell receives independent Poisson spike trains filtered through external or recurrent AMPA receptors, but in the present reduced models (as in that of [49]), additive Gaussian noise is superimposed on the external input current to each population. In Appendix B we derive an explicit description of these inputs by converting a Poisson spike train of mean frequency *f* to fluctuations in a single synaptic variable *S* about its average value *f τ*, where *τ* is the synaptic decay time constant. (See (43): the derivation is done for the distribution of fluctuations in the quantity *gS*, but the standard deviation is found to scale linearly with *g*, so it may be factored out as described in Appendix B.) The distribution of interspike fluctuations is not normal (see Figure 16(b)) but has zero mean and variance

$$\mathrm{Var}\phantom{\rule{0.1em}{0ex}}(\mathrm{\Delta}S)=\frac{\mathit{f\; \tau}}{\mathit{f\; \tau}+2}.$$

(27)

(*a*) Each incoming spike, with interspike intervals exponentially distributed with average *1*/f, raises S(t) by g, after which it decays with time constant τ. (*b*) Probability distributions of ξ (left) and η (right); note flipped **...**

We then appeal to the Central Limit Theorem to argue that the summed effect of single spikes in single neurons (of *N* total) on the population-averaged synaptic variable is approximately Gaussian with zero mean and variance

$$\frac{\mathrm{Var}\phantom{\rule{0.1em}{0ex}}(\eta )}{{N}^{2}}=\frac{\mathit{f\; \tau}}{{N}^{2}(\mathit{f\; \tau}+2)}.$$

(28)

We can now calculate *I _{noise,k}* for inclusion in (15), which expresses

$$d{I}_{\mathit{noise},k}=-{I}_{\mathit{noise},k}\frac{\mathit{dt}}{{T}_{\mathit{AMPA}}}+{J}_{\mathit{AMPA},\mathit{ext},k}\sqrt{\frac{{f}^{2}\tau}{N(\mathit{f\; \tau}+2)}}\mathit{dW}.$$

(29)

Equipped with the f-I curves and Gaussian noise approximation described above, we can now implement neuromodulation in the four-population model by scaling all glutamatergic synaptic currents by *γ _{E}* and all GABA-ergic currents by

We next consider the model with f-I curve *ϕ _{E,a}*(

The match in reward rate contours can be further improved by using the function *ϕ _{E,b}*(

To better understand the behaviors seen in Figures 1 and and2,2, we now study bifurcations of the four-population models as the average stimulus current *μ*_{0} and coherence *E* vary. Computations are done without noise, using the software XPPAUT [20], and to reveal the system’s full structure we explore *μ*_{0} values well beyond the physiologically realistic range, including “negative stimuli” (*μ*_{0} < 0). Figures 5--88 show branches of equilibria for characteristic points on the neuromodulation plane in terms of the synaptic variables *S*_{NMDA,1} and *S*_{NMDA,2}, written here and henceforth as *S*_{1} and *S*_{2} for brevity.

Bifurcation diagrams for the four-population model with f-I curves ϕ_{E,a}(I_{syn,j}) for E = *0* at (γ_{E}, γ_{I}) = (*1, 1*) (*a*) and (*2, 2*) (*b*); cf. Figure 4(*a*). Stable branches shown bold; branches for S_{2} are identical due to symmetry.

Bifurcation diagrams for the four-population model with f-I curves ϕ_{E,b}(I_{syn,p}) and ϕ_{I} (I_{syn,I}) and E = *0.128* at (γ_{E}, γ_{I}) = (*1, 1*); cf. Figure 4(*b*). (*a*) S_{1} versus μ_{0}; (*b*) S_{2} versus μ_{0}. Stable branches are **...**

We start with the f-I curve *ϕ _{E,a}* of Figure 4(a), with rescaled glutamatergic currents as described in section 3.4, and at first set

The bifurcation diagram for (*γ _{E}, γ_{I}*) = (2, 2) is similar to that for (1, 1), but the two pitchforks move to

These results help explain the network dynamics responsible for the behavior mapped in Figure 1. Prior to the start of a trial, *μ*_{0} = 0 and the state remains near the low-low attractor, if it exists. In both panels of Figure 5 the pitchfork indicating loss of stability of this attractor lies in the range 20–35 Hz, so that, when *μ*_{0} jumps to 40 Hz at stimulus onset, the system state diverges toward a choice attractor. In [49], the third high-high attractor does not appear until much larger *μ*_{0}, creating a pure decision dynamic in which the high-low attractors are the only stable states available. Due to the changes in synaptic currents made to match the region of good performance in the present model, this high-high state is also stable for all *μ*_{0} (0, 40) Hz, but its basin of attraction does not include the physiologically relevant region in which both populations lie below threshold at stimulus onset, so the decision dynamics and effective tristability remain intact. This is illustrated by phase portraits of the two-population system in section 4.2.3.

The bifurcation diagrams also explain false alarm rates seen in the behavioral simulations. At *μ* = 0 noise can take the place of stimulus by carrying the state outside the low-low attractor’s basin. After crossing its boundary, the state moves toward a decision attractor, causing more impulsive choices. This is representative of the upper range of good performance in the peak discriminability region of Figure 2.

The diagrams of Figure 5 may seem confusing, because in the right panel, for (*γ _{E}, γ_{I}*) = (2, 2), the loss of stability occurs at a higher value of

For (*γ _{E}, γ_{I}*) = (0.5, 1) and (2.5, 0.25), in the blue and red regions of Figure 2, the system has unique monotonic branches of fixed points at low and high

We now turn to an asymmetric case. Figures 6(a) and 6(b) show fixed point branches for *S*_{1} and *S*_{2} at (*γ _{E}, γ_{I}*) = (1, 1) and

Bifurcation diagrams for the four-population model with f-I curves ϕ_{E,a}(I_{syn,j}) and E = *0.128* at (γ_{E}, γ_{I}) = (*1, 1*); cf. Figure 4(*a*). (*a*) S_{1} versus μ_{0}; (*b*) S_{2} versus μ_{0}. Stable branches are shown in bold.

We now consider the second four-population model that uses *ϕ _{E,b}* of (24) for pyramidal cells and the threshold-linear f-I curve (25) for interneurons, again focusing on behavior in the peak discriminability region of Figure 2(a). This model recreates the neuromodulation plane even better than that with

Bifurcation diagrams for the four-population model with f-I curves ϕ_{E,b}(I_{syn,p}) and ϕ_{I} (I_{syn,I}) and E = *0* at (γ_{E}, γ_{I}) = (*1, 1) (a*) and (γ_{E}, γ_{I}) = (*2, 2) (b*); cf. Figure 4(*b*). Stable branches are shown in **...**

We again start with symmetric cases. Comparing Figures 7(a) and 5(a), we see that the branch structures at (*γ _{E}, γ_{I}*) = (1, 1) are similar but that the higher

Figure 8 shows an asymmetric case at (*γ _{E}, γ_{I}*) = (1, 1) with

Figure 9 shows reaction time (*RT*) distributions for the asymmetric (*E* = 0.128) condition for both four-population models at (*γ _{E}, γ_{I}*) = (1, 1) and (2, 2). The reaction times measured in behavioral experiments are sums of decision times

The bifurcation studies of section 3.5 show that the four-population, 11-dimensional system (17)-(20) preserves key behaviors of the spiking model, but it is still too complex for full analytical study. Moreover, we wish to relate the original spiking network to one- and two-dimensional OU and leaky accumulator models [46, 8]. We therefore perform a further reduction to a two-population model, closely following the methods of [49], but relaxing assumptions made in that work in order to accommodate the wider range of synaptic currents due to NE modulation. This requires retaining two additional dynamic equations, and our model is therefore four- instead of two-dimensional as in [49]. The two-population model runs faster than the four-population models by a factor of four and allows analysis in a lower-dimensional space.

Reduction to two populations is based on separation of time scales [30, 27]. The synaptic time constants for *AMPA* and *GABA* are fast (*T _{AMPA}* = 2 ms;

We next recognize that firing rates also rapidly approach the values set by the f-I curves in (20) since they too are dominated by *T _{AMPA}*. We may therefore set

$$\frac{d{S}_{1}}{\mathit{dt}}=-\frac{{S}_{1}}{{T}_{\mathit{NMDA}}}+0.641(1-{S}_{1})\frac{{\nu}_{1}}{1000},$$

(30)

$$\frac{d{S}_{2}}{\mathit{dt}}=-\frac{{S}_{2}}{{T}_{\mathit{NMDA}}}+0.641(1-{S}_{2})\frac{{\nu}_{2}}{1000},$$

(31)

$$\frac{d{\nu}_{1}}{\mathit{dt}}=-\frac{{\nu}_{1}-{\varphi}_{E,b}({I}_{\mathit{syn},1})}{{T}_{\mathit{pop}2}},$$

(32)

$$\frac{d{\nu}_{2}}{\mathit{dt}}=-\frac{{\nu}_{2}-{\varphi}_{E,b}({I}_{\mathit{syn},2})}{{T}_{\mathit{pop}2}}$$

(33)

(here, as in section 3.5, we drop the subscript *NMDA* in *S _{j}*). However, input currents, which also include contributions due to nonselective and inhibitory neurons that no longer enter the dynamical equations, must be approximated in a self-consistent manner. Specifically, (30)-(33) are complicated by the fact that

$${I}_{\mathit{syn},1}={\alpha}_{1}{S}_{1}+{\alpha}_{2}{S}_{2}+{\beta}_{1}{\nu}_{1}+{\beta}_{2}{\nu}_{2}+{I}_{\mathit{const},1},$$

(34)

$${I}_{\mathit{syn},2}={\alpha}_{2}{S}_{1}+{\alpha}_{1}{S}_{2}+{\beta}_{2}{\nu}_{1}+{\beta}_{1}{\nu}_{2}+{I}_{\mathit{const},2},$$

(35)

as in [49], where the parameter pairs *α _{j}, β_{j}* reflect the 1 ↔ 2 symmetry of the network. Accounting for all the components of the currents in (34)-(35) is difficult. The firing rates for the selective populations cannot be directly replaced by their steady states

The derivation begins with the current equations (11)-(15) for the four-population model, which include synaptic variables for populations 3 and I that do not appear explicitly in (30)-(33). To account for these we break (11) into currents for each synapse type and population, expressing them as linear functions of *S _{j}, ν_{j}*, and constants, and gather like terms to form (34)-(35). In this section we calculate the coefficients of the input current to population 1

$$\begin{array}{lll}{I}_{\mathit{syn},1}& =& {I}_{\mathit{NMDA},1,1}+{I}_{\mathit{NMDA},2,1}+{I}_{\mathit{NMDA},3,1}+{I}_{\mathit{AMPA},1,1}+{I}_{\mathit{AMPA},2,1}\\ & +& {I}_{\mathit{AMPA},3,1}+{I}_{\mathit{GABA},I,1}+{I}_{\mathit{AMPA},\mathit{ext},1}+{I}_{\mathit{stim},1}+{I}_{\mathit{noise},1},\end{array}$$

(36)

where *I _{type,j,k}* denotes the current from population

$$\begin{array}{lll}{I}_{\mathit{NMDA},1,1}& =& {S}_{1}{N}_{1}{\omega}_{+}\phantom{\rule{0.2em}{0ex}}{J}_{\mathit{NMDA},p};\phantom{\rule{0.2em}{0ex}}{I}_{\mathit{NMDA},2,1}={S}_{2}{N}_{2}{\omega}_{-}\phantom{\rule{0.2em}{0ex}}{J}_{\mathit{NMDA},p};\\ {I}_{\mathit{NMDA},3,1}& =& {\psi}_{3}{N}_{3}{\omega}_{-}\phantom{\rule{0.2em}{0ex}}{J}_{\mathit{NMDA},p};\phantom{\rule{0.2em}{0ex}}{I}_{\mathit{AMPA},1,1}={\nu}_{1}\frac{{T}_{\mathit{AMPA}}}{1000}{N}_{1}{\omega}_{+}\phantom{\rule{0.2em}{0ex}}{J}_{\mathit{AMPA},p};\\ {I}_{\mathit{AMPA},2,1}& =& {\nu}_{2}\frac{{T}_{\mathit{AMPA}}}{1000}{N}_{2}{\omega}_{-}\phantom{\rule{0.2em}{0ex}}{J}_{\mathit{AMPA},p};\phantom{\rule{0.2em}{0ex}}{I}_{\mathit{AMPA},3,1}={\varphi}_{3}\frac{{T}_{\mathit{AMPA}}}{1000}{N}_{3}{\omega}_{-}\phantom{\rule{0.2em}{0ex}}{J}_{\mathit{AMPA},p};\\ {I}_{\mathit{GABA},I,1}& =& {\varphi}_{I}\frac{{T}_{\mathit{GABA}}}{1000}{N}_{I}\phantom{\rule{0.2em}{0ex}}{J}_{\mathit{GABA},p};\phantom{\rule{0.2em}{0ex}}{I}_{\mathit{AMPA},\mathit{ext},1}={N}_{\mathit{ext}}{\varphi}_{\mathit{ext}}\frac{{T}_{\mathit{AMPA}}}{1000}{J}_{\mathit{AMPA},p,\mathit{ext}};\\ {I}_{\mathit{stim},1}& =& {J}_{\mathit{AMPA},\mathit{ext},p}\phantom{\rule{0.1em}{0ex}}{\mu}_{0}(1+E)\frac{{T}_{\mathit{AMPA}}}{1000},\end{array}$$

(37)

and *I*_{noise,1} is determined by (29) of section 3.3. As in the four-population model of section 3, the parameters *J _{type,k}* denote peak synaptic currents, and they retain the same values. The factor

$${\psi}_{3}=\frac{0.641{\varphi}_{3}{T}_{\mathit{NMDA}}}{1000+0.641{\varphi}_{3}{T}_{\mathit{NMDA}}}.$$

(38)

Finally, *ϕ*_{3} and *ϕ _{I}* are the steady-state firing rates of the rapidly equilibrating populations 3 and I, which are calculated as described below and in Appendix C.

The currents of (37) are problematic in that they contain the quantities *ϕ*_{3}, *ψ*_{3}(*ϕ*_{3}), and *ϕ _{I}*, which do not appear explicitly in (34)-(35). These quantities must be expressed in terms of

The main idea is to show that the firing rate *ϕ*_{3} of population 3 can be approximated at a constant value *ϕ*_{3,0} throughout the region of good performance. Self-consistency between populations 3 and *I* is then relatively easy to establish, and the coefficients *α _{j}* and

In Appendix C.1 we linearize *ψ*_{3}(*ϕ*_{3}) about an appropriate baseline value *ϕ*_{3,0} and calculate consistent baseline firing rates
${\varphi}_{3}^{\ast}$ and
${\varphi}_{I}^{\ast}$ for populations 3 and *I*. In Appendices C.2–C.3 we consider the effects of *S _{j}* and

$$\begin{array}{lll}{\alpha}_{1}& =& {N}_{j}{\omega}_{+}\phantom{\rule{0.2em}{0ex}}{J}_{\mathit{NMDA},p}+{N}_{I}\phantom{\rule{0.2em}{0ex}}{J}_{\mathit{GABA},p}\phantom{\rule{0.2em}{0ex}}\left(\frac{{T}_{\mathit{GABA}}}{1000}\right)\phantom{\rule{0.2em}{0ex}}\frac{{c}_{I}\phantom{\rule{0.2em}{0ex}}{N}_{1}\phantom{\rule{0.2em}{0ex}}{J}_{\mathit{NMDA},I}}{{\mathrm{\Gamma}}_{I}},\\ {\alpha}_{2}& =& {N}_{j}{\omega}_{-}\phantom{\rule{0.2em}{0ex}}{J}_{\mathit{NMDA},p}+{N}_{I}\phantom{\rule{0.2em}{0ex}}{J}_{\mathit{GABA},p}\phantom{\rule{0.2em}{0ex}}\left(\frac{{T}_{\mathit{GABA}}}{1000}\right)\phantom{\rule{0.2em}{0ex}}\frac{{c}_{I}\phantom{\rule{0.2em}{0ex}}{N}_{1}\phantom{\rule{0.2em}{0ex}}{J}_{\mathit{NMDA},I}}{{\mathrm{\Gamma}}_{I}},\\ {\beta}_{1}& =& {N}_{j}{\omega}_{+}\phantom{\rule{0.2em}{0ex}}{J}_{\mathit{AMPA},p}\frac{{T}_{\mathit{AMPA}}}{1000}+{N}_{I}\phantom{\rule{0.2em}{0ex}}{J}_{\mathit{GABA},p}\phantom{\rule{0.2em}{0ex}}\left(\frac{{T}_{\mathit{GABA}}}{1000}\right)\phantom{\rule{0.2em}{0ex}}\frac{{c}_{I}\phantom{\rule{0.2em}{0ex}}{N}_{1}\phantom{\rule{0.2em}{0ex}}{J}_{\mathit{AMPA},I}\phantom{\rule{0.1em}{0ex}}\frac{{T}_{\mathit{AMPA}}}{1000}}{{\mathrm{\Gamma}}_{I}},\\ {\beta}_{2}& =& {N}_{j}{\omega}_{-}\phantom{\rule{0.2em}{0ex}}{J}_{\mathit{AMPA},p}\frac{{T}_{\mathit{AMPA}}}{1000}+{N}_{I}\phantom{\rule{0.2em}{0ex}}{J}_{\mathit{GABA},p}\phantom{\rule{0.2em}{0ex}}\left(\frac{{T}_{\mathit{GABA}}}{1000}\right)\phantom{\rule{0.2em}{0ex}}\frac{{c}_{I}\phantom{\rule{0.2em}{0ex}}{N}_{1}\phantom{\rule{0.2em}{0ex}}{J}_{\mathit{AMPA},I}\phantom{\rule{0.2em}{0ex}}\frac{{T}_{\mathit{AMPA}}}{1000}}{{\mathrm{\Gamma}}_{I}},\end{array}$$

(39)

with ${\mathrm{\Gamma}}_{I}=1-{c}_{I}\phantom{\rule{0.2em}{0ex}}{N}_{I}\phantom{\rule{0.2em}{0ex}}{J}_{\mathit{GABA},I}\phantom{\rule{0.2em}{0ex}}\frac{{T}_{\mathit{GABA}}}{1000}$. The constant currents are (cf. (85))

$$\begin{array}{ccc}\hfill {I}_{\mathit{const},j}& =& {I}_{\mathit{AMPA},\mathit{ext},j}+{I}_{\mathit{stim},j}+{I}_{\mathit{Noise},j}+{I}_{\mathit{NMDA},3,j}^{\ast}+{I}_{\mathit{AMPA},3,j}^{\ast}+{I}_{\mathit{GABA},I,j}^{\ast},\hfill \\ \hfill {I}_{\mathit{AMPA},\mathit{ext},j}& =& {N}_{\mathit{ext}}{\varphi}_{\mathit{ext}}\frac{{T}_{\mathit{AMPA}}}{1000}{J}_{\mathit{AMPA},\mathit{ext},p},\hfill \\ \hfill {I}_{\mathit{stim},1}& =& {J}_{\mathit{AMPA},\mathit{ext},p}\phantom{\rule{0.1em}{0ex}}{\mu}_{0}(1+E/100)\frac{{T}_{\mathit{AMPA}}}{1000},\hfill \\ \hfill {I}_{\mathit{stim},2}& =& {J}_{\mathit{AMPA},\mathit{ext},p}\phantom{\rule{0.1em}{0ex}}{\mu}_{0}(1-E/100)\frac{{T}_{\mathit{AMPA}}}{1000},\hfill \\ \hfill {I}_{\mathit{NMDA},3,j}^{\ast}& =& {N}_{3}{\omega}_{-}\phantom{\rule{0.2em}{0ex}}{J}_{\mathit{NMDA},p}({K}_{\mathit{NMDA}}\phantom{\rule{0.1em}{0ex}}({\varphi}_{3}^{\ast}-{\varphi}_{3,0})+\psi ({\varphi}_{3,0})),\hfill \\ \hfill {I}_{\mathit{AMPA},3,j}^{\ast}& =& {N}_{3}{\omega}_{-}\phantom{\rule{0.2em}{0ex}}{J}_{\mathit{AMPA},p}\frac{{T}_{\mathit{AMPA}}}{1000}{\varphi}_{3}^{\ast},\hfill \\ \hfill {I}_{\mathit{GABA},3,I}^{\ast}& =& {N}_{I}\phantom{\rule{0.2em}{0ex}}{J}_{\mathit{GABA},p}\frac{{T}_{\mathit{GABA}}}{1000}{\varphi}_{I}^{\ast}.\hfill \end{array}$$

(40)

The noise currents are the same as those derived in section 3.3. The baseline firing rates are ${\varphi}_{3}^{\ast}={\varphi}_{3,0}=1\phantom{\rule{0.2em}{0ex}}\mathrm{Hz}$, and ${\varphi}_{I}^{\ast}$ is calculated in section C.5 as

$${\varphi}_{I}^{\ast}=\frac{{\varphi}_{I,0}+{c}_{I}\phantom{\rule{0.2em}{0ex}}({I}_{\mathit{AMPA},\mathit{ext},I}+{I}_{\mathit{NMDA},3,I}+{I}_{\mathit{AMPA},3,I}-{I}_{\mathit{thresh},I})}{{\mathrm{\Gamma}}_{I}}.$$

(41)

Finally, recall that the gains *γ _{E}* and

We now explore the effects of neuromodulation on the two-population model, ending with a bifurcation analysis under various gain conditions to compare with the spiking network and four-population models.

We numerically integrate the system defined in (30)-(33), with currents as in (34)-(35), using a timestep Δ*t* = 0.2 ms. The system behaves poorly with *T*_{pop2} = *T _{AMPA}* in (32)-(33), because inhibition due to the interneurons changes with each timestep and Δ

Contour map of reward rate over the glutamatergic-GABA-ergic (γ_{E}, γ_{I}) gain modulation plane for the two-population model with f-I curve ϕ_{E,b} for selective populations and ϕ_{3} = const over part of the region and linearly **...**

Simulating the system as described above, we explore the neuromodulation plane using the same parameters as the four-population model with the f-I curve *ϕ _{E,b}* of (24), obtaining the results shown in Figure 10. The location and curvature of the high reward rate ridge compare well with both Figures 1(a) and 4(b): that this match with the spiking network results is achieved with the same parameters for the two- and four-population models is a particularly good feature of

Figures 11--1212 show bifurcation diagrams for (*γ _{E}, γ_{I}*) = (1, 1) and (2, 2) for both symmetric (

Since the dynamics of the firing rates (*ν*_{1}, *ν*_{2}) are much faster than the NMDA synaptic variables (*S*_{1}, *S*_{2}), the two-population system can be studied through its dynamics projected onto the (*S*_{1}, *S*_{2})-plane. Nullclines for *S*_{1} and *S*_{2} can be determined by fixing *S*_{1} and *S*_{2}, allowing the firing rates *ν*_{1} and *ν*_{2} to rapidly equilibrate, and finding conditions for which *dS*_{1}/*dt* = 0 and *dS*_{2}/*dt* = 0 in (30)-(31). Figure 13 shows these nullclines, the fixed points with their stability types, and sample noisy trajectories of the four-dimensional flow projected on the (*S*_{1}, *S*_{2})-plane. The two components of each nullcline—a line parallel to an axis and an approximate parabola—derive from the threshold in the f-I curve *ϕ _{E,b}*: the latter result from substitution of (34)-(35) into (30)-(33), and the lines correspond to states with synaptic currents that remain below threshold. Note that, for sufficiently large

Dynamics of two-population model projected onto the (S_{1}, S_{2})-plane, with nullclines for S_{1} shown in orange and nullclines for S_{2} shown in green. Stable points are shown as filled circles, saddles with one-dimensional stable manifolds as open triangles, **...**

We consider the standard gains (*γ _{E}, γ_{I}*) = (1, 1) without stimulus (top left), with symmetric stimuli (top right), and with asymmetric stimuli with

Impulsive trials occur without stimulus for (*γ _{E}, γ_{I}*) = (2, 1.2), as in the lower right panel. The two low saddles and the unstable fixed point again lie near the low-low attractor, decreasing its basin of attraction, and since noise is twice as large as before, essentially all trials are impulsive. Similar phase-plane structures for two-dimensional neural systems appear in [21].

The model of neuromodulation developed in [19] and reduced to lower dimensions above can implement a speed-accuracy tradeoff to optimize performance in 2AFC decision tasks. Reward rates can be maximized by appropriately balancing speed with accuracy [44], and subjects can modulate this balance in response to specific task instructions [33]. Changing task conditions such as stimulus strength or *RSI* can change the speed-accuracy tradeoff that maximizes the reward rate (see (1)), e.g., decreasing *RSI* shifts optima toward faster, less accurate performance. The models for speed-accuracy tradeoff in [39, 45, 24, 41, 42] modulate the gain of an input-output relationship or the drift rate of a DD model, and others [44, 8] implement threshold modulation, although these strategies are interchangeable in some models and both can simulate behavior very well. From a basic physiological perspective, however, identical motor responses should correspond to a similar level of firing in motor cortex, and thus threshold modulation may not be plausible.

Tonic (steady) locus coeruleus (LC) activity is known to affect cognitive performance via NE levels [5]. Experiments in a target-detection task found an inverse relationship between interstimulus interval and tonic LC firing rate [34], the latter changing from 3 to 5 Hz as the stimulus presentation frequency increased from 0.6/s to 1.0/s. Motivated by this, here we examine the ability of tonic NE levels, as represented in the four-population model of section 3 with f-I curve *ϕ _{E,b}*, to implement a physiologically derived speed-accuracy tradeoff. Such an investigation with a spiking model would be prohibitively time-consuming, and here the reduced model demonstrates its usefulness.

Simulations were done for gains lying along four line segments in a box of side length 0.6 around the point (*γ _{E}, γ_{I}*) = (1, 1) in the neuromodulation plane (see Figure 4(b)), and locations that maximize reward rate were determined for

Results for linear slices through neuromodulation plane, for fixed gains (dotted black), *2Δ*γ_{E} = *Δ*γ_{I} (magenta), *Δ*γ_{E} = *Δ*γ_{I} (black), *Δ*γ_{E} = *4Δ*γ_{I} (green), and glutamatergic **...**

(*a*) Optimal reward rates as a function of RSI for each set of combinations of (γ_{E}, γ_{I}) shown in Figure 4(*b*). (*b*) Optimal gain path in neuromodulation plane; arrows denote decreasing RSIs.

The results for fixed gains at (*γ _{E}, γ_{I}*) = (1, 1) without neuromodulation are shown in dashed black at the bottom of the reward rate plot. Not surprisingly, this produces the worst overall performance, although decreasing

Performances for *γ _{E}* =

The ability of NE modulation in the present models to implement a speed-accuracy tradeoff as seen in experiments [34] suggests new experiments relating tonic LC firing rate to changes in *RSI*, with more *RSI* values tested than the two of [34]. With a linear mapping from tonic LC to NE levels [6], and assuming a linear mapping from NE to synaptic gains, our models could be fit to this data. This would inform the mapping from LC to synaptic gains, and factoring out the tonic LC to NE part would produce estimates of NE effects on synaptic gains in the full spiking model.

The optimal speed-accuracy tradeoff strategy in these models corresponds to a neuromodulation ratio with the highest slope that preserves robustness but still achieves fast, impulsive performance before the phasic response is disrupted [5]. For the limited gain range explored here, the slice Δ*γ _{E}* = 4 Δ

The actual path of NE effects through the neuromodulation plane may vary across brain regions, possibly due to different distributions of NE receptors. For instance, effects of NE have been found to be similar but with different frequencies of occurrence in thalamus and whisker barrel cortex of rats [16]. Effects of NE may also depend on phasic LC and the time course of NE concentration changes. Extending our models to different tasks and brain regions may require modifying basic model parameters such as the relative synaptic strength *ω*_{+}. These parameters can be fit to behavioral data under fixed task conditions, and task conditions can then be varied with these parameters fixed and changes in performance determined in terms of *γ _{E}* and

This study of physiologically realistic models of neuromodulation builds on the spiking neuronal network simulations of [19], extending from linear slices to the two-dimensional neuromodulation plane mapped in section 2.2. Many different system behaviors are observed in the spiking model, but it is difficult to understand underlying mechanisms in such a high-dimensional setting. The low-dimensional models developed in sections 3–4 facilitate this investigation: their bifurcation diagrams, in particular, are helpful in interpreting the results of Figure 2 in terms of population dynamics of cells correlated with the two choices. In the region where the full model exhibits failure of pooled inhibition, only a high-high attractor exists without stable states to which either selective population could be suppressed. In the below threshold region only low-low stable states exist and decision thresholds are never reached. In the peak discriminability region, stimulus-induced loss of stability of the low-low state can force decisions. Impulsive choices can also occur due to noise perturbing the state across a basin boundary toward a choice attractor, and higher noise levels at high excitatory gains make such jumps more likely. At still higher gains, the low-low state is unstable in the absence of stimuli, and all trials are impulsive. However, firing rates of the selective populations continue to split until the high-low choice attractors disappear in the no discriminability region.

Unlike earlier models that were tuned to fit reaction time distributions for a fixed set of synaptic gains [49], the present low-dimensional models were developed to match the global neuromodulation results of Figures 1--2.2. Such models cannot accurately describe all aspects of neural systems, but they do allow the study of specific behaviors in analytically tractable settings. Indeed, knowledge of the underlying neurophysiology is currently insufficient to formulate detailed models, as in the present case of NE modulation. Low-dimensional models should then be judged by their ability to increase the understanding of basic system phenomena and to suggest further experiments. In this regard our two- and four-population models are ready to incorporate physiological mappings from LC activity via NE concentration to synaptic gains, once such mappings are experimentally determined.

The main purpose of these models is to link physiological effects of NE and other neuromodulators to behaviors in decision-making tasks, as illustrated in section 5. Different locations in the neuromodulation plane correspond to different levels of subject performance. High-dimensional, cellular-level models obviously permit better representations of neurophysiology, but they are difficult to parameterize and computationally intensive. As demonstrated in [47, 19], such models can include the cellular and synaptic detail needed to relate neuromodulator effects to behavioral outcomes, but many trials must be simulated to amass reliable data and fit multiple experimental conditions, which demands extensive computational resources.

In contrast, one-dimensional DD and OU models and two-dimensional leaky accumulators, with appropriate parameterization, can describe behavioral aspects of two-alternative decisions remarkably well [46, 25, 26, 8, 18]. They are not only fast and simple to simulate, but they admit analytical investigation. Our reduction from 9200 to 11 and finally 4 variables provides a bridge across this dimensional divide that promises to better relate neuromodulation at the cellular scale to bifurcation of attractors and low-dimensional models in cognitive psychology. For the computational cost of simulating 500 trials in a single gain condition in the spiking model, it is possible to simulate 500 four-population or 2000 two-population trials over a grid of 900 gain conditions. This enables studies such as that of section 5, which are not computationally feasible in the spiking model. Moreover, as more neurophysiological data become available, the low-dimensional models can be improved.

In summary, this set of physiologically inspired models provides a bridge from neuromodulation at the single-cell level to decision-making behavior and thereby allows improvements in knowledge in each area to increase understanding in the other. The spiking model of neuromodulation provides a map from neurophysiology to behavior, but many of the neurophysiological details have yet to be experimentally determined, especially for NE dynamics under phasic LC activity. Model simulations raise important questions and suggest new experiments, as described here and in [19], but there is a limit to how far the spiking model can be taken at this time. In contrast, the reduced models are computationally efficient enough to allow data-fitting in LC-recording experiments with varying *RSI*s, creating a map from behavior back to neurophysiology.

*E* = 0.128, *μ*_{0} = 40 Hz, threshold = 20 Hz.

*V _{thresh}* = −50 mV,

*T _{AMPA}* = 2ms,

*g _{AMPA,ext,p}* = 2.1 nS,

*V _{E,AMPA}* = 0 mV,

*J _{AMPA,ext,p}* = 0.11025 nA,

*c _{p}* = 310 Hz/nA,

*ϕ _{p,max}* = 100 Hz,

In order to calculate the effect of neuromodulation on the noise current, we include the synaptic conductance *g _{AMPA,ext,k}* in the synaptic variables

$${S}_{\mathit{peak}}\phantom{\rule{0.1em}{0ex}}\mathrm{exp}(-1/\mathit{f\; \tau})={S}_{\mathit{peak}}-g\phantom{\rule{1em}{0ex}}\mathrm{or}\phantom{\rule{1em}{0ex}}{S}_{\mathit{peak}}=\frac{g}{1-\mathrm{exp}(-1/\mathit{f\; \tau})}.$$

(42)

The average is computed by integrating over the interspike interval:

$${S}_{\mathit{avg}}=f{\mathit{\int}}_{0}^{1/f}{S}_{\mathit{peak}}{e}^{-t/\tau}\mathit{dt}=\mathit{f\; \tau \; g}.$$

(43)

For a Poisson spike train of mean frequency *f* interspike intervals *ξ* are exponentially distributed with density

$$f{e}^{-\mathit{f\xi}}\mathit{d\xi},$$

(44)

so that the jumps in *S _{type,k}* are also exponentially distributed; see Figure 16(a). If

$$\eta =g-(\mathit{g\; f\; \tau}+g)(1-{e}^{-\xi /\tau}).$$

(45)

These fluctuations can be positive or negative, depending upon the interspike interval.

Solving (45) for *ξ* in terms of *η* and differentiating, we obtain

$$\xi =-\tau \phantom{\rule{0.2em}{0ex}}\mathrm{log}\phantom{\rule{0.2em}{0ex}}\left(1-\frac{g-\eta}{\mathit{g\; f\; \tau}+g}\right)\phantom{\rule{1em}{0ex}}\mathrm{and}\phantom{\rule{1em}{0ex}}\mathit{d\xi}=\frac{-\mathit{\tau \; d\eta}}{(1-\frac{g-\eta}{\mathit{g\; f\; \tau}+g})(\mathit{g\; f\; \tau}+g)},$$

(46)

and substituting (45) into (44) yields the distribution for the fluctuations in *η*:

$$\frac{\mathit{f\; \tau}}{\mathit{g\; f\; \tau}+g}\phantom{\rule{0.2em}{0ex}}{\left(\frac{\mathit{g\; f\; \tau}+\eta}{\mathit{g\; f\; \tau}+g}\right)}^{\mathit{f\; \tau}-1}\mathit{d\eta}.$$

(47)

Here *η* varies from *g* to −*g f τ* as *ξ* varies in (0,∞); see Figure 16(b). Reversing the limits of integration so that *η* varies from −*g f τ* to *g* adds the extra negative sign which makes (47) a proper nonnegative probability distribution. At the end of this appendix we check that this distribution has zero mean and show that its variance is

$$\mathrm{Var}(\eta )=\frac{f{g}^{2}\tau}{\mathit{f\; \tau}+2}.$$

(48)

Hence, as noted in section 3.3, the standard deviation varies linearly with *g*, and *g* may be factored out in (27)-(29). The fluctuations in current input to neuron *j* in the spiking neuronal model are equal to Δ(*gS*)(*V _{j} − V_{E}*) (cf. (9)), and once

The Central Limit Theorem [22] implies that the distribution of a sum of *n* independent, identically distributed random variables, with mean *a* and variance *b*, converges to a normal distribution with mean *na* and variance *nb* as *n* → ∞. We now appeal to this to argue that the Poisson spike trains entering individual cells in the spiking model can be replaced by additive Gaussian noise in the averaged synaptic variable ODEs (17)-(19) of the population model.

For a single cell with input frequency *f* = 2.4 kHz and timesteps Δ*t* = 0.1 ms (appropriate for numerical integration of (17)-(19)), *n* = *f* Δ*t* = 0.24 is small, and so the argument does not apply. However, the variables *S _{type,j/I}* in (17)-(19) are obtained by averaging over populations of

$$\frac{n\phantom{\rule{0.1em}{0ex}}\mathrm{Var}(\eta )}{n{N}^{2}}=\frac{\mathrm{Var}\phantom{\rule{0.1em}{0ex}}(\eta )}{{N}^{2}}.$$

(49)

For large *N* this is small, so that *S _{avg}* for a single cell will be close to

We now consider changes *S*(*t* + *s*) − *S*(*s*) due to fluctuations from time *s* to *s* + *t* (*t* ≥ 0). These have expected value 0 since all increments have expected value zero, and since ≈ *N f t* incoming spikes occur in (*s, s* + *t*), the variance of their summed effect is

$$\mathrm{Var}\phantom{\rule{0.1em}{0ex}}(S(t+s)-S(s))=\frac{N\phantom{\rule{0.2em}{0ex}}\mathit{f\; t}\mathrm{Var}\phantom{\rule{0.1em}{0ex}}(\eta )}{{N}^{2}}=\frac{\mathit{f\; t}\mathrm{Var}\phantom{\rule{0.1em}{0ex}}(\eta )}{N}\stackrel{\mathrm{def}}{=}{\sigma}^{2}t.$$

(50)

From the previous paragraph, this summed noise is normally distributed, so it satisfies the requirements of a Wiener process [23]. We may therefore model the noise current *I* by a stochastic differential equation with additive Gaussian noise, as in (29).

Finally we check that the mean and variance of the distribution (47) are as stated above. The mean is given by

$$\mu (\eta )={{\mathit{\int}}_{-\mathit{g\; f\; \tau}}^{g}\frac{\mathit{f\; \tau}}{g+\mathit{g\; f\; \tau}}\phantom{\rule{0.1em}{0ex}}\left(\frac{\mathit{g\; f\; \tau}+\eta}{\mathit{g\; f\; \tau}+g}\right)}^{\mathit{f\; \tau}-1}\mathit{\eta \; d\eta},$$

(51)

which may be integrated by parts to give

$${\left[\eta \phantom{\rule{0.2em}{0ex}}{\left(\frac{\mathit{g\; f\; \tau}+\eta}{\mathit{g\; f\; \tau}+g}\right)}^{\mathit{f\; \tau}}\right]}_{-\mathit{g\; f\; \tau}}^{g}-{{\mathit{\int}}_{-\mathit{g\; f\; \tau}}^{g}\phantom{\rule{0.2em}{0ex}}\left(\frac{\mathit{g\; f\; \tau}+\eta}{\mathit{g\; f\; \tau}+g}\right)}^{\mathit{f\; \tau}}\mathit{d\eta}=g-{\left[{\left(\frac{\mathit{g\; f\; \tau}+\eta}{\mathit{g\; f\; \tau}+g}\right)}^{\mathit{f\; \tau}+1}\frac{\mathit{g\; f\; \tau}+g}{\mathit{f\; \tau}+1}\right]}_{-\mathit{g\; f\; \tau}}^{g}=0.$$

(52)

The expression for variance,

$${{\mathit{\int}}_{-\mathit{g\; f\; \tau}}^{g}\frac{\mathit{f\; \tau}}{g+\mathit{g\; f\; \tau}}\phantom{\rule{0.1em}{0ex}}\left(\frac{\mathit{g\; f\; \tau}+\eta}{\mathit{g\; f\; \tau}+g}\right)}^{\mathit{f\; \tau}-1}{\eta}^{2}\phantom{\rule{0.1em}{0ex}}\mathit{d\eta},$$

(53)

is integrated by parts twice, first to yield

$$\mathrm{Var}\phantom{\rule{0.1em}{0ex}}(\eta )={g}^{2}-\frac{2}{{(\mathit{g\; f\; \tau}+g)}^{\mathit{f\; \tau}}}{\mathit{\int}}_{-\mathit{g\; f\; \tau}}^{g}\eta {(\mathit{g\; f\; \tau}+\eta )}^{\mathit{f\; \tau}}\phantom{\rule{0.1em}{0ex}}\mathit{d\eta},$$

and then finally

$$\begin{array}{l}\mathrm{Var}\phantom{\rule{0.1em}{0ex}}(\eta )={g}^{2}-\frac{2}{{(\mathit{g\; f\; \tau}+g)}^{\mathit{f\; \tau}}}\phantom{\rule{0.2em}{0ex}}{\left[\eta \frac{{(\mathit{g\; f\; \tau}+\eta )}^{\mathit{f\; \tau}+1}}{\mathit{f\; \tau}+1}\right]}_{-\mathit{g\; f\; \tau}}^{g}-\frac{1}{\mathit{f\; \tau}+1}{\mathit{\int}}_{-\mathit{g\; f\; \tau}}^{g}\phantom{\rule{0.2em}{0ex}}{(\mathit{g\; f\; \tau}+\eta )}^{\mathit{f\; \tau}+1}\mathit{d\eta}\\ \phantom{\rule{1em}{0ex}}={g}^{2}-2{g}^{2}+\frac{2{(\mathit{g\; f\; \tau}+g)}^{\mathit{f\; \tau}+2}}{(\mathit{f\; \tau}+1)\phantom{\rule{0.1em}{0ex}}(\mathit{f\; \tau}+2)\phantom{\rule{0.1em}{0ex}}{(\mathit{g\; f\; \tau}+g)}^{\mathit{f\; \tau}}}=\frac{f{g}^{2}\tau}{\mathit{f\; \tau}+2}.\end{array}$$

(54)

In Appendices C.1–C.5 we provide details of the self-consistency calculations, sketched in section 4.1.2, that close the ODEs (34)-(35).

The firing rates of populations 3 and I must be expressed as functions of *S _{j}, ν_{j}*, and constant input currents. This is done by separating the inputs to these populations, calculating the output corresponding to each input, and summing them. For such a linear superposition to work, the f-I curves must also be linear (or piecewise-linear with fluctuations confined to a single linear piece). We use the form

$${\varphi}_{j}({I}_{\mathit{syn},j})={\varphi}_{j,0}+{c}_{j}\theta ({I}_{\mathit{syn},j}-{I}_{\mathit{thresh},j}),\phantom{\rule{1em}{0ex}}j=3,I,$$

(55)

where *ϕ*_{j,0} is the minimum firing rate for *I _{syn,j}* <

To match the f-I curve *ϕ _{E,b}* in the critical region below the decision threshold (20 Hz), we set

$${\psi}_{3}({\varphi}_{3})\approx {\psi}_{3}^{\ast}+{K}_{\mathit{KMDA}}{\varphi}_{3},$$

(56)

where

$${\psi}_{3}^{\ast}=-{K}_{\mathit{NMDA}}{\varphi}_{3,0}+{\psi}_{3}({\varphi}_{3,0})\phantom{\rule{1em}{0ex}}\mathrm{and}\phantom{\rule{1em}{0ex}}{K}_{\mathit{NMDA}}=\frac{0.641{T}_{\mathit{NMDA}}}{{(1000+0.641{T}_{\mathit{NMDA}}{\varphi}_{3,0})}^{2}}.$$

(57)

(Here we use *ϕ*_{3,0} = 1 Hz as for *ϕ _{E,b}* in the four-population model, but the derivation is general.) If

$$\mathrm{\Delta}{\varphi}_{3}={c}_{3}\mathrm{\Delta}{I}_{\mathit{syn},3}\phantom{\rule{1em}{0ex}}\mathrm{and}\phantom{\rule{1em}{0ex}}\mathrm{\Delta}{\varphi}_{I}={c}_{I}\phantom{\rule{0.1em}{0ex}}\mathrm{\Delta}{I}_{\mathit{syn},I}.$$

(58)

To begin to separate *ϕ*_{3} and *ϕ _{I}*, we ignore inputs from populations 1 and 2 and focus on the external currents to populations 3 and I and recurrent connections in populations 3 and I. For

$${\varphi}_{j}({I}_{\mathit{syn},j})={\varphi}_{j,0}+{c}_{j}({I}_{\mathit{syn},j}-{I}_{\mathit{thresh},j}),\phantom{\rule{1em}{0ex}}j=3,I.$$

(59)

Define
${\varphi}_{3}^{\ast}$ and
${\varphi}_{I}^{\ast}$ as the solutions to this linear system. If
${\varphi}_{3}^{\ast}>{\varphi}_{3,0}$ and
${\varphi}_{I}^{\ast}>{\varphi}_{I,0}$, then the system is in the purely linear range and (58) applies. If either condition fails, *ϕ*_{3} and *ϕ _{I}* must be treated otherwise, as described in sections C.4–C.5.

We solve for the baseline firing rates
${\varphi}_{3}^{\ast}$ and
${\varphi}_{I}^{\ast}$ from (55), using the components of the currents *I*_{syn,3} and *I _{syn,I}*, except for the inputs from populations 1 and 2:

$${\varphi}_{3}^{\ast}={\varphi}_{3,0}+{c}_{3}\phantom{\rule{0.1em}{0ex}}({I}_{\mathit{AMPA},\mathit{ext},3}+{I}_{\mathit{NMDA},3,3}+{I}_{\mathit{AMPA},3,3}+{I}_{\mathit{GABA},I,3}-{I}_{\mathit{thresh},p}),$$

(60)

$${\varphi}_{I}^{\ast}={\varphi}_{I,0}+{c}_{I}\phantom{\rule{0.1em}{0ex}}({I}_{\mathit{AMPA},\mathit{ext},I}+{I}_{\mathit{NMDA},3,I}+{I}_{\mathit{AMPA},3,I}+{I}_{\mathit{GABA},I,I}-{I}_{\mathit{thresh},I}),$$

(61)

in which

$$\begin{array}{l}{I}_{\mathit{NMDA},3,3}={N}_{3}{J}_{\mathit{NMDA},p}{\psi}_{3}({\varphi}_{3}^{\ast});\phantom{\rule{0.2em}{0ex}}{I}_{\mathit{NMDA},3,I}={N}_{3}\phantom{\rule{0.1em}{0ex}}{J}_{\mathit{NMDA},I}{\psi}_{3}({\varphi}_{3}^{\ast});\\ {I}_{\mathit{AMPA},3,3}={N}_{3}\phantom{\rule{0.1em}{0ex}}{J}_{\mathit{AMPA},P}\frac{{T}_{\mathit{AMPA}}}{1000}{\varphi}_{3}^{\ast};\phantom{\rule{0.2em}{0ex}}{I}_{\mathit{AMPA},3,I}={N}_{3}\phantom{\rule{0.1em}{0ex}}{J}_{\mathit{AMPA},I}\frac{{T}_{\mathit{AMPA}}}{1000}{\varphi}_{3}^{\ast};\\ {I}_{\mathit{GABA},I,3}={N}_{I}\phantom{\rule{0.1em}{0ex}}{J}_{\mathit{GABA},p}\frac{{T}_{\mathit{GABA}}}{1000}{\varphi}_{I}^{\ast};\phantom{\rule{0.2em}{0ex}}{I}_{\mathit{GABA},I,I}={N}_{I}\phantom{\rule{0.1em}{0ex}}{J}_{\mathit{GABA},I}\frac{{T}_{\mathit{GABA}}}{1000}{\varphi}_{I}^{\ast},\\ {I}_{\mathit{AMPA},\mathit{ext},3}={N}_{\mathit{ext}}{\varphi}_{\mathit{ext}}\frac{{T}_{\mathit{AMPA}}}{1000}{J}_{\mathit{AMPA},\mathit{ext},p,}\\ {I}_{\mathit{AMPA},\mathit{ext},I}={N}_{\mathit{ext}}{\varphi}_{\mathit{ext}}\frac{{T}_{\mathit{AMPA}}}{1000}{J}_{\mathit{AMPA},\mathit{ext},I.}\end{array}$$

(62)

Substituting the expressions (62) into (60)-(61), using (56), and collecting terms yields a linear system of the form

$${a}_{1}{\varphi}_{3}^{\ast}+{b}_{1}{\varphi}_{I}^{\ast}+{d}_{1}=0,$$

(63)

$${a}_{2}{\varphi}_{3}^{\ast}+{b}_{2}{\varphi}_{I}^{\ast}+{d}_{2}=0,$$

(64)

where

$$\begin{array}{lll}{a}_{1}& =& {c}_{3}{N}_{3}\phantom{\rule{0.2em}{0ex}}\left({J}_{\mathit{AMPA},p}\frac{{T}_{\mathit{AMPA}}}{1000}+{J}_{\mathit{NMDA},p}{K}_{\mathit{NMDA}}\right)-1,\\ {b}_{1}& =& {c}_{3}\phantom{\rule{0.1em}{0ex}}{N}_{I}\phantom{\rule{0.1em}{0ex}}{J}_{\mathit{GABA},p}\frac{{T}_{\mathit{GABA}}}{1000},\\ {d}_{1}& =& {\varphi}_{3,0}+{c}_{3}\phantom{\rule{0.1em}{0ex}}[{I}_{\mathit{AMPA},\mathit{ext},3}+{N}_{3}\phantom{\rule{0.1em}{0ex}}{J}_{\mathit{NMDA},p}(-{K}_{\mathit{NMDA}}{\varphi}_{3,0}+{\psi}_{3,0})-{I}_{\mathit{thresh},p}],\\ {a}_{2}& =& {c}_{I}{N}_{3}\phantom{\rule{0.2em}{0ex}}\left({J}_{\mathit{AMPA},I}\frac{{T}_{\mathit{AMPA}}}{1000}+{J}_{\mathit{NMDA},I}{K}_{\mathit{NMDA}}\right),\\ {b}_{2}& =& {c}_{I}\phantom{\rule{0.1em}{0ex}}{N}_{I}\phantom{\rule{0.1em}{0ex}}{J}_{\mathit{GABA},I}\frac{{T}_{\mathit{GABA}}}{1000}-1,\\ {d}_{2}& =& {\varphi}_{I,0}+{c}_{I}\phantom{\rule{0.1em}{0ex}}[{I}_{\mathit{AMPA},\mathit{ext},I}+{N}_{3}\phantom{\rule{0.1em}{0ex}}{J}_{\mathit{NMDA},I}\phantom{\rule{0.1em}{0ex}}(-{K}_{\mathit{NMDA}}{\varphi}_{3,0}+{\psi}_{3,0})-{I}_{\mathit{thresh},i}].\end{array}$$

(65)

Equations (63)-(64) have the solution

$${\varphi}_{3}^{\ast}=\frac{-{d}_{1}{b}_{2}+{d}_{2}{b}_{1}}{\mathrm{\Gamma}},\phantom{\rule{1em}{0ex}}{\varphi}_{I}^{\ast}=\frac{{d}_{1}{a}_{2}-{d}_{2}{a}_{1}}{\mathrm{\Gamma}},$$

(66)

where Γ = *a*_{1}*b*_{2} − *a*_{2}*b*_{1}. If (58) holds, the effects of *S _{j}* and

$${I}_{\mathit{NMDA},3,j}^{\ast}={N}_{3}{\omega}_{-}\phantom{\rule{0.1em}{0ex}}{J}_{\mathit{NMDA},p}({K}_{\mathit{NMDA}}({\varphi}_{3}^{\ast}-{\varphi}_{3,0})+\psi ({\varphi}_{3,0})),$$

(67)

$${I}_{\mathit{AMPA},3,j}^{\ast}={N}_{3}{\omega}_{-}\phantom{\rule{0.1em}{0ex}}{J}_{\mathit{AMPA},p}\frac{{T}_{\mathit{AMPA}}}{1000}{\varphi}_{3}^{\ast},$$

(68)

$${I}_{\mathit{GABA},3,I}^{\ast}={N}_{I}\phantom{\rule{0.1em}{0ex}}{J}_{\mathit{GABA},p}\frac{{T}_{\mathit{GABA}}}{1000}{\varphi}_{I}^{\ast}.$$

(69)

These currents (67)-(69) become the part of the constant terms of (34) due to populations 3 and *I*. The currents into populations 3 and *I* due to *S _{j}* and

Each of the four state variables that remain in the ODEs (34)-(35) affects the four-population system through a similar set of currents: the current back to its own population, the current to the other selective population, the current to the nonselective population, and the current to the interneuron population. Those to populations 3 and I are especially difficult to estimate because increasing the input to each population changes the currents of populations 3 and I among themselves. To show how these currents affect the system, we first describe the effects of *S*_{1}, which influence the four populations through the following *NMDA* currents:

$${I}_{\mathit{NMDA},1,1}={N}_{1}{\omega}_{+}\phantom{\rule{0.1em}{0ex}}{J}_{\mathit{NMDA},p}{S}_{1},\phantom{\rule{1em}{0ex}}{I}_{\mathit{NMDA},1,2}={N}_{1}{\omega}_{-}\phantom{\rule{0.1em}{0ex}}{J}_{\mathit{NMDA},p}{S}_{1};$$

(70)

$${I}_{\mathit{NMDA},1,3}={N}_{1}1{J}_{\mathit{NMDA},p}{S}_{1},\phantom{\rule{1em}{0ex}}{I}_{\mathit{NMDA},1,I}={N}_{1}1{J}_{\mathit{NMDA},I}{S}_{1}.$$

(71)

Since populations 1 and 2 are explicitly tracked in (30)-(33), we need only calculate the changes in *ϕ*_{3} and *ϕ _{I}* due to

$$\Delta {\varphi}_{3}\left({S}_{1}\right)={c}_{3}\Delta {I}_{\mathit{syn},3}\left({S}_{1}\right)\phantom{\rule{1em}{0ex}}\mathrm{and}\phantom{\rule{1em}{0ex}}\Delta {\varphi}_{I}\phantom{\rule{0.1em}{0ex}}\left({S}_{1}\right)={c}_{I}\Delta {I}_{\mathit{syn},I}\phantom{\rule{0.1em}{0ex}}\left({S}_{1}\right),$$

(72)

which imply that

$$\begin{array}{lll}\mathrm{\Delta}{\varphi}_{3}({S}_{1})& =& {c}_{3}\phantom{\rule{0.2em}{0ex}}[{I}_{\mathit{NMDA},1,3}+{N}_{3}\phantom{\rule{0.2em}{0ex}}\left(\frac{{T}_{\mathit{AMPA}}}{1000}{J}_{\mathit{AMPA},p}+{K}_{\mathit{NMDA}}\phantom{\rule{0.1em}{0ex}}{J}_{\mathit{NMDA},p}\right)\phantom{\rule{0.2em}{0ex}}\mathrm{\Delta}{\varphi}_{3}\\ & +& {N}_{I}\frac{{T}_{\mathit{GABA}}}{1000}{J}_{\mathit{GABA},p}\mathrm{\Delta}{\varphi}_{I}\phantom{\rule{0.2em}{0ex}}],\end{array}$$

(73)

$$\begin{array}{lll}\mathrm{\Delta}{\varphi}_{I}({S}_{1})& =& {c}_{I}\phantom{\rule{0.2em}{0ex}}[{I}_{\mathit{NMDA},1,I}+{N}_{3}\phantom{\rule{0.2em}{0ex}}\left(\frac{{T}_{\mathit{AMPA}}}{1000}{J}_{\mathit{AMPA},I}+{K}_{\mathit{NMDA}}\phantom{\rule{0.1em}{0ex}}{J}_{\mathit{NMDA},I}\right)\phantom{\rule{0.2em}{0ex}}\mathrm{\Delta}{\varphi}_{3}\\ & +& {N}_{I}\frac{{T}_{\mathit{GABA}}}{1000}{J}_{\mathit{GABA},I}\mathrm{\Delta}{\varphi}_{I}].\end{array}$$

(74)

Collecting terms, (73)-(74) have the same form as (63)-(64), with *a*_{1}, *b*_{1}, *a*_{2}, and *b*_{2} as in (65) and

$${d}_{1}={c}_{3}{I}_{\mathit{NMDA},1,3},\phantom{\rule{0.4em}{0ex}}{d}_{2}={c}_{I}{I}_{\mathit{NMDA},1,I},$$

(75)

so that the solution to (73)-(74) is

$$\mathrm{\Delta}{\varphi}_{3}({S}_{1})=\frac{-{c}_{3}{b}_{2}{N}_{1}{J}_{\mathit{NMDA},p}+{c}_{I}{b}_{1}{N}_{1}{J}_{\mathit{NMDA},I}}{\mathrm{\Gamma}}{S}_{1},$$

(76)

$$\mathrm{\Delta}{\varphi}_{I}({S}_{1}\phantom{\rule{0.1em}{0ex}})=\frac{-{c}_{3}{a}_{2}{N}_{1}{J}_{\mathit{NMDA},p}-{c}_{I}{a}_{1}{N}_{1}{J}_{\mathit{NMDA},I}}{\mathrm{\Gamma}}{S}_{1},$$

(77)

where Γ = *a*_{1}*b*_{2} − *a*_{2}*b*_{1}. The effects on population 1 are now calculated as follows, starting with *I*_{NMDA,3,1}:

$$\begin{array}{ccc}\hfill {I}_{\mathit{NMDA},3,1}& =& {I}_{\mathit{NMDA},3,1}^{\ast}+\mathrm{\Delta}{I}_{\mathit{NMDA},3,1},\hfill \\ \hfill \mathrm{\Delta}{I}_{\mathit{NMDA},3,1}\phantom{\rule{0.1em}{0ex}}({S}_{1})& =& {N}_{3}{\omega}_{-}\phantom{\rule{0.1em}{0ex}}{J}_{\mathit{NMDA},p}{K}_{\mathit{NMDA}}\mathrm{\Delta}{\varphi}_{3}({S}_{1})\hfill \\ & =& {N}_{3}{\omega}_{-}\phantom{\rule{0.1em}{0ex}}{J}_{\mathit{NMDA},p}{K}_{\mathit{NMDA}}\frac{-{c}_{3}{b}_{2}{N}_{1}{J}_{\mathit{NMDA},p}+{c}_{I}{b}_{1}{N}_{1}\phantom{\rule{0.1em}{0ex}}{J}_{\mathit{NMDA},I}}{\mathrm{\Gamma}}{S}_{1}.\hfill \end{array}$$

(78)

The other currents due to *S*_{1} are calculated similarly, yielding
${I}_{\mathit{AMPA},3,1}={I}_{\mathit{AMPA},3,1}^{\ast}+\mathrm{\Delta}{I}_{\mathit{AMPA},3,1}$ and
${I}_{\mathit{GABA},I,1}={I}_{\mathit{GABA},I,1}^{\ast}+\mathrm{\Delta}{I}_{\mathit{GABA},I,1}$, with

$$\begin{array}{lll}\mathrm{\Delta}{I}_{\mathit{AMPA},3,1}({S}_{1})& =& {N}_{3}{\omega}_{-}\phantom{\rule{0.1em}{0ex}}{J}_{\mathit{AMPA},p}\phantom{\rule{0.2em}{0ex}}\left(\frac{{T}_{\mathit{AMPA}}}{1000}\right)\phantom{\rule{0.2em}{0ex}}\frac{-{c}_{3}{b}_{2}{N}_{1}\phantom{\rule{0.1em}{0ex}}{J}_{\mathit{NMDA},p}+{c}_{I}{b}_{1}{N}_{1}\phantom{\rule{0.1em}{0ex}}{J}_{\mathit{NMDA},I}}{\mathrm{\Gamma}}{S}_{1},\\ \mathrm{\Delta}{I}_{\mathit{GABA},I,1}({S}_{1})& =& {N}_{1}\phantom{\rule{0.1em}{0ex}}{J}_{\mathit{GABA},p}\frac{{T}_{\mathit{GABA}}}{1000}\frac{{c}_{3}{a}_{2}{N}_{1}\phantom{\rule{0.1em}{0ex}}{J}_{\mathit{NMDA},p}-{c}_{I}{a}_{1}{N}_{1}\phantom{\rule{0.1em}{0ex}}{J}_{\mathit{NMDA},I}}{\mathrm{\Gamma}}{S}_{1}.\end{array}$$

Due to symmetry, the currents to population 2 from populations 3 and I are identical to those to population 1 and the effect of *S*_{2} can be determined by interchanging subscripts 1 and 2 in the above formulae. The effects of *ν*_{1} and *ν*_{2} collected in the *β*’s are similar to the effects of *S _{j}* in the

We are now in a position to compute the coefficients *α _{j}, β_{j}*, and

$$\begin{array}{lll}{\alpha}_{1}& =& {\alpha}_{1a}+{\alpha}_{1b}+{\alpha}_{1c},\phantom{\rule{0.4em}{0ex}}\mathrm{where}\phantom{\rule{0.4em}{0ex}}{\alpha}_{1a}=\frac{{I}_{\mathit{NMDA},j,j}}{{S}_{j}},\\ {\alpha}_{1b}& =& \frac{\mathrm{\Delta}{I}_{\mathit{NMDA},3,j}({S}_{j})+\mathrm{\Delta}{I}_{\mathit{AMPA},3,j}({S}_{j})}{{S}_{j}},\phantom{\rule{0.4em}{0ex}}\mathrm{and}\phantom{\rule{0.4em}{0ex}}{\alpha}_{1c}=\frac{\mathrm{\Delta}{I}_{\mathit{GABA},I,j}({S}_{j})}{{S}_{j}},\end{array}$$

(79)

which imply

$$\begin{array}{lll}{\alpha}_{1a}& =& {N}_{j}{\omega}_{+}\phantom{\rule{0.1em}{0ex}}{J}_{\mathit{NMDA},p,}\\ {\alpha}_{1b}& =& {N}_{3}{\omega}_{-}\phantom{\rule{0.2em}{0ex}}\left(\frac{{T}_{\mathit{AMPA}}}{1000}{J}_{\mathit{AMPA},p}+{K}_{\mathit{NMDA}}\phantom{\rule{0.2em}{0ex}}{J}_{\mathit{NMDA},p}\right)\phantom{\rule{0.2em}{0ex}}\left(\frac{{c}_{3}{b}_{2}{N}_{j}\phantom{\rule{0.2em}{0ex}}{J}_{\mathit{NMDA},p}-{c}_{I}{b}_{1}{N}_{j}\phantom{\rule{0.2em}{0ex}}{J}_{\mathit{NMDA},I}}{\mathrm{\Gamma}}\right),\\ {\alpha}_{1c}& =& {N}_{I}\frac{{T}_{\mathit{GABA}}}{1000}{J}_{\mathit{GABA},p}\phantom{\rule{0.2em}{0ex}}\left(\frac{{c}_{3}{b}_{2}{N}_{j}\phantom{\rule{0.2em}{0ex}}{J}_{\mathit{NMDA},p}-{c}_{I}{b}_{1}{N}_{j}\phantom{\rule{0.2em}{0ex}}{J}_{\mathit{NMDA},I}}{\mathrm{\Gamma}}\right).\end{array}$$

(80)

The coefficient *α*_{2} is computed in the same manner:

$$\begin{array}{lll}{\alpha}_{2}& =& {\alpha}_{2a}+{\alpha}_{2b}+{\alpha}_{2c},\phantom{\rule{0.2em}{0ex}}\mathrm{where}\phantom{\rule{0.2em}{0ex}}{\alpha}_{2a}=\frac{{I}_{\mathit{NMDA},2,1}}{{S}_{2}}={N}_{2}{\omega}_{-}\phantom{\rule{0.1em}{0ex}}{J}_{\mathit{NMDA},p},\\ {\alpha}_{2b}& =& \frac{\mathrm{\Delta}{I}_{\mathit{NMDA},3,1}\left({S}_{2}\right)+\mathrm{\Delta}{I}_{\mathit{AMPA},3,1}\left({S}_{2}\right)}{{S}_{2}}={\alpha}_{1b},\phantom{\rule{0.2em}{0ex}}\mathrm{and}\phantom{\rule{0.2em}{0ex}}{\alpha}_{2c}=\frac{\mathrm{\Delta}{I}_{\mathit{GABA},I,1}\left({S}_{2}\right)}{{S}_{2}}={\alpha}_{1c}.\end{array}$$

(81)

The last two equalities hold because the currents to populations 1 and 2 from populations 3 and I are the same. We have already noted the symmetry *S*_{1} ↔ *S*_{2} inherent in the reduced system.

We next calculate the *β _{j}*’s, which capture the effect of the currents mediated by the AMPA synapses of populations 1 and 2 through their firing rates

$$\begin{array}{lll}{\beta}_{1}& =& {\beta}_{1a}+{\beta}_{1b}+{\beta}_{1c},\phantom{\rule{0.2em}{0ex}}\mathrm{where}\phantom{\rule{0.2em}{0ex}}{\beta}_{1a}=\frac{{I}_{\mathit{AMPA},j,j}}{{\nu}_{j}},\\ {\beta}_{1b}& =& \frac{\mathrm{\Delta}{I}_{\mathit{NMDA},3,j}\phantom{\rule{0.2em}{0ex}}({\nu}_{j})+\mathrm{\Delta}{I}_{\mathit{AMPA},3,j}\phantom{\rule{0.2em}{0ex}}({\nu}_{j})}{{\nu}_{j}},\phantom{\rule{0.2em}{0ex}}\mathrm{and}\phantom{\rule{0.2em}{0ex}}{\beta}_{1c}=\frac{\mathrm{\Delta}{I}_{\mathit{GABA},I,j}\phantom{\rule{0.2em}{0ex}}({\nu}_{j})}{{\nu}_{j}},\end{array}$$

(82)

which imply

$$\begin{array}{l}{\beta}_{1a}={N}_{j}{\omega}_{+}\phantom{\rule{0.1em}{0ex}}{J}_{\mathit{AMPA},p}\frac{{T}_{\mathit{AMPA}}}{1000},\\ {\beta}_{1b}={N}_{3}{\omega}_{-}\phantom{\rule{0.2em}{0ex}}\left(\frac{{T}_{\mathit{AMPA}}}{1000}{J}_{\mathit{AMPA},p}+{K}_{\mathit{NMDA}}\phantom{\rule{0.1em}{0ex}}{J}_{\mathit{NMDA},p}\right)\phantom{\rule{0.2em}{0ex}}{N}_{j}{T}_{\mathit{AMPA}}\phantom{\rule{0.2em}{0ex}}\left(\frac{{c}_{3}{b}_{2}\phantom{\rule{0.2em}{0ex}}{J}_{\mathit{AMPA},p}-{c}_{I}{b}_{1}\phantom{\rule{0.2em}{0ex}}{J}_{\mathit{AMPA},I}}{1000\mathrm{\Gamma}}\right),\\ {\beta}_{1c}={N}_{I}\frac{{T}_{\mathit{GABA}}}{1000}{J}_{\mathit{GABA},p}\phantom{\rule{0.2em}{0ex}}\left(\frac{{c}_{3}{b}_{2}{N}_{j}\phantom{\rule{0.2em}{0ex}}{J}_{\mathit{AMPA},p}{T}_{\mathit{AMPA}}-{c}_{I}{b}_{1}\phantom{\rule{0.2em}{0ex}}{N}_{j}\phantom{\rule{0.2em}{0ex}}{J}_{\mathit{AMPA},I}{T}_{\mathit{AMPA}}}{1000\mathrm{\Gamma}}\right);\end{array}$$

(83)

and

$$\begin{array}{lll}{\beta}_{2}& =& {\beta}_{2a}+{\beta}_{2b}+{\beta}_{2c},\phantom{\rule{0.2em}{0ex}}\mathrm{where}\phantom{\rule{0.2em}{0ex}}{\beta}_{2a}=\frac{{I}_{\mathit{AMPA},2,1}}{{\nu}_{2}}={N}_{j}{\omega}_{-}\phantom{\rule{0.1em}{0ex}}{J}_{\mathit{AMPA},p}\frac{{T}_{\mathit{AMPA}}}{1000},\\ {\beta}_{2b}& =& \frac{\mathrm{\Delta}{I}_{\mathit{NMDA},3,1}({\nu}_{2})+\mathrm{\Delta}{I}_{\mathit{AMPA},3,1}({\nu}_{2})}{{\nu}_{2}}={\beta}_{1b},\phantom{\rule{0.2em}{0ex}}\mathrm{and}\phantom{\rule{0.4em}{0ex}}{\beta}_{2c}=\frac{\mathrm{\Delta}{I}_{\mathit{GABA},I,1}({\nu}_{2})}{{\nu}_{2}}={\beta}_{1c}.\end{array}$$

(84)

Finally, we enumerate the constant currents in (34)-(35):

$$\begin{array}{ccc}\hfill {I}_{\mathit{const},j}& =& {I}_{\mathit{AMPA},\mathit{ext},j}+{I}_{\mathit{stim},j}+{I}_{\mathit{Noise},j}+{I}_{\mathit{NMDA},3,j}^{\ast}+{I}_{\mathit{AMPA},3,j}^{\ast}+{I}_{\mathit{GABA},I,j}^{\ast},\hfill \\ \hfill {I}_{\mathit{AMPA},\mathit{ext},j}& =& {N}_{\mathit{ext}}{\varphi}_{\mathit{ext}}\frac{{T}_{\mathit{AMPA}}}{1000}{J}_{p,\mathit{AMPA},\mathit{ext}},\hfill \\ \hfill {I}_{\mathit{stim},1}& =& {J}_{\mathit{AMPA},\mathit{ext},p}{\mu}_{0}\left(1+E/100\right)\frac{{T}_{\mathit{AMPA}}}{1000},\hfill \\ \hfill {I}_{\mathit{stim},2}& =& {J}_{\mathit{AMPA},\mathit{ext},p}{\mu}_{0}\left(1-E/100\right)\frac{{T}_{\mathit{AMPA}}}{1000},\hfill \\ \hfill {I}_{\mathit{NMDA},3,j}^{\ast}& =& {N}_{3}{\omega}_{-}\phantom{\rule{0.1em}{0ex}}{J}_{\mathit{NMDA},p}({K}_{\mathit{NMDA}}({\varphi}_{3}^{\ast}-{\varphi}_{3,0})+\psi ({\varphi}_{3,0})),\hfill \\ \hfill {I}_{\mathit{AMPA},3,j}^{\ast}& =& {N}_{3}{\omega}_{-}\phantom{\rule{0.1em}{0ex}}{J}_{\mathit{AMPA},p}\frac{{T}_{\mathit{AMPA}}}{1000}{\varphi}_{3}^{\ast},\hfill \\ \hfill {I}_{\mathit{GABA},3,I}^{\ast}& =& {N}_{I}\phantom{\rule{0.1em}{0ex}}{J}_{\mathit{GABA},p}\frac{{T}_{\mathit{GABA}}}{1000}{\varphi}_{I}^{\ast}.\hfill \end{array}$$

(85)

The noise currents are the same as in the four-population model.

Currents of the form in (34)-(35) can be calculated provided that *ϕ _{I}* and

$${\varphi}_{3}^{\ast}>{\varphi}_{3,0},\phantom{\rule{0.4em}{0ex}}{\varphi}_{I}^{\ast}>{\varphi}_{I,0},\phantom{\rule{0.4em}{0ex}}\mathrm{\Gamma}\ne 0$$

(86)

(the latter must hold if *ϕ*_{3} and *ϕ _{I}* are to be well defined). The values for
${\varphi}_{I}^{\ast}$ and
${\varphi}_{3}^{\ast}$ are calculated in section C.1 as they vary with

(Left panel) Two-population linearity conditions: the red curve marks the boundary of
${\mathrm{\varphi}}_{\mathrm{I}}^{\ast}>{\mathrm{\varphi}}_{\mathrm{I},\mathit{0}}$ for case (*a*), which fails to its upper left, the thick blue curve shows *Γ* = *0*, and the black curve shows that I_{syn,} **...**

Above the black curve,
${\varphi}_{3}^{\ast}$ is fixed at *ϕ*_{3,0} and
${\varphi}_{I}^{\ast}$ is recalculated as described in section C.5, and the results are examined to verify if case (B) applies. The thick magenta curve denotes the new linear regime for
${\varphi}_{I}^{\ast}$, which now includes the region of good performance. The *α*’s, *β*’s, and constant currents for fixed *ϕ*_{3} can therefore be recalculated as described in section C.5. This requires *I*_{syn,3} to remain below threshold (region above black curve) for all time in the presence of inputs from populations 1 and 2. This assumption is now tested.

Above the black curve
${\varphi}_{3}^{\ast}={\varphi}_{3,0}$, but if *I*_{syn,3} were to increase above *I _{thresh,p}* due to inputs from populations 1 and 2, the assumption of constant

$$\frac{d{\varphi}_{3}}{\mathit{dI}}=\frac{d}{\mathit{dI}}\phantom{\rule{0.2em}{0ex}}\left(\frac{{c}_{3}{b}_{2}I-{c}_{I}{b}_{1}I/1.3}{\mathrm{\Gamma}}\right)=\frac{{c}_{3}{b}_{2}-{c}_{I}{b}_{1}/1.3}{\mathrm{\Gamma}};$$

(87)

cf. (76). The right panel of Figure 17 shows that *dϕ*_{3}/*dI* < 0 for *γ _{E}* >≈ 0.9, which covers much of the region of interest. If

For the small portion of the region of good performance in which *dϕ*_{3}/*dI* > 0, *I*_{syn,3} is far enough below threshold that although *I*_{syn,3} increases with the net effect of populations 1 and 2, their limited firing rates in the green region of Figure 2 do not allow *I*_{syn,3} to reach *I _{thresh,p}*. Applying case (B) parameters to the entire region between the magenta and black curves in the left panel has the effect of eliminating the loss of discriminability region from Figure 2 and extending “peak discriminability region” performance all the way to the boundary of the no discriminability region. This is a shortcoming of this assumption but does not critically affect the model.

Below the black curve in the left panel of Figure 17, the case (A) model parameters from section C.3 apply. To the left of the magenta curve, which is fully in the blue region, there are insufficient external inputs to any population, all firing rates remain at their minimum values (case C), and the model is trivial. Case (D) never occurs because *I*_{syn,3} > *I _{thresh,p}*

We now solve for
${\varphi}_{I}^{\ast}$, the baseline value of *ϕ _{I}*, in case (B) of section C.4 (

$${\varphi}_{I}^{\ast}={\varphi}_{I,0}+{c}_{I}\phantom{\rule{0.2em}{0ex}}\left({I}_{\mathit{AMPA},\mathit{ext},I}+{I}_{\mathit{NMDA},3,I}+{I}_{\mathit{AMPA},3,I}+{N}_{I}\phantom{\rule{0.2em}{0ex}}{J}_{\mathit{GABA},I}\frac{{T}_{\mathit{GABA}}}{1000}{\varphi}_{I}^{\ast}-{I}_{\mathit{thresh},i}\right),$$

(88)

where

$${I}_{\mathit{NMDA},3,I}={N}_{3}\phantom{\rule{0.2em}{0ex}}{J}_{\mathit{NMDA},p}\psi ({\varphi}_{3,0}),\phantom{\rule{1em}{0ex}}{I}_{\mathit{AMPA},3,I}={N}_{3}\phantom{\rule{0.2em}{0ex}}{J}_{\mathit{AMPA},p}\frac{{T}_{\mathit{AMPA}}}{1000}{\varphi}_{3,0},$$

(89)

and *I _{AMPA,ext,I}* is the same as in the last equation of (62). Setting
${\mathrm{\Gamma}}_{I}=1-{c}_{I}\phantom{\rule{0.1em}{0ex}}{N}_{I}\phantom{\rule{0.1em}{0ex}}{J}_{\mathit{GABA},I}\phantom{\rule{0.1em}{0ex}}\frac{{T}_{\mathit{GABA}}}{1000}$, we obtain

$${\varphi}_{I}^{\ast}=\frac{{\varphi}_{I,0}+{c}_{I}\phantom{\rule{0.2em}{0ex}}({I}_{\mathit{AMPA},\mathit{ext},I}+{I}_{\mathit{NMDA},3,I}+{I}_{\mathit{AMPA},3,I}-{I}_{\mathit{thresh},I})}{{\mathrm{\Gamma}}_{I}}.$$

(90)

Here *α*_{1a}, *α*_{2a}, *β*_{1a}, and *β*_{2a} are the same as under case (A) in section C.3, (79)-(84), and since *ϕ*_{3} is constant, *α*_{1b} = *α*_{2b} = *β*_{1b} = *β*_{2b} = 0.

The effect of population I can be calculated as follows:

$$\begin{array}{l}\mathrm{\Delta}{\varphi}_{I}=\frac{{c}_{I}\phantom{\rule{0.2em}{0ex}}{I}_{\mathit{syn},I}}{{\mathrm{\Gamma}}_{I}},\phantom{\rule{0.4em}{0ex}}\mathrm{\Delta}{I}_{\mathit{GABA},I,j}={N}_{I}\phantom{\rule{0.2em}{0ex}}{J}_{\mathit{GABA},p}\phantom{\rule{0.2em}{0ex}}\left(\frac{{T}_{\mathit{GABA}}}{1000}\right)\phantom{\rule{0.2em}{0ex}}\left(\frac{{c}_{I}}{{\mathrm{\Gamma}}_{I}}\right)\phantom{\rule{0.2em}{0ex}}\mathrm{\Delta}{I}_{\mathit{syn},I},\\ \mathrm{and}\phantom{\rule{0.4em}{0ex}}{\alpha}_{1c}={N}_{I}\phantom{\rule{0.2em}{0ex}}{J}_{\mathit{GABA},p}\phantom{\rule{0.2em}{0ex}}\left(\frac{{T}_{\mathit{GABA}}}{1000}\right)\phantom{\rule{0.2em}{0ex}}\frac{{c}_{I}\phantom{\rule{0.2em}{0ex}}{N}_{1}\phantom{\rule{0.2em}{0ex}}{J}_{\mathit{NMDA},I}}{{\mathrm{\Gamma}}_{I}},\phantom{\rule{1em}{0ex}}{\alpha}_{2c}={\alpha}_{1c},\\ \phantom{\rule{2em}{0ex}}{\beta}_{1c}={N}_{I}\phantom{\rule{0.2em}{0ex}}{J}_{\mathit{GABA},p}\phantom{\rule{0.2em}{0ex}}\left(\frac{{T}_{\mathit{GABA}}}{1000}\right)\phantom{\rule{0.2em}{0ex}}\frac{{c}_{I}\phantom{\rule{0.2em}{0ex}}{N}_{1}\phantom{\rule{0.2em}{0ex}}{J}_{\mathit{AMPA},I}\frac{{T}_{\mathit{AMPA}}}{1000}}{{\mathrm{\Gamma}}_{I}},\phantom{\rule{1em}{0ex}}{\beta}_{2c}={\beta}_{1c}.\end{array}$$

(91)

The constant currents are the same as in (85) but with ${\varphi}_{3}^{\ast}={\varphi}_{3,0}$ and ${\varphi}_{I}^{\ast}$ as calculated in (90). These new parameter values then define the model for the case (B) condition.

To the left of the magenta curve in Figure 17, the firing rates trivially remain at the values *ϕ*_{1} = *ϕ*_{2} = *ϕ*_{3} = *ϕ*_{p,0} and *ϕ _{I}* =

^{*}This work was supported by the National Institute of Mental Health (MH62196, Cognitive and Neural Mechanisms of Conflict and Control, Silvio M. Conte Center), and by the Air Force Research Laboratory (FA9550-07-1-0537).

1. Abbott LF, Chance ES. Drivers and modulators from push-pull and balanced synaptic input. Prog Brain Res. 2005;149:147–155. [PubMed]

2. Amit DJ, Brunel N. Dynamics of a recurrent network of spiking neurons before and following learning. Netw Comput Neural Syst. 1997;8:373–404.

3. Amit DJ, Brunel N. Model of global spontaneous activity and local structured activity during delay periods in the cerebral cortex. Cereb Cortex. 1997;7:237–274. [PubMed]

4. Amit DJ, Tsodyks MV. Quantitative study of attractor neural network retrieving at low spike rates I: Substrate-spikes, rates and neuronal gain. Network. 1991;2:259–274.

5. Aston-Jones G, Cohen JD. An integrative theory of locus coeruleus-norepinephrine function: Adaptive gain and optimal performance. Annu Rev Neurosci. 2005;28:403–450. [PubMed]

6. Berridge CW, Abercrombie ED. Relationship between locus coeruleus discharge rates and rates of norepinephrine release within neocortex as assessed by in vivo microdialysis. Neuroscience. 1999;93:1263–1270. [PubMed]

7. Berridge CW, Waterhouse BD. The locus coeruleus-noradrenergic system: Modulation of behavioral state and state-dependent cognitive processes. Brain Res Rev. 2003;42:33–84. [PubMed]

8. Bogacz R, Brown E, Moehlis J, Holmes P, Cohen JD. The physics of optimal decision making: A formal analysis of models of performance in two alternative forced choice tasks. Psychological Review. 2006;113:700–765. [PubMed]

9. Bouret S, Sara SJ. Network reset: A simplified overarching theory of locus coeruleus noradrenaline function. Trends Neurosci. 2005;28:574–582. [PubMed]

10. Britten KH, Shadlen MN, Newsome WT, Movshon JA. The analysis of visual motion: A comparison of neuronal and psychophysical performance. J Neurosci. 1992;12:4745–4765. [PubMed]

11. Britten KH, Shadlen MN, Newsome WT, Movshon JA. Response of neurons in macaque MT to stochastic motion signals. Vis Neurosci. 1993;10:1157–1169. [PubMed]

12. Brunel N, Chance FS, Fourcaud N, Abbott LF. Effects of synaptic noise and filtering on the frequency response of spiking neurons. Phys Rev Lett. 2001;86:2186–2189. [PubMed]

13. Brunel N, Wang X-J. Effects of neurmodulation in a cortical network model of object working memory dominated by recurrent inhibition. J Comput Neurosci. 2001;11:63–85. [PubMed]

14. Dayan P, Abbott LF. Theoretical Neuroscience: Computational and Mathematical Modeling of Neural Systems. MIT Press; Cambridge, MA: 2001.

15. Deco G, Rolls ET. Decision-making and Weber’s law: A neurophysiological model. Eur J Neurosci. 2006;24:901–906. [PubMed]

16. Devilbiss DM, Waterhouse BD. The effects of tonic locus ceruleus output on sensory-evoked responses of ventral posterior medial thalamic and barrel field cortical neurons in the awake rat. J Neurosci. 2004;24:10773–10785. [PubMed]

17. Doya K. Modulators of decision making. Nat Neurosci. 2008;11:410–416. [PubMed]

18. Eckhoff P, Holmes P, Law C, Connolly PM, Gold JI. On diffusion processes with variable drift rates as models for decision making during learning. New J Phys. 2008;10 015006. [PMC free article] [PubMed]

19. Eckhoff P, Wong K-F, Holmes P. Optimality and robustness of a biophysical decision-making model under norepinephrine modulation. J Neurosci. 2009;29:4301–4311. [PMC free article] [PubMed]

20. Ermentrout B. Software Environ Tools. Vol. 14. SIAM; Philadelphia: 2002. Simulating, Analyzing, and Animating Dynamical Systems: A Guide to XPPAUT for Researchers and Students.

21. Ermentrout B. Phase-plane analysis of neural nets. In: Arbib MA, editor. The Handbook of Brain Theory and Neural Networks. MIT Press; Cambridge, MA: 2003. pp. 881–885.

22. Feller W. An Introduction to Probability Theory and Its Applications. 2. Wiley; New York: 1957.

23. Gardiner CW. Handbook of Stochastic Methods. 2. Springer; New York: 1985.

24. Gilzenrat MG, Holmes BD, Rajkowski J, Aston-Jones G, Cohen JD. Simplified dynamics in a model of noradrenergic modulation of cognitive performance. Neural Networks. 2002;15:647–663. [PubMed]

25. Gold JI, Shadlen MN. Neural computations that underlie decisions about sensory stimuli. Trends in Cognitive Science. 2001;5:10–16. [PubMed]

26. Gold JI, Shadlen MN. Banburismus and the brain: Decoding the relationship between sensory stimuli, decisions, and reward. Neuron. 2002;36:299–308. [PubMed]

27. Guckenheimer J, Holmes PJ. Nonlinear Oscillations, Dynamical Systems and Bifurcations of Vector Fields. Springer-Verlag; New York: 1983.

28. Horwitz GD, Newsome WT. Separate signals for target selection and movement specification in the superior colliculus. Science. 1999;284:1158–1161. [PubMed]

29. Horwitz GD, Newsome WT. Target selection for saccadic eye movements: Prelude activity in the superior colliculus during a direction-discrimination task. J Neurophysiol. 2001;86:2543–2558. [PubMed]

30. Jones CKRT. Lecture Notes in Math 1609. Springer; Berlin: 1995. Geometric singular perturbation theory, in Dynamical Systems (Montecatini Terme, 1994) pp. 44–118.

31. Kim JN, Shadlen MN. Neural correlates of a decision in the dorsolateral prefrontal cortex. Nat Neurosci. 1999;2:176–185. [PubMed]

32. Marti D, Deco G, Mattia M, Gigante G, Del Giudice P. A fluctuation-driven mechanism for slow decision processes in reverberant networks. PLoS ONE. 2008;3:e2534. [PMC free article] [PubMed]

33. Palmer J, Huk AC, Shadlen MN. The effect of stimulus strength on the speed and accuracy of a perceptual decision. J Vis. 2005;5:376–404. [PubMed]

34. Rajkowski J, Majczynski H, Clayton E, Aston-Jones G. Activation of monkey locus coeruleus neurons varies with difficulty and performance in a target detection task. J Neurophysiol. 2004;92:361–371. [PubMed]

35. Renart A, Brunel N, Wang X-J. Mean-field theory of irregularly spiking neuronal populations and working memory in recurrent cortical networks. In: Feng J, editor. Computational Neuroscience: A Comprehensive Approach. Chapman and Hall/CRC; Boca Raton, FL: 2004. pp. 431–490.

36. Ricciardi LM. Diffusion Processes and Related Topics in Biology. Springer; Berlin: 1977.

37. Roitman J, Shadlen M. Response of neurons in the lateral intraparietal area during a combined visual discrimination reaction time task. J Neurosci. 2002;22:9475–9489. [PubMed]

38. Sara SJ. The locus coeruleus and noradrenergic modulation of cognition. Nat Rev Neuro. 2009;10:211–223. [PubMed]

39. Servan-Schreiber D, Printz H, Cohen JD. A network model of catecholamine effects: Gain, signal-to-noise ratio, and behavior. Science. 1990;249:892–895. [PubMed]

40. Shadlen MN, Newsome WT. Neural basis of a perceptual decision in the parietal cortex (area LIP) of the rhesus monkey. J Neurophysiol. 2001;86:1916–1936. [PubMed]

41. Shea-Brown E, Gao J, Holmes P, Bogacz R, Gilzenrat MS, Cohen JD. Simple neural networks that optimize decisions. Internat J Bifur Chaos Appl Sci Engrg. 2005;15:803–826.

42. Shea-Brown E, Gilzenrat M, Cohen JD. Optimization of decision making in multilayer networks: The role of locus coeruleus. Neural Comput. 2008;20:2863–2894. [PubMed]

43. Shea-Brown E, Moehlis J, Holmes P, Clayton E, Rajkowski J, Aston-Jones G. The influence of spike rate and stimulus duration on noradrenergic neurons. J Comput Neurosci. 2004;17:13–29. [PubMed]

44. Simen P, Cohen JD, Holmes P. Rapid decision threshold modulation by reward rate in a neural network. Neural Networks. 2006;19:1013–1026. [PMC free article] [PubMed]

45. Usher M, Cohen JD, Servan-Schreiber D, Rajkowski J, Aston-Jones G. The role of locus coeruleus in the regulation of cognitive performance. Science. 1999;283:549–554. [PubMed]

46. Usher M, McClelland JL. On the time course of perceptual choice: The leaky competing accumulator model. Psych Rev. 2001;108:550–592. [PubMed]

47. Wang X-J. Probabilistic decision making by slow reverberation in cortical circuits. Neuron. 2002;36:955–968. [PubMed]

48. Waterhouse BD, Devilbiss D, Fleischer D, Sessler FM, Simpson KL. New perspectives on the functional organization and postsynaptic influences of the locus ceruleus efferent projection system. Adv Pharmacol. 1998;42:749–754. [PubMed]

49. Wong KF, Wang X-J. A recurrent network mechanism of time integration in perceptual decisions. J Neurosci. 2006;26:1314–1328. [PubMed]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |