|Home | About | Journals | Submit | Contact Us | Français|
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 , 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 , 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 .
In  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 , 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 , 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 . 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(103) 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  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 , 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 . 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
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  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 N1 = N2 = 240 pyramidal cells, a nonselective population of N3 = 1120 pyramidal cells, and a population of NI = 400 inhibitory interneurons. The state variables are the (subthreshold) transmembrane voltages Vj(t) (mV) of each cell j, the network’s internal gating variables SAMPA,j(t), SNMDA,j(t), and SGABA,j(t), and gating variables SAMPA,ext,j(t) external to the local circuit, which include noisy inputs from MT neurons as well as other brain areas, that enter all 2000 cells. These Stype,j’s describe integrated presynaptic activities which are summed and then multiplied by connection weights and synaptic conductances to yield an effective conductance in units of nS: hence only NI SGABA,j’s are required for the interneurons, and NE glutamatergic synaptic variables for the pyramidal cells. Connection weights ωj,k among the cells depend on the identity of pre- and postsynaptic populations and are supposed to have formed in a Hebbian context, with connections within each selective population enhanced and connections to other excitatory populations suppressed [2, 3].
The ODEs describing this leaky integrate-and-fire neuronal model take the form
where type = AMPA or GABA, Ttype is the time constant for that synapse type, and the summation of delta functions in (3) is over spikes in the presynaptic neuron. Here j [1, 2000] identifies presynaptic cells with j [1, 1600] for pyramidal cells and j [1601, 2000] for interneurons, with different cell types tracking their appropriate synaptic variables.
The slower NMDA synapses require two ODEs for each j [0, 1600] with separate time constants describing the rise and fall of SNMDA,j:
When Vj(t) crosses a threshold Vthresh at time the cell emits a delta function , after which Vj is instantaneously reset and held at Vreset for an absolute refractory period τref.
Spikes from cells external to the local network, with a combined presynaptic firing rate fext = 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 , which would normally be taken as independently generated Poisson processes, are approximated by drawing from independent Gaussian distributions with mean and variance
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 SAMPA,ext,j(t) are calculated as
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
Here the subscript rec (henceforth omitted) indicates recurrent connections within the network, VE and VI are the glutamatergic and GABA-ergic synaptic reversal potentials, gtype,k are synaptic conductances, and gL in (2) denotes the leak conductance. For each synapse type there are two conductance values: gtype,p for postsynaptic pyramidal neurons and gtype,I for postsynaptic interneurons. Stimuli are included by adding terms of the form μ0(1 ± E) to fext in (6) defining μS,AMPA,ext,k for the selective cells, where μ0 is a background input and E [−1, 1] is the coherence defined above. These will be used as bifurcation parameters in section 3.5. All parameter values for the model are listed in Appendix A.
Equations (2)-(9) define a 9200-dimensional dynamical system. Specifically, k ranges from 1 to 2000 for voltages Vk and external inputs SAMPA,ext,k, but j ranges only from 1 to NE = 1600 for the excitatory variables SAMPA,j, SNMDA,j, and xj, and for SGABA,j.
Neuromodulation can be modeled by changing the leak gL and synaptic conductances gtype,k. In , we showed that modulating gL has similar effects to gtype,k, and we shall not vary gL here. Experiments have yet to distinguish any differential effects of NE on NMDA and AMPA (glutamatergic) synapses, so we introduce two parameters, γE which scales gNMDA and gAMPA equally, and γI which scales gGABA. These are allowed to vary separately, since there is insufficient physiological data to determine them from NE concentrations [NE]. In  both tonic and phasic (dynamic) NE release were modeled (cf. [45, 43]), but here we consider only tonic NE levels. When a mapping from [NE] to (γE, γI) is determined experimentally, it can be implemented in our model. Our approach may also be generalized to include other neuromodulators which affect the same synapse types.
In  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 (γE, γI) [0, 3] × [0, 3], plotting values of the performance measures that enter (1).
Figure 1(a) shows that no trials are rewarded for γE < 0.65, because low excitation cannot overcome leak conductance, even without inhibition. For γE > 0.65 good performance exists between the no-choice region on the left and the impulsive region on the right, with reward rates gradually dropping along a line of slope γI/γE = 1. A line of slope 1 can remain in the good region far to the right but quickly encounters no-choices to the left. This failure mode is reversed for steeper slopes as no-choices appear for high γE, and impulsive choices as γE drops below 1 (cf. the 3:1 ratio case ), and a slope γI : γE = 2 keeps performance above the impulsive regime longest before reaching γE = 0.65.
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 γE to escape this regime, which extends rightward at higher γI.
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 ω+ = 1.7 within each selective population. This region ends with loss of stability of the low-low state in the absence of stimulus, so that every trial becomes impulsive.
The loss of discriminability region marks a transition from splitting to race model behavior , 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  to incorporate neuromodulation. This yields 11-dimensional systems of stochastic ODEs in which populations become homogeneous units described by average firing rates νj(t), j = 1, 2, 3, and νI(t) (four variables), the inhibitory population carries one population-averaged synaptic variable SGABA(t), and each excitatory population has two such variables SAMPA,j(t) and SNMDA,j(t) (six in all). The synaptic variables will be determined by firing rates and firing rates by net input currents to each population. We start by calculating these currents.
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 = (Vreset+Vthresh)/2 to estimate synaptic currents that enter each of the four cell populations as follows:
(multiplication of gtype,k in (nS) and ( − Vtype) in (mV) requires division by 1000 to convert Jtype,k to nA, the unit of choice for currents). These effective currents to populations replace the currents into individual cells in (9). Different Jtype,k’s are used for pyramidal populations (j = 1, 2, 3, also written Jtype,p) and interneurons (Jtype,I), because gtype,p ≠ gtype,I. Since ( − VGABA < 0) the inhibitory currents are negative, as expected. We include in JNMDA,k the factor [1 + (1/3.57) exp(− 0.062)]−1 ≈ 0.12 to incorporate the effects of Mg2+ block (cf. (9)). Neuromodulation is now accomplished by γE and γI scaling the Jtype,k’s which contain the gtype,k’s that were previously scaled in section 2.2.
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:
Here ωj,k is the weight from population j to population k, being ω+ for recurrent connections within a selective population, ω− from one selective population to another excitatory population, and ωj,k = 1 otherwise. Nj is the number of excitatory cells in population j and NI the number of interneurons. External inputs IAMPA,ext,k are given by
where Next is the number of connections per cell from presynaptic neurons external to the network (here 800), νext is their average firing rate (3 Hz), and conversion between time constants (ms) and firing rates (Hz) requires division by 1000. Inoise,k is calculated below in section 3.3 and Appendix B. Finally, Istim,k = JAMPA,ext,p μ0 (1 ± E) for the selective populations when stimuli are present.
Following [13, 49], we account for the rise in SNMDA,j due to presynaptic spikes by including a term F(ψ(νj)) = ψj/[TNMDA(1 − ψj)], where ψj is the steady state of SNMDA,j. For a Poisson spike train, this is well approximated by
, which yields F (ψ (νj)) 0.641νj/1000.
and the firing rates obey
Here j = 1, 2, 3, I, and the input-output functions or f-I (firing rate-current) curves ϕj(Isyn,j) model the relationships between input currents and quasi–steady-state firing rates for each population, which all four νj’s approach with time constant TAMPA = 2 ms . We next describe several different f-I curves, distinguishing pyramidal cells and interneurons in each case.
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  such as (2), the expected average interspike interval for solutions to travel from reset voltage Vreset to threshold Vthresh may be calculated analytically, given a mean input current Isyn with white noise fluctuations [36, 4, 35]. The expected firing rate is the reciprocal of this interval plus the absolute refractory period τref:
At low firing rates ϕLIF (Isyn) approaches the reciprocal of the integral in (21) alone and can be approximated by a simpler expression , except at asymptotically low input currents, where the denominator changes sign:
This is plotted as the blue curve labeled Abbott–Chance. Here cp/I is the asymptotic slope, Ithresh,p/I is the threshold for input current, and gp/I is the noise factor that creates nonzero spontaneous activity; values for these parameters are listed in Appendix A. Equation (22) can be improved to account for higher firing rates by including τref, as follows:
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 Jtype,k, they can reasonably approximate the region of good performance in the neuromodulation plane, but steady-state firing rates are too high because ϕAC(Isyn) from (22) exceeds the limit 1/τref = 500 Hz in the upper-right part of the high reward rate region (left part of the peak discriminability region in Figure 2), and ϕE,a only reduces firing rates of the winning selective population in this region to 200 Hz. In addition, with constant synaptic coupling these modifications can capture only a region near (γE, γI) = (1, 1). The full region [0, 3] × [0, 3] of Figure 1 was reproduced only after rescaling the effective synaptic currents Jtype derived from the spiking model (by just under 33% for glutamatergic 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:
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 Isyn 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., ϕ(Isyn,p) ≤ 20 Hz. Several features of ϕE,b accomplish this. First, the spontaneous firing activity ϕp,0 is set to 1 Hz instead of 0. Second, Ithresh,p is reduced from 0.403 nA to 0.384 nA to start the linear rise at a lower input current. In addition, cp is increased from 310 Hz/nA to 352 Hz/nA to produce a steeper initial slope. Noting that the parameter gp/I in (22)-(23) sets the sharpness of the transition from zero firing rate to the rising portion, we set gp = 1 for ϕE,b to create a sharper threshold. Finally, ϕE,b saturates at ϕp,0 + ϕp,max = 101 Hz as Isyn → ∞; cf. the red curve of Figure 3(a). This not only matches ϕE,a of (23) in the relevant region below 20 Hz, but as shown in Figure 4(b), it also fits the spiking network neuromodulation data well, with all parameters within 6% of the values derived from that model. The function ϕE,b allows us to capture global neuromodulatory behavior while retaining the synaptic couplings derived in section 3.1.
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(Isyn,I) with a spontaneous firing rate of 3 Hz, slope of 600 Hz/nA, and cutoff at Ithresh,i = 0.29 nA achieves this; see Figure 3(b) (black):
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  by a piecewise-linear function with the same current threshold and initial slope as the selective populations, which will continue to use ϕE,b:
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 ), 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
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
We can now calculate Inoise,k for inclusion in (15), which expresses IAMPA,ext,k as the sum of a constant average current and Inoise,k (of mean zero). The fluctuations in S are multiplied by JAMPA,ext,k to get the fluctuations in Inoise,k, and this relationship is maintained for varying gains since the standard deviation of fluctuations in gS scales linearly with g (see Appendix B). This is a population average of external inputs in the full model (cf. (8)). Summarizing, each Inoise,k is a sample path of the OU process
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 γI. The model of  was fitted for a single point (γE, γI) = (1, 1) in the neuromodulation plane by changing peak synaptic currents and weights from those of the spiking network as follows: JGABA,p = 1.2JGABA,I instead of the full model’s 1.3JGABA,I and ω+ = 1.8 instead of the full model’s 1.7. Here we wish to obtain a good fit over a broad range of (γE, γI). Computing reward rates over [0, 3] × [0, 3] as in Figure 1, we found that this model is overdominated by excitation: the ridge of optimal performance analogous to that in Figure 1(a) rises too steeply in the (γE, γI) plane. Moreover, as noted in section 3.2, population firing rates are too high in the upper part of this ridge (results not shown).
We next consider the model with f-I curve ϕE,a(Isyn,j) that saturates at 1/τref and with Gaussian noise and synaptic currents derived above (results not shown). The slope of the ridge is reduced but is still too high. In addition, the behavior at (γE, γI) = (1, 1) is characterized by no-choice trials. This may be corrected by setting JGABA,p = 1.393JGABA,I to adjust the slope and scaling glutamatergic currents by a factor in the range [1.3, 1.33]. (This adjustment is perhaps required by our assumption of uniform = −52.5 mV (cf. (10)), which most strongly affects GABA currents, since VE,GABA = −70 mV and VE,AMPA = VE,NMDA = 0 mV.) This shifts the good performance region to include (1, 1), resulting in the reward rate contours of Figure 4(a), which compare reasonably well with those of Figure 1(a). Steady-state population firing rates can exceed those of the spiking network by 100%, but these are less critical for our purposes, since the 2AFC protocol focuses on trajectories up to the 20 Hz threshold and the model performs well in this region.
The match in reward rate contours can be further improved by using the function ϕE,b(Isyn,p) of (24) for the three pyramidal populations and the piecewise-linear function ϕI,b(Isyn,I) of (25) for the interneurons, as shown in Figure 4(b). For this final model, external AMPA, recurrent AMPA and NMDA currents, and JGABA,I are exactly as derived from the spiking network, and JGABA,p = 1.367 JGABA,I is increased by slightly over 5% compared with the spiking network value of 1.3 JGABA,I. These modified f-I curves properly situate the region of good behavior, moreover, and this model captures the concave curvature of the high reward rate ridge. Moreover, it allows simulation of 500 trials over 900 gain conditions in the same time as 500 trials of the spiking model for a single gain condition: this is a speed-up of O(103).
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 , 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 SNMDA,1 and SNMDA,2, written here and henceforth as S1 and S2 for brevity.
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 E = 0, so that the vectorfield of (17)-(20) is reflection-symmetric about the Stype,1 = Stype,2 and ν1 = ν2 hyperplanes. Figure 5(a) shows branches of fixed points at (γE, γI) = (1, 1) in the lower left of the peak discriminability region in Figure 2. For μ0 = 0 the system is tristable: a state with the activities of both selective populations low coexists with a symmetric pair of “memory” attractors [47, 49] corresponding to choices, each with one population high and the other suppressed to zero. The low-low state loses stability in a subcritical pitchfork bifurcation at μ0 ≈ 25, passes through two saddle-node bifurcations, and restabilizes as a high-high state at μ0 ≈ −5, while the stable high-low states survive until they vanish in saddle-nodes at μ0 ≈ 154. The low-low state persists for all μ0 < 0, and the high-low states vanish in saddle-nodes at μ ≈ −118, leaving the low-low state as the sole attractor. Cycling μ0 therefore allows memories of choices to be set and cleared in successive trials (e.g., by a strong inhibitory input).
The bifurcation diagram for (γE, γI) = (2, 2) is similar to that for (1, 1), but the two pitchforks move to μ0 = 35 (low) and μ0 = −219 (high): see Figure 5(b). This implies that both the low-low and high-high states coexist with the symmetric pair of high-low states over the range μ0 (−219, 35), and that all four are stable (in this range there are nine fixed points).
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 , 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 μ0 than in the left panel, but the system exhibits more impulsive trials for (γE, γI) = (2, 2) than for (1, 1). This occurs due to the linear scaling of noise with γE, which doubles its magnitude at (2, 2).
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 Sj, respectively (not shown here). In the latter case μ0 and E have little influence and both Sj’s remain near 1, because both firing rates lie in the saturated part of the f-I curve where changes in input have little effect relative to the overall strong excitation within the network. No decisions are made in the blue region because there are no choice attractors, and splitting does not occur in the no discriminability region because there is a single high-high attractor.
We now turn to an asymmetric case. Figures 6(a) and 6(b) show fixed point branches for S1 and S2 at (γE, γI) = (1, 1) and E = 0.128 for comparison with Figure 5. The symmetric branch is now perturbed and the pitchforks become saddle-node bifurcations at μ0 ≈ 20 and μ0 ≈ 28, thus forcing decisions for μ0 > 20 as solutions approach one of the choice attractors. The asymmetry due to nonzero coherence favors the state with larger S1 (choice 1 correct) over that with larger S2 when noise is present, since the former has a larger basin of attraction. These memory states persist for μ0 (−105, 183) and μ0 (−136, 133), respectively, terminating in saddle-node bifurcations much as in Figure 5. For μ0 > 296 a unique high-high stable state exists, analogous to the one in Figure 5(a). We also computed bifurcation diagrams for (γE, γI) = (2, 2) and (0.5, 1), as for the E = 0 symmetric case above. At (2, 2) the result was generally similar to Figure 6, with some reconnections among branches typical of codimension 2 bifurcations ; a single monotonic branch with S1 > S2 was evident at (0.5, 1).
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 ϕE,a (Figure 4(b)), and the bifurcation diagrams, shown in Figures 7 and and8,8, remain similar to those of Figures 5 and and66.
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 Sj values, and hence steady-state firing rates, are substantially reduced. This is as expected, since ϕE,b saturates at a lower value than ϕE,a (cf. Figure 3(a)). At (γE, γI) = (2, 2) (Figure 7(b)) the high-low states extend over a wider range in μ0 and have higher firing rates, but the branches are much like those for ϕE,a shown in Figure 5(b). For (γE, γI) = (0.5, 1) and (2.5, 0.25) this system also has unique fixed point branches (not shown here), as in the ϕE,a case.
Figure 8 shows an asymmetric case at (γE, γI) = (1, 1) with E = 0.128. Fixed point branches are again similar to those of the analogous ϕE,a case (Figure 6), but the saddle-node in which the low-low state loses stability now occurs at μ0 ≈ 44, higher than the 40 Hz value used in spiking model simulations. The high-high state appears in a saddle-node at μ0 ≈ 20, so that, as for the E = 0 cases of Figure 7, the low-low, high-high, and two high-low states are all stable over a short range of μ0. Although the low-low state loses stability above the stimulus-on condition μ0 = 40, decisions still occur because the noise suffices to perturb states across the low-low basin boundary.
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 DT and nondecision latencies (NDL = 250 ms; cf. (1)). The standard gain cases on the left both produce typical RT distributions with longer tails to the right, while the high-gain conditions shown at right have tighter, faster RT distributions. The lower-right panel (ϕE,b, (γE, γI) = (2, 2)) also has a long tail of impulsive trials at the left, corresponding to the low-stimulus long exponential distribution noted in . In this condition, firing rates of the selective populations tend to start nearer threshold and therefore to cross it sooner and with much less variance following stimulus onset. There are more impulsive choices than for (γE, γI) = (1, 1), and there is no long tail to the right. Changes in the RT distributions along with the corresponding changes in accuracy can implement a speed-accuracy tradeoff, as described in section 5.
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 , 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 . 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 (TAMPA = 2 ms; TGABA = 5 ms), while that for NMDA decay is relatively slow (TNMDA = 100 ms). SAMPA,j and SGABA,j therefore converge rapidly to quasi-steady states that closely track the population firing rates in (18)-(19) so that Stype,j ≈ νjTtype/1000. This eliminates three ODEs for the excitatory populations and one for the inhibitory population, giving a seven-dimensional system.
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 TAMPA. We may therefore set νj(t) ≈ ϕj(Isyn,j(t)) for the nonselective and interneuron populations and remove two firing rate ODEs. Now only three ODEs describing NMDA dynamics in the excitatory populations remain, along with the firing rate ODEs for populations 1 and 2. Finally, when stimuli are on, the nonselective population j = 3 typically exhibits a lower and much less variable firing rate than the selective populations so that SNMDA,3 can be replaced by its steady-state value and its ODE may also be removed, leaving four ODEs for SNMDA,1, SNMDA,2, ν1, and ν2:
(here, as in section 3.5, we drop the subscript NMDA in Sj). 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 Isyn,j contains terms which depend on Sj and ϕE,b(Isyn,j), so that its vectorfield is defined recursively. Ideally, we seek linear relationships of the form
as in , 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 νj = ϕE,b(Isyn,j) to yield a two-dimensional model, as in , since the nonlinearity of ϕE,b(Isyn,j) makes the equations infinitely recursive. To close the system, (32)-(33) are therefore retained to track the firing rates as they approach their steady states. The calculations are summarized below and detailed in Appendix C. Given the other approximations, the rate of approach of νj to its quasi-steady state ϕ(Isyn,j) also must be changed from TAMPA to Tpop2, as noted below in section 4.2.1.
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 Sj, ν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 Isyn,1; those for population 2 follow by interchanging subscripts 1 ↔ 2 and replacing (1 + E) by (1 − E). From (11)-(15) we have
where Itype,j,k denotes the current from population j to population k due to the given synapse,
and Inoise,1 is determined by (29) of section 3.3. As in the four-population model of section 3, the parameters Jtype,k denote peak synaptic currents, and they retain the same values. The factor ψ3 in INMDA,3,1 is the steady-state value of S3 from (16), but with ν3 replaced by ϕ3:
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 S1, S2, ν1, ν2, and constant currents Iconst,j. The calculation employs the piecewise-linear f-I functions of (25) and (26), and one must determine regions of the parameter space in which the input currents exceed threshold, so that these functions are purely linear.
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 βj and constant currents Iconst,j in (34)-(35) can be computed. This calculation is carried out in several stages, as detailed in Appendix C.
In Appendix C.1 we linearize ψ3(ϕ3) about an appropriate baseline value ϕ3,0 and calculate consistent baseline firing rates and for populations 3 and I. In Appendices C.2–C.3 we consider the effects of Sj and νj on all four populations and calculate expressions for αj and βj, assuming pure linearity. We then check, in Appendix C.4, if the input currents do in fact remain above threshold throughout the relevant part of the (γE, γI)-neuromodulation plane. We find that the current due to the nonselective cells (population 3) does not satisfy this condition. We then use the self-consistency results again to show that the relationship between the currents of populations 3 and I implies that ϕ3 can be held constant at . This necessitates resolving for and recomputing some components of the coefficients αj and βj, which is done in Appendix C.5. The resulting values for the constant coefficients in (34)-(35) are as follows (cf. (79)-(85) and (91)):
with . The constant currents are (cf. (85))
The noise currents are the same as those derived in section 3.3. The baseline firing rates are , and is calculated in section C.5 as
Finally, recall that the gains γE and γI enter the glutamatergic- and GABA-ergic currents as multiplicative factors in JNMDA,p/I, JAMPA,p/I, and JGABA,p/I, respectively.
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 Tpop2 = TAMPA in (32)-(33), because inhibition due to the interneurons changes with each timestep and Δt TAMPA = 2 ms. Recurrent excitation due to AMPA currents depending on νj therefore takes longer to reach quasi-equilibrium, while the inhibitory current increases immediately. This affects the quasi-equilibrium values for νj, and the system is incorrectly dominated by inhibition. If we instead set Tpop2 = Δt, the solution remains numerically stable without overshoot, and the neuromodulation results of Figure 10 closely match those for the corresponding four-population model. In this condition, the model is effectively two-dimensional.
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 ϕE,b. We henceforth consider only ϕE,b. We have already noted that ϕE,a gives unrealistically high firing rates in the four-population model; used in the two-population model it fails catastrophically at high gains (results not shown).
Figures 11--1212 show bifurcation diagrams for (γE, γI) = (1, 1) and (2, 2) for both symmetric (E = 0) and asymmetric (E = 0.128) cases, done by determining the coefficients in the current equations from αj’s and βj’s appropriate to these specific locations in the modulation plane. (Note that the change in timescale for Tpop2 is irrelevant in determining fixed points of the noise-free system.) Fixed point branches are very similar to those of the corresponding four-population models (Figures 5--8):8): all bifurcations are of the same types, and occur within μ0 = 1 – 2 Hz of their locations in that system. Interpretation of system behaviors and dynamics is therefore the same as in section 3.5, and we find that the two-population reduction successfully captures behavior over a substantial part of the neuromodulation plane.
Since the dynamics of the firing rates (ν1, ν2) are much faster than the NMDA synaptic variables (S1, S2), the two-population system can be studied through its dynamics projected onto the (S1, S2)-plane. Nullclines for S1 and S2 can be determined by fixing S1 and S2, allowing the firing rates ν1 and ν2 to rapidly equilibrate, and finding conditions for which dS1/dt = 0 and dS2/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 (S1, S2)-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 gp, νj = ϕE,b(Isyn,j) ≈ ϕp,0 for Isyn,j < Ithresh,p, so that (30)-(31) become linear. Above threshold ϕE,b rises monotonically, behaving qualitatively like a threshold-linear function and producing the parabola-like curves when νj = ϕE,b(Isyn,j) is solved for νj in terms of Sj and substituted into (30)-(31).
We consider the standard gains (γE, γI) = (1, 1) without stimulus (top left), with symmetric stimuli (top right), and with asymmetric stimuli with E = 0.128 (lower left). Without stimulus, all trials remain near the low-low fixed point, even in the presence of noise, and both choice attractors persist, representing memory states. With stimulus (top right and bottom left), the low saddles and unstable fixed point approach the low-low stable fixed point, decreasing its basin of attraction and allowing more noisy trials to escape. Although the low-low attractor persists in this model due to the form of the modified f-I curve, it has only a small basin of attraction and the choice attractors dominate the dynamics, as in . The nullclines in the lower-left panel illustrate the effects of asymmetric stimuli, and it is also clear that the high-high state does not significantly influence the dynamics.
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 .
The model of neuromodulation developed in  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 , and subjects can modulate this balance in response to specific task instructions . 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 . Experiments in a target-detection task found an inverse relationship between interstimulus interval and tonic LC firing rate , 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 RSIs in a given range. Figures 14 and and1515 show the results. Each curve in the upper-left panel of Figure 14 plots the reward rate at the optimal gain location (along its segment) for the RSIs specified along the horizontal axis. Continuing clockwise in Figure 14, we show the optimal γE’s, RT’s, and accuracies for the same range of 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 RSI does increase reward rate (see (1)). The most robust case of , ΔγI = 2 ΔγE (magenta curve), does only marginally better, showing that robustness and the ability of neuromodulation to adapt performance under changing conditions are inversely related. Moreover, optimal γ’s (and therefore optimal LC firing rates) decrease with decreasing RSI for this gain modulation ratio (Figure 14, top right), opposite to what is seen experimentally. The other ratios, including those exhibiting the best performance, match the experimental results of .
Performances for γE = γI (solid black) and ΔγE = 4 ΔγI (green) successively improve on both the previous cases. The blue curve (overlapping the green) shows the result for glutamatergic modulation alone, with γI 1. Both excitation-dominated slices (blue and green curves) reach the fast impulsive regime at low RSIs and can therefore achieve the sharp upturn in reward rate that more robust ratios cannot. Finally, when neuromodulation conditions are allowed to vary over the entire region γE, γI [0.8, 1.4], further small improvements are obtained over the best “linear slice” results (red curve in Figure 15(a)). The optimal path through the (γE, γI)-neuromodulation plane which produces the red reward rate curve as RSI changes is shown in Figure 15(b). It rises with slope less than 1 over most of the range before plunging toward (γE, γI) = (1.4, 0.8) as RSIs fall below 300 ms. This sharply negative slope may not occur physiologically due to NE, but other cognitive mechanisms such as a switch to automatic responding can achieve similar effects (cf. the right-hand end of the optimal performance curve in ).
The ability of NE modulation in the present models to implement a speed-accuracy tradeoff as seen in experiments  suggests new experiments relating tonic LC firing rate to changes in RSI, with more RSI values tested than the two of . With a linear mapping from tonic LC to NE levels , 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 . For the limited gain range explored here, the slice ΔγE = 4 ΔγI (green line in Figure 4) achieves this. The more GABA-ergic modulation, the more stable the system to changes in NE concentration, and adding modulation of GABA-ergic synapses does not prevent optimal performance as long as the ratio is dominated by excitation. These results may differ in multilayer decision models , but that case is not examined here.
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 . 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 γI.
This study of physiologically realistic models of neuromodulation builds on the spiking neuronal network simulations of , 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 , 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 , 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 RSIs, creating a map from behavior back to neurophysiology.
E = 0.128, μ0 = 40 Hz, threshold = 20 Hz.
Vthresh = −50 mV, Vreset = −55 mV, VL = −70 mV, gL,p = 25 nS, gL,I = 20 nS, Cp = 0.5 nF, CI = 0.2 nF, τm,p = 20 ms, τm,I = 10 ms, τref,p = 2 ms, τref,I = 1 ms.
TAMPA = 2ms, TNMDA = 100 ms, TNMDA,rise = 2ms, TGABA = 5 ms.
gAMPA,ext,p = 2.1 nS, gAMPA,ext,I = 1.65 nS, gAMPA,p = 0.05 nS, gAMPA,I = 0.04 nS, gNMDA,p = 0.165 nS, gNMDA,I = 0.13 nS, gGABA,p = 1.0 * 1.3 nS, gGABA,I = 1.0 nS.
JAMPA,ext,p = 0.11025 nA, JAMPA,ext,I = 0.08505 nA, JAMPA,p = 0.002625 nA, JAMPA,I = 0.0021, JNMDA,p = 0.0010487, JNMDA,I = 0.0008262, JGABA,p = −0.0175 * 1.3, JGABA,I = −0.0175.
cp = 310 Hz/nA, Ithresh,p = 0.403 nA, gp = 0.16; cI = 615 Hz/nA, Ithresh,I = 0.289 nA, gI = 0.087.
ϕp,max = 100 Hz, ϕp,0 = 1 Hz, gp = 1, cp = 352 Hz/nA, Ithresh,p = 0.384 nA; ϕI,0 = 3 Hz, Ithresh,I = 0.29 nA, cI = 600 Hz/nA.
In order to calculate the effect of neuromodulation on the noise current, we include the synaptic conductance gAMPA,ext,k in the synaptic variables SAMPA,ext,k for this derivation. We first consider a simple case of periodic spike inputs. The values of SNMDA, SAMPA, SAMPA,ext, and SGABA for each population are averages of the Stype,k values for each neuron in the population. For a single cell with an input of frequency f, this approaches a periodic oscillation. Each spike increases Stype,k by the peak synaptic conductance g; it then decays exponentially with synaptic time constant τ, so periodicity is achieved when the decay between input spikes equals g. Letting Speak denote the value immediately after the spike input, we therefore require
The average is computed by integrating over the interspike interval:
For a Poisson spike train of mean frequency f interspike intervals ξ are exponentially distributed with density
so that the jumps in Stype,k are also exponentially distributed; see Figure 16(a). If S > g f τ, it decays faster than at its average value, and if S < g f τ, the decay is slower, and the average is stable. To get the true diffusive noise for use in (29), we consider fluctuations from the average. Setting Speak ≈ Savg + g, we compute the change in S from immediately before one input spike to immediately before the next as
These fluctuations can be positive or negative, depending upon the interspike interval.
Solving (45) for ξ in terms of η and differentiating, we obtain
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
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)(Vj − VE) (cf. (9)), and once g is factored out of Δ(gS) it is multiplied by (Vj − VE). Once fluctuations are averaged over a population, this product must use the population average voltage , and Δ(S) is multiplied by gAMPA,ext,k( − VE) = JAMPA,ext,k (cf. (10)) as seen in (29).
The Central Limit Theorem  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 Stype,j/I in (17)-(19) are obtained by averaging over populations of Nj and NI cells. The result of n fluctuations drawn from the zero-mean distribution of Figure 16(b) (right) entering each cell has variance n Var(η), so that, summing over N cells, the population average changes by a random quantity with zero mean and variance nVar(η)/N2. Since n = N f Δt, the Gaussian approximation is much better justified, and we may model the effect of a single fluctuation on the population-averaged synaptic variable as a zero mean Gaussian with variance
For large N this is small, so that Savg for a single cell will be close to g f τ, justifying our assumption prior to (45).
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
From the previous paragraph, this summed noise is normally distributed, so it satisfies the requirements of a Wiener process . 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
which may be integrated by parts to give
The expression for variance,
is integrated by parts twice, first to yield
and then finally
The firing rates of populations 3 and I must be expressed as functions of Sj, ν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
where ϕj,0 is the minimum firing rate for Isyn,j < Ithresh,j and Ithresh,j is the threshold current for population j.
To match the f-I curve ϕE,b in the critical region below the decision threshold (20 Hz), we set Ithresh,p = 0.384 nA for pyramidal cells and Ithresh,I = 0.29 nA for interneurons (see black curves in Figures 3(a) and 3(b)). We approximate ψ3(ϕ3) by linearizing (38) about a value ϕ3,0:
(Here we use ϕ3,0 = 1 Hz as for ϕE,b in the four-population model, but the derivation is general.) If Isyn,3 > Ithresh,p, then ϕ3 and the approximation (56) both vary linearly with changes in Isyn,3, and the same is true for ϕI, so if both input currents are above their respective thresholds, increments in current and firing rate are related by
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 Isyn,j > Ithresh,j, the step function θ disappears and (55) simplifies to
Define and as the solutions to this linear system. If and , 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 and from (55), using the components of the currents Isyn,3 and Isyn,I, except for the inputs from populations 1 and 2:
where Γ = a1b2 − a2b1. If (58) holds, the effects of Sj and νj can be calculated separately and then summed. The base currents to populations 1 and 2 due to and are then constant and can be calculated as
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 Sj and νj affect the currents from populations 3 and I to each other and back to populations 1 and 2 as calculated in the next section.
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 S1, which influence the four populations through the following NMDA currents:
Since populations 1 and 2 are explicitly tracked in (30)-(33), we need only calculate the changes in ϕ3 and ϕI due to S1 and the resulting currents back to populations 1 and 2. Setting Δϕ3(S1) to be the change in ϕ3 as a function of S1 and similarly for ϕI, we get
which imply that
where Γ = a1b2 − a2b1. The effects on population 1 are now calculated as follows, starting with INMDA,3,1:
The other currents due to S1 are calculated similarly, yielding and , with
Due to symmetry, the currents to population 2 from populations 3 and I are identical to those to population 1 and the effect of S2 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 Sj in the α’s, but with JNMDA,p/IS1 replaced by .
We are now in a position to compute the coefficients αj, βj, and Iconst,j in (34)-(35). We start by writing α1 = α1a + α1b + α1c, in which α1a is the effect of Si directly on itself, α1b is the effect of Si on itself through the change in the AMPA and NMDA currents of population 3, and α1c is the effect of S1 on itself through the GABA-mediation of population I.
The coefficient α2 is computed in the same manner:
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 S1 ↔ S2 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 νj. These are similar to the αj’s, but with JNMDA,p/I replaced by JAMPA,p/ITAMPA/1000 where appropriate. We find
Currents of the form in (34)-(35) can be calculated provided that ϕI and ϕ3 each stay on a single side of the thresholds in their f-I curves. For example, provided Isyn,3 < Ithresh,p for all time (at a given location in the neuromodulation plane), regardless of inputs from populations 1 and 2, ϕ3 remains in its constant regime and can be fixed at ϕ3,0. If Isyn,3 > Ithresh,p for all time, regardless of inputs from populations 1 and 2, ϕ3 remains in its purely linear regime and the approximations of sections C.1–C.3 hold. Similar conditions apply for ϕI, and there are then four possibilities for which the parameters of (34)-(35) can be calculated: (A) Isyn,I and Isyn,3 above threshold for all time, (B) Isyn,I above threshold and Isyn,3 below threshold for all time, (C) Isyn,I and Isyn,3 below threshold for all time, and (D) Isyn,I below threshold and Isyn,3 above threshold for all time (this turns out to be impossible in the present network structure). The derivation in sections C.1–C.3 assumes case (A), and three requirements must be met for this assumption to hold, which then implies linearity in (60)-(61):
(the latter must hold if ϕ3 and ϕI are to be well defined). The values for and are calculated in section C.1 as they vary with γE and γI scaling the effective currents Jtype,k, and the left panel of Figure 17 shows the regions in the neuromodulation plane over which these inequalities are satisfied. The red curve defines the boundary to the right of which crosses ϕI,0 and enters its linear regime, and enters its linear regime below the black curve. These constraints are stricter than Γ ≠ 0, which is defined by the thick blue curve. The condition (black curve) is the strictest constraint, and case (A) applies only below this curve. The region of validity for case (A), then, does not overlap with any part of the region of interest in which RR > 0, and the α’s and β’s derived in section C.3 do not apply in the region of good behavior of the higher-dimensional models.
Above the black curve, is fixed at ϕ3,0 and 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 , 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 Isyn,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 , but if Isyn,3 were to increase above Ithresh,p due to inputs from populations 1 and 2, the assumption of constant ϕ3 corresponding to case (B) would fail. An input to population 3 from the selective populations corresponds to an input to the interneurons in the ratio JNMDA,p : JNMDA,I ≈ 1.3 : 1, and we can therefore estimate their effects via the derivative
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 dϕ3/dI < 0 and , ϕ3 will remain at this value, regardless of the inputs from the selective populations, which have a net suppressive effect, and case (B) holds for all time. Thus, in the region above the black curve in the right panel of Figure 17 and right of the magenta curve in the left panel, which contains most of the region of good performance, case (B) applies and we may repeat the two-population model derivation with ϕ3 = ϕ3,0 held constant, as described in section C.5.
For the small portion of the region of good performance in which dϕ3/dI > 0, Isyn,3 is far enough below threshold that although Isyn,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 Isyn,3 to reach Ithresh,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 Isyn,3 > Ithresh,p Isyn,I > Ithresh,I, as Ithresh,p > 1.32Ithresh,I and Isyn,3 < 1.32Isyn,I due to network parameters. The lower-right transition from to < 0 occurs as Γ = 0. This zero of the determinant Γ is responsible for the chain of erratic model behaviors seen in the lower right region of Figure 10.
We now solve for , the baseline value of ϕI, in case (B) of section C.4 (Isyn,I > Ithresh,I and Isyn,3 < Ithresh,p). This is done by computing the firing rate given the constant inputs, including the now-constant inputs from population 3:
and IAMPA,ext,I is the same as in the last equation of (62). Setting , we obtain
The effect of population I can be calculated as follows:
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 = ϕI,0 for all time, and every trial is a no-choice trial.
*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).