Modeling neural activity at the ensemble level
Abstract
Here we demonstrate that the activity of neural ensembles can be quantitatively modeled. We first show that an ensemble dynamical model (EDM) accurately approximates the distribution of voltages and average firing rate per neuron of a population of simulated integrateandfire neurons. EDMs are highdimensional nonlinear dynamical models. To faciliate the estimation of their parameters we present a dimensionality reduction method and study its performance with simulated data. We then introduce and evaluate a maximumlikelihood method to estimate connectivity parameters in networks of EDMS. Finally, we show that this model an methods accurately approximate the highgamma power evoked by pure tones in the auditory cortex of rodents. Overall, this article demonstrates that quantitatively modeling brain activity at the ensemble level is indeed possible, and opens the way to understanding the computations performed by neural ensembles, which could revolutionarize our understanding of brain function.
Contents
1 Introduction
If we observe a fluid at the molecular level we see random motions, but if we look at it macroscopically we may see a smooth flow. An intriguing possibility is that by analyzing brain activity at a macroscopic level, i.e., at the level of neural ensembles, we may discover patterns not apparent at the singleneuron level, that are as useful as velocity or temperature are to understand, and predict, the motion of fluids.
Technology frequently drives science. For instance, thanks to the development of microelectordes in the 1930’s, we now know with exquiste detail computations performed by single neurons. We are now experiencing a dramatic increase in our capacity to monitor the activity of larger and larger populations of neurons with higher and higher spatial and temporal resolution. These new ensemble recordings may soon allow us to uncover crucial computations performed by neural ensembles.
Here we present results of developing and evaluating a mathematical model and estimation methods to characterize the activity of ensembles of neurons from electrophysiological data.
Section 2 reports the evaluation of an ensemble dyanmical model (EDMs). And EDMs is a highdimensional nonlinear dynamical models. To estimate its parameters it is convenient to reduce the number of parameters. We escribe a dimensionality reduction method for EDMs in Section 3. Section 4 presents a maximumlikelihood method to estimate parameters of EDMs. Finally, in Section 5 we show that EDMs can accurately approximate highgamma electroencephalographic activity evoked by pure tones in the auditory cortex of rodents.
2 Building EDMs
We wanted to learn how to build Ensemble Density Models (EDMs), dynamical models of the state variables (e.g., transmembrane potential, time since last spike, etc.) of a population of identical neurons, starting from a dynamical model of the state variables of a single neuron. For Integrate and Fire (IF) model of single neurons, an EDM should provide the ensemble probability density function (pdf) , from which to compute the probability of finding a neuron in the ensemble with with a given transmembrane voltage at time (i.e., ). It should also provide the average firing rate per neuron in the population, . To construct EDMs we chose the methodology described in Omurtag et al. (2000).
To evaluate the EDM we compared its outputs ( and ) with those derived from the direct simulation of a population of 9,000 IF neurons. The value of the density function derived from the direct simulation was the proportion of IF neurons having a voltage at time , and the average firing rate per neuron derived from the the direct simulation was the proportion of cells in the population at time with voltage at threshold.
An EDM is driven by an external excitatory current and modulated by an external inhibitory current. In addition, every cell in the EDM receives inputs from other neurons in the population. A fraction of these inputs is excitatory, and the remainder are inhibitory. These intrapopulation inputs act as feedback mechanisms to the EDM. The mathematical representation of an EDM for a population of IF neurons is given in Equations 1–3, modified from Equations 26, 39 and 47 in Omurtag et al. (2000). External excitatory and inhibitory currents appears as and , respectively, in Equation 3, and excitatory and inhibitory feedback are given by the terms and , respectively, in Equation 3.
(1)  
(2)  
(3)  
We first verified that for a single population the outputs of the EDM matched those of the direct simulations (Section 2.1). We next built a network of excitatory and inhibitory populations, and again compared the outputs of the EDM and those of the direct stimulation (Section 2.2).
2.1 One Population
2.1.1 One Population of Independent Neurons
The top panel in Figure 1 shows the ensemble probability density function, , calculated by integrating the differential equation of an EDM, Equation 1. The bottom panel shows and approximation to this pdf obtained from the histogram of voltages of a direct simulation of a population of 9,000 IF neurons. Note the large similarity of the pdfs at all time points.
The black line in Figure 2 shows the average firing rate per neuron, , calculated using the EDM, Equation 2. The grey line shows the average firing rate per neuron calculated from a direct simulation of a population of 9,000 IF neurons, as the proportion of cells with voltage at threshold. Note the almost perfect match between these firing rates.
Figure 3 is as Figure 2 but for a step input current that jumps from 0 to 800 impulses per second at time . Note that by 0.4 seconds after the step in the input current the average firing rate per neuron has reached a new steady state around 10 impulses per second, as revealed by the EDM (black line in Figure 3) and by the direct simulation of a population of 9,000 IF neurons (grey line, in Figure 3)). Figure 4 shows the pdf of voltages at this new steady state calculated by the EDM, Equation 2 (black line in Figure 4) and approximated using the histogram of voltages at 0.4 ms of a direct simulation of 9,000 IF neurons (grey line in Figure 4).
2.1.2 One Population with Feedback
The previous Figures showed results from the simulation of a single population of independent neurons. Figure 5 is as Figure 2, but shows the average firing rate per neuron in a population where each neuron receives excitatory inputs from ten other neurons in the population (i.e., G=10 and f=1 in Equation 3). For comparison, the dashed line shows the average firing rate per neuron from the population without feedback.
2.2 A Network of Populations
The previous sections evaluated EDMs in a single population of neurons. Here we evaluate EDMs of excitatory and inhibitory populations combined in the network of populations shown in Figure 7. The network is driven by an excitatory input to the excitatory population. This population has excitatory feedback (i.e., G=5 and f=1 in Equation 3), and its average firing rate per neuron output, scaled by a constant , drives the inhibitory population. The inhibitory population has inhibitory feedback (i.e., G=5 and f=0 in Equation 3), and its average firing rate per neuron output , scaled by a constant , modulates the activity of the excitatory population.
Figures 8 and 9 show the activity in the excitatory and inhibitory populations, respectively in the network of Figure 7. The upper panels in these figures show the population activities computed by the EDMs and the bottom panels the activity derived from the direct simulation of 9,000 IF neurons. The blue curves, scaled along the left axis, show the average spike rate per neuron in the population, the magenta and yellow curves, scaled along the right axis, show the excitatory and inhibitory external inputs to the population, respectively, and the magenta and yellow dashed curves, also scaled along the right axis, show the excitatory and inhibitory feedback inputs to the population. We should have performed the direct simulations with a larger number of IF neurons to obtain smoother spike rates and currents. Nevertheless, the activities computed by the EDMs are in close agreement with those derived from the direct simulation.
2.3 Partial Conclusions
Given a dynamical model of a neuron, we now know how to derive an EDM for a population of such neurons. For an IF model of a neuron, here we have shown that EDMs accurately approximate population activity (i.e., the pdf of the transmembrane voltage, , and the average firing rate per neuron, ). The next step in this sub project is to estimate connectivity parameters (e.g., and in Figure 7) from simulated data.
3 Reducing dimensionality in EDMs
In the previous section we showed that Ensemble Density Models (EDMs) accurately approximated the average firing rate per neuron and the probability density function (pdf) of direct simulations of ensembles of integrateandfire (IF) neurons. In networks of EDMs, we want to estimate connectivity parameters and state variables (i.e., the pdfs of the different ensembles) from recorded ensemble firing rates. The state space of each previously reported EDM contained 210 variables. To make the estimation of parameters and state variables in networks of EDMs feasible/efficient it would be helpful to find lowdimensional approximations of EDMs. Here we report the approximation power of one such lowdimensional approximation method. This method was inspired by the moving basis technique in Knight (2000).
3.1 Method to find lowdimensional approximations of EDMs
The evolution of the ensemble pdf, , is given by:
(4) 
where is a differential operator that depends on the stimulus . The normalized voltage in Equation 4 ranges in the unit interval. To numerically solve this equation, we discretize , and , , giving a discretization of Equation 4:
(5) 
where is the matrix representation of the differential operator . Equation 5 is a system of differential equations. The objective of the dimensionality reduction method described below is to approximate the evolution of using a system of differential equations, where . For this, the ensemble pdf is represented in a basis of eigenvectors of the differential matrix , where many coefficients in the new representation can be discarded without much loss in approximation power.
3.1.1 Representing the ensemble pdf in a new basis
Let and be the eigenvectors and eigenvalues, respectively, of :
(6) 
or in matrix notation:
(7) 
where is the nth column of the matrix and is the nth diagonal element of the diagonal matrix .
Assuming the eigenvectors are linearly independent, we represent as:
(8) 
or in matrix notation:
(9) 
(10) 
Using a backward difference to approximate the time derivative of :
(11) 
in Equation 5 we obtain:
(12) 
(13) 
or:
(14) 
Equation 14 provides the evolution of the coefficients in Equation 9. It just expresses the evolution of the ensemble pdf in Equation 5 in another basis. It is an exact formula (i.e., it is not an approximation) and has the same dimensionality as Equation 5. The importance of Equation 14 is that one can approximate the evolution of the ensemble pdf by discarding many components of , as we explain next.
3.1.2 Reducing dimensionality in the new basis
To understand why one can discard many coefficients in the representation of the ensemble pdf of Equation 9 without much loss in approximation power, consider the case where the stimulus, , is constant. In such a case, the eigenvectors and eigenvalues will neither depend on the stimulus; i.e., and . Then Equations 5 and 9 reduce to:
(15) 
and:
(16) 
Taking derivatives in Equation 16 we obtain:
(17) 
(18) 
(19) 
or:
(20) 
with solution:
(21) 
Thus, the absolute value of a low dimensional coefficient evolves as:
(22) 
Because all eigenvalues have negative real part (with the exception of the zero eigenvalue), the coefficients associated with nonzero eigenvalues will decay to zero, and the speed of this decay will be proportional to the absolute value of the real part of the corresponding eigenvalue. Therefore, to achieve dimensionality reduction in Equation 14 we may discard those coefficients associated with eigenvalues with larger absolute value of their real part, since these coefficients will rapidly decay to zero.
3.2 Evaluation of the method
We first study how well lowdimensional EDMs approximate the average firing rate per neuron in direct stimulations (Section 3.2.1) and then how they approximate the ensemble pdf (Section 3.2.2).
These studies were peformed in data simulated from the network of EDMs illustrated in Figure 7. The network was driven by an excitatory sinusoidal input to the excitatory population (shown by the dotted curve and scaled along the right axis in the top panel of Figure 10). The average firing per neuron in this population, scaled by a constant , drove the inhibitory population (shown by the dotted curve and scaled along the right axis in the top panel of Figure 11). In turn, the average firing rate per neuron in the inhibitory population, scaled by a constant , inhibited the excitatory population. Both populations contained feedback (i.e., each neuron received 10 inputs from ten other cells in the same population and 80% of these inputs were inhibitory).
3.2.1 Firing Rates
The top panels of Figures 10 and 11 show the average firing rate per neuron obtained from direct simulation (grey curve), from a fulldimensional EDM (red curve), and from lowdimensional EDMs (the blue, cyan, and red curves correspond to 17, 5, and 1 moving basis, respectively). The full dimensional EDM and its low dimensinal approximations with 17 and 5 moving basis almost perfectly approximate the average firing rate per neuron of the direct simulation. The approximation power of the EDM with only one moving basis is not as good, but it look reasonable.
3.2.2 Ensemble Proability Density Functions
We compared the normalized histogram of number of directly simulated neuron per voltage bin (i.e., the ensemble pdf from direct simulation) with the pdfs calculated with EDMs (i.e., Equation 5). For this we computed at every time step the KullbackLeibler (KL) divergence (in bits) between the pdfs obtained by direct stimulation and those obtained from EDMs. These KL divergences are shown in the bottom panels of Figures 10 and 11.
We see that the pdfs obtained from EDMs were good approximation of the pdfs from the direct simulation at times of large average firing rate per neuron (between 0.73 and 0.85 seconds and between 1.03 and 1.2 seconds in the bottom panel of Figure 10, and between 0.72 and 0.82 seconds and between 1.03 and 1.18 seconds in the bottom panel of Figure 11).
At times of low averaged firing rate per neuron the difference between the pdfs obtained from direct simulation and those obtained from EDMs were an order of magnitude larger for the inhibitory than for the excitatory ensemble. This is probably because the excitatory ensemble was driven by a large and smooth sinusoidal input while the inhibitory ensemble was driven by the weaker and nonsmooth average firing rate per neuron of the excitatory ensemble. Also, at most time points, the larger the number of moving basis in lowdimensional EDMs, the better the EDM pdf approximated the pdf obtained from direct simulation, as can be seen more clearly in the bottom panel of Figure 11.
To try to understand why the pdfs obtained from direct simulation were different from those obtained from lowdimensional EDMs, we plotted these pdfs for the inhibitory ensemble in an interval of low average firing rate per neuron (from 0.81 seconds in Figure 12 to 1.09 seconds in Figure 26). In the transition between weak to zero average firing rate per neuron (from 0.81 seconds in Figure 12 to 0.93 seconds in Figure 18) the lowdimensional pdfs moved faster towards lowers voltage than the pdfs from the fulldimensional EDM and those from direct simulation. Similarly, in the transition between zero to weak average firing rate per neuron (from 0.95 seconds in Figure 19 to 1.09 seconds in Figure 26) the lowdimensional pdfs moved faster towards higher voltages than the pdfs from the fulldimensional EDM and those from direct simulation. This suggests that the moving basis discarded in the lowdimensional approximations of EDMs help prevent the EDM pdf to transition too fast to and away from the pdf corresponding to large average firing rate per neuron.
3.3 Partial conclusions
We have developed a method to reduce the dimensionality of the state space of EDMs. We observed that the pdfs from lowdimensional EDMs move faster to and away from the highaveragefiringrateperneuron pdf than the pdf derived from direct simulation or that from full EDMS. However, these differences did not have a large impact on the average firing rate per neuron produced by lowdimensional EDMs, which almost perfectly approximates the average firing rate per neuron computed by direct simulation.
Our next step is to estimate connectivity parameters and state space variables in networks of EDMs from recorded spike rates. We will first do this with simulated data.
4 Estimating parameters of EDMs
Here we evaluate a maximum likelihood method to estimate the connectivity parameters in the EDMs network in Figure 27. We seek the connectivity parameters for which a set of pairs measurements, , are most probable; ie:
The pair of measurements at time , comprises measurements from the excitatory and inhibitory ensembles; i.e., .
4.1 Noise Model
As a first approximation we make the following three assumptions on the noise of the measurements. These are the assumptions made in previous work by our collaborators (Kostuk et al., 2012).

Noise is independent in time; i.e., .

Noise is independent in the different populations; i.e., .

Gaussian noise with a known variance, , in each population; i.e., , where is the activity generated by the ensemble at time .
With these assumptions the log likelihood function of the model parameters reduces to:
The red lines in Figure 28 plots the activity generated by the excitatory and inhibitory ensemble with connectivity parameters and . The blue lines plot the noisy measurements () that we use below for parameter estimation.
4.2 Optimization Surface
The levelplot in Figure 29 shows the optimization surface for the set of parameters show in the axes. We see that it has a convex shape with peak at the true parameters. Thus, an iterative gradientascent optimization procedure should climb to the maximumlikelihood parameters from any starting set of parameters.
4.3 Gradient of LogLikelihood Function
To compute the gradient of the loglikelihood function at a given set of parameters , one needs to integrate the EDMs to generate at each time step the activities of both ensembles, and , as well as the ensembles pdfs, and . With these quantities in hand, the gradient of the loglikelihood function can be computed recursively in a second integration step (Equation 23).
(23) 
The white arrows in Figure 29 point in the direction of the gradient, and are scaled according to the gradient magnitude. Note that, as expected, the gradient direction is perpendicular to the level lines.
With a convex loglikelihood function, for which we can compute its gradient, we can now use an iterative gradient ascent procedure to maximize this function. The green line in Figure 29 shows a gradientascent trajectory that in only three steps accurately approximated the maximum likelihood parameters.
5 Modeling evoked auditory activity in rodents with EDMs
Here we show that the average firingrate per neuron predicted by EDMs well approximates the evoked highgamma power (HGP) from auditory neurons in response to stimulation with puretones.
5.1 Recordings
We used very highresolution simultaneous surface and laminar recordings (Figure 30) in anesthetized rodents stimulated with pure tones at different frequencies and amplitudes. The blue trace in Figure 31 shows the highgamma power (HGP) evoked by the presentation of tones in one sample electrode of the array. Each vertical dotted color line marks the onset of a tone of a corresponding frequency.
We see that these evoked waveforms are very stereotypical. Some waveforms have a large peak followed by a bump. Other waveforms have two peaks, with a smaller peak preceding a larger one. We also see inf Figure 32 waveforms that have more than one peak preceding a larger one and waveforms having a larger peak followed by a smaller one.
5.2 Qualitative analysis
5.3 Quantitative analysis
Using a similar method as described above to learn connectivity parameters in networks of EDMs, we etimated parameters so that an EDM approximated as close as possible the HGP evoked by the presentation of pure tones to rodents. In total we estimated 13 parameters: three parameters of an EDM, two parameters for its initial condition, four parameters for the sinusoidal excitatory input, and four parameters for the inhibitory input. The recorded and approximated HGP are given by the blue and red solid curves, respectively, in Figure 37.
5.4 Partial Conclusions
These preliminary results show that EDMs, besides accurately approximating the average firing rate of ensembles of simulated IF neurons, well approximate the HGP in the auditory cortex of rats evoked by auditory tones.
6 Conclusions
We have shown that an ensemble model accurately reproduce the probability density function of the transmembrane voltage, as well as the averagefiring rate per neurona, in a large ensemble of integrateandfire simulated neurons (Section 2). We developed and evaluated methods to reduce the dimensionality (Section 3) and estimate parameters (Section 4). Finally we demonstrated the feasibility of EDMs to model the highgamma power evoked by pure tones in the auditory cortex of rodents.
The possiblity of quantitatively model the activity of ensemble of neurons may allow us to uncover fundamental computations performed by neural ensembles.
References
 Knight [2000] B.W. Knight. Dynamics of encoding in neuron populations: some general mathematical features. Neural Computation, 12:473–518, 2000.
 Kostuk et al. [2012] M. Kostuk, B.A. Toth, C.D. Meliza, D. Margoliash, and H.D.I. Abarbanel. Dynamical estimation of neuron and network properties II: path integral monte carlo methods. Biological Cybernetics, 106:155–167, 2012.
 Omurtag et al. [2000] A. Omurtag, B.W. Knight, and L. Sirovich. On the simulation of large populations of neurons. Journal of Computational Neuroscience, 8:51–63, 2000.