WIXLIB

I&F point neurons. In ref. 12 we described an asymptoticreduction of this large-scale point-neuron network to a CGkinetic representation that is well suited for scale-up. We showedthat this kinetic theory is dynamically accurate and numericallyefficient, even when the original point-neuron network is fluctuation-dominant. However, the coarse-graining procedure inthe derivation of the kinetic theory includes a local average overtime, thus removing detailed spike-timing information. Here, toregain this spike information, we take the hybrid approach,identifying a distinguished subnetwork of point neurons that willbe retained and replacing the other neurons with a CG kineticdescription.Now, we turn to the description of the essential steps used inthis hybridization procedure. Partition the original network intotwo subpopulations, a ‘‘minority’’ population comprising thedistinguished subnetwork of embedded neurons and a ‘‘majority’’ population that is to be coarse-grained and constitutes a CGbackground patch. The ith neuron in this majority populationsatisfies the I&F dynamics, d G G i i i BD,dt i 共v, gE, gI; t兲dgEdgI ,0 E共v; t兲 gE 共gE, gI兩v; t兲dgEdgI ,0 I共v; t兲 gI 共gE, gI兩v; t兲dgEdgI ,0where (g E, g I兩v; t) is the conditional probability densityfunction, i.e., (v E, g E, g I; t) (g E, g I兩v; t) (v; t). Uponclosure assumptions (12), our CG kinetic theory of the dynamics(1) reduces to closed equations for (v; t), E(v; t), and I(v; t): t v兵 共v, E, I兲 其, t E [1a][2a] E1 关共v V E兲 兴共 E g E共t兲兲 E 共v兲 v v E ,[2b]2[1b]together with the spike-to-reset dynamics. Here, i describesthe external drive, i describes the cortico-cortical couplings among neurons in the CG background, and BDdescribes thefeedback from the remaining, embedded neurons. , E, Ilabel the excitation and inhibition; , E, and I denote themembrane, excitatory, and inhibitory timescales, respectively; Vis the membrane potential; G E and G I are excitatory andinhibitory conductance, respectively, V R, V E, and V I are theresting, excitatory, and inhibitory reversal potentials, respectively. The external input Ei E f E兺 (t T Ei ) uses aPoisson spike train {T Ei } with strength f E and a prescribed ratev 0E(t). S, C, labels simple and complex cells, respectively.Here, simple and complex cells are modeled by the followingEnetwork architecture: ES 1, C 0, i.e., simple cells are drivenexternally with excitatory inputs in addition to their inputs fromother cells, whereas complex cells are not driven externally andreceive inputs only from other cortical cells.¶ In our model, forthe sake of presentation, half of the neurons are simple, the other half are complex, unless otherwise stated. BDare the inputs(feedback) to the neuron of type ( , ) in the background fromthe neurons in the embedded subnetwork, which is symbolizedby D.Next, we replace these ‘‘background neurons’’ of Eq. 1 with aCG kinetic representation. As described in ref. 12, the CGkinetic theory studies the evolution of theprobability densityN function: (v, g E, g I; t) [(1兾N )兺 i 1{ [v V i(t)] [g E G Ei (t)] [g I G Ii (t)]}], where the expectation is taken overall realizations of incoming Poisson spike trains {T Ei }, and overpossible random initial conditions. N is the number of neuronsin the background population of ( , )-type. Define the marginaland conditional moments:¶Note that, in this article, we merely use the terms ‘‘simple’’ and ‘‘complex’’ to refer to thisparticular network architecture. Neurons in the visual cortex are classified (16) ‘‘simple’’ or‘‘complex,’’ with simple cells responding to visual stimulation in an essentially linearfashion (for example, responding to sinusoidally modulated standing gratings at thefundamental frequency, with the magnitude of response sensitive to the spatial phase ofthe grating pattern), and with complex cells that respond nonlinearly (with a significantsecond harmonic) in a phase-insensitive manner. See refs. 14 and 15 for a model of simpleand complex cells which captures these physiological features.Cai et al. 2d V 共 V i V R兲 G Ei 共t兲共 V i V E兲dt i G Ii 共t兲共 V i V I兲,冕冕冕 t I I1 关共v V I兲 兴共 I g I共t兲兲 I 共v兲 v v I ,[2c]where (v, E, I) 1[(v V R) (v V E) E (v V I) I] and g 共t兲 0 共t兲f S 冘 pm 共t兲 S BD 2 冋1 2 f 2 共t兲 共S 兲2 0 2共S BD兲2 冘 pm̃ 共t兲Ñ ,冘pm̃ 共t兲, 冘 pm 共t兲N 册[3a][3b]where 0 (t) is the (temporally modulated) rate for the externalPoisson spike train for -population, and m (t) is the averagepopulation firing rate per neuron of ( , )-type. m̃ (t) capturesthe feedback from the embedded point neurons (see SupportingText, which is published on the PNAS web site). Ñ is the totalnumber of embedded neurons of ( , )-type in the subnetwork.(Note that all the symbols with a tilde stand for the corresponding variables for the subnetwork of embedded point neurons.) denotes the cortico-cortical synaptic connection strengthS describesfrom -population to -population of -cell type. S BDthe coupling strength from the embedded subnetwork of population to the -population in the background. The stochastic nature of synaptic connection (17) is modeled by a Bernoullirandom variable with release probability p (or synaptic failureprobability 1 p). (For details, see Supporting Text.) Thesekinetic equations are posed for V I ⱕ v ⱕ V T, under two-pointboundary conditions in v, which are derived from the fact thatthe total flux of neurons across firing threshold V T is equal to thetotal flux at reset V R (12), together with vanishing flux at v V I.The dynamics of the remaining, i.e., the embedded, neurons isnow described byPNAS 兩 September 28, 2004 兩 vol. 101 兩 no. 39 兩 14289NEUROSCIENCE 共v; t兲

d Ṽ 共Ṽ i V R兲 G̃ Ei 共t兲共Ṽ i V E兲dt i G̃ Ii 共t兲共Ṽ i V I兲,d G̃ i G̃ i i i DB,dt[4a][4b] i and i describe the external drive to the embeddedwhere neurons and their cortico-cortical couplings, respectively, and DB S DB兺 兺 p̂ (t t̂ ) describes the feedback from theCG populations to the embedded neurons (Supporting Text). Instead of using merely the average contribution to this feedback asdescribed by the mean rate pN m (t), we invoke the Poissonassumption and reconstruct stochastic spiking times {t̂ E } and {t̂ I }by using the total population rates pN E m E(t) and pN I m I (t) as themeans for these Poisson spike trains, to capture spike fluctuations in our hybrid approach. Note that this Poisson spikereconstruction (PSR) automatically accounts for the spike fluctuation in the input via possibly sparse connections from thebackground, even if the number of neurons in the background isvery large. We will demonstrate the significance of this PSRbelow.Eqs. 2 and 4 constitute a network of point neurons (Eq. 4)embedded within, and fully interacting with, a dynamical CGbackground represented by kinetic theory (Eq. 2). The computation effort of this network involves only the simulation of a fewembedded I&F neurons plus three (1 1)-D partial differentialequations for simple and complex cells in the CG background.For statistical information, this approach allows us to drasticallyreduce computational cost otherwise needed in the simulation ofthe original, full I&F network (12) while retaining spike-timinginformation within the CG model.ResultsWe have verified our hybrid approach by comparing the embedded network dynamics with the original full I&F networkwith the same underlying network architecture and parameters.First, we addressed accuracy and present an excellent comparison for systems in which different parts of the total network arereplaced by the kinetic theory description.Fig. 1 shows raster plots (spike times for each neuron) for fivestimulus periods for a single CG patch containing 400 neurons.The raster plot for the full network of I&F neurons is shown inFig. 1a, with the results for three embedded networks shown inFig. 1 b–d. In these embedded networks, the eliminated pointneurons are replaced by the kinetic theory CG background (Eq.2). Shown in Fig. 1 b–d are the spike rasters for the remaining,i.e., embedded I&F neurons, for comparison with the sameneurons in the full network (Fig. 1a). In Fig. 1 b and c,respectively, half complex (simple) excitatory cells are describedby the kinetic theory, whereas in Fig. 1d, half complex and halfsimple cells, again all excitatory, are described by the kinetictheory. Even in this case, the firing patterns of the embeddedneurons (Fig. 1d) agree well with the full I&F network (Fig. 1a).The networks shown in Fig. 1 are operating in a fluctuationdominant regime (12). This regime can be a natural operatingpoint for cortical networks (2, 18–20), whether there are externalstimuli or not. The kinetic theory gives a valid, effective reduceddescription of such situations (12). Fig. 2 accentuates this pointby comparing results for the same full network of 400 neuronsas in Fig. 1, with those of the CG background represented by aCG mean firing-rate representation, m (t), replacing kinetictheory. This mean firing-rate background is obtained by solving14290 兩 01Fig. 1. Comparisons of raster plots of firing dynamics of the original I&Fnetwork, with networks with different components replaced by the kinetictheory description, for a single CG patch containing 400 neurons, 200 of whichare simple (150 excitatory and 50 inhibitory) and 200 complex¶ (150 excitatoryand 50 inhibitory) with network connection probability p 0.25. (The modelnetworks in all figures always have 75% excitatory neurons and 25% inhibitory neurons.) External drive is modeled by Poisson-spike trains with the ratef 0E(t) 20. (1 0.5 sin(2 f t)) with f 10 Hz and f 0.01, five stimulusperiods are shown. Network parameters E 5 ms, I 10 ms, VR 70 mV,VT 55 mV, VE 0 mV, VI 80 mV (the same for all figuresEI IE IIthroughout) and SS (SEES , SS , SS , SS ) (0.2, 0.4, 0.2, 0.4), and SC (0.2, 0.4,0.2, 0.4). [The order convention of (SEE, SEI, SIE, SII) will be used throughout.]Throughout, neurons indexed from 1 to 50 are simple excitatory, and neuronsindexed from 51 to 100 are complex excitatory (inhibitory cells not shown). (a)The full network of I&F neurons. (b) Embedded networks with half of thecomplex cells coarse-grained. (c) Embedded network with half of the simplecells coarse-grained. (d) The embedded network with half of the simple andhalf of the complex cells coarse-grained. All CG populations are described bythe kinetic theory (Eq. 2) with PSR (see text).m 共t兲 log兩1 U 共兵m 其兲兩for U 1;m 共t兲 0, otherwise,[5]for , E, I, S, C, where 1(1 g E g I),U ({m }) (( V T V R )兾( V S V R )), and V S ( ) 1(V R g EV E g IV I), where g are defined in Eq. 3a.Eq. 5 can be derived (2) from the kinetic theory (Eq. 2) by takingthe limit N 3 , and f 3 0, f 0 finite). Fig. 2 b–d clearlydisplay the inaccuracy of the firing patterns of the remainingembedded neurons when compared with the same neurons in thefull I&F network (Fig. 2a). In particular, note the completefailure of this mean firing-rate embedding (Fig. 2 b–d) to capturethe firing of the complex cells in this fluctuation-driven case,whereas the kinetic theory embedding provides a faithful representation of the firing of these complex cells as shown in Fig.1. Fig. 3 further corroborates this point. Shown are dynamicfiring-rate curves for four representations in comparison withthe full I&F network simulations of only simple cells. As clearlyshown in Fig. 3, either purely mean-rate representation (as in theapproach of ref. 11), or the hybrid approach using mean-rate onlyfor the CG background, does not yield a quantitatively accuratedescription. We note that results, as shown in Figs. 1–3, can alsobe viewed as the validation of assumptions used in our CGkinetic theory. We will use only kinetic theory (Eq. 2) below todescribe the CG background, instead of Eq. 5.Our hybrid approach contains a distinct step of capturing theeffects of fluctuations in the feedback to the embedded networkCai et al.

from the CG background, which is accomplished by the Poissonspike reconstruction as described above. Here, we contrast thisstep to a feedback of only the mean rate computed using thegeneral kinetic theory (Eq. 2) (not Eq. 5!) without the PSR; i.e., DBis replaced by S DB兺 pN m (t) in Eq. 4b. One mightbelieve that a coupling of the kinetic theory to embedded-pointFig. 3. Dynamic firing-rate curves for excitatory cells, for four representations in comparison with the full 400 I&F simple cell network simulations (thicksolid line) with network connection probability p 0.25. (i) All excitatory andinhibitory cells are coarse-grained using our kinetic theory (Eq. 2) (thin solidline). (ii) All excitatory and inhibitory cells are coarse-grained using the meanfiring-rate representation (Eq. 5) (dotted line). (iii) All inhibitory cells arecoarse-grained with the kinetic theory and the excitatory cells are modeled asembedded I&F neurons (thick dashed line). (iv) All inhibitory cells are coarsegrained using the mean firing-rate representation (Eq. 5) and the excitatorycells are modeled as embedded I&F neurons (thin dashed line). The externaldrive is modeled by Poisson spike trains with the rate f 0E(t) 13. (1 0.25 sin(2 f t), f 10 Hz, and f 0.01.Cai et al.neurons through the direct use of this modulated mean firingrate would be just as effective. The results shown in Fig. 4establish that it is not, i.e., it does not capture the fluctuationswell enough for an accurate description. Fig. 4a shows firing-ratecurves for excitatory complex neurons as a function of theaverage input conductance G input fE 0E for steady inputs,whereas Fig. 4b shows cycle-averaged firing rate under temporally periodically modulated inputs. In both the steady state anddynamically modulated cases (Fig. 4), the firing rates of embedded neurons compare quantitatively very well with those of thefull I&F network simulation when the PSR is used in combination of kinetic theory, whereas the modeling error in the rate ofthe embedded approach can be rather substantial without theproper fluctuation effect in the feedback from the CG background. Indeed, PSR is needed in fluctuation-dominated systems; thus, elsewhere in this article, we always use PSR.Having addressed the modeling accuracy issues of the hybridapproach, we turn to the demonstration of its full power inmodeling situations when there is a distinguished subnetwork ofneurons which are sparse yet strongly coupled. Here, we consideran idealized network containing such a subnetwork, both toillustrate the necessity of the embedded network to model thissituation efficiently and to contrast with another modelingsituation in which ‘‘test neurons’’ within a kinetic theory canprovide an adequate description.The idealized network constitutes a four-population networkof 400 neurons; each neuron is either simple or complex,excitatory or inhibitory. There are only eight complex cells in thisnetwork, forming a distinct subnetwork that is very stronglycoupled. [Synaptic coupling strengths for these strongly coupledcomplex cells are S̃ C (0.4, 0.4, 0.8, 0.4) in contrast to S S (0.1, 0.2, 0.1, 0.2) for the simple cells in the background andS CDB (0.25, 0.2, 0.25, 0.2) and S SBD (0.025, 0.05, 0.025,0.05) for the couplings between the subnetwork and the CGbackground.] Fig. 5 displays the dynamics of this networkcharacterized by (i) raster plots, (ii) cross-correlation as quantified by the conditional firing rate, i.e., the probability per unittime that the jth neuron fires, given the condition that the ithneuron fires, averaged over all i and j, and (iii) ISI distributions.These characteristics from the simulation of the full I&F network (Lower) are successfully reproduced by the hybrid approach (Upper), when all simple cells are coarse-grained usingPNAS 兩 September 28, 2004 兩 vol. 101 兩 no. 39 兩 14291NEUROSCIENCEFig. 2. Comparisons of raster plots of firing dynamics of the original I&Fnetwork with different components of the total network coarse-grained bythe mean firing-rate description (Eq. 5), for a single CG patch with the samenetwork architecture and parameters as in Fig. 1. Throughout, neurons indexed from 1 to 50 are simple excitatory, and neurons indexed from 51 to 100are complex excitatory (inhibitory cells not shown). (a) The full network of I&Fneurons. (b) Embedded networks with half of the complex cells coarsegrained. (c) Embedded network with half of the simple cells coarse-grained.(d) The embedded network with half of the simple and half of the complexcells coarse-grained. All coarse-grained populations are described by Eq. 5,rather than by kinetic theory (Eq. 2).Fig. 4. Steady-state firing-rate curves as a function of the average inputconductance (a) and cycle-averaged firing-rate curves for the embedded,excitatory complex cells¶ in a four-population model with both excitatory andinhibitory cells of simple and complex types in comparison with the full I&Fsimulations (thick solid line) of 400 neurons with network connection probability p 0.25 (b). Half of the simple cells are coarse-grained using thekinetic theory (Eq. 2) with the Poisson spike reconstruction in the feedback tothe embedded point neurons (thin solid line), in contrast to where only themean rate from the kinetic theory (Eq. 2) is used in the feedback (dashed line).The external drive is modeled by Poisson spike trains with the rate f 0E(t) 13. (1 0.5 sin(2 f t), f 10 Hz and f 0.02. Cycle-averaged firingrate is the rate averaged many periods of the stimulus at 10 Hz by m(t) 兺nNt 1 m(t nT)兾Nt for 0 t T, where T is the period.

Fig. 5. Dynamics of strongly coupled embedded subnetwork on a background under an external input rate that is constant in time. Plotted are: (i)Raster plots, only the eight complex cells¶ in the distinct subpopulation areshown (neurons labeled 1– 6 are excitatory, and those labeled 7– 8 are inhibitory); (ii) cross-correlation (see text); and (iii) ISI distributions. (Lower) Simulation of the full I&F network; (Upper) Hybrid approach when all cells in thebackground are coarse-grained by using the kinetic theory, except for theeight distinct complex cells that are treated as an embedded network of pointneurons.the kinetic theory, and the eight distinct neurons are treated asan embedded network of point neurons. These results demonstrate that our approach can successfully capture high-orderstatistics between two different neurons (e.g., cross-correlation)as well as for the same neuron (e.g., ISI). Instead of this constantexternal stimuli, a more useful and interesting situation fordemonstrating the power of the hybrid approach is when thebackground actually has its own response dynamics to timedependent stimuli. Here, the hybrid network really shows itsability to capture the distinct dynamics of the strongly interactingembedded subnetwork as demonstrated in the Fig. 6 Lower.Clearly, in addition to the cross-correlation induced by theexternal periodic drive (as captured by the test neurons in Fig.6 Upper), the strongly interacting embedded neurons exhibit astronger cross-correlation with more fine structure, which reflects relatively coherent recurrent dynamics of their own subnetwork. Furthermore, had we used merely test neurons toextract spike information of the embedded network, we wouldhave obtained a singled-peaked ISI distribution instead of adouble-peaked ISI distribution for the interacting embeddedneurons (Fig. 6).If one merely desires firing statistics of individual neurons inthe background, test neurons in the embedded subnetwork willextract that information from a background whose full dynamicsis modeled by kinetic theory. Fig. 7 shows excellent agreementbetween the ISI distribution of firing statistics of a sampleneuron in a network, obtained by the full I&F network simulation, and that of a test neuron in the CG background descriptionby kinetic theory.To further illustrate the usefulness of the embedded neurons(actually in the sense of test neurons), we use our hybridapproach to construct a RTC (Fig. 8). RTC methods enable theexperimentalist to probe the dynamical response of a corticalnetwork from a system-analysis viewpoint. Based on ‘‘spiketriggered counting,’’ these methods are most directly modeled byusing detailed spike-timing information and thus are natural forthe applications of embedded point-neuron networks. In one14292 兩 01Fig. 6. Dynamics of strongly coupled embedded subnetwork on a background driven by a time-dependent stimulus. Plotted are: (i) Raster plots, onlythe eight complex cells¶ in the distinct subpopulation are shown (neuronslabeled 1– 6 are excitatory, and those labeled 7– 8 are inhibitory); (ii) crosscorrelation; and (iii) ISI distributions. (Lower) The hybrid approach when allcells in the background are coarse-grained by using the kinetic theory, and theeight distinct complex cells are treated as an embedded network of pointneurons. (Upper) A subnetwork of eight test neurons in interaction with thebackground but without interaction among themselves. The external drive ismodeled by Poisson spike trains with the rate f 0E(t) 20. (1 0.25sin(2 f t), f 1 Hz, and f 0.01.version of these methods (21), the experimentalist sequentiallypresents an anesthetized monkey with a grating stimulus, chosenrandomly at every refresh time T refresh from a collection ofgrating patterns with varying orientations and spatial phases, andmeasures the response of an individual cell in V1 extracellularly.For each spike, one asks at time earlier, what was theorientation and grating angle of the stimulus pattern presented?Running over all spikes, one measures the probability P( , ; )that, at time before each spike, the stimulus had orientation and spatial phase . P( , ; ) provides an estimate of thefirst-order linear response function of the cortical system.Fig. 7. Extraction of ISI distributions of neurons in the background, testneuron method. (a) Test neuron driven by a CG neuronal patch described bythe kinetic theory. (b) Sample neuron from a full network of I&F neurons. Wenote that the firing rate of the neuron in the full I&F network is 46 spikes persec, whereas the test neuron driven by the CG background has 45 spikes persec, with only a 2% of relative error.Cai et al.

We have constructed the RTC for neurons on a ring thatmodels the orientation-tuning dynamics of V1 neurons (9, 12, 22

An embedded network approach for scale-up of fluctuation-driven systems with preservation of spike information David Cai*†, Louis Tao‡, and David W. McLaughlin*§ *Courant Institute of Mathematical Sciences and §Center for Neural Science, New York University, New York, NY 10012; and ‡Department of Mathematical Sciences, New Jersey Institute of Technology, Newark, NJ 07102