JN AJP: Gastrointestinal and Liver Physiology
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
 QUICK SEARCH:   [advanced]


     


J Neurophysiol 95: 1099-1114, 2006. First published October 19, 2005; doi:10.1152/jn.00491.2005
0022-3077/06 $8.00
This Article
Right arrow Abstract Freely available
Right arrow Full Text (PDF)
Right arrow All Versions of this Article:
95/2/1099    most recent
00491.2005v1
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Right arrow Citation Map
Services
Right arrow Email this article to a friend
Right arrow Similar articles in this journal
Right arrow Similar articles in ISI Web of Science
Right arrow Similar articles in PubMed
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via ISI Web of Science (4)
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Miller, P.
Right arrow Articles by Wang, X.-J.
Right arrow Search for Related Content
PubMed
Right arrow PubMed Citation
Right arrow Articles by Miller, P.
Right arrow Articles by Wang, X.-J.

Power-Law Neuronal Fluctuations in a Recurrent Network Model of Parametric Working Memory

Paul Miller and Xiao-Jing Wang

Volen Center for Complex Systems, Brandeis University, Waltham, Massachusetts

Submitted 11 May 2005; accepted in final form 14 October 2005


    ABSTRACT
 TOP
 ABSTRACT
 INTRODUCTION
 METHODS
 RESULTS
 DISCUSSION
 APPENDIX
 GRANTS
 ACKNOWLEDGMENTS
 REFERENCES
 
In a working memory system, persistent activity maintains information in the absence of external stimulation, therefore the time scale and structure of correlated neural fluctuations reflect the intrinsic microcircuit dynamics rather than direct responses to sensory inputs. Here we show that a parametric working memory model capable of graded persistent activity is characterized by arbitrarily long correlation times, with Fano factors and power spectra of neural activity described by the power laws of a random walk. Collective drifts of the mnemonic firing pattern induce long-term noise correlations between pairs of cells, with the sign (positive or negative) and amplitude proportional to the product of the gradients of their tuning curves. None of the power-law behavior was observed in a variant of the model endowed with discrete bistable neural groups, where noise fluctuations were unable to cause long-term changes in rate. Therefore such behavior can serve as a probe for a quasi-continuous attractor. We propose that the unusual correlated fluctuations have important implications for neural coding in parametric working memory circuits.


    INTRODUCTION
 TOP
 ABSTRACT
 INTRODUCTION
 METHODS
 RESULTS
 DISCUSSION
 APPENDIX
 GRANTS
 ACKNOWLEDGMENTS
 REFERENCES
 
Neuronal activity in the cerebral cortex is highly irregular (Buzsaki 2004Go; Compte et al. 2003Go; Shadlen and Newsome 1994Go Softky and Koch 1993Go). Fluctuations in spike discharges are often correlated between simultaneously recorded cells on a trial-by-trial basis (Bair et al. 2001Go; Gawne and Richmond 1993Go; Lee et al. 1998Go; Zohary et al. 1994Go); this covariance greatly impacts the extent to which noise can be averaged out over a neuronal pool and so determines the accuracy with which a stimulus feature can be extracted from a large population of neurons (Abbott and Dayan 1999Go; Averbeck and Lee 2004Go; Shadlen and Newsome 1996Go; Shamir and Sompolinsky 2004Go; Sompolinsky et al. 2002Go). Correlations arise from common inputs caused by both overlapping afferents and local synaptic interconnections. For example, cells in the primary visual cortex become correlated because of shared thalamic inputs as well as lateral connections within the cortex. Exactly how correlated noise enhances or decreases the coding efficiency for an oriented visual stimulus critically depends on the (feedforward vs. recurrent) network architecture (Seriès et al. 2004Go).

Correlated fluctuations provide a valuable probe into the dynamic nature of a neural network. For responses of sensory cells to external stimuli, the effects of direct sensory inputs and local connectivity are confounded. In contrast, neurons in working memory systems exhibit persistent activity in the absence of external stimulation, when the subject must hold information transiently in the mind. Hence fluctuations of persistent activity (Compte et al. 2003Go; Constantinidis and Goldman-Rakic 2002Go) are believed to result from the intrinsic dynamics of a working memory microcircuit. The cellular mechanisms of mnemonic persistent activity represent a topic of intense current interest (Brody et al. 2003Go; Major and Tank 2004Go; Wang 2001Go). Of special interest are working memory circuits that encode an analog quantity, such as those integrating a velocity signal into a persistent change in position or direction (Aksay et al. 2000Go; Taube and Bassett 2003Go) or storing in active short-term memory the spatial location or temporal frequency of a sensory stimulus (Funahashi et al. 1989Go; Romo et al. 1999Go). After a transient input, neurons in such a network exhibit sustained changes of activity that are tuned with a stimulus feature in a graded manner. This suggests the system could possess a continuous family of stable self-sustained activity states. Neuronal tuning curves that arise from such activity states are often bell-shaped functions or monotonic functions of the encoded quantity.

For a parametric working memory system, if the delay firing rates of two monotonically tuned neurons are plotted against each other, a curved line results, as shown in Fig. 1A. Such a curved line is referred to as a line attractor (Seung 1996Go). The realization of a continuum of monotonically tuned persistent firing states requires fine-tuning of network parameters (Miller et al. 2003Go; Seung et al. 2000bGo). In an alternative scenario, the network could display multiple persistent states, not truly a continuum, by virtue of a set of discrete stable activity patterns that can be robustly realized without fine-tuning of parameters (Goldman et al. 2003Go; Koulakov et al. 2002Go). Different stimuli cause the network to change discontinuously from one pattern of activity to another, whereas background noise does not cause a switch in the state of persistent activity. Hence the discrete network model has a limited sensitivity to small differences in stimuli, but is robust to noise in the circuit. Whether the continuous or discrete model better describes neural integrators remains an open question (Major and Tank 2004Go).


Figure 1
View larger version (26K):
[in this window]
[in a new window]
 
FIG. 1. Network properties. A: schematic representation of a line attractor. Top: tuning curves of 3 different neurons. Bottom: representation of the continuous attractor in the planes of pairs of neurons. Arrow indicates direction in which fluctuations are not rapidly damped in time, because they shift the system along the continuous attractor of the network. B: schematic model architecture. There are 2 networks of positively and negatively monotonic neurons, respectively. Each network has 12 excitatory pyramidal cell populations (squares) and 12 inhibitory interneuron populations (circles). Synaptic connections are stronger within the same population than between populations. Connectivity is asymmetrical, so that the activation threshold, {Theta}, by stimulus is the lowest for neural population 1 and highest for neural population 12. The 2 networks interact through pyramid-to-interneuron connections, resulting in cross-inhibition. Stimulus inputs are given equally to all excitatory cells within a network, but they differ across the 2 oppositely tuned networks.

 
Thus far, most experimental and computational studies have been concerned with the trial-averaged firing rates. Fluctuation analysis offers a different approach to characterize and potentially test distinct working memory models. In a continuous attractor network, noise is integrated in time, so the mnemonic activity could drift in the same manner as a random walk. Because correlation functions for a random walk do not decay exponentially with time (Mandelbrot and Ness 1968Go), we reasoned that unusual fluctuation properties should be a conspicuous feature of continuous "parametric" memory networks, in which mnemonic neural firing is monotonically tuned to an analog quantity (Aksay et al. 2000Go; Miller et al. 2003Go; Romo et al. 1999Go; Seung 1996Go; Seung et al. 2000bGo). A more robust discrete integrator network (Goldman et al. 2003Go; Koulakov et al. 2002Go) would not show the same random drifts in response to noise, unless the noise were strong enough to cause transitions between the multiple discrete states (Miller 2006Go). If the scale of noise fluctuations sets the coarseness of the system, a system is considered as quasi-continuous whenever the difference in discrete rates is less than the noise-induced variation in rate during a trial.

The purpose of this study is to analyze temporal fluctuations in firing patterns of working memory models with monotonically tuned persistent states. We show that power-law fluctuations are a salient feature of continuous attractor networks.


    METHODS
 TOP
 ABSTRACT
 INTRODUCTION
 METHODS
 RESULTS
 DISCUSSION
 APPENDIX
 GRANTS
 ACKNOWLEDGMENTS
 REFERENCES
 
Computational network model

We developed a cortical microcircuit model for the task of parametric working memory, as published previously (Miller et al. 2003Go). The task requires macaque monkeys to encode and maintain in memory a vibrational frequency, before comparison with a second vibrational frequency. In the secondary somatosensory cortex of a monkey, the spiking response of a neuron to the vibrotactile stimuli varies either positively monotonically or negatively monotonically with the stimulus frequency. Neurons in the prefrontal and premotor cortices, which receive inputs from the secondary somatosensory cortex, show tuned mnemonic activity that can be sustained for up to 6 s during a delay after the end of the initial stimulus (before a 2nd comparison stimulus). Our model network is intended to reflect the activity of prefrontal cortical neurons. We included two sets of populations, corresponding to sets of neurons that receive inputs separately from the secondary somatosensory cortex. One set of inputs is positively monotonically tuned while the other set is negatively monotonically tuned to the vibrational stimulus frequency.

We used interconnected integrate-and-fire neurons, grouped into populations that had strong intrapopulation connections and weaker connections between populations. The connection strengths decrease exponentially with the difference between population numbers (Fig. 1B). The decay in connection strength was more gradual from high-to-low population numbers, so that the neurons in populations with low index received more excitatory current from other cells than those with higher index. This led to a range of thresholds for responses to external input, with the low-index populations being the most excitable and having the highest spontaneous rates.

We couple the two sets of oppositely tuned populations by cross-inhibition. The positively monotonic populations are labeled 1 for the most excitable to 12 for the least excitable, and the negatively monotonic neurons are labeled 1* for the most excitable to 12* for the least excitable. Cross-inhibition is not uniform, but strongest between populations labeled by i and (14 – i)*, where i runs from 2 to 12. Hence the increase in activity of a highly excitable positively monotonic population, after a weak stimulus, helps to reduce the activity of a less excitable negatively monotonic population. A larger stimulus results in a strong increase in activity for the less excitable positively monotonic populations and a concurrent decrease in activity of the more excitable negatively monotonic populations.

We realized a continuous attractor, with moderate recurrent excitation, where the total synaptic weight from all other excitatory populations is of the same order as the self-excitation within a population. For each population, the background external input amplitude and intrapopulation recurrent excitation strength were adjusted to a particular combination, corresponding to a point called a "cusp" in the two-dimensional parameter space, indicated by C in Fig. 2A. For a strength of recurrent excitation greater than at the cusp, bistability can exist between distinct up (persistently active) and down (resting) states, with a gap between them (A in Fig. 2A). For weaker recurrent excitation than at the cusp, the population firing rate varies smoothly with the excitatory drive (B in Fig. 2A). At the cusp, there is a possibility of an approximately vertical line segment, indicating a range of stable firing rates for the population with a fixed background excitatory input (position C, Fig. 2A, right).


Figure 2
View larger version (33K):
[in this window]
[in a new window]
 
FIG. 2. Continuous or discrete attractor networks. A: schematic phase diagram showing possible stable states for a single population as a function of excitatory drive (x-axis) and strength of intrinsic excitatory feedback (y-axis). Boundary of a bistable region, where 2 stable states exist with a gap between them (A) is given by a cusp (C). For recurrent excitation weaker than at the cusp, the system varies continuously from low to high activity (B). For recurrent excitation greater than at the cusp, bistability exists in a range of external drive. B: trial-averaged neural firing rate for a positively monotonic excitatory cell in the network with quasi-continuous states. Solid bar represents cue duration of 1 s. Colors from blue through green to red represent increasing stimulus frequency. C: trial-averaged neural firing rate for a positively monotonic excitatory readout cell in the network with discrete states. Same representation as B. D: tuning curves showing average firing rate during the delay for 3 neurons from different populations in the quasi-continuous model. The points X, Y, and Z mark the stimuli used in Fig. 3B. E: tuning curves with average firing rate in the delay for 3 neurons from different bistable populations in the network with discrete states (open symbols) and readout population (filled symbols). Only the readout cells are designed to have realistic properties that can be compared with experimental data.

 
We generated a set of discrete attractors for robust working memory (Koulakov et al. 2002Go) by creating many bistable populations with a range of thresholds. We achieved this by increasing the recurrent connection strength within each population while reducing the connection strengths between populations (see point A in Fig. 2A) and increasing the average leak conductances. The bistable populations each converge onto a readout population of 400 neurons, whose firing rates encode in a step-like manner the number of active bistable populations. We use the statistical properties of the readout cells for a fair comparison between the discrete model and the continuous model. The set of bistable populations should just be considered as one particular mechanism for generating the more realistic properties of the readout cells, which are suitable for experimental comparison.

We chose stimulus strengths for the discrete attractor, so that each stimulus would cause the discrete system to be in the same state for all trials with that stimulus. In this way, we can use the discrete system as a control to show that typical networks with a single stable state after a stimulus do not exhibit the same type of fluctuations seen in the continuous attractor network.

Single neuron parameters

Individual neurons are simulated using the single-compartment leaky-integrate-and-fire model (Tuckwell 1988Go), such that the membrane potential, Vi, of cell, i, follows the current-balance equation

Formula 1(1)
where Cm is the total membrane capacitance, gL is the leak conductance, VL is the leak potential, gE and VE are the conductance and reversal potential for excitatory channels and gI and VI are the conductance and reversal potential for inhibitory channels, respectively. gext and gcue are the fixed conductances for background noisy input and applied, stimulus-dependent input, respectively, whereas sext and scue are the corresponding time-dependent gating variables (see Eqs. 5 and 6). When the membrane potential reaches a threshold, Vthr, the neuron spikes, and the membrane potential is reset at Vreset for an absolute refractory period, {tau}ref, before continuing to follow Eq. 1.

The total synaptic drive for excitation or inhibition (SE or SI) is the sum of synaptic inputs from all presynaptic neurons j

Formula 2(2)
where Wj->i is the relative synaptic weight from cell j to cell i, and sj is the synaptic current gating variable activated by the presynaptic neuron j firing spikes at times tspike,j. Specifically, for excitatory synapses, we have

Formula 3(3)
and for inhibitory synapses

Formula 4(4)
with synaptic time constants {tau}s. The probability of vesicular release, <PR(t)>, is described in the next subsection.

Background noisy input to all neurons is simulated using uncorrelated Poisson spike trains at a rate, rext, through nonsaturating synapses, of conductance gext, which are gated according to

Formula 5(5)
with synaptic time constant {tau}ext following spikes at times, tspike,ext.

Similarly, during the stimulus, Poisson spike trains of rate, {lambda}, generate additional excitation through {alpha}-amino-3-hydroxy-5-methylisoxazole-4-propionic acid receptor (AMPAR)-mediated synapses of conductance, gcue, multiplied by a gating variable, scue, which follows

Formula 6(6)
The rate, {lambda}, represents the combined spike rate of multiple afferents from secondary somatosensory cortex, and increases or decreases linearly with stimulus frequency for the positively monotonic or negatively monotonic populations, respectively.

In the network models presented here, background and stimulus inputs are mediated by {alpha}-amino-3-hydroxy-5-methylisoxazole-4-propionic acid (AMPA) receptors, with {tau}ext = 2ms, recurrent excitation through N-methyl-D-aspartate (NMDA) receptors (Wang 1999Go, 2001Go) with {tau}s = 100 ms and VE = 0 mV, and inhibition through GABAA receptors with {tau}s = 10 ms and VI = –70 mV. In the continuous attractor network, the cellular parameters are as follows for excitatory cells: Cm = 0.5 nF, gL = 38.4 nS, VL = –70 mV,Vreset = –60 mV, Vthr = –45 mV, {tau}ref = 2 ms, gext = 6 nS, rext = 1.2 kHz, gE = 36 nS, gI = 12 nS. For inhibitory cells, they are as follows: Cm = 0.2 nF, gL = 17.6 nS, VL = –70 mV, Vreset = –60 mV, Vthr = –50 mV, {tau}ref = 1 ms, gext = 1.6 nS, rext = 1.8 kHz,gE = gcue = 36 nS, gI = 12 nS. The discrete integrator network has a range of leak conductances for the excitatory cells, equally spaced from gL = 30.4 nS for cells in population-1 to gL = 40 nS for cells in population-12. The inhibitory neurons have gL = 20 nS, gext = 3 nS, and rext = 1.0 kHz. Otherwise single-cell properties are the same as in the continuous attractor network.

Short-term plasticity of excitatory synapses

All excitatory synapses exhibit short-term presynaptic facilitation and depression (Hempel et al. 2000Go; Varela et al. 1997Go). We implement the scheme described by Matveev and Wang (2000)Go, which assumes a docked pool of vesicles containing neurotransmitter, where each released vesicle is replaced with a time constant, {tau}d. The finite pool of vesicles leads to synaptic saturation, because when the presynaptic neuron fires more rapidly than vesicles are replaced, no extra excitatory transmission is possible. Such synaptic depression contributes to stabilizing persistent activity at relatively low rates.

We assume that there is at most one vesicle release per spike; hence the release probability at any individual synapse, PR(t), is

Formula 7(7)
where pv(t) is the release probability for an individual vesicle and N(t) is the number of docked vesicles (smaller than a maximum N0). We make the simplification that there are many synapses between each pair of connected neurons, such that the average release probability per synapse, <PR(t)>, simply scales the amplitude of synaptic transmission, as shown in Eq. 3. Similarly, we do not keep track of a discrete N(t) for every individual synapse, but assume that the several synapses between two neurons have a binomial distribution with an average number of docked vesicles, <N(t)> (where the brackets <> represent the average over this binomial distribution, with mean <N(t)>) and maximum N0. Hence <N(t)> is a continuous variable obeying

Formula 8(8)
decreasing by <PR(t)> after a spike at time tspike. By averaging over the binomial distribution, we have

Formula 9(9)

Formula 10(10)
and this value of <PR(t)> is used in Eq. 3.

The vesicular release probability is given by the product of three gating variables, pv(t) = O1(t)O2(t)O3(t). A gating variable Oi(t) (i = 1,2,3) increases because of calcium influx triggered by an action potential, followed by a decay with time constant {tau}fi between spikes. Specifically, the following simple update rule is used: a gating variable Oi(t) (i = 1,2,3) follows

Formula 11(11)

Our simulations use the following values for the parameters in the continuous network: N0 = 16, {tau}d = 0.5 s, Cf1 = 0.45, {tau}f1 = 50 ms,Cf2 = 0.75, {tau}f2 = 200 ms, Cf3 = 0.9, {tau}f3 = 2 s. Excitatory neurons in the discrete attractor network are identical except with {tau}d = 0.1 s.

Network connectivity

Connection strengths between neurons depend only on their group numbers and are all-to-all. All weights are normalized by (i.e., divided by) the number of neurons in the presynaptic group, so that average network properties should be independent of the system size.

The set of weights WEI, WIE, WII all follow the same form

Formula 12(12)
where Ngrps = 12 is the total number of groups used, WmaxEI is the maximum connection strength between groups of the same label (i = j), and {sigma}EI determines the breadth of connections to other groups. WmaxIE, {sigma}IE, WmaxII, and {sigma}II have similar definitions.

The recurrent excitation has a slightly different form. First, the connections within the same group are significantly stronger than those between groups, so we define a separate set of parameters for the Wi->iEE =Wi. Second, the connection strengths between different groups i and j are asymmetric

Formula 13(13)
for i > j and

Formula 14(14)
if i < j, where AEE is an asymmetry factor. The continuous attractor network has AEE = 1.5 so that connections are stronger from higher to lower threshold neurons.

The cross-inhibition, WcrossEI, is the strength of connection from each inhibitory population, labeled by i to the excitatory population of opposite tuning labeled by 14 – i for 2 ≤ i ≤ 12.

The full set of parameters are as follows for the continuous network:

W0EE = 0.16, W1 = 0.244, W2 = 0.239, W3 = 0.237, W4 = 0.238, W5 = 0.239, W6 = 0.24, W7 = 0.241, W8 = 0.242, W9 = 0.243,W10 = 0.244, W11 = 0.245, W12 = 0.246; {sigma}1 = 0.5, {sigma}2 = 0.4, {sigma}3 = 0.39, {sigma}4 = 0.385, {sigma}5 = 0.385, {sigma}6 = 0.388, {sigma}7 = 0.392, {sigma}8 = 0.397, {sigma}9 = 0.402, {sigma}10 = 0.408, {sigma}11 = 0.414, {sigma}12 = 0.42; AEE = 1.5;WmaxEI = 1.65, {sigma}EI = 0.25, WmaxIE = 0.5, {sigma}IE = 0.2, WmaxII = 2.0, {sigma}II = 0.5. WcrossEI = 0.25.

For the discrete network: W0EE = 0.14, W1 = 0.35, W2 = 0.365,W3 = 0.378, W4 = 0.39, W5 = 0.401, W6 = 0.412, W7 = 0.423,W8 = 0.434, W9 = 0.445, W10 = 0.455, W11 = 0.465, W12 = 0.475; {sigma}1 = 10, {sigma}2 = 10, {sigma}3 = 10, {sigma}4 = 10, {sigma}5 = 10, {sigma}6 = 10, {sigma}7 = 10, {sigma}8 = 10, {sigma}9 = 10, {sigma}10 = 10, {sigma}11 = 10, {sigma}12 = 10; AEE = 1.0;WmaxEI = 0.3, {sigma}EI = 0.4, WmaxIE = 0.3, {sigma}IE = 0.4, WmaxII = 0.5, {sigma}II = 0.5; WcrossEI = 0.25.

The readout cells of the discrete network are excited by cells from the populations, i, with weights WiER as follows: W1ER = 0.45, W2ER = 0.4, W3ER = 0.35, W4ER = 0.4, W5ER = 0.25, W6ER = 0.2, W7ER = 0.2, W8ER = 0.2, W9ER = 0.2, W10ER = 0.2, W11ER = 0.2, W12ER = 0.2.

Covariance functions

In the following sections, a key quantity that enters the calculations and affects all the statistics is the covariance in spike rates, covij(r) (t1, t2) between neurons i and j (which can refer to the same neuron with i = j) at times t1 and t2. The special case covii(r) (t1, t2) is the variance in spike rates across trials for a neuron i at time t1. When analyzing data, we count the number of spikes by neuron, i, in trial, n, in specific bins of width {delta}t (where {delta}t = 25 ms unless otherwise stated) between t = (k – 1/2){delta}t and t = (k + 1/2){delta}t as Nin(k). Hence the relevant covariance in spike count becomes

Formula 15(15)
where we have substituted for the trial-averaged count, Formula 15, explicitly in the final line. We achieve the average across different trials, n, in our computer simulations by using different initializations of the random number generator.

Correlation functions and correlograms

We calculate the unnormalized cross-correlogram by averaging the covariance function over the measurement interval according to the formula (Bair et al. 2001Go; Brody 1998Go, 1999Go)

Formula 16(16)
The time lag is {tau} and the total temporal window of data measurement is T. Both T and t are integral multiples of {delta}{tau}. The integration window for detection of spikes offset by {tau} is limited to T{tau}. The above definition is for positive t. For negative {tau}, the summation limits are from k = –{tau}/{delta}{tau} to T/{delta}{tau} and prefactor is 1/(T + {tau}) such that Cij(–{tau}) = Cji({tau}).

We normalize the covariance functions by dividing by the geometric mean of the variance at the two time-points (Brody 1999Go). This leads to the correlation coefficient, rij(t,t + {tau}), between spike times t and t + {tau} of neurons i and j

Formula 17(17)

We average over time, t

Formula 18(18)
to produce the correlograms in Fig. 3, A and C.


Figure 3
View larger version (28K):
[in this window]
[in a new window]
 
FIG. 3. Long-time correlograms. A: autocorrelogram measured during delay period for 3 populations of the quasi-continuous network (black, blue, and green) and 1 from the discrete network (red). Central dip is caused by the refractory period. Note that the neurons in the network with quasi-continuous states have a non-zero correlation for large time-lags, decaying slowly, approximately linearly with time. Inset: autocorrelogram reaches 0 after a few tens of milliseconds for a neuron in the discrete network. B: size of offset in the unnormalized autocorrelogram is cue-dependent and is greatest at the steepest parts of the tuning curve for any particular neuron. Stimuli correspond to those marked X, Y, and Z in Fig. 2D. C: cross-correlograms between neurons of the same population (2–2), two closely tuned positively monotonic populations (2–3), a low-threshold and high-threshold population (2–12) and between a positively and negatively tuned population (2–12*). The lack of time-symmetry between oppositely tuned populations-2 and 12* is a sign that in this case, the firing rate of the negatively tuned population (12*) drifts for several seconds before affecting the firing rate of the positively tuned population (2).

 
Figure 4B is calculated using a variant of Eq. 16 with T {tau} replaced by T' to remove any artificial dependence on {tau} because of its inclusion in the integration limit. Different values of T' produce the four different curves in Fig. 4B.


Figure 4
View larger version (30K):
[in this window]
[in a new window]
 
FIG. 4. Random walk behavior after the 14-Hz cue. A: trial-to-trial variance of average firing rate in population-2 of the continuous attractor network increases approximately linearly with time (black), but with an offset (at t = 0) indicating significant trial-to-trial variation in the responses to the stimulus by the end of the cue. Mean rate for the population following the cue varies between 8 and 11 Hz. Readout population from the discrete attractor shows a low, constant variance (red). Inset: each thin line represents the population-average firing rate (with trial-averaged mean rate subtracted) for an individual trial after the end of the stimulus (10 trials are shown of 50 used to calculate the variance for the population in the continuous attractor network). B: dependence of the unnormalized autocorrelogram on the measurement interval, T'. We vary the integration window for spikes as a function of {tau} to maintain a constant range, T', for the first spike in the pair. Hence the window for a pair of spikes is equal to T' + {tau}. The 4 curves in ascending order are for values of T' = 4, 8, 12, and 16 s such that the autocorrelation increases monotonically with T'. For a perfect random walk, curves would be independent of {tau} and proportional to T'. Results are average from 100 neurons (25 in each of 4 populations) over 50 trials for a single cue (black) and from 25 cells in the readout population of the discrete network (red). C: log-log plot of the power spectrum for cells from population-2 of the continuous attractor network, shows decay at low frequencies close to a power-law (close to linear on the log-log plot). Filled circles for frequencies, f = {omega}n/2{pi}, with odd n, and open circles for even n ({omega}n = n{pi}/T with T = 10 s). Data are calculated from the spike trains of 25 cells of each population to produce the figure. Inset: power spectrum showing low-frequency behavior for cells from 1 population (black circles), with the readout population of the discrete attractor as comparison (red crosses). D: Fano factor as a function of time. Dark lines, for neurons in the quasi-continuous network the variance in spike count normalized by the mean spike count (Fano factor) increases with time during the delay after cues that give a large offset in the autocorrelogram; red, neurons in the discrete network have approximately constant Fano factors over the same time interval, with smaller values for cues leading to higher firing rates and more regular interspike intervals. Analysis predicts that a random-walk of the rate leads to a Fano factor that increases quadratically with time (dashed gray line).

 
Power spectrum

The appropriate power spectrum for a nonstationary process is the Wigner-Ville spectrum (Flandrin 1989Go), which is well established in signal processing and physics (Mallat 1999Go). It is defined as a function of frequency, {omega}, and time, t, as

Formula 19(19)
where cov(t1,t2) is the covariance function (between 2 times, t1 and t2) described in an earlier subsection. The integration over {tau}, the time-difference ({tau} = t2 t1) leads to a result independent of t, the mean time [t = (t1 = t2)/2] for a stationary process.

For a nonstationary process, we can calculate a power spectrum in frequency alone, P({omega}), independent of time, t, by averaging the Wigner-Ville spectrum across the measurement period (Flandrin 1989Go) as

Formula 20(20)
In practice, to obtain good statistics from noisy spike trains, such a temporal average is essential.

Calculation of power spectrum from spike trains

In Fig. 4C we plot the temporally averaged Wigner-Ville spectrum, P({omega}), for neurons in our simulations, using a measurement interval of T = 10 s. Specifically, we calculated the average power spectrum by combining spikes from different neurons in a population and using Eqs. 15, 19, and 20 to give

Formula 21(21)
where {omega}n = (n{tau})/T, and the sums are over all binned spikes [where k is the time index for the bin that contains spikes at times t such that k – ({delta}t)/2 < t ≤ k + ({delta}t)/2)] from Npop = 25 different neurons from the same population (with identical network inputs, but different background noise) in trial-{lambda} out of Ntrials. In the log-log plot of Fig. 4C, we show values for low-n up until statistical noise causes the power to be negative for some data points (in which case we cannot take the logarithm) and estimate the power coefficient, {alpha}, from a straight line fit of log-power versus log-frequency through the points shown.

Fano factor

The Fano factor, F(T), is a measurement of trial-to-trial variability of the spike count of an individual neuron. It is obtained from the variance in spike count during a measurement window of temporal length, T, and is expressed as a function of time

Formula 22(22)

Noise correlation

We evaluate the noise correlation by counting the total number of spikes during a delay of 6 s after the cue offset, for a given neuron, i, in the nth trial, as Ni,n. The noise correlation, Xij, between two neurons, i and j, is defined as (Lee et al. 1998Go)

Formula 23(23)
where <Ni> = ({Sigma}nNi,n)/Ntrial is the trial-averaged number of spikes for the ith neuron. Clearly Xii = 1, the maximum possible correlation when i = j, and in general, –1 ≤ Xij ≤ 1. The noise correlation is calculated for each cue and averaged across cues in Fig. 5, A and B.


Figure 5
View larger version (27K):
[in this window]
[in a new window]
 
FIG. 5. Noise correlation of persistent activity. A: normalized correlations in the spike counts for delay activity after all cues. Correlations of 3 positively monotonic excitatory neurons with 2 neurons from each population are shown. Note that correlations are greatest between cells of the same population (points neighboring the same-cell peak of magnitude 1) or between cells in populations with similar tuning curves, and are negative between oppositely tuned cells. B: solid or dashed black curves: distributions of noise correlations in the continuous attractor network between cells of which the 2 tuning curves have the same signs of slope (predominantly positive correlations) or the opposite signs of slope (negative correlations), respectively. Gray: 2 distributions of noise correlations for cells in the discrete attractor network, with means of +0.001 and –0.001 that are barely distinguishable. C: covariance in total spike count of pairs of neurons, as a function of the product of the slopes of their tuning curves. Data are binned along the x-axis, with mean and error bars of each bin plotted. Data are taken during the delay period, for each cue, and between all pairs of neurons with 2 representatives from each population. We fitted a tanh function to the tuning curve of average rate as a function of stimulus and used the gradient of this curve at each stimulus as df/ds. Theory predicts that the 2 quantities are proportional to each other, such that the points would fall about a straight line through the origin. Steeper gray line is a fit through data in the positive quadrant (for neurons with the same sign of tuning curve). Notably, the 2nd gray line, which fits the data in the negative quadrant, is less steep.

 
The total spike count covariance, plotted in Fig. 5C, is simply the numerator in the preceding function, that is

Formula 23


    RESULTS
 TOP
 ABSTRACT
 INTRODUCTION
 METHODS
 RESULTS
 DISCUSSION
 APPENDIX
 GRANTS
 ACKNOWLEDGMENTS
 REFERENCES
 
We examined spike-time correlations in a recurrent network model of noisy spiking neurons for somatosensory parametric working memory (Miller et al. 2003Go). The model was designed to simulate prefrontal cortical neurons in monkeys during a somatosensory delayed discrimination task (Romo et al. 1999Go). In this task, the monkeys were trained to discriminate the frequencies of two vibrotactile stimuli, presented before and after a delay period of 3–6 s, so that the behavioral response (a lever press to indicate whether the 1st or 2nd stimulus was at a higher frequency) depended on remembering the first stimulus frequency across the delay period. It was discovered that neurons in prefrontal and premotor cortices exhibited persistent delay activity that was tuned monotonically to the initial frequency. Such neurons maintain the necessary information for the animal to perform the delayed comparison task.

The model network contains 12,000 neurons that fire spontaneously and stochastically (with a CV averaging >1 at firing rates <10 Hz, dropping to 0.5 as firing rates increase to 40 Hz) because of noisy excitatory inputs. The neurons are arranged in two sets of 12 neuronal populations. In line with the experimental data (Romo et al. 1999Go), the activity of one set of populations increases with the stimulus frequency, whereas the activity of the other set decreases with increasing stimulus frequency. Hence the tuning of populations is positively monotonic or negatively monotonic, respectively. Recurrent excitation dominates within each set of populations, whereas the two oppositely tuned sets are connected by cross-inhibition (see Fig. 1B). The results comparing the mnemonic activity of our model network with the delay activity of monkey neurons are published elsewhere (Miller et al. 2003Go).

In our model, after a transient stimulus, neurons are able to fire persistently, with slow drift, over a range of firing rates. The feedback to each population is tuned (schematically to position C in Fig. 2A) so that each population is on the cusp of bistability, enabling the system to be close to possessing a continuous attractor (Aksay et al. 2000Go; Durstewitz 2003Go; Loewenstein and Sompolinsky 2003Go; Miller et al. 2003Go; Seung 1996Go; Seung et al. 2000bGo). The delay activity after different strengths of stimulus is shown in Fig. 2, B and D. Neurons are able to fire persistently over a wide range of rates. Initial stimuli of different intensities cause a transient increase in firing rates, before the network settles into a state of approximately constant firing rate, but with fluctuations caused by noise. We label the system in this case as quasi-continuous.

Alternatively, we can adjust parameters so that each population provides strong excitatory feedback to itself and receives less excitation from other sources (position A in Fig. 2A). In such a case, the population can be strongly bistable, often with a large gap in firing rates between a down state and an up state. If groups of neurons are bistable, but have differing excitabilities or thresholds, the network can possess a number of robust, discrete stable states (Goldman et al. 2003Go; Koulakov et al. 2002Go), with the number of states approximately equal to the number of bistable groups. Figure 2, C and E, shows the persistent activity for such a multistable network. The possible firing rates of an individual neuron have a large gap, between states where the rate is almost zero, and states where the firing rate is high (Fig. 2E, open symbols). Such a large gap in firing rates results from the strong recurrent feedback within a population, so that once the neurons in a population are able to fire a little, they excite each other with strong synaptic input, causing increasing activity in an escalating manner. Synaptic depression and receptor saturation limit the rise in firing rate to a persistent level of 40–50 Hz in our model. Any transition between the multiple discrete states in the network involves a discrete jump in the rate of at least one population of neurons, from their spontaneous value to their lowest persistent rate (across the gap in Fig. 2A). Such a large jump in firing rates of neurons in a population between stable states makes the discrete system more stable to noise fluctuations (Koulakov et al. 2002Go).

To produce more realistic cells in the discrete attractor network, to be compared with those in the continuous attractor network, we added a readout population to the discrete network. The readout population receives excitatory input that is a weighted sum from all populations of the same sign of tuning in the network. The tuning curve of the readout population (Fig. 2E, solid symbols) more closely resembles those of neurons in the continuous network and experimental data (Romo et al. 1999Go). Cells from the readout population also have a similar, rate-dependent CV (although about 10% lower) to neurons in the continuous attractor network. In all statistical analyses, we use such cells from the readout population of the discrete network.

Hence using the same framework of interconnected populations of neurons, we compared the noise and correlations of these two types of memory system. Our main results are independent of details of the model. However, they rely on noise moving the system along a continuous range of states but being unable to cause large jumps in activity between discrete states. A set of discrete attractors with small gaps in rate between states, or with large noise fluctuations, so that noise could cause spontaneous transitions between the discrete states, would show the same behavior we describe here for a continuous attractor (Miller 2006Go).

Random walk of rate

If a continuous attractor underlies our parametric working memory model, fluctuations in the state of the network are described by a random walk. To show this conceptual point, consider a simplified, linear model of a typical cell in a neural group whose firing rate r(t) obeys the following dynamic equation

Formula 24(24)
where GI is the mean firing rate produced by the input current I, {eta}(t) is an uncorrelated noise, Formula 24 is the noise amplitude, and {tau}s is a synaptic or membrane time constant (a few to tens of milliseconds). Noise leads to an autocorrelation which decays exponentially

Formula 24
Parametric working memory requires neurons to display sustained changes of firing activity by feedback mechanisms either within a cell (Camperi and Wang 1998Go; Egorov et al. 2002Go; Goldman et al. 2003Go; Loewenstein and Sompolinsky 2003Go) or through a reverberatory network (Koulakov et al. 2002Go; Miller et al. 2003Go; Seung 1996Go; Seung et al. 2000bGo). In either case, a positive feedback implies a component of the input current that increases with its own firing rate, so making the simplification of linear feedback, I = I0 + Wr, the above equation becomes

Formula 25(25)
where {tau}eff(W) = {tau}s/(1 – GW) > {tau}s. If there is no feedback, W = 0, and in time {tau}s, the firing rate reaches its steady-state value, GI0, as in Eq. 24. With strong feedback (GW -> 1), one has {tau}eff >> {tau}s. In the limit {tau}eff -> {infty}, the first term in the right hand side of the equation vanishes. The firing rate integrates the input I as well as noise over time, in the sense of calculus: Formula 25. Moreover, after the input is withdrawn, the firing rate is maintained stationary at any level (within an operating range) except for noise-induced drifts; hence the system behaves like a line attractor. The average rate is given by < r > = GI0 tI/{tau}s, where tI is the time when the input is withdrawn, whereas the trial-to-trial variance in rate increases linearly with time: <(r<r>)2> = At.

The analogy between the random walk model and our parametric working memory model lies in the fact that the continuum of persistent activity states form a one-dimensional curved line in the space of population firing rates, as schematically shown in Fig. 1A. If the firing rates are temporarily perturbed by noise to move off this line, the network's dynamics cause the firing rates to return back to the line, but do not prevent any drifts along the continuous attractor. Such noise-induced drifts along the attractor are similar to a random walk. To assess whether the simple random walk model can quantitatively capture the behavior of our large-scale and highly nonlinear network model, we analyze the fluctuation properties of population firing rates and spike counts.

Random walk analysis of a line attractor

A closer analysis of the effects of noise on a system with a continuous attractor reveals that the correlation functions are determined by the products of gradients of the tuning curves of two neurons (Ben-Yishai et al. 1995Go; Pouget et al. 1998Go). The argument is shown in Fig. 1A, where the top three figures depict the tuning curves (average firing rate as a function of stimulus) for three neurons. The first two are positively monotonic; the third is negatively monotonic. When the system has a continuous attractor, the only fluctuations that are not rapidly damped out are those along the attractor. The attractor can be seen as a curve in a plot of the firing rate of one neuron versus another, with each point in the curve corresponding to the two firing rates of persistent activity after a particular stimulus. Noise has the same effect as the stimulus in shifting the firing rates along the attractor. Two positively monotonic neurons [with firing rates r1,r2 and tuning curves f1(s),f2(s)] either both increase or decrease their firing rates together, after a larger or smaller stimulus, and after noise fluctuations. Hence their noise correlations are expected to be positive.

If we consider the firing rate of mnemonic persistent activity as a function of stimulus, s, and noise, {eta}(t), expanded to first order in the noise, we have

Formula 26(26)

Formula 27(27)
The noise correlation and unnormalized cross-correlation functions are proportional to the covariance of firing rates, Formula 27, which is given by the product of the noise terms in Eqs. 26 and 27. Hence the cross-correlation is proportional to the product of gradients of the two tuning curves [(df1/ds). (df2/ds)] (Ben-Yishai et al. 1995Go; Pouget et al. 1998Go). If the gradients of tuning curves of two neurons have opposite sign [as in Fig. 1A; (df1/ds)(df3/ds) < 0], the noise correlation is negative between the corresponding neurons. That is, a fluctuation along the attractor simultaneously causes the positively monotonic neurons to increase their firing rates and negatively monotonic neurons to decrease their rates or vice versa. Similarly, for an individual neuron, the covariance function is proportional to (dfi/ds)2, so all effects of noise are greater after cues that fall on a steep part of a neuron's tuning curve.

Noise during the stimulus presentation leads to variation in the initial encoded values of s and has a correlated effect on the firing rate of neurons in the same way as noise during the delay. The main effect of both stimulus noise and drift noise in a line attractor is encoded in the one-dimensional variable, s. The firing rate of each neuron responds to the network noise by an amount proportional to the square of the gradient of its tuning curve. In the preceding analysis, we assumed that fluctuations are small compared with changes in the gradient of the tuning curve, so dfi/ds remains constant for each neuron for a complete trial, allowing Eqs. 26 and 27 to be expanded only to first order. If the first-order approximation is valid, the proportionality constant, A, is the same for both stimulus noise and drift noise. That allows us to write the total noise at time t as

Formula 28(28)
for any neuron, i, where t0 quantifies noise during the stimulus and is a constant for all neurons that form the line attractor and Ai {propto} (dfi/ds)2.

Other results for a random walk are described in detail in the APPENDIX. We summarize the key features that can be compared with our network results here. The covariance function, cov(t1,t2), for a random walk starting at t = 0 is independent of the time difference, |t2t1|, and instead equals the variance of the process given above in Eq. 28 where t is the earlier of t1 or t2 (Gillespie 1992Go) (cf. Eqs. A4 and A5).

The other statistical measures presented in this paper are derived from the covariance or variance functions. The time-average of the covariance function is proportional to the integration limits used in the time-averaging (see Eqs. A6A9). The average (Eq. 16), leads to a linear decay in |{tau}|

Formula 29(29)
for the unnormalized correlation function of a random walk, which will be compared with Fig. 3A–C . An alternative method results in curves independent of {tau}, but linearly increasing with T'

Formula 30(30)
for the unnormalized correlation function of a random walk, which will be compared with Fig. 4B.

A hallmark of random walk behavior is an inverse-square power law for the Wigner-Ville power spectrum (Eq. 19) (Flandrin 1989Go). When averaged over a finite time interval, T, specific frequencies, {omega}n = n{pi}/T should be used, so that the inverse-square power law is apparent in the averaged spectrum. We find that initial variance in the random walk leads to oscillations in the power spectrum between even and odd frequencies (Eqs. A11 and A12), such that

Formula 31(31)
for even n and

Formula 32(32)
for odd n.

The Fano factor, F(T), for a random walk, also shows power-law behavior, increasing quadratically with time. Initial variance in the starting points contributes a linear term, such that for a Poisson-like renewal process

Formula 33(33)
In the following subsections, we present the results of our spiking network and compare them with those expected for a random walk in a line attractor.

Correlation function

Fluctuations are commonly quantified by autocorrelation functions of single neurons and cross-correlation functions between them (Bair et al. 2001Go; Brody 1999Go; Perkel et al. 1967aGo,bGo; Singer and Gray 1995Go). The autocorrelation function is a function of two times (1 for each spike) and should only be written as a function of the time-difference, or time-lag, {tau}, if the underlying process is stationary (Gardiner 1985Go). Many processes measured in neuroscience are not stationary, so the shuffle-correction is used to remove the principle effects of systematic variations in the average rate. We use the term correlogram to denote the shuffle-corrected covariance function, normalized by the product of SD in spike count (producing a correlation coefficient), averaged across the measurement interval (Brody 1999Go). We use the term unnormalized correlogram when referring to a temporal average of the covariance function.

Trial-to-trial differences in the underlying rate have been shown to give rise to artifacts in cross-correlograms at short time scales (Brody 1998Go, 1999Go). In contrast, we report that the correlation functions are affected on long time scales by the linear increase in variance of firing rate as a function of time that arises from the temporal integration of noise in a quasi-continuous attractor.

The autocorrelograms for three cells in our graded memory network are shown in black, green, and blue in Fig. 3A. Strikingly, temporal correlations between spikes persist for many seconds and seem to show a gradual, linear decrease with the time lag. In contrast, a network with robust, discrete states does not show any long-term correlations (Fig. 3A, red curve). In fact, the autocorrelograms for all neurons fall to zero during a time scale on the order of the synaptic time constant (100 ms for NMDA receptors in this work; Fig. 3A, inset).

The autocorrelogram in Fig. 3A was calculated by averaging over memory activity after different stimuli. If we compute the correlations separately after each stimulus, we find that only a small subset of stimuli result in such large, long-term persistence in the correlations. This is because, if the stimulus lies on a relatively flat part of a neuron's tuning curve, small network fluctuations do not cause significant changes in the firing rate of that neuron. Tuning curves of neurons in our network are not linear in general, but resemble sigmoidal functions, as seen for the majority of neurons in the experimental data (Miller et al. 2003Go; Romo et al. 1999Go). It is the steepest part of the tuning curve that is most sensitive to fluctuations. If the tuning curve of a neuron is linear (unlike those seen in our model), the effect of fluctuations during the delay should be independent of stimulus.

In Fig. 3B, we calculated the unnormalized autocorrelogram separately for delay activity after three stimuli, corresponding to the positions X, Y, and Z in Fig. 2D. To obtain better statistics, we averaged cross-correlations between 25 neurons in a single population. Because the population fires asynchronously and all neurons within one population receive the same network input, their autocorrelation functions are essentially the same. The average cross-correlation between all pairs of neurons in the same population is identical to the average autocorrelation of those neurons, except at very short time scales (less than {approx}25 ms). It is only the {delta}-function peak at {tau} = 0 and the dip in the autocorrelation at short {tau} (because of the refractory period and recovery from reset; Fig. 3A) that distinguishes autocorrelations from cross-correlations between cells within a population of our network.

Figure 3C shows that the cross-correlogram is greater between neurons with similar tuning curves (such as 2 neurons in the same population, 2–2, or neighboring population, 2–3) than between neurons with very different thresholds (2–12). The cross-correlogram is negative between neurons of opposite tuning (2–12*) and is temporally asymmetric. The asymmetry is a sign that firing rates of one population can drift for several seconds before the change in rate is strong enough to affect the populations with opposite sign of tuning. This is because cross-inhibition is not the dominant feature in our network and suggests that, on a time scale of a few seconds, our system resembles more closely two weakly coupled continuous attractors (where drift in 1 attractor provides input to the other) than a single set of stable points.

In all cases, autocorrelograms and cross-correlograms do not decay exponentially as exp(–t/{tau}) with a characteristic time constant, {tau} (see the text after Eq. 24), but decay more linearly, as is typical for a random walk.

Increasing variance of rate

We found that the variance of single-trial population firing rates (calculated with 100 trials) increases approximately linearly from the stimulus offset for up to 10 s of the mnemonic period (Fig. 4A) . The nonzero intercept at t = 0 reflects trial-to-trial noise in the system's response to the stimulus.

The linear time course with nonzero intercept for the firing rate variance is consistent with a random walk (see Eqs. A4 and A5 and Noise during the stimulus), with an initial distribution of starting points caused by variability of neuronal response during the stimulus presentation (Eq. A2). Note that the variance of neurons in the discrete system remains constant and low (Fig. 4A, red curve).

Covariance independent of time lag

Such a random walk also provides an explanation for the observed near-linear behavior of the correlograms reported in Fig. 3. Because random walks are nonstationary processes (the variance increases with time), correlation functions depend separately on the two times of measurement (t1,t1 + {tau}) and not just on their difference, {tau} (see METHODS). Hence when evaluating correlations as a function of the time lag, {tau}, by integrating over t1, the precise measurement interval affects the result. It can be shown mathematically that a linear increase in the variance of rates as a function of time leads to a linear increase in the unnormalized correlogram as a function of the integration interval (Miller 2006Go; Saleh 1978Go). In the standard formulation, the integration interval is (T – |{tau}|), where {tau} is the time-lag and T is the total time of measurement. With fixed T, this leads to a linear decrease with {tau}, as seen in Fig. 3B.

In Fig. 4B, we tested whether the unnormalized autocorrelogram is truly a function of the measurement interval using an alternative integration window to calculate the temporal average (the sum over k in Eq. 16). We calculate the unnormalized autocorrelogram by including spikes over a range of T' + {tau} for each value of {tau} so that the integral over t1 is over a fixed value, T' = T {tau}max, rather than T{tau}. Hence the measured correlation interval is independent of {tau} and proportional to T'. We plotted curves for the average unnormalized autocorrelogram of 100 neurons (25 neurons in each of 4 populations) using increasing values of T' (4, 8, 12, and 16 s for the ascending curves). For a perfect random walk, the curves in Fig. 4B would be constant with {tau} and equally spaced along the y-axis (see METHODS). The curves do have a y-offset that increases monotonically with T', whereas they vary little with {tau}, so the network's activity is behaving qualitatively like a random walk during the delay.

Power law of power spectra

Furthermore, for a random walk, the power spectrum, P({omega}n), of the spike train should include a delta-function at the origin and at low frequencies an inverse-square law decay, such that P({omega}n) {propto} {omega}n{alpha} where a equals 2 (Flandrin 1989Go) and {omega}n = n{pi}/T. A pure power law decay should appear as a straight line on a log-log plot, with the slope of the line giving the exponent (Fig. 4C). We estimated the exponent, {alpha}, for several cells by fitting the linear region of log-log power spectra. We fitted curves separately through odd and even data points, because these were offset due to noise during the stimulus (see Eq. A12 and Noise during the stimulus). The fits yielded values ranging from {alpha} = 1.6 to {alpha} = 2.1 ({alpha} = 1.77 and 1.66 for the odd and even frequencies, respectively, of population-2, shown in Fig. 4C). The fact