Searching for optimal stimuli: ascending a neuron’s response function

Abstract

Many methods used to analyze neuronal response assume that neuronal activity has a fundamentally linear relationship to the stimulus. However, some neurons are strongly sensitive to multiple directions in stimulus space and have a highly nonlinear response. It can be difficult to find optimal stimuli for these neurons. We demonstrate how successive linear approximations of neuronal response can effectively carry out gradient ascent and move through stimulus space towards local maxima of the response. We demonstrate search results for a simple model neuron and two models of a highly selective neuron.

Appendices

Appendix A: Linear kernel is gradient of smoothed function

This appendix provides a proof of the following: if the noise Z is Gaussian noise, then the expected value of the linear kernel $\mathbf{h}(\hat{\mathbf{x}})$ will point in the direction of the gradient of $F_{\text{ave}}(\mathbf{X})=E(F(\mathbf{X}+\mathbf{Z}))$, the smoothed version of F(X). As a consequence, the expected value of the linear kernel will point in the direction in which the function $F_{\text{ave}}$ increases most rapidly.

If one adds zero mean noise Z with covariance matrix C _Z to the original stimulus $\hat{\mathbf{x}}$ and presents the combined stimulus $\mathbf{X}=\hat{\mathbf{x}}+\mathbf{Z}$ to a neuron with spiking probability given by Eq. (1), the expected value of the linear kernel (Eq. (4)) would be

$$ E(\mathbf{h}(\hat{\mathbf{x}})) = C_{\mathbf{Z}}^{-1}E((\mathbf{X}-\hat{\mathbf{x}})(R - \mu_R(\hat{\mathbf{x}}))) = C_{\mathbf{Z}}^{-1}E(\mathbf{Z}R), $$

(8)

where the μ _R drops out become Z is mean zero. Because $E(R|\mathbf{Z})=F(\hat{\mathbf{x}}+\mathbf{Z})$, we can rewrite this as

$$ \begin{array}{rll} E(\mathbf{h}(\hat{\mathbf{x}})) &=& C_{\mathbf{Z}}^{-1}E(\mathbf{Z}F(\hat{\mathbf{x}}+\mathbf{Z})) \\ &=& \int_{\mathbf{R}^n} \rho_{\mathbf{Z}}(\mathbf{Z}) C_{\mathbf{Z}}^{-1}\mathbf{Z} F(\hat{\mathbf{x}}+\mathbf{Z}) d\mathbf{Z}, \end{array} $$

(9)

where ρ _Z is the probability distribution of the noise.

If the noise is Gaussian, then its probability density function is that of the multivariate normal distribution

$$ \rho_{\mathbf{Z}}(\mathbf{Z}) = \frac{1}{(2\pi)^n | C_{\mathbf{Z}}|^{1/2}} \exp\left(-\frac{1}{2} \mathbf{Z}^T C_{\mathbf{Z}}^{-1} \mathbf{Z}\right). $$

We will use this to rewrite the integrand in Eq. (9). Differentiating ρ(Z) by the ith component of Z, we calculate

$$ \begin{array}{rll} {\frac{\partial {\rho_{\mathbf{Z}}(\mathbf{Z})}}{\partial {z_i}}} &=& \frac{1}{(2\pi)^n | C_{\mathbf{Z}}|^{1/2}}\frac{\partial}{\partial {z_i}}\exp\left(-\frac{1}{2} \mathbf{Z}^T C_{\mathbf{Z}}^{-1} \mathbf{Z}\right) \\ &=& \rho_{\mathbf{Z}}(\mathbf{Z}) \frac{\partial}{\partial{z_i}}\left(-\frac{1}{2} \mathbf{Z}^T C_{\mathbf{Z}}^{-1} \mathbf{Z}\right) \\ &=& - \rho_{\mathbf{Z}}(\mathbf{Z}) \left(\frac{C_{\mathbf{Z}}^{-1} + (C_{\mathbf{Z}}^{-1})^T}{2} \mathbf{Z}\right)_i = - \rho_{\mathbf{Z}}(\mathbf{Z}) \left(C_{\mathbf{Z}}^{-1} \mathbf{Z}\right)_i, \end{array} $$

where the subscript i indicates the ith component of the vector in parentheses, and we used the fact that the covariance matrix is symmetric. Combining this for all components, we find that

$$ \rho_{\mathbf{Z}}(\mathbf{Z}) C_{\mathbf{Z}}^{-1} \mathbf{Z} = -\nabla\rho_{\mathbf{Z}}(\mathbf{Z}) $$

so that Eq. (9) becomes

$$ E(\mathbf{h}(\hat{\mathbf{x}})) = -\int_{\mathbf{R}^n} \nabla \rho_{\mathbf{z}}(\mathbf{z}) F(\hat{\mathbf{x}}+\mathbf{z}) d\mathbf{Z}. $$

(10)

Integrate each component by parts and transfer the derivative onto F. The boundary terms disappear since ρ _z (z) goes to zero for large values of z. We end up transferring the gradient onto F and find that the expected value of the linear kernel is the expected value of the gradient of F, averaged over the noise,

$$ E(\mathbf{h}(\hat{\mathbf{X}})) = \int_{\mathbf{R}^n} \rho_{\mathbf{Z}}(\mathbf{Z}) \nabla F(\hat{\mathbf{X}}+\mathbf{Z}) d\mathbf{Z} = E(\nabla F(\hat{\mathbf{X}}+\mathbf{Z})). $$

(11)

This result is related to Stein’s Lemma (Stein 1981).

If we define $F_{\text{ave}}(\mathbf{X})$ to be a version of F(X) that is smoothed by the noise,

$$ F_{\text{ave}}(\mathbf{X}) = \int_{\mathbf{R}^n} \rho_{\mathbf{z}}(\mathbf{z}) F(\mathbf{X}+\mathbf{z}) d\mathbf{z} = E(F(\mathbf{X}+\mathbf{Z})), $$

(12)

then $F_{\text{ave}}(\hat{\mathbf{z}}) = E(F(\hat{\mathbf{x}}+\mathbf{Z})) = E(R) = E(\mu_R(\hat{\mathbf{X}}))$ and $\nabla F_{\text{ave}}(\hat{\mathbf{x}}) = E(\nabla F(\hat{\mathbf{x}}+\mathbf{Z}))= E(\mathbf{h}(\hat{\mathbf{x}}))$. We see that the expected value of the linear kernel is the gradient of $F_{\text{ave}}$ and that the expected value of the linear approximation (3) of F is the differential linear approximation (2) of $F_{\text{ave}}$:

$$ F_{\text{ave}}(\mathbf{X}) \approx F_{\text{ave}}(\hat{\mathbf{x}}) + \nabla F_{\text{ave}}(\hat{\mathbf{x}}) \cdot (\mathbf{X}-\hat{\mathbf{x}}). $$

In particular, this means the expected value of the linear kernel $\mathbf{h}(\hat{\mathbf{x}})$ defined by Eq. (4) will point in the direction in which the smoothed version (Eq. (12)) of F increases most rapidly.

Appendix B: Calculating noise-averaged response along a line

For the line search algorithm, the goal is to locate along the line $\mathbf{x}_s = \hat{\mathbf{x}}_i + s \mathbf{h}(\hat{\mathbf{x}}_i)$ specified by the linear kernel the maximum of

$$ F_{\text{ave}}(\mathbf{X}) =E(F(\mathbf{X}+\mathbf{Z}))= \int_{\mathbf{R}^n} \rho_{\mathbf{Z}}(\mathbf{z}) F(\mathbf{X}+\mathbf{z}) d\mathbf{z}, $$

(13)

where ρ _z (z) is the probability density of a single realization of the noise Z. This search can require many data points. In this appendix, we explain how to reduce the number of data points needed.

The simplest idea would be to estimate $F_{\text{ave}}(\mathbf{x}_s)$ at points along the line by generating multiple realizations of the noise Z and taking the average number of spikes elicited by the stimulus x _s  + Z. In principle, this method would work, but it would require many realizations of the noise for each point, so the line search would be very expensive in terms of the amount of data that would need to be collected from a neuron.

Another way to estimate the integral for $F_{\text{ave}}(\mathbf{x}_s)$, at least in principle, would be to sample Z uniformly over a bounded subset of R ⁿ (that contained most of the mass of ρ _z (z)), and then take the average of the number of spikes elicited by X + Z, weighted by ρ _z (Z). This method would be a grossly inefficient method to estimate the whole integral, as one would waste effort collecting data from unlikely Z and multiplying the result by a very small ρ _z (Z).

We demonstrate how a combination of the two approaches can be used to efficiently estimate a series of values of $F_{\text{ave}}(\mathbf{x}_s)$ along the line $\mathbf{x}_s= \hat{\mathbf{x}}_i + s \mathbf{h}(\hat{\mathbf{x}}_i)$. The first step is to separate integral (13) into two parts: the one dimensional integral in the direction of h and the n − 1 dimensional integral in the orthogonal complement. Let z _h  = z · h be the magnitude of z in the direction of the linear kernel and let z _⊥ = z − z _h h denote the perpendicular components of z. Write the probability distribution ρ _z (z) = ρ _h (z _h ) ρ _⊥(z _⊥| z _h ), where ρ _h (z _h ) is the marginal distribution of z _h and ρ _⊥(Z _⊥| z _h ) is the distribution of z _⊥ conditioned on a value of z _h . The integral (13) evaluated at $\mathbf{X}=\mathbf{x}_s=\hat{\mathbf{x}}_i + s \mathbf{h}$ can hence be rewritten as

$$ \begin{array}{rll} F_{\text{ave}}(\mathbf{x}_s) &=& \int_{-\infty}^\infty dz_h \rho_h(z_h) \int_{\mathbf{R}^{n-1}} d\mathbf{z}_\perp \rho_\perp(\mathbf{z}_\perp | z_h) \\ && \times F(\hat{\mathbf{x}}_i+(s+z_h)\mathbf{h} + \mathbf{z}_\perp). \end{array} $$

(14)

Next, approximate the conditional distribution ρ _⊥(z _⊥| z _h ) as though it did not depend on z _h , replacing it with the marginal distribution ρ _⊥(z _⊥). This approximation will allow us to rewrite the n − 1 dimensional integral so that it does not explicitly depend on s. We simply change variables to u = z _h  + s in the one dimensional integral, obtaining

$$ \begin{array}{rll} F_{\text{ave}}(\mathbf{x}_s) &\approx& \int_{-\infty}^\infty du\, \rho_h(u-s) \int_{\mathbf{R}^{n-1}} d\mathbf{z}_\perp \rho_\perp(\mathbf{z}_\perp )\\ && \times F(\hat{\mathbf{x}}_i+u\mathbf{h} + \mathbf{z}_\perp). \end{array}$$

(15)

Given a fixed $\hat{\mathbf{x}}_i$ and h, the n − 1 dimensional integral is just a function of the parameter u,

$$ I(u) = \int_{\mathbf{R}^{n-1}} \rho_\perp(\mathbf{z}_\perp ) F(\hat{\mathbf{x}}_i+u\mathbf{h} + \mathbf{z}_\perp) d\mathbf{z}_\perp, $$

(16)

and the approximation of $F_{\text{ave}}(\mathbf{x}_s)$ is just a convolution of I(u) with the marginal density ρ _h

$$ F_{\text{ave}}(\mathbf{x}_s) \approx \int_{-\infty}^\infty \rho_h(u-s) I(u)du.$$

(17)

To minimize the amount of data needed to simultaneously estimate $F_{\text{ave}}(\mathbf{x}_s)$ for many s, we collect data to estimate I(u) using Eq. (16) rather than estimating Eq. (13) directly. We sample data along a uniformly spaced grid of u and weight the results proportional to ρ _h (u − s), according to Eq. (17). This procedure decreases the efficiency for calculating a single $F_{\text{ave}}(\mathbf{x}_s)$, as some samples will be multiplied by a small ρ _h (u − s). However, since we can reuse the same data points to estimate $F_{\text{ave}}(\mathbf{x}_s)$ at many different points x _s , the procedure increases the efficiency for estimating a whole series of $F_{\text{ave}}(\mathbf{x}_s)$.

To estimate $F_{\text{ave}}(\mathbf{x}_s)$ via Eq. (17), we need to approximate I(u). We use the following algorithm to sample data whose expected value is I(u). We generate one independent realization Z of the noise and subtract off its projection onto h, creating Z _⊥ = Z − (Z ·h) h. Then, we present the stimulus $\mathbf{X} = \hat{\mathbf{x}}_i+u\mathbf{h} + \mathbf{Z}_\perp$ to the neuron and record the number of spikes elicited, denoting the result by S(u). In this way, E(S(u)) = I(u). Although we would have to repeat this step many times and average S(u) to obtain a reasonable estimate of I(u), we sample only one data point per u because we aren’t interested in I(u) itself.

The estimate of $F_{\text{ave}}(\mathbf{x}_s)$ via Eq. (17) also requires the calculation of the marginal distribution ρ _h (z _h ) of z _h  = Z ·h, which could depend on the direction h. Depending on the distribution ρ(z) of the noise, the marginal distribution could be calculated analytically or estimated via Monte Carlo simulation using many realizations of the noise Z. (Since this sampling doesn’t involve the neuron, the concerns for data efficiency don’t apply here.) However, as discussed in Section 2.2.1, using non-Gaussian noise Z for $F_{\text{ave}}(\mathbf{X}_s)$ may result in artifacts. One might still use non-Gaussian Z if the neuron didn’t respond well to x _s  + Z for Gaussian Z. Even so, the calculation at this point does not involve stimulating the neuron, which means the argument of obtaining a robust neuronal response does not apply to choice of ρ _h . Hence, we recommend approximating ρ _h using a normal distribution with standard deviation set to the σ of the noise magnitude.

Combining these calculations, we estimate $F_{\text{ave}}(\mathbf{x}_s)$ along the line via Eq. (17) as follows. First, for a given step size Δs, we sample I(u) at grid points kΔs, obtaining S(kΔs). Each S(kΔs) is the number of spikes elicited from the neuron from a single presentation of the stimulus $\mathbf{X} =\hat{\mathbf{x}}_i + k\Delta s \mathbf{h} + \mathbf{Z}_\perp$, where for each k, we generate an independent realization of Z to calculate Z _⊥ = Z − (Z·h)h. Second, we numerically sample ρ _h (jΔs) for integer j between − J and J, where J = ⌈2σ/Δs ⌉, normalizing the result to vector components $\hat{\rho}_j$ that sum to 1. This estimate of $\hat{\rho}_j$ does not involve any presentation of the stimulus to the neuron. We obtain the estimate of $F_{\text{ave}}(\mathbf{x}_s)$, for $\mathbf{x}_s = \hat{\mathbf{x}}_i + s\mathbf{h}$:

$$ F_{\text{ave}}(\mathbf{x}_s) \approx T(s)=T(k\Delta s)= \sum\limits_{j=-J}^J \hat{\rho}_jS((k-j)\Delta s). $$

(18)

As the spike counts S(jΔs) appear in estimates for many T(kΔs), Eq. (18) can be used to efficiently estimate $F_{\text{ave}}$ along the line x _s .

T(s) is the weighted average of 2J + 1 samples of spike counts S(·). The differing weights mean that T(s) is not invariant under permutation of the S, making it trickier to estimate the variance in T via resampling. We estimated confidence intervals of T via an approximate jackknife algorithm (Efron 1982). We divided the samples into ten parts, where each part is composed of every tenth S(jΔs). Although the dependence of T on those ten parts was not exactly invariant under permutation, each part sampled the range of weights $\hat{\rho}_j$. We obtained a jackknife estimate of the standard deviation δT(s) of T(s) based on calculations of T that omit one of those parts.

We calculate T(s) initially for 0 ≤ s ≤ 4σ. Then, we continue to calculate T(s) for larger s and compare the new value T(s) with those from a previous window of (s − 4σ,s). Denoting the minimum and maximum value of T in the window by m and M, respectively, we stop calculating additional T(s) if T(s) is smaller than in the window (T(s) ≤ m) or if T(s) hasn’t changed much in the window and doesn’t appear to be rising (M − m ≤ mean δT in window and T(s) < maximum T − δT in window).

When either criterion is met, we stop the line search and find $s_{\text{max}}$, defined as the value of s that maximizes T(s) over all values calculated. If $T(s_{\text{max}}) - \delta T(s_{\text{max}})$ is greater than our previous maximum estimate of $F_{\text{ave}}(\mathbf{X})$, then we conclude we have found a new stimulus $\mathbf{x}_{s_{\text{max}}} = \hat{\mathbf{x}}_i + s_{\text{max}} \mathbf{h}(\hat{\mathbf{x}}_i)$ that appears to elicit a higher firing rate in the neuron and set $\hat{\mathbf{x}}_{i+1} = \mathbf{x}_{s_{\text{max}}}$.