text:
# 1 Introduction

## 1 Introduction

On short timescales, one can conceive of a single neuron as a computational device that maps inputs at its synapses into a sequence of action potentials or spikes. To a good approximation, the dynamics of this mapping are determined by the kinetic properties of ion channels in the neuron’s membrane. In the fifty years since the pioneering work of Hodgkin and Huxley, we have seen the evolution of an ever more detailed description of channel kinetics, making it plausible that the short time dynamics of almost any neuron we encounter will be understandable in terms of interactions among a mixture of diverse but known channel types \[Hille, 1992, Koch, 1999\]. The existence of so nearly complete a microscopic picture of single neuron dynamics brings into focus a very different question: What does the neuron compute? Although models in the Hodgkin–Huxley tradition define a dynamical system which will reproduce the behavior of the neuron, this description in terms of differential equations is far from our intuition about—or the formal description of—computation.

The problem of what neurons compute is one instance of a more general problem in modern quantitative biology and biophysics: Given a progressively more complete microscopic description of proteins and their interactions, how do we understand the emergence of function? In the case of neurons, the proteins are the ion channels, and the interactions are very simple—current flows through open channels, charging the cell’s capacitance, and all channels experience the resulting voltage. Arguably there is no other network of interacting proteins for which the relevant equations are known in such detail; indeed some efforts to understand function and computation in other networks of proteins make use of analogies to neural systems \[Bray, 1995\]. Despite the relative completeness of our microscopic picture for neurons, there remains a huge gap between the description of molecular kinetics and the understanding of function. Given some complex dynamic input to a neuron, we might be able to simulate the spike train that will result, but we are hard pressed to look at the equations for channel kinetics and say that this transformation from inputs to spikes is equivalent to some simple (or perhaps not so simple) computation such as filtering, thresholding, coincidence detection or feature extraction.

Perhaps the problem of understanding computational function in a model of ion channel dynamics is a symptom of a much deeper mathematical difficulty. Despite the fact that all computers are dynamical systems, the natural mathematical objects in dynamical systems theory are very different from those in the theory of computation, and it is not clear how to connect these different formal schemes. Finding a general mapping from dynamical systems to their equivalent computational functions is a grand challenge, but we will take a more modest approach.

We believe that a key intuition for understanding neural computation is the concept of feature selectivity: While the space of inputs to a neuron—whether we think of inputs as arriving at the synapses or being driven by sensory signals outside the brain—is vast, individual neurons are sensitive only to some restricted set of features in this vast space. The most general way to formalize this intuition is to say that we can compress (in the information theoretic sense) our description of the inputs without losing any information about the neural output \[Tishby et al., 1999\]. We might hope that this selective compression of the input data has a simple geometric description, so that the relevant bits about the input correspond to coordinates along some restricted set of relevant dimensions in the space of inputs. If this is the case, feature selectivity should be formalized as a reduction of dimensionality \[de Ruyter van Steveninck and Bialek, 1988\], and this is the approach we follow here. Closely related work on the use of dimensionality reduction to analyze neural feature selectivity has been described in recent papers \[Bialek and de Ruyter van Steveninck, 2002, Sharpee et al., 2002a\].

Here we develop the idea of dimensionality reduction as a tool for analysis of neural computation, and apply these tools to the Hodgkin–Huxley model. While our initial goal was to test new analysis methods in the context of a presumably simple and well understood model, we have found that the Hodgkin–Huxley neuron performs a computation of surprising richness. Preliminary accounts of these results have appeared \[Agüera y Arcas, 1998, Agüera y Arcas et al., 2001\].

## 2 Dimensionality reduction

Neurons take input signals at their synapses and give as output a sequences of spikes. To characterize a neuron completely is to identify the mapping between neuronal input and the spike train the neuron produces in response. In the absence of any simplifying assumptions, this requires probing the system with every possible input. Most often, these inputs are spikes from other neurons; each neuron typically has of order $`N10^3`$ presynaptic connections. If the system operates at 1 msec resolution and the time window of relevant inputs is 40 msec, then we can think of a single neuron as having an input described by a $`4\times 10^4`$ bit word—the presence or absence of a spike in each 1 msec bin for each presynaptic cell—which is then mapped to a one (spike) or zero (no spike). More realistically, if average spike rates are $`10\mathrm{s}^1`$, the input words can be compressed by a factor of ten. In this picture, a neuron computes a Boolean function over roughly 4000 variables. Clearly one cannot sample every one of the $`2^{4000}`$ inputs to identify the neural computation. Progress requires making some simplifying assumption about the function computed by the neuron so that we can vastly reduce the space of possibilities over which to search. We use the idea of dimensionality reduction in this spirit, as a simplifying assumption which allows us to make progress but which also must be tested directly.

The ideas of feature selectivity and dimensionality reduction have a long history in neurobiology. The idea of receptive fields as formulated by Hartline, Kuffler, and Barlow for the visual system gave a picture of neurons as having a template against which images would be correlated \[Hartline, 1940, Kuffler, 1953, Barlow, 1953\]. If we think of images as vectors in a high dimensional space, with coordinates determined by the intensities of each pixel, then the simplest receptive field models describe the neuron as sensitive to only one direction or projection in this high dimensional space. This picture of projection followed by thresholding or some other nonlinearity to determine the probability of spike generation was formalized in the linear perceptron \[Rosenblatt, 1958, Rosenblatt, 1962\]. In subsequent work Barlow and others \[Barlow et al., 1964\] characterized neurons in which the receptive field has subregions in space and time such that summation is at least approximately linear in each subregion but these summed signals interact nonlinearly, for example to generate direction selectivity and motion sensitivity. We can think of Hubel and Wiesel’s description of complex and hypercomplex cells \[Hubel and Wiesel, 1962\] again as a picture of approximately linear summation within subregions followed by nonlinear operations on these multiple summed signals. More formally, the proper combination of linear summation and nonlinear or logical operations may provide a useful bridge from receptive field properties to proper geometric primitives in visual computation \[Iverson and Zucker, 1995\]. In the same way that a single receptive field or perceptron model has one relevant dimension in the space of visual stimuli, these more complex cells have as many relevant dimensions as there are independent subregions of the receptive field. While this number is larger than one, it still is much smaller than the full dimensionality of the possible spatiotemporal variations in visual inputs.

The idea that neurons in the auditory system might be described by a filter followed by a nonlinear transformation to determine the probability of spike generation was the inspiration for de Boer’s development \[de Boer and Kuyper, 1968\] of triggered or reverse correlation. Modern uses of reverse correlation to characterize the filtering or receptive field properties of a neuron often emphasize that this approach provides a “linear approximation” to the input/output properties of the cell, but the original idea was almost the opposite: Neurons clearly are nonlinear devices, but this is separate from the question of whether the probability of generating a spike is determined by a simple projection of the sensory input onto a single filter or template. In fact, as explained by Rieke et al. (1997), linearity is seldom a good approximation for the neural input/output relation, but if there is one relevant dimension then (provided that input signals are chosen with suitable statistics) the reverse correlation method is guaranteed to find this one special direction in the space of inputs to which the neuron is sensitive. While the reverse correlation method is guaranteed to find the one relevant dimension if it exists, the method does not include any way of testing for other relevant dimensions, or more generally for measuring the dimensionality of the relevant subspace.

The idea of characterizing neural responses directly as the reduction of dimensionality emerged from studies \[de Ruyter van Steveninck and Bialek, 1988\] of a motion sensitive neuron in the fly visual system. In particular, this work suggested that it is possible to estimate the dimensionality of the relevant subspace, rather than just assuming that it is small (or equal to one). More recent work on the fly visual system has exploited the idea of dimensionality reduction to probe both the structure and adaptation of the neural code \[Brenner et al., 2000a, Fairhall et al., 2001\] and the nature of the computation which extracts the motion signal from the spatiotemporal array of photoreceptor inputs \[Bialek and de Ruyter van Steveninck, 2002\]. Here we review the ideas of dimensionality reduction from previous work; extensions of these ideas begin in Section 3.

In the spirit of neural network models we will simplify away the spatial structure of neurons and consider time dependent currents $`I(t)`$ injected into a point–like neuron. While this misses much of the complexity of real cells, we will find that even this system is highly nontrivial. If the input is an injected current, then the neuron maps the history of this current, $`I(t<t_0)`$, into the presence or absence of a spike at time $`t_0`$. More generally, we might imagine that the cell (or our description) is noisy, so that there is a probability of spiking $`P[\mathrm{spike}\mathrm{at}t_0|I(t<t_0)]`$ which depends on the current history. The dependence on the history of the current means that the input signal still is high dimensional, even without spatial dependence. Working at time resolution $`\mathrm{\Delta }t`$ and assuming that currents in a window of size $`T`$ are relevant to the decision to spike, the input space is of dimension $`D=T/\mathrm{\Delta }t`$, where $`D`$ is often of order 100.

The idea of dimensionality reduction is that the probability of spike generation is sensitive only to some limited number of dimensions $`K`$ within the $`D`$ dimensional space of inputs. We begin our analysis by searching for linear subspaces, that is a set of signals $`s_1,s_2,\mathrm{},s_K`$ which can be constructed by filtering the current,

$$s_\mu =_0^{\mathrm{}}𝑑tf_\mu (t)I(t_0t),$$
(1)
so that the probability of spiking depends only on this small set of signals,

$$P[\mathrm{spike}\mathrm{at}t_0|I(t<t_0)]=P[\mathrm{spike}\mathrm{at}t_0]g(s_1,s_2,\mathrm{},s_K),$$
(2)
where the inclusion of the average probability of spiking, $`P[\mathrm{spike}\mathrm{at}t_0]`$, leaves $`g`$ dimensionless. If we think of the current $`I(t_0T<t<t_0)`$ as a $`D`$ dimensional vector, with one dimension for each discrete sample at spacing $`\mathrm{\Delta }t`$, then the filtered signals $`s_i`$ are linear projections of this vector. In this formulation, characterizing the computation done by a neuron involves three steps:

1. Estimate the number of relevant stimulus dimensions $`K`$, hopefully much less than the original dimensionality $`D`$.
2. Identify a set of filters which project into this relevant subspace.
3. Characterize the nonlinear function $`g(\stackrel{}{s})`$.

The classical perceptron–like cell of neural network theory would have only one relevant dimension, given by the vector of weights, and a simple form for $`g`$, typically a sigmoid.

Rather than trying to look directly at the distribution of spikes given stimuli, we follow de Ruyter van Steveninck and Bialek (1988) and consider the distribution of signals conditional on the response, $`P[I(t<t_0)|\mathrm{spike}\mathrm{at}t_0]`$, also called the response conditional ensemble (RCE); these are related by Bayes’ rule,

$$\frac{P[\mathrm{spike}\mathrm{at}t_0|I(t<t_0)]}{P[\mathrm{spike}\mathrm{at}t_0]}=\frac{P[I(t<t_0)|\mathrm{spike}\mathrm{at}t_0]}{P[I(t<t_0)]}.$$
(3)
We can now compute various moments of the RCE. The first moment is the spike triggered average stimulus (STA),

$$\mathrm{STA}(\tau )=[dI]P[I(t<t_0)|\mathrm{spike}\mathrm{at}t_0]I(t_0\tau ),$$
(4)
which is the object that one computes in reverse correlation \[de Boer and Kuyper, 1968, Rieke et al., 1997\]. If we choose the distribution of input stimuli $`P[I(t<t_0)]`$ to be Gaussian white noise, then for a perceptron–like neuron sensitive to only one direction in stimulus space, it can be shown that the STA or first moment of the response conditional ensemble is proportional to the vector or filter $`f(\tau )`$ that defines this direction \[Rieke et al., 1997\].

Although it is a theorem that the STA is proportional to the relevant filter $`f(\tau )`$, in principle it is possible that the proportionality constant is zero, most plausibly if the neuron’s response has some symmetry, such as phase invariance in the response of high frequency auditory neurons. It also is worth noting that what is really important in this analysis is the Gaussian distribution of the stimuli, not the “whiteness” of the spectrum. For nonwhite but Gaussian inputs the STA measures the relevant filter blurred by the correlation function of the inputs and hence the true filter can be recovered (at least in principle) by deconvolution. For nonGaussian signals and nonlinear neurons there is no corresponding guarantee that the selectivity of the neuron can be separated from correlations in the stimulus \[Sharpee et al., 2002a\].

To obtain more than one relevant direction (or to reveal relevant directions when symmetries cause the STA to vanish), we proceed to second order and compute the covariance matrix of fluctuations around the spike triggered average,

$$C_{\mathrm{spike}}(\tau ,\tau ^{})=[dI]P[I(t<t_0)|\mathrm{spike}\mathrm{at}t_0]I(t_0\tau )I(t_0\tau ^{})\mathrm{STA}(\tau )\mathrm{STA}(\tau ^{}).$$
(5)
In the same way that we compare the spike triggered average to some constant average level of the signal in the whole experiment, we compare the covariance matrix $`C_{\mathrm{spike}}`$ with the covariance of the signal averaged over the whole experiment,

$$C_{\mathrm{prior}}(\tau ,\tau ^{})=[dI]P[I(t<t_0)]I(t_0\tau )I(t_0\tau ^{}),$$
(6)
to construct the change in the covariance matrix

$$\mathrm{\Delta }C=C_{\mathrm{spike}}C_{\mathrm{prior}}.$$
(7)
With time resolution $`\mathrm{\Delta }t`$ in a window of duration $`T`$ as above, all of these covariances are $`D\times D`$ matrices. In the same way that the spike triggered average has the clearest interpretation when we choose inputs from a Gaussian distribution, $`\mathrm{\Delta }C`$ also has the clearest interpretation in this case. Specifically, if inputs are drawn from a Gaussian distribution then it can be shown that \[Bialek and de Ruyter van Steveninck, 2002\]:

1. If the neuron is sensitive to a limited set of $`K`$ input dimensions as in Eq. (2), then $`\mathrm{\Delta }C`$ will have only $`K`$ nonzero eigenvalues.111As with the STA, it is in principle possible that symmetries or accidental features of the function $`g(\stackrel{}{s})`$ would cause some of the $`K`$ eigenvalues to vanish, but this is very unlikely. In this way we can measure directly the dimensionality $`K`$ of the relevant subspace.
2. If the distribution of inputs are both Gaussian and white, then the eigenvectors associated with the nonzero eigenvalues span the same space as that spanned by the filters $`\{f_\mu (\tau )\}`$.
3. For nonwhite (correlated) but still Gaussian inputs, the eigenvectors span the space of the filters $`\{f_\mu (\tau )\}`$ blurred by convolution with the correlation function of the inputs.

Thus the analysis of $`\mathrm{\Delta }C`$ for neurons responding to Gaussian inputs should allow us to identify the subspace of inputs of relevance and to test specifically the hypothesis that this subspace is of low dimension.

Several points are worth noting. First, except in special cases, the eigenvectors of $`\mathrm{\Delta }C`$ and the filters $`\{f_\mu (\tau )\}`$ are not the principal components of the response conditional ensemble, and hence this analysis of $`\mathrm{\Delta }C`$ is not a principal component analysis. Second, the nonzero eigenvalues of $`\mathrm{\Delta }C`$ can be either positive or negative, depending on whether the variance of inputs along that particular direction is larger or smaller in the neighborhood of a spike. Third, although the eigenvectors span the relevant subspace, these eigenvectors do not form a preferred coordinate system within this subspace. Finally, we emphasize that dimensionality reduction—identification of the relevant subspace—is only the first step in our analysis of the computation done by a neuron.

## 3 Measuring the success of dimensionality reduction

The claim that certain stimulus features are most relevant is in effect a model for the neuron, so the next question is how to measure the effectiveness or accuracy of this model. Several different ideas have been suggested in the literature as ways of testing models based on linear receptive fields in the visual system \[Stanley et al., 1999, Keat et al., 2001\] or linear spectrotemporal receptive fields in the auditory system \[Theunissen et al., 2000\]. These methods have in common that they introduce a metric to measure performance—for example, mean square error in predicting the firing rate as averaged over some window of time. Ideally we would like to have a performance measure which avoids any arbitrariness in the choice of metric, and such metric free measures are provided uniquely by information theory \[Shannon, 1948, Cover and Thomas, 1991\].

Observing the arrival time $`t_0`$ of a single spike provides a certain amount of information about the input signals. Since information is mutual, we can also say that knowing the input signal trajectory $`I(t<t_0)`$ provides information about the arrival time of the spike. If “details are irrelevant” then we should be able to discard these details from our description of the stimulus and yet preserve the mutual information between the stimulus and spike arrival times; for an abstract discussion of such selective compression see \[Tishby et al., 1999\]. In constructing our low dimensional model, we represent the complete ($`D`$ dimensional) stimulus $`I(t<t_0)`$ by a smaller number ($`K<D`$) of dimensions $`\stackrel{}{s}=(s_1,s_2,\mathrm{},s_K)`$.

The mutual information $`[I(t<t_0);t_0]`$ is a property of the neuron itself, while the mutual information $`[\stackrel{}{s};t_0]`$ characterizes how much our reduced description of the stimulus can tell us about when spikes will occur. Necessarily our reduction of dimensionality causes a loss of information, so that

$$[\stackrel{}{s};t_0][I(t<t_0);t_0],$$
(8)
but if our reduced description really captures the computation done by the neuron then the two information measures will be very close. In particular if the neuron were described exactly by a lower dimensional model—as for a linear perceptron, or for an integrate and fire neuron \[Agüera y Arcas and Fairhall, 2002\]—then the two information measures would be equal. More generally, the ratio $`[\stackrel{}{s};t_0]/[I(t<t_0);t_0]`$ quantifies the efficiency of the low dimensional model, measuring the fraction of information about spike arrival times that our $`K`$ dimensions capture from the full signal $`I(t<t_0)`$.

As shown by Brenner et al. (2000), the arrival time of a single spike provides an information

$$[I(t<t_0);t_0]I_{\mathrm{one}\mathrm{spike}}=\frac{1}{T}_0^T𝑑t\frac{r(t)}{\overline{r}}\mathrm{log}_2\left[\frac{r(t)}{\overline{r}}\right],$$
(9)
where $`r(t)`$ is the time dependent spike rate, $`\overline{r}`$ is the average spike rate, and $`\mathrm{}`$ denotes an average over time. In principle, information should be calculated as an average over the distribution of stimuli, but the ergodicity of the stimulus justifies replacing this ensemble average with a time average. For a deterministic system like the Hodgkin–Huxley equations, the spike rate is a singular function of time: given the inputs $`I(t)`$, spikes occur at definite times with no randomness or irreproducibility. If we observe these responses with a time resolution $`\mathrm{\Delta }t`$, then for $`\mathrm{\Delta }t`$ sufficiently small the rate $`r(t)`$ at any time $`t`$ either is zero or corresponds to a single spike occurring in one bin of size $`\mathrm{\Delta }t`$, that is $`r=1/\mathrm{\Delta }t`$. Thus the information carried by a single spike is

$$I_{\mathrm{one}\mathrm{spike}}=\mathrm{log}_2\overline{r}\mathrm{\Delta }t.$$
(10)
On the other hand, if the probability of spiking really depends only on the stimulus dimensions $`s_1,s_2,\mathrm{},s_K`$, we can substitute

$$\frac{r(t)}{\overline{r}}\frac{P(\stackrel{}{s}|\mathrm{spike}\mathrm{at}t)}{P(\stackrel{}{s})}.$$
(11)
Replacing the time averages in Eq. (9) with ensemble averages, we find

$$[\stackrel{}{s};t_0]I_{\mathrm{one}\mathrm{spike}}^\stackrel{}{s}=d^KsP(\stackrel{}{s}|\mathrm{spike}\mathrm{at}t)\mathrm{log}_2\left[\frac{P(\stackrel{}{s}|\mathrm{spike}\mathrm{at}t)}{P(\stackrel{}{s})}\right];$$
(12)
for details of these arguments see \[Brenner et al., 2000b\]. This allows us to compare the information captured by the $`K`$ dimensional reduced model with the true information carried by single spikes in the spike train.

For reasons that we will discuss in the following section, and as was pointed out in \[Agüera y Arcas et al., 2001, Agüera y Arcas and Fairhall, 2002\], we will be considering isolated spikes, i.e. spikes separated from previous spikes by a period of silence. This has important consequences for our analysis. Most significantly, as we will be considering spikes that occur on a background of silence, the relevant stimulus ensemble, conditioned on the silence, is no longer Gaussian. Further, we will need to refine our information estimate.

The derivation of Eq. (9) makes clear that a similar formula must determine the information carried by the occurrence time of any event, not just single spikes: we can define an event rate in place of the spike rate, and then calculate the information carried by these events \[Brenner et al., 2000b\]. In the present case we wish to compute the information obtained by observing an isolated spike, or equivalently by the event silence+spike. This is straightforward—we replace the spike rate by the rate of isolated spikes and Eq. (9) will give us the information carried by the arrival time of a single isolated spike. The problem is that this information includes both the information carried by the occurrence of the spike and the information conveyed in the condition that there were no spikes in the preceding $`t_{\mathrm{silence}}`$ msec \[for an early discussion of the information carried by silence see de Ruyter van Steveninck and Bialek (1988)\]. We would like to separate these contributions, since our idea of dimensionality reduction applies only the triggering of a spike, not to the temporally extended condition of nonspiking.

To separate the information carried by the isolated spike itself, we have to ask how much information we gain by seeing an isolated spike given that the condition for isolation has already been met. As discussed by Brenner et al. (2000b), we can compute this information by thinking about the distribution of times at which the isolated spike can occur. Given that we know the input stimulus, the distribution of times at which a single isolated spike will be observed is proportional to $`r_{\mathrm{iso}}(t)`$, the time dependent rate or peri–stimulus time histogram for isolated spikes. With proper normalization we have

$$P_{\mathrm{iso}}(t|\mathrm{inputs})=\frac{1}{T}\frac{1}{\overline{r}_{\mathrm{iso}}}r_{\mathrm{iso}}(t),$$
(13)
where $`T`$ is duration of the (long) window in which we can look for the spike and $`\overline{r}_{\mathrm{iso}}`$ is the average rate of isolated spikes. This distribution has an entropy

$`S_{\mathrm{iso}}(t|\mathrm{inputs})`$ $`=`$ $`{\displaystyle _0^T}𝑑tP_{\mathrm{iso}}(t|\mathrm{inputs})\mathrm{log}_2P_{\mathrm{iso}}(t|\mathrm{inputs})`$ (14)
$`=`$ $`{\displaystyle \frac{1}{T}}{\displaystyle _0^T}𝑑t{\displaystyle \frac{r_{\mathrm{iso}}(t)}{\overline{r}_{\mathrm{iso}}}}\mathrm{log}_2\left[{\displaystyle \frac{1}{T}}{\displaystyle \frac{r_{\mathrm{iso}}(t)}{\overline{r}_{\mathrm{iso}}}}\right]`$ (15)
$`=`$ $`\mathrm{log}_2(T\overline{r}_{\mathrm{iso}}\mathrm{\Delta }t)\mathrm{bits},`$ (16)
where again we use the fact that for a deterministic system the time dependent rate must be either zero or the maximum allowed by our time resolution $`\mathrm{\Delta }t`$. To compute the information carried by a single spike we need to compare this entropy with the total entropy possible when we don’t know the inputs.

It is tempting to think that without knowledge of the inputs, an isolated spike is equally likely to occur anywhere in the window of size $`T`$, which leads us back to Eq. (10) with $`\overline{r}`$ replaced by $`\overline{r}_{\mathrm{iso}}`$. In this case, however, we are assuming that the condition for isolation has already been met. Thus, even without observing the inputs, we know that isolated spikes can occur only in windows of time whose total length is $`T_{\mathrm{silence}}=TP_{\mathrm{silence}}`$, where $`P_{\mathrm{silence}}`$ is the probability that any moment in time is at least $`t_{\mathrm{silence}}`$ after the most recent spike. Thus the total entropy of isolated spike arrival times (given that the condition for silence has been met) is reduced from $`\mathrm{log}_2T`$ to

$$S_{\mathrm{iso}}(t|\mathrm{silence})=\mathrm{log}_2(TP_{\mathrm{silence}}),$$
(17)
and the information that the spike carries beyond what we know from the silence itself is

$`\mathrm{\Delta }I_{\mathrm{iso}\mathrm{spike}}`$ $`=`$ $`S_{\mathrm{iso}}(t|\mathrm{silence})S_{\mathrm{iso}}(t|\mathrm{inputs})`$ (18)
$`=`$ $`{\displaystyle \frac{1}{T}}{\displaystyle _0^T}𝑑t{\displaystyle \frac{r_{\mathrm{iso}}(t)}{\overline{r}_{\mathrm{iso}}}}\mathrm{log}_2\left[{\displaystyle \frac{r_{\mathrm{iso}}(t)}{\overline{r}_{\mathrm{iso}}}}P_{\mathrm{silence}}\right]`$ (19)
$`=`$ $`\mathrm{log}_2\left(\overline{r}_{\mathrm{iso}}\mathrm{\Delta }t\right)+\mathrm{log}_2P_{\mathrm{silence}}\mathrm{bits}.`$ (20)
This information, which is defined independently of any model for the feature selectivity of the neuron, provides the benchmark against which our reduction of dimensionality will be measured. To make the comparison, however, we need the analog of Eq. (12).

Equation (19) provides us with an expression for the information conveyed by isolated spikes in terms of the probability that these spikes occur at particular times; this is analogous to Eq. (9) for single (nonisolated) spikes. If we follow a path analogous to that which leads from Eq. (9) to Eq. (12) we find an expression for the information which an isolated spike provides about the $`K`$ stimulus dimensions $`\stackrel{}{s}`$:

$`\mathrm{\Delta }I_{\mathrm{iso}\mathrm{spike}}^\stackrel{}{s}`$ $`=`$ $`{\displaystyle 𝑑\stackrel{}{s}P(\stackrel{}{s}|\mathrm{iso}\mathrm{spike}\mathrm{at}t)\mathrm{log}_2\left[\frac{P(\stackrel{}{s}|\mathrm{iso}\mathrm{spike}\mathrm{at}t)}{P(\stackrel{}{s}|\mathrm{silence})}\right]}`$ (21)
$`+\mathrm{log}_2P(\mathrm{silence}|\stackrel{}{s}),`$
where the prior is now also conditioned on silence: $`P(\stackrel{}{s}|\mathrm{silence})`$, is the distribution of $`\stackrel{}{s}`$ given that $`\stackrel{}{s}`$ is preceded by a silence of at least $`t_{\mathrm{silence}}`$. Notice that this silence conditioned distribution is not knowable a priori, and in particular it is not Gaussian; $`P(\stackrel{}{s}|\mathrm{silence})`$ must be sampled from data.

The last term in Eq. (21) is the entropy of a binary variable which indicates whether particular moments in time are silent given knowledge of the stimulus. Again, since the HH model is deterministic, this conditional entropy should be zero if we keep a complete description of the stimulus. In fact we’re not interested in describing those features of the stimulus which lead to silence, nor is it fair (as we will see) to judge the success of dimensionality reduction by looking at the prediction of silence which necessarily involves multiple dimensions. To make a meaningful comparison, then, we will assume that there is a perfect description of the stimulus conditions leading to silence, and focus on the stimulus features that trigger the isolated spike. When we approximate these features by the $`K`$ dimensional space $`\stackrel{}{s}`$, we capture an amount of information

$$\mathrm{\Delta }I_{\mathrm{iso}\mathrm{spike}}^\stackrel{}{s}=𝑑\stackrel{}{s}P(\stackrel{}{s}|\mathrm{iso}\mathrm{spike}\mathrm{at}t)\mathrm{log}_2\left[\frac{P(\stackrel{}{s}|\mathrm{iso}\mathrm{spike}\mathrm{at}t)}{P(\stackrel{}{s}|\mathrm{silence})}\right].$$
(22)
This is the information which we can compare with $`\mathrm{\Delta }I_{\mathrm{iso}\mathrm{spike}}`$ in Eq. (20) to determine the efficiency of our dimensionality reduction.

## 4 Characterising the Hodgkin-Huxley neuron

For completeness we begin with a brief review of the dynamics of the space–clamped Hodgkin–Huxley neuron \[Hodgkin and Huxley, 1952\]. Hodgkin and Huxley modelled the dynamics of the current through a patch of membrane with ion–specific conductances:

$$C\frac{dV}{dt}=I(t)\overline{g}_\mathrm{K}n^4\left(VV_\mathrm{K}\right)\overline{g}_{\mathrm{Na}}m^3h\left(VV_{\mathrm{Na}}\right)\overline{g}_l\left(VV_l\right),$$
(23)
where $`I(t)`$ is injected current, $`\mathrm{K}`$ and $`\mathrm{Na}`$ subscripts denote potassium– and sodium–related variables, respectively, and $`l`$ (for “leakage”) terms include all other ion conductances with slower dynamics. $`C`$ is the membrane capacitance. $`V_\mathrm{K}`$ and $`V_{\mathrm{Na}}`$ are ion-specific reversal potentials, and $`V_l`$ is defined such that the total voltage $`V`$ is exactly zero when the membrane is at rest. $`\overline{g}_\mathrm{K}`$, $`\overline{g}_{\mathrm{Na}}`$ and $`\overline{g}_l`$ are empirically determined maximal conductances for the different ion species, and the gating variables $`n`$, $`m`$ and $`h`$ (on the interval $`[0,1]`$) have their own voltage dependent dynamics:

$`dn/dt`$ $`=`$ $`(0.01V+0.1)(1n)\mathrm{exp}(0.1V)0.125n\mathrm{exp}(V/80)`$
$`dm/dt`$ $`=`$ $`(0.1V+2.5)(1m)\mathrm{exp}(0.1V1.5)4m\mathrm{exp}(V/18)`$
$`dh/dt`$ $`=`$ $`0.07(1h)\mathrm{exp}(0.05V)h\mathrm{exp}(0.1V4).`$ (24)
We have used the original values for these parameters, except for changing the signs of the voltages to correspond to the modern sign convention: $`C=1\mu \mathrm{F}/\mathrm{cm}^2`$, $`\overline{g}_\mathrm{K}=36\mathrm{mS}/\mathrm{cm}^2`$, $`\overline{g}_{\mathrm{Na}}=120\mathrm{mS}/\mathrm{cm}^2`$, $`\overline{g}_l=0.3\mathrm{mS}/\mathrm{cm}^2`$, $`V_\mathrm{K}=12\mathrm{mV}`$, $`V_{\mathrm{Na}}=+115\mathrm{mV}`$, $`V_l=+10.613\mathrm{mV}`$. We have taken our system to be a $`\pi \times 30^2\mu \mathrm{m}^2`$ patch of membrane. We solve these equations numerically using fourth order Runge–Kutta integration.

The system is driven with a Gaussian random noise current $`I(t)`$, generated by smoothing a Gaussian random number stream with an exponential filter to generate a correlation time $`\tau `$. It is convenient to choose $`\tau `$ to be longer than the time steps of numerical integration, since this guarantees that all functions are smooth on the scale of single time steps. Here we will always use $`\tau =0.2`$ msec, a value that is both less than the timescale over which we discretize the stimulus for analysis, and far less than the neuron’s capacitative smoothing timescale $`RC3`$ msec. $`I(t)`$ has a standard deviation $`\sigma `$, but since the correlation time is short the relevant parameter usually is the spectral density $`S=\sigma ^2\tau `$; we also add a DC offset $`I_0`$. In the following, we will consider two parameter regimes, $`I_0=0`$, and $`I_0`$ a finite value, which leads to more periodic firing.

The integration step size is fixed at $`0.05\mathrm{msec}`$. The key numerical experiments were repeated at a step size of $`0.01\mathrm{msec}`$ with identical results. The time of a spike is defined as the moment of maximum voltage, for voltages exceeding a threshold (see Fig. 1), estimated to subsample precision by quadratic interpolation. As spikes are both very stereotyped and very large compared to subspiking fluctuations, the precise value of this threshold is unimportant; we have used $`+20`$ mV.

### 4.1 Qualitative description of spiking

The first step in our analysis is to use reverse correlation, Eq. (4), to determine the average stimulus feature preceding a spike, the spike triggered average (STA). In Fig. 1(a) we display the STA in a regime where the spectral density of the input current is $`6.5\times 10^4\mathrm{nA}^2\mathrm{msec}`$. The spike triggered averages of the gating terms $`n^4`$ (proportion of open potassium channels) and $`m^3h`$ (proportion of open sodium channels) and the membrane voltage $`V`$ are plotted in parts (b) and (c). The errorbars mark the standard deviation of the trajectories of these variables.

As expected, the voltage and gating variables follow highly stereotyped trajectories during the $``$5 msec surrounding a spike: First, the rapid opening of the sodium channels causes a sharp membrane depolarization (or rise in $`V`$); the slower potassium channels then open and repolarize the membrane, leaving it at a slightly lower potential than rest. The potassium channels close gradually, but meanwhile the membrane remains hyperpolarized and, due to its increased permeability to potassium ions, at lower resistance. These effects make it difficult to induce a second spike during this $``$15 msec “refractory period.” Away from spikes, the resting levels and fluctuations of the voltage and gating variables are quite small. The larger values evident in Fig. 1(b) and (c) by $`\pm 15`$ msec are due to the summed contributions of nearby spikes.

The spike triggered average current has a largely transient form, so that spikes are on average preceded by an upward swing in current. On the other hand, there is no obvious “bottleneck” in the current trajectories, so that the current variance is almost constant throughout the spike. This is qualitatively consistent with the idea of dimensionality reduction: If the neuron ignores most of the dimensions along which the current can vary, then the variance—which is shared almost equally among all dimensions for this near white noise—can change only by a small amount.

### 4.2 Interspike interaction

Although the STA has the form of a differentiating kernel, suggesting that the neuron detects edge–like events in the current vs. time, there must be a DC component to the cell’s response. We recall that for constant inputs the Hodgkin–Huxley model undergoes a bifurcation to constant frequency spiking, where the frequency is a function of the value of the input above onset. Correspondingly the STA does not sum precisely to zero; one might think of it as having a small integrating component that allows the system to spike under DC stimulation, albeit only above a threshold.

The system’s tendency to periodic spiking under DC current input also is felt under dynamic stimulus conditions, and can be thought of as a strong interaction between successive spikes. We illustrate this by considering a different parameter regime with a small DC current and some added noise ($`I_0=0.11\mathrm{nA}`$ and $`S=0.8\times 10^4\mathrm{nA}^2\mathrm{sec}`$). Note that the DC component puts the neuron in the metastable region of its $`fI`$ curve, Fig. 2. In this regime the neuron tends to fire quasi-regular trains of spikes intermittently, as shown in Fig. 3. We will refer to these quasi–regular spike sequences as “bursts” (note that this term is often used to refer to compound spikes in neurons with additional channels; such events do not occur in the Hodgkin–Huxley model).

Spikes can be classified into three types—those initiating a spike burst; those within a burst; and those ending a burst. The minimum length of the silence between bursts is taken in this case to be $`70\mathrm{msec}`$. Taking these three categories of spike as different “symbols” \[de Ruyter van Steveninck and Bialek, 1988\], we can determine the average stimulus for each. These are shown in Fig. 4 with the spike at $`t=0`$.

In this regime, the initial spike of a burst is preceded by a rapid oscillation in the current. Spikes within a burst are affected much less by the current; the feature immediately preceding such spikes is similar in shape to a single “wavelength” of the leading spike feature, but is of much smaller amplitude, and is temporally compressed into the interspike interval. Hence, although it is clear that the timing of a spike within a burst is determined largely by the timing of the previous spike, the current plays some role in affecting the precise placement. This also demonstrates that the shape of the STA is not the same for all spikes; it depends strongly and nontrivially on the time to the previous spike, and this is related to the observation that subtly different patterns of two or three spikes correspond to very different average stimuli \[de Ruyter van Steveninck and Bialek, 1988\]. For a reader of the spike code, a spike within a burst conveys a different message about the input than the spike at the onset of the burst. Finally, the feature ending a burst has a very similar form to the onset feature, but reversed in time. Thus, to a good approximation, the absence of a spike at the end of a burst can be read as the opposite of the onset of the burst.

In summary, this regime of the HH neuron is similar to a “flip-flop”, or 1-bit memory. Like its electronic analogue, the neuron’s memory is preserved by a feedback loop, here implemented by the interspike interaction. Large fluctuations in the input current at a certain frequency “flip” or “flop” the neuron between its silent and spiking states. However, while the neuron is spiking, further details of the input signal are transmitted by precise spike timing within a burst. If we calculate the spike triggered average of all spikes for this regime, without regard to their position within a burst, then as shown in Fig. 5 the relatively well localized leading spike oscillation of Fig. 4 is replaced by a long-lived oscillating function resulting from the spike periodicity. This is shown explicitly by comparing the overall STA with the spike autocorrelation, also shown in Fig. 5. This same effect is seen in the STA of the burst spikes, which in fact dominates the overall average. Prediction of spike timing using such an STA would be computationally difficult, due to its extension in time, but, more seriously, unsuccessful, as most of the function is an artifact of the spike history rather than the effect of the stimulus.

While the effects of spike interaction are interesting, and should be included in a complete model for spike generation, we wish here to consider only the current’s role in initiating spikes. Therefore, as we have argued elsewhere, we limit ourselves initially to the cases in which interspike interaction plays no role \[Agüera y Arcas et al., 2001, Agüera y Arcas and Fairhall, 2002\]. These “isolated” spikes can be defined as spikes preceded by a silent period $`t_{\mathrm{silence}}`$ long enough to ensure decoupling from the timing of the previous spike. A reasonable choice for $`t_{\mathrm{silence}}`$ can be inferred directly from the interspike interval distribution $`P(\mathrm{\Delta }t)`$, illustrated in Fig. 6. For the HH model, as in simpler models and many real neurons \[Brenner et al., 1998\], the form of $`P(\mathrm{\Delta }t)`$ has three noteworthy features: a refractory “hole” during which another spike is unlikely to occur, a strong mode at the preferred firing frequency, and an exponentially decaying, or Poisson tail. The details of all three of these features are functions of the parameters of the stimulus, and certain regimes may be dominated by only one or two features. The emergence of Poisson statistics in the tail of the distribution implies that these events are independent, so we can infer that the system has lost memory of the previous spike. We will therefore take isolated spikes to be those preceded by a silent interval $`\mathrm{\Delta }tt_{\mathrm{silence}}`$, where $`t_{\mathrm{silence}}`$ is well into the Poisson regime. Note that the burst onset spikes of Fig. 4 are isolated spikes by this definition.

## 5 Isolated spike analysis

Focusing now on isolated spikes, we proceed to a second–order analysis of the current fluctuations around the isolated spike triggered average, Fig. 7. We consider the response of the HH neuron to currents $`I(t)`$ with mean $`I_0=0`$ and spectral density of $`S=6.5\times 10^4\mathrm{nA}^2\mathrm{sec}`$. Isolated spikes in this regime are defined by $`t_{\mathrm{silence}}=60`$ msec.

### 5.1 How many dimensions?

As explained in Section 2, our path to dimensionality reduction begins with the computation of covariance matrices for stimulus fluctuations surrounding a spike. The matrices are accumulated from stimulus segments 200 samples in length, roughly corresponding to sampling at the correlation timescale for sufficiently long to capture the relevant features. Thus we begin in a 200 dimensional space. We emphasize that the theorem which connects eigenvalues of the matrix $`\mathrm{\Delta }C`$ to the number of relevant dimensions is valid only for truly Gaussian distributions of inputs, and that by focusing on isolated spikes we are essentially creating a nonGaussian stimulus ensemble, namely those stimuli which generate the silence out of which the isolated spike can appear. Thus we expect that the covariance matrix approach will give us a heuristic guide to our search for lower dimensional descriptions, but we should proceed with caution.

The “raw” isolated spike triggered covariance $`C_{\mathrm{iso}\mathrm{spike}}`$ and the corresponding covariance difference $`\mathrm{\Delta }C`$, Eq. (7), are shown in Fig. 8. The matrix shows the effect of the silence as an approximately translationally invariant band preceding the spike, the second order analogue of the constant negative bias in the isolated spike STA, Fig.7. The spike itself is associated with features localized to $`\pm 15`$msec. In Fig. 9 we show the spectrum of eigenvalues of $`\mathrm{\Delta }C`$ computed using a sample of 80,000 spikes. Before calculating the spectrum, we multiply $`\mathrm{\Delta }C`$ by $`C_{\mathrm{prior}}^1`$. This has the effect of giving us eigenvalues scaled in units of the input standard deviation along each dimension. Because the correlation time is short, $`C_{\mathrm{prior}}`$ is nearly diagonal.

While the eigenvalues decay rapidly, there is no obvious set of outstanding eigenvalues. To verify that this is not an effect of finite sampling, Fig. 10 shows the spectrum of eigenvalue magnitudes as a function of sample size $`N`$. Eigenvalues which are truly zero up to the noise floor determined by sampling decrease like $`\sqrt{N}`$. We find that a sequence of eigenvalues emerges stably from the noise.

These results do not, however, imply that a low dimensional approximation cannot be identified. The extended structure in the covariance matrix induced by the silence requirement is responsible for the apparent high dimensionality. In fact, as has been shown in \[Agüera y Arcas and Fairhall, 2002\], the covariance eigensystem includes modes that are local and spike–associated, and others that are extended and silence-associated, and thus irrelevant to a causal model of spike timing prediction. Fortunately, because extended silences and spikes are (by definition) statistically independent, there is no mixing between the two types of modes. To identify the spike–associated modes, we follow the diagnostic of \[Agüera y Arcas and Fairhall, 2002\], computing the fraction of the energy of each mode concentrated in the period of silence, which we take to be $`60t40`$ msec. The energy of a spike–associated mode in the silent period is due entirely to noise, and will therefore decrease like $`1/n_{\mathrm{spikes}}`$ with increasing sample size, while this energy remains of order unity for silence modes. Carrying out the test on the covariance modes, we obtain Fig. 11, which shows that the first and fourth modes rapidly emerge as spike–associated. Two further spike–associated modes appear over the sample shown, with the suggestion of other, weaker modes yet to emerge. The two leading silence modes are shown in Fig. 13. Those shown are typical; most modes resemble Fourier modes, as the silence condition is close to time translationally invariant.

Examining the eigenvectors corresponding to the two leading spike–associated eigenvalues, which for convenience we will denote $`s_1`$ and $`s_2`$ (although they are not the leading modes of the matrix), we find (Fig. 12) that the first mode closely resembles the isolated spike STA, and the second is close to the derivative of the first. Both modes approximate differentiating operators; there is no linear combination of these modes which would produce an integrator.

If the neuron filtered its input and generated a spike when the output of the filter crosses threshold, we would find two significant dimensions associated with a spike. The first dimension would correspond simply to the filter, as the variance in this dimension is reduced to zero (for a noiseless system) at the occurrence of a spike. As the threshold is always crossed from below, the stimulus projection onto the filter’s derivative must be positive, again resulting in a reduced variance. It is tempting to suggest, then, that filtered threshold crossing is a good approximation to the HH model, but we will see that this is not correct.

### 5.2 Evaluating the nonlinearity

At each instant of time we can find the projections of the stimulus along the leading spike–associated dimensions $`s_1`$ and $`s_2`$. By construction, the distribution of these signals over the whole experiment, $`P(s_1,s_2)`$, is Gaussian. The appropriate prior for the isolation condition, $`P(s_1,s_2|\mathrm{silence})`$, differs only subtly from the Gaussian prior. On the other hand, for each spike we obtain a sample from the distribution $`P(s_1,s_2|\mathrm{iso}\mathrm{spike}\mathrm{at}t_0)`$, leading to the picture in Fig. 14. The prior and spike conditional distributions are clearly better separated in two dimensions than in one, which means that the two dimensional description captures more information than projection onto the spike triggered average alone. Surprisingly, the spike conditional distribution is curved, unlike what we would expect for a simple thresholding device. Furthermore, the eigenvalue of $`\mathrm{\Delta }C`$ which we associate with the direction of threshold crossing (plotted on the $`y`$-axis in Fig. 14) is positive, indicating increased rather than decreased variance in this direction. As we see, projections onto this mode are almost equally likely to be positive or negative, ruling out the threshold crossing interpretation.

Combining eqs. (2) and (3) for isolated spikes, we have

$$g(s_1,s_2)=\frac{P(s_1,s_2|\mathrm{iso}\mathrm{spike}\mathrm{at}t_0)}{P(s_1,s_2|\mathrm{silence})},$$
(25)
so that these two distributions determine the input/output relation of the neuron in this 2D space \[Brenner et al., 2000a\]. Recall that although the subspace is linear, $`g`$ can have arbitrary nonlinearity. Fig. 14 shows that this input/output relation has clear structure, but also some fuzziness. As the HH model is deterministic, the input/output relation should be a singular function in the continuous space of inputs—spikes occur only when certain exact conditions are met. Of course, finite time resolution introduces some blurring, and so we need to understand whether the blurring of the input/output relation in Fig. 14 is an effect of finite time resolution or a real limitation of the 2D description.

### 5.3 Information captured in two dimensions

We will measure the effectiveness of our description by computing the information in the 2D approximation, according to the methods described in Sect. 3. If the two dimensional approximation were exact we would find that $`I_{\mathrm{iso}\mathrm{spike}}^{s_1,s_2}=I_{\mathrm{iso}\mathrm{spike}}`$; more generally one finds $`I_{\mathrm{iso}\mathrm{spike}}^{s_1,s_2}I_{\mathrm{iso}\mathrm{spike}}`$, and the fraction of the information captured measures the quality of the approximation. This fraction is plotted in Fig. 15 as a function of time resolution. For comparison, we also show the information captured in the one dimensional case, considering only the stimulus projection along the STA.

We find that our low dimensional model captures a substantial fraction of the total information available in spike timing in a HH neuron over a range of time resolutions. The approximation is best near $`\mathrm{\Delta }t=3\mathrm{msec}`$, reaching 75%. Thus the complex nonlinear dynamics of the Hodgkin–Huxley model can be approximated by saying that the neuron is sensitive to a two dimensional linear subspace in the high dimensional space of input signals, and this approximate description captures up to 75% of the mutual information between input currents and (isolated) spike arrival times.

The dependence of information on time resolution, Fig. 15, shows that the absolute information captured saturates for both the 1D and 2D cases, at $`3.2`$ and $`4.1`$ bits respectively. Hence, for smaller $`\mathrm{\Delta }t`$, the information fraction captured drops. The model provides, at its best, a time resolution of 3 msec, so that information carried by more precise spike timing is lost in our low dimensional projection. Might this missing information be important for a real neuron? Stochastic HH simulations with realistic channel densities suggest that the timing of spikes in response to white noise stimuli is reproducible to within 1–2 msec \[Schneidman et al., 1998\], a figure which is comparable to what is observed for pyramidal cells in vitro \[Mainen and Sejnowski, 1995\], as well in vivo in the fly’s visual system \[de Ruyter van Steveninck et al., 1997, Lewen et al., 2001\], the vertebrate retina \[Berry II et al., 1997\], the cat LGN \[Reinagel and Reid, 2000\] and the bat auditory cortex \[Dear et al., 1993\]. This suggests that such timing details may indeed be important. We must therefore ask why our approximation seems to carry an inherent time resolution limitation, and why, even at its optimal resolution, the full information in the spike is not recovered.

For many purposes, recovering 75% of the information at $`3`$ msec resolution might be considered a resounding success. On the other hand, with such a simple underlying model we would hope for a more compelling conclusion. From a methodological point of view it behooves us to ask what we are missing in our 2D model, and perhaps the methods we use in finding the missing information in the present case will prove applicable more generally.

## 6 What is missing?

The obvious first approach to improving the 2D approximation is to add more dimensions. Let us consider the neglected modes. We recall from Fig. 11 that, in simulations with very large numbers of spikes, we can isolate at least two more modes that have significant eigenvalues and are associated with the isolated spike rather than the preceding silence; these are shown in Fig. 16. We see that these modes look like higher order derivatives, which makes some sense since we are missing information at high time resolution. On the other hand, if all we are doing is expanding in a basis of increasingly higher order derivatives it is not clear that we will do qualitatively better by including one or two more terms.

Our original model attempted to approximate the conditional response ensemble as lying within a $`K`$–dimensional linear subspace of the original $`D`$–dimensional input space. The obvious generalization of this idea is to regard the subspace as a curved, but still low dimensional manifold. Several methods have been proposed for the general problem of identifying low dimensional nonlinear manifolds \[Cottrell et al., 1988, Boser et al., 1992, Guyon et al., 1993, Oja and Karhunen, 1995, Roweis and Saul, 2000\], but these various approaches share the disadvantage that the manifold or, equivalently, the relevant set of features remains implicit. Our hope is to understand the behaviour of the neuron explicitly; we therefore wish to obtain an explicit representation of the spiking surface.

A first approach to representing the curved manifold is to try to form a locally linear tiling. Beginning with first order statistics, there is only one natural direction along which to parametrize—the spike triggered average. We sort the stimulus histories according to their projection onto the spike triggered average, divide them into bins, and recompute the averages over the bins individually. If these conditional averages have a component orthogonal to the overall STA, then this component is a locally meaningful basis vector for the manifold at that “slice.”. This procedure results in the family of curves orthogonal to the STA shown in Fig. 17.

The analysis shown here was carried out with 80,000 isolated spikes; note that a similar number of spikes cannot resolve more than two spike–associated covariance modes in the covariance matrix analysis. The spikes were divided into eight bins according to their projection onto the STA, with 10,000 spikes per bin. Applying singular value decomposition to the family of curves shows that there are at least four significant independent directions in stimulus space apart from the STA. This gives us a lower bound on the embedding dimension of the manifold.

Computing the information as a function of $`\mathrm{\Delta }t`$ using this locally linear model, we obtain the curve shown in Fig.18, where the results can be compared against the information found from the STA alone and from the covariance modes. The information from the new model captures a maximum of 4.8 bits, recovering $`90\%`$ of the information at a time resolution of approximately 1 msec.

## 7 Discussion

The Hodgkin–Huxley equations describe the dynamics of four degrees of freedom, and almost since these equations first were written down there have been attempts to find simplifications or reductions. FitzHugh and Nagumo proposed a two dimensional system of equations that approximate the HH model \[Fitzhugh, 1961, Nagumo et al., 1962\], and this has the advantage that one can visualize the trajectories directly in the plane and thus achieve an intuitive graphical understanding of the dynamics and its dependence on parameters. The need for reduction in the sense pioneered by FitzHugh and by Nagumo et al. has become only more urgent with the growing use of increasingly complex HH–style model neurons with many different channel types. With this problem in mind, Kepler, Abbott and Marder have introduced reduction methods which are more systematic, making use of the difference in timescales among the gating variables \[Kepler et al., 1992, Abbott and Kepler, 1990\].

In the presence of constant current inputs, it makes sense to describe the Hodgkin–Huxley equations as a four dimensional autonomous dynamical system; by well known methods in dynamical systems theory, one could consider periodic input currents by adding an extra dimension. The question asked by FitzHugh and Nagumo was whether this four or five dimensional description could be reduced to two or three dimensions.

Closer in spirit to our approach is the work by Kistler et al., who focused in particular on the interaction among successive action potentials \[Kistler et al., 1997\]. They argued that one could approximate the Hodgkin–Huxley model by a nearly linear dynamical system with a threshold, identifying threshold crossing with spiking, provided that each spike generated either a change in threshold or an effective input current which influences the generation of subsequent spikes.

The notion of model dimensionality considered here is distinct from the dynamical systems perspective in which one simply counts the system’s degrees of freedom. Here we are attempting to find a description of the dynamics which is essentially functional or computational. We have identified the output of the system as spike times, and our aim is to construct as complete a description as possible of the mapping between input and output. The dimensionality of our model is that of the space of inputs relevant for this mapping. There is no necessary relationship between these two notions of dimensionality. For example, in a neural network with two attractors, a system described by a potentially large number of variables, there might be a simple rule (perhaps even a linear filter) which allows us to look at the inputs to the network and determine the times at which the switching events will occur. Conversely, once we leave the simplified world of constant or periodic inputs, even the small number of differential equations describing a neuron’s channel dynamics could in principle be equivalent to a very complicated set of rules for mapping inputs into spike times.

In our context, simplicity is (roughly) feature selectivity: the mapping is simple if spiking is determined by a small number of features in the complex history of inputs. Following the ideas which emerged in the analysis of motion sensitive neurons in the fly \[de Ruyter van Steveninck and Bialek, 1988, Bialek and de Ruyter van Steveninck, 2002\], we have identified “features” with “dimensions” and searched for low dimensional descriptions of the input history which preserve the mutual information between inputs and outputs (spike times). We have considered only the generation of isolated spikes, leaving aside the question of how spikes interact with one another as considered by Kistler et al.. For these isolated spikes we began by searching for projections onto a low dimensional linear subspace of the originally $`200`$ dimensional stimulus space, and we found that a substantial fraction of the mutual information could be preserved in a model with just two dimensions. Searching for the information that is missing from this model, we found that rather than adding more (Euclidean) dimensions we could capture $`90\%`$ of the information at high time resolution by keeping a two dimensional description but allowing these dimensions to curve, so that the neuron is sensitive to stimulus features which lie on a two dimensional manifold.

The geometrical picture of neurons as being sensitive to features that are defined by a low dimensional manifold in stimulus space is attractive and, as noted in the Introduction, corresponds to a widely shared intuition about the nature of neuronal feature selectivity. While curved manifolds often appear as the targets for learning in complex neural computations such as invariant object recognition, the idea that such manifolds appear already in the description of single neuron computation we believe to be novel.

While we have exploited the fact that long simulations of the Hodgkin–Huxley model are quite tractable to generate large amounts of “data” for our analysis, it is important that, in the end, our construction of a curved manifold as the relevant stimulus subspace involves a series of computations which are just simple generalizations of the conventional reverse correlation or spike triggered average. This suggests that our approach can be applied to real neurons without requiring qualitatively larger data sets than might have been needed for a careful reverse correlation analysis. In the same spirit, recent work has shown how covariance matrix analysis of the fly’s motion sensitive neurons can reveal nonlinear computations in a four dimensional subspace using data sets of fewer than $`10^4`$ spikes \[Bialek and de Ruyter van Steveninck, 2002\]. Low dimensional (linear) subspaces can be found even in the response of model neurons to naturalistic inputs if one searches directly for dimensions which capture the largest fraction of the mutual information between inputs and spikes \[Sharpee et al., 2002a\], and again the errors involved in identifying the relevant dimensions are comparable to the errors in reverse correlation \[Sharpee et al., 2002b\]. All of these results point to the practical feasibility of describing real neurons in terms of nonlinear computation on low dimensional relevant subspaces in a high dimensional stimulus space.

Our reduced model of the Hodgkin–Huxley neuron both illustrates a novel approach to dimensional reduction and gives new insight into the computation performed by the neuron. The reduced model is essentially that of an edge detector for current trajectories, but is sensitive to a further stimulus parameter, producing a curved manifold. An interpretation of this curvature will be presented in a forthcoming manuscript. This curved representation is able to capture almost all information that isolated spikes convey about the stimulus, or conversely, allow us to predict isolated spike times with high temporal precision from the stimulus. The emergence of a low dimensional curved manifold in a model as simple as the Hodgkin–Huxley neuron suggests that such a description may be also appropriate for biological neurons.

Our approach is limited in that we address only isolated spikes. This restricted class of spikes nonetheless has biological relevance; for example, in vertebrate retinal ganglion cells \[Berry II and Meister, 1999\] and in rat somatosensory cortex \[Panzeri et al., 2001\], the first spike of a burst has been shown to convey distinct (and the majority of the) information. However, a clear next step in this program is to extend our formalism to take into account interspike interaction. For neurons or models with explicit long timescales, adaptation induces very long range history dependence which complicates the issue of spike interactions considerably. A full understanding of the interaction between stimulus and spike history will therefore, in general, involve understanding the meanings of spike patterns \[de Ruyter van Steveninck and Bialek, 1988, Brenner et al., 2000b\] and the influence of the larger statistical context \[Fairhall et al., 2001\]. Our results point to the need for a more parsimonious description of self–excitation, even for the simple case of dependence only on the last spike time.

We would like to close by reminding the reader of the more ambitious goal of building bridges between the burgeoning molecular level description of neurons and the functional or computational level. Armed with a description of spike generation as a nonlinear operation on a low dimensional, curved manifold in the space of inputs, it is natural to ask how the details of this computational picture are related to molecular mechanisms. Are neurons with more different types of ion channels sensitive to more stimulus dimensions, or do they implement more complex nonlinearities in a low dimensional space? Are adaptation and modulation mechanisms that change the nonlinearity separable from those which change the dimensions to which the cell is sensitive? Finally, while we have shown how a low dimensional description can be constructed numerically from observations of the input/output properties of the neuron, one would like to understand analytically why such a description emerges and whether it emerges universally from the combinations of channel dynamics selected by real neurons.

## Acknowledgments

We thank N. Brenner for discussions at the start of this work, and M. Berry for comments on the manuscript.


format: text

text:
# Noncommutative fields in curved space

## I Introduction

Usual noncommutative theories Seiberg  are based on a constant antisymmetric quantity of rank two. Consequently, they do not exihibt the Lorentz symmetry for spacetime dimensions higher than two. In this way, it does not make much sense to look for noncommutative theories invariant by general coordinate transformations or even noncommutative quantum fields in curved background.

The question of Lorentz invariance in noncommutative theories has been recently considered Carlson , where the antissymetric tensor is taken as an independent antisymmetric tensor operator. Our purpose in the present paper is to consider the main aspects of noncommutative fields in curved space. Our approach is such that the noncommutativity does not affect the metric but it does on the geommetry. For simplicity, we deal with the particular case of $`D=2`$ where there is a constant flat limit for the noncommutative antisymmetric tensor. The only difference for dimensions higher than two is that this limit should be spacetime dependent. We deal with real scalar fields, but what will be presented here can be naturally extended for any kind of fields.

Our paper is organized as follows. In Sec. II we present the main aspects of the problem where we see that it is necessity to adopt a more general definition for the Moyal product. In Sec. III we discuss in details the property of commutativity for quadractic terms in the action which is not kept in curved space, and in Sec. IV we consider the question of diffeomorphism symmetry. We left Sec. V for some concluding remarks.

## II Presenting the problem

Let us first consider the usual action for real scalar fields in curved spacetime at $`D=2`$,

$$S=\frac{1}{2}d^2x\sqrt{g}g^{\mu \nu }_\mu \varphi _\nu \varphi .$$
(2.1)
Of course, this theory is invariant under general coordinate transformation. The origin of this symmetry is related fo the (first-class) constraints Dirac

$$\pi ^{\mu \nu }=\frac{}{\dot{g}_{\mu \nu }}=0,$$
(2.2)
where $``$ is the Lagrangian density corresponding to the action (2.1). The components of $`\pi ^{\mu \nu }`$ given by Eq. (2.2) are three (primary) constraints. Constructing the Dirac Hamiltonian and imposing that these constraints do not evolve in time, one obtain two secondary constraints. Repeating this procedure, no more constraints are obtained and all of them are first-class Dirac . According to Castellani Castellani , primary first-class constraints are related to symmetries of the theory. In the present case, two of them concern to the general coordinate transformation (embodied in the diffeomorphism algebra generated by the two secondary constraints) and the last primary constrains is related to the conformal (Weyl) symmetry.

If one direct follows the general rule of transforming usual theories in noncommutative ones by replacing product of fields by the Moyal product Seiberg , we would have

$$S=\frac{1}{2}d^2x\sqrt{g}g^{\mu \nu }_\mu \varphi _\nu \varphi ,$$
(2.3)
where $``$ means the Moyal product, which for any two fields $`\varphi _1`$ and $`\varphi _2`$, it is defined as

$$\varphi _1\varphi _2=\mathrm{exp}\left(\frac{ik}{2}ϵ^{\mu \nu }_\mu ^x_\nu ^y\right)\varphi _1(x)\varphi _2(y)|_{x=y},$$
(2.4)
and $`k`$ is a constant with mass dimension minus two and $`ϵ^{\mu \nu }`$ is the usual Levi-Civita tensor (in this paper, we adopt the convention $`ϵ^{01}=1`$).

The action (2.3) does not exhibit the diffeomorphism symmetry because the canonical momentum conjugate to $`g_{\mu \nu }`$ is not a constraint anymore. Further, notice that the Moyal product must also acts on the elements of $`g=detg_{\mu \nu }`$, making the relation (2.3) meaningless for any consistent use.

The natural alternative is to consider the derivative operators of the Moyal product in curved space as covariant ones. In this way, they do not affect the spacetime geommetry and the inconsistencies mentioned above disappear. So, for the same two fields $`\varphi _1`$ and $`\varphi _2`$, we make

$$\varphi _1\varphi _2=\mathrm{exp}\left(\frac{ik}{2\sqrt{g}}ϵ^{\mu \nu }D_\mu ^xD_\nu ^y\right)\varphi _1(x)\varphi _2(y)|_{x=y},$$
(2.5)
where the symbol $``$ is to distinguish from the Moyal product in the flat case. It is important to emphasize that there is no problem related to the commutativity of the covariant derivatives because they act in different spacetime coordinates. Notice that in curved space, even in two dimensions, the characterist antisymmetric tensor $`\theta _{\mu \nu }`$ is spacetime dependent, $`\theta _{\mu \nu }(x)=\sqrt{g}ϵ_{\mu \nu }`$ (for the contravariant notation, $`\theta ^{\mu \nu }(x)=ϵ^{\mu \nu }/\sqrt{g}`$).

Considering the Moyal product given by (2.5), it does not affect products involving the metric tensor. Consequently, instead of the action (2.3) we simply have

$$S=\frac{1}{2}d^2x\sqrt{g}g^{\mu \nu }_\mu \varphi _\nu \varphi .$$
(2.6)
However, the price paid is that some usual properties of the Moyal product are lost. For example, it affects quadratic terms of the action (we are going to analyse this point with details in the next section). It is easy to see that the same occurs for cyclic properties.

## III A fundamental characteristic of the new Moyal product

The expansion of the exponential operator in (2.5) leads to

$`\varphi _1\varphi _2=\varphi _1\varphi _2+{\displaystyle \frac{ik}{2\sqrt{g}}}ϵ^{\mu \nu }_\mu \varphi _1_\nu \varphi _2`$
$`+{\displaystyle \frac{1}{2!}}\left({\displaystyle \frac{ik}{2\sqrt{g}}}\right)^2ϵ^{\mu \nu }ϵ^{\rho \lambda }D_\rho _\mu \varphi _1D_\lambda _\nu \varphi _2`$
$`+\mathrm{}`$ (3.1)
Under integration over $`d^2x\sqrt{g}`$ we directly see that the second term is zero, but for the third one we have (we drop the initial factor $`(ik/2)^2/2!`$)

$`{\displaystyle d^2x\sqrt{g}\frac{ϵ^{\mu \nu }}{\sqrt{g}}\frac{ϵ^{\rho \lambda }}{\sqrt{g}}D_\rho _\mu \varphi _1D_\lambda _\nu \varphi _2}`$
$`={\displaystyle d^2x\sqrt{g}\left(g^{\mu \lambda }g^{\nu \rho }g^{\mu \rho }g^{\nu \lambda }\right)D_\rho _\mu \varphi _1D_\lambda _\nu \varphi _2},`$
$`={\displaystyle d^2x\sqrt{g}[D_\mu ,D_\nu ]^\mu \varphi _1^\nu \varphi _2},`$
$`={\displaystyle d^2x\sqrt{g}g^{\mu \rho }R_{\lambda \rho \mu \nu }^\lambda \varphi _1^\nu \varphi _2},`$
$`={\displaystyle d^2x\sqrt{g}Rg^{\mu \nu }_\mu \varphi _1_\nu \varphi _2},`$ (3.2)
where we have used some identitities listed in the Appendix A. We observe that by virtue of the curvature this term is not zero. The same occurs for the next terms of the expansion (3.1). Let us write some of them below (also dropping the factors involving $`k`$)

$$d^2xϵ^{\mu \nu }\frac{ϵ^{\rho \lambda }}{\sqrt{g}}\frac{ϵ^{\eta \xi }}{\sqrt{g}}D_\eta D_\rho _\mu \varphi _1D_\xi D_\lambda _\nu \varphi _2=\frac{1}{2}d^2xϵ^{\mu \nu }\left(_\mu RD_\nu _\lambda \varphi _1^\lambda \varphi _2+\frac{1}{2}R^2_\mu \varphi _1_\nu \varphi _2\right),$$
(3.3)
$`{\displaystyle d^2xϵ^{\mu \nu }\frac{ϵ^{\rho \lambda }}{\sqrt{g}}\frac{ϵ^{\eta \xi }}{\sqrt{g}}\frac{ϵ^{\theta \gamma }}{\sqrt{g}}D_\theta D_\eta D_\rho _\mu \varphi _1D_\gamma D_\xi D_\lambda _\nu \varphi _2}`$
$`={\displaystyle \frac{1}{2}}{\displaystyle }d^2x\sqrt{g}R(3D_\mu \mathrm{}\varphi _1D^\mu \mathrm{}\varphi _23D_\alpha D_\mu _\nu \varphi _1D^\alpha D^\mu ^\nu \varphi _2`$
$`RD_\mu \mathrm{}\varphi _1^\mu \varphi _2R_\mu \varphi _1D^\mu \mathrm{}\varphi _2+R^2_\mu \varphi _1^\mu \varphi _2),`$ (3.4)
where $`\mathrm{}`$ is the Laplace-Beltrami operator $`\mathrm{}=g^{\mu \nu }D_\mu D_\nu `$.

## IV Diffeomorphism symmetry

An interesting question concerns to symmetries. It is well-known that noncommutative theories exhibit the unitarity problem Gomis . This is so due to the presence of higher derivatives. But the problem here is more subtle because de number of higher derivatives increases indefinitely. This means that the unitarity problem is accompanied by a more fundamental question related to the dynamics of these theories. This is so because dynamics is given by the last order of its time derivative (in lower orders take place momenta and constraints). If there is an infinite number of time derivatives, there is an ambiguity in the limit of these two regions.

A way of understanding this problem in a deeper way consists in analyzing separately the terms of the expansion (isolately, there is no problem related to the ambiguity of dynamics and constraints) Barc1 . Since there is dynamics for the terms of (3.1), that is simpler than the corresponding ones of (2.6), let us consider the second term of the expansion of expression (3.1). Denoting it by $`S^{(2)}`$, and taking $`\varphi _1=\varphi _2=\varphi `$, we have

$$S^{(2)}=\frac{1}{2}d^2xϵ^{\mu \nu }\frac{ϵ^{\rho \lambda }}{\sqrt{g}}D_\rho _\mu \varphi D_\lambda _\nu \varphi .$$
(4.1)
Of course, expression (4.1) can be written in terms of the scalar curvature as given in (3.2). However, for our purposes here, it is better to keep it in the way it is. Let us calculate the momenta. For convenience, we take them as surface terms lying on a hypersurface orthogonal to the time direction Landau

$$\overline{\delta }S^{(2)}=_{t_0}^t𝑑t𝑑xϵ^{\mu \nu }\frac{ϵ^{\rho \lambda }}{\sqrt{g}}D_\rho _\mu \varphi \overline{\delta }\left(D_\lambda _\nu \varphi \right),$$
(4.2)
where $`\overline{\delta }`$ means variation of fields under derivatives (because just these terms will contribute to the momenta). In this way, we have for delta term which appears in the equation above

$`\overline{\delta }\left(D_\lambda _\nu \varphi \right)`$ $`=`$ $`D_\lambda _\nu \overline{\delta }\varphi \overline{\delta }\mathrm{\Gamma }_{\lambda \nu }^\alpha _\alpha \varphi `$
$`=`$ $`D_\lambda _\nu \delta \varphi g^{\alpha \beta }\left(_\nu \overline{\delta }g_{\lambda \beta }{\displaystyle \frac{1}{2}}_\beta \overline{\delta }g_{\lambda \nu }\right)_\alpha \varphi ,`$
where in the last step it was used the symmetry between $`\lambda `$ and $`\nu `$ indices, given in (4.2). Developing in a separate way the terms which will appear after replacing (LABEL:4.3) into (4.2), we have

$`{\displaystyle _{t_0}^t}𝑑t{\displaystyle 𝑑xϵ^{\mu \nu }\frac{ϵ^{\rho \lambda }}{\sqrt{g}}D_\rho _\mu \varphi D_\lambda _\nu \overline{\delta }\varphi }`$
$`={\displaystyle _{t_0}^t}dt{\displaystyle }dx\{_\lambda \left(ϵ^{\mu \nu }{\displaystyle \frac{ϵ^{\rho \lambda }}{\sqrt{g}}}D_\rho _\mu \varphi _\nu \overline{\delta }\varphi \right)`$
$`_\nu \left(ϵ^{\mu \nu }{\displaystyle \frac{ϵ^{\rho \lambda }}{\sqrt{g}}}D_\lambda D_\rho _\mu \varphi \overline{\delta }\varphi \right)],`$
$`={\displaystyle }dx\{[_1\left({\displaystyle \frac{1}{\sqrt{g}}}D_1_0\varphi \right)+{\displaystyle \frac{ϵ^{\rho \lambda }}{\sqrt{g}}}D_\lambda D_\rho _1\varphi ]\overline{\delta }\varphi `$
$`+{\displaystyle \frac{1}{\sqrt{g}}}D_1_1\varphi \overline{\delta }\dot{\varphi }\}.`$ (4.4)
For the first term involving $`\overline{\delta }g_{\alpha \beta }`$, we obtain

$`{\displaystyle _{t_0}^t}𝑑t{\displaystyle 𝑑x\frac{ϵ^{\mu \nu }ϵ^{\rho \lambda }}{\sqrt{g}}D_\rho _\mu \varphi g^{\alpha \beta }_\nu \overline{\delta }g_{\lambda \beta }_\alpha \varphi }`$
$`={\displaystyle _{t_0}^t}𝑑t{\displaystyle 𝑑x_\nu \left(\frac{ϵ^{\mu \nu }ϵ^{\rho \lambda }}{\sqrt{g}}D_\rho _\mu \varphi g^{\alpha \beta }_\alpha \varphi \overline{\delta }g_{\lambda \beta }\right)},`$
$`={\displaystyle 𝑑x\frac{ϵ^{\rho \mu }}{\sqrt{g}}D_\rho _1\varphi g^{\alpha \nu }_\alpha \varphi \overline{\delta }g_{\mu \nu }}.`$ (4.5)
Similarly, for the last term the result is

$`{\displaystyle _{t_0}^t}𝑑t{\displaystyle 𝑑x\frac{ϵ^{\mu \nu }ϵ^{\rho \lambda }}{\sqrt{g}}D_\rho _\mu \varphi g^{\alpha \beta }_\alpha \varphi _\beta \overline{\delta }g_{\lambda \nu }}`$
$`={\displaystyle 𝑑x\frac{ϵ^{\lambda \nu }ϵ^{\rho \mu }}{\sqrt{g}}D_\rho _\lambda \varphi g^{\alpha 0}_\alpha \varphi \overline{\delta }g_{\mu \nu }}.`$ (4.6)
Replacing (4.4)-(4.6) into (4.2) we identify the momenta conjugate to $`\varphi `$, $`\dot{\varphi }`$ and $`g_{\mu \nu }`$ (which are the coefficients of the variations of the corresponding fields). Denoting them respectively by $`p`$, $`p^{(1)}`$ and $`\pi ^{\mu \nu }`$, we have

$`p=_1\left({\displaystyle \frac{1}{\sqrt{g}}}D_1_0\varphi \right){\displaystyle \frac{ϵ^{\mu \nu }}{\sqrt{g}}}D_\mu D_\nu _1\varphi ,`$ (4.7)
$`p^{(1)}={\displaystyle \frac{1}{\sqrt{g}}}D_1_1\varphi ,`$ (4.8)
$`\pi ^{\mu \nu }={\displaystyle \frac{1}{2\sqrt{g}}}_\alpha \varphi (ϵ^{\rho \mu }g^{\alpha \nu }D_\rho _1\varphi `$
$`+ϵ^{\rho \nu }g^{\alpha \mu }D_\rho _1\varphi +ϵ^{\lambda \nu }ϵ^{\rho \mu }g^{\alpha 0}D_\rho _\lambda \varphi ).`$ (4.9)
Notice that $`\pi ^{\mu \nu }`$ is symmetric in $`\mu `$ and $`\nu `$ indices, as it should be. However, the identification of the primary first-class constraints (responsible for the symmetries of the theory) is not so apparent as it was in the case of the action (2.1), where the constraints are simple given by (2.2). Let us explicitly write down the components of $`\pi ^{\mu \nu }`$.

$`\pi ^{01}`$ $`=`$ $`{\displaystyle \frac{1}{2\sqrt{g}}}g^{1\alpha }_\alpha \varphi D_1_1\varphi ,`$ (4.10)
$`\pi ^{00}`$ $`=`$ $`{\displaystyle \frac{1}{2\sqrt{g}}}g^{0\alpha }_\alpha \varphi D_1_1\varphi ,`$ (4.11)
$`\pi ^{11}`$ $`=`$ $`{\displaystyle \frac{1}{2\sqrt{g}}}_\alpha \varphi \left(2g^{\alpha 1}D_0_1\varphi +g^{\alpha 0}D_0_0\varphi \right).`$ (4.12)
If one combine (4.8), (4.10) and (4.11), we observe that $`\pi ^{01}`$ and $`\pi ^{00}`$ are actually constraints ($`\pi ^{11}`$ is not). Denoting these constraints by $`\psi `$ and $`\chi `$ we have

$`\psi =\pi ^{01}+{\displaystyle \frac{1}{2}}g^{1\alpha }_\alpha \varphi p^{(1)}0,`$ (4.13)
$`\chi =\pi ^{00}+{\displaystyle \frac{1}{2}}g^{0\alpha }_\alpha \varphi p^{(1)}0,`$ (4.14)
where the symbol $``$ means weakly zero Dirac . Using the fundamental Poisson brackets

$`\{\dot{\varphi }(x),p^{(1)}(y)\}=\delta (xy),`$ (4.15)
$`\{g_{\mu \nu }(x),\pi ^{\rho \lambda }\}={\displaystyle \frac{1}{2}}\left(\delta _\mu ^\rho \delta _\nu ^\lambda +\delta _\mu ^\lambda \delta _\nu ^\rho \right)\delta (xy)`$ (4.16)
and also

$$\{g^{\mu \nu }(x),\pi ^{\rho \lambda }\}=\frac{1}{2}\left(g^{\mu \rho }g^{\nu \lambda }+g^{\mu \lambda }g^{\nu \rho }\right)\delta (xy),$$
(4.17)
obtained from (4.16), we can actually show that $`\psi `$ and $`\chi `$ are primary first-class constraints, which are related to the diffeomorphism symmetry. In the case of action (2.1), there was one extra primary first-class constraints responsible for the conformal symmetry. In the present case, this symmetry is missing.

The next step of the present procedure would be the obtainment of secondary constraints (what has to be done in a Hamiltonian formalism keeping velocities Barc2 ) and the construction of the diffeomorphism algebra. This is just a question of a hard algebraic work, but we belive that nothing new will appear and the diffeomorphism algebra should be verified, since the theory is invariant by general coordinate transformation.

## V Conclusion

In this paper we have studied noncommutative fields in curved space. We consider a formulation where the noncommutativity does not affect the metric but it does on the geommetry. This procedure might be taken as an intermediary step to obtain a noncommutative gravitational theory without inconsistencies on the flat space.

We have deal with a spacetime dimensions $`D=2`$, where the flat limit of the noncommutative parameter is a constant tensor. Of course, this could have been presented for higher dimensions. We envisage two independent formulations for such theories. One of them is to consider the parameter $`\theta _{\mu \nu }`$ depending on $`x`$. In this way, it would be some antisymmetric tensor field. Another possibility, followed by Carlson et al. Carlson , is to consider $`\theta _{\mu \nu }`$ as an independent spacetime coordinate Barc3 .

###### Acknowledgements.

This work is supported in part by Conselho Nacional de Desenvolvimento Científico e Tecnológico - CNPq (Brazilian Research Agency). One of us, J.B.-N. has also the support of PRONEX 66.2002/1998-9. I am in debt with Professor R. Amorim for many helpful discussions.

## Appendix A Useful relations in curved space

For some vector and antisymmetric tensor field ($`V^\mu `$ and $`F^{\mu \nu }`$), we have

$$D_\mu \left(\sqrt{g}V^\mu \right)=_\mu \left(\sqrt{g}V^\mu \right),$$
(A.1)
$$D_\nu \left(\sqrt{g}F^{\mu \nu }\right)=_\nu \left(\sqrt{g}F^{\mu \nu }\right).$$
(A.2)
In the particular case of $`D=2`$, the following relations are true

$$\frac{ϵ^{\mu \nu }}{\sqrt{g}}\frac{ϵ^{\rho \lambda }}{\sqrt{g}}=g^{\mu \lambda }g^{\nu \rho }g^{\mu \rho }g^{\nu \lambda },$$
(A.3)
$$R_{\mu \nu \rho \lambda }=\left(g_{\mu \rho }g_{\nu \lambda }g_{\mu \lambda }g_{\nu \rho }\right)\frac{R}{2},$$
(A.4)
$$R_{\mu \nu }=g_{\mu \nu }\frac{R}{2},$$
(A.5)