How do we find a good estimate for the parameters $\theta$ for a chosen model class? The general recipe is as follows. The first step is to introduce a model that makes sense biophysically, and incorporates our prior knowledge in a tractable manner. Next we write down the likelihood of the observed data given the model parameters, along with a prior distribution that encodes our prior beliefs about the model parameters. Finally, we compute the posterior distribution of the model parameters given the observed data, using Bayes’ rule, which states that

the left-hand side is the desired posterior distribution, while the right-hand side is just the product of the likelihood and the prior.

In the current setting, we need to write down the likelihood $p(D|X,\mbox{\boldmath\(k\)})$ of the observed spike data $D$ given the model parameter $k$ and the observed set of stimuli summarized in the matrix $X$, and then we may employ standard likelihood optimization methods to obtain the maximum likelihood (ML) or maximum a posteriori (MAP) solutions for $k$, defined by

where the maximization runs over all possible parameter choices.

We assume that spike counts per bin follow a conditional Poisson distribution, given $\rho(t)$:

cf. text and exercises of Chapter 7. For example, with the rate parameter of the Poisson distribution given by a GLM or SRM model $\rho(t)=f(\mbox{\boldmath\(k\)}\cdot\mbox{\boldmath\(x\)}_{t})$, we have

Here $\Pi_{t}$ denotes the product over all time steps. We recall that, by definition, $n_{t}!=1$ for $n_{t}=0$.

Our aim is to optimize the parameter $k$. For a given observed spike train, the spike count numbers $n_{t}$ are fixed. In this case, $({\text{d}}t)^{n_{t}}/(n_{t})!$ is a constant which is irrelevant for the parameter optimization. If we work with a fixed time step ${\text{d}}t$ and drop all units, we can therefore reshuffle the terms and consider the logarithm of the above likelihood

If we choose the time step ${\text{d}}t$ shorter than the half-width of an action potential, say 1 ms, the spike count variable $n_{t}$ can only take the values zero or one. For small ${\text{d}}t$, the likelihood of Eq. (10.22) is then identical to that of the SRM with escape noise, defined in Chapter 9; cf. Eqs. (9.10) and (9.15) - (9.17).

We don’t have an analytical expression for the maximum of the likelihood defined in Eq. (10.22), but nonetheless we can numerically optimize this function quite easily if we are willing to make two assumptions about the nonlinear function $f(.)$. More precisely, if we assume that

(i) $f(u)$ is a convex (upward-curving) function of its scalar argument $u$, and

(ii) $\log f(u)$ is concave (downward-curving) in $u$,

then the loglikelihood above is guaranteed to be a concave function of the parameter $k$, since in this case the loglikelihood is just a sum of concave functions of $k$ (383).

This implies that the likelihood has no non-global maximum (also called local maximum). Therefore the maximum likelihood parameter $\hat{\mbox{\boldmath\(k\)}}_{\rm ML}$ may be found by numerical ascent techniques; cf. Fig. 10.2. Functions $f(.)$ satisfying these two constraints are easy to think of: for example, the standard linear rectifier and the exponential function both work.

Fitting model parameters proceeds as follows: we form the (augmented) matrix $X$ where each row is now

Similarly, the parameter vector is in analogy to Eq. (10.10)

here $b=u_{\rm rest}-\vartheta$ is a constant offset term which we want to optimize.

We then calculate the log-likelihood

and compute the ML or maximum a posteriori (MAP) solution for the model parameters $k$ by a concave optimization algorithm. Note that, while we still assume that the conditional spike count $n_{t}$ within a given short time bin is drawn from a one-dimensional $Poiss(\rho(t){\text{d}}t)$ distribution given $\rho(t)$, the resulting model displays strong history effects (since $\rho(t)$ depends on the past spike trains) and therefore the output of the model, considered as a vector of counts $D=\{n_{t}\}$, is no longer a Poisson process, unless $\mbox{\boldmath\(\eta\)}=0$. Importantly, because of the refractory effects incorporated by a strong negative $\eta$ at small times, the spike count variable $n_{t}$ cannot take a value larger than $1$ if ${\text{d}}t$ is in the range of one or a few milliseconds. Therefore, we can expect that interspike intervals are correctly reproduced in the model; see Fig. 10.5 for an example application.

Finally, we may expand the definition of $X$ to include observations of other spike trains, and therefore develop GLMs not just of single spike trains, but network models of how populations of neurons encode information jointly (95; 379; 521; 399). The resulting model is summarized as follows: Spike counts are conditionally Poisson distributed given $\rho_{i}(t)$ $n_{i,t}\sim Poiss(\rho_{i}(t){\text{d}}t)$ with a firing rate

Here, $\rho_{i}(t)$ denotes the instantaneous firing rate of the $i$-th cell at time $t$ and $\mbox{\boldmath\(k\)}_{i}$ is the cell’s linear receptive field including spike-history effects; cf. Eq. (10.25). The net effect of a spike of neuron $i^{\prime}$ onto the membrane potential of neuron $i$ is summarized by $\epsilon_{i^{\prime},j}$; these terms are summed over all past spike activity $n_{i^{\prime},t-j}$ in the population of cells from which we are recording simultaneously. In the special case that we record from all neurons in the population, $\epsilon_{i^{\prime},j}$ can be interpreted as the excitatory or inhibitory postsynaptic potential caused by a spike of neuron $i^{\prime}$ a time $j\,{\text{d}}t$ earlier.

In the linear regression case it is well-known that estimates of the components of the parameter vector $k$ can be quite noisy when the dimension of $k$ is large. The noisiness of the estimate $\hat{\mbox{\boldmath\(k\)}}_{ML}$ is roughly proportional to the dimensionality of $k$ (the number of parameters in $k$ that we need to estimate from data) divided by the total number of observed samples (382). The same ‘overfitting’ phenomenon (estimator variability increasing with number of parameters) occurs in the GLM context. A variety of methods have been introduced to ‘regularize’ the estimated $k$, to incorporate prior knowledge about the shape and/or magnitude of the true $k$ to reduce the noise in $\hat{\mbox{\boldmath\(k\)}}_{ML}$. One basic idea is to restrict $k$ to lie within a lower-dimensional subspace; we then employ the same fitting procedure to estimate the coefficients of $k$ within this lower-dimensional basis (model selection procedures may be employed to choose the dimensionality of this subspace (521)).

A slightly less restrictive approach is to maximize the posterior

(with $k$ allowed to take values in the full original basis), instead of the likelihood $p(D|X,\mbox{\boldmath\(k\)})$; here $p(\mbox{\boldmath\(k\)})$ encodes our a priori beliefs about the true underlying $k$.

It is easy to incorporate this maxima a posteriori idea in the GLM context (383): we simply maximize

Whenever $\log p(\mbox{\boldmath\(k\)})$ is a concave function of $k$, this ‘penalized’ likelihood (where $\log p(\mbox{\boldmath\(k\)})$ acts to penalize improbable values of $k$) is a concave function of $k$, and ascent-based maximization may proceed (with no local maximum) as before; cf. Fig. 10.4.

Suppose an experimentalist has injected a time-dependent current $I(t)$ into a neuron and has recorded with a second electrode the voltage $u^{\rm exp}(t)$ of the same neuron. The voltage trajectory contains spikes $S(t)=\sum_{f}\delta(t-t^{(f)})$ with firing times $t^{(1)},t^{(2)},...t^{(N)}$. The natural approach would be to write down the joint likelihood of observing both the spike times and the subthreshold membrane potential (381). A simpler approach would be to maximize the likelihood of observing the spike separately from the likelihood of observing the membrane potential.

Given the input $I(t)$ and the spike train $S(t)$, the voltage of a SRM is given by Eq. (10.9) and we can adjust the parameters of the filter $\kappa$ and $\eta$ so as to minimize the squared error Eq. (10.5). We now fix the parameters for the membrane potential and maximize the likelihood of observing the spike times given our model voltage trajectory $u(t)$. We insert $u(t)$ into the escape rate function $\rho(t)=f(u(t)-\vartheta(t))$ which contains the parameters of the threshold

We then calculate the log-likelihood

and compute the ML or maximum a posteriori (MAP) solution for the model parameters $k$ (which are the parameters of the threshold - the subthreshold voltage parameters are already fixed) by an optimization algorithm for concave functions.

Fig. 10.5 shows an example application. Both voltage in the subthreshold regime and spike times are nicely reproduced. Therefore, we can expect that interspike intervals are correctly reproduced as well. In order to quantify the performance of neuron models, we need the develop criteria of ‘goodness of fit’ for subthreshold membrane potential, spike timings, and possibly higher-order spike train statistics. This is the topic of section 10.3; we will return to similar applications in the next chapter.

The GLM encoding framework described here can be extended in a number of important directions. We briefly describe two such directions here.

First, as we have described the GLM above, it may appear that the model is restricted to including only linear dependencies on the stimulus $\mbox{\boldmath\(x\)}_{t}$, through the $\mbox{\boldmath\(k\)}\cdot\mbox{\boldmath\(x\)}_{t}$ term. However, if we modify our input matrix $X$ once again, to include nonlinear transformations $\mathcal{F}_{j}(\mbox{\boldmath\(x\)})$ of the stimulus $x$, we may fit nonlinear models of the form

efficiently by maximizing the loglikelihood $\log p(D|X,\mbox{\boldmath\(a\)})$ with respect to the weight parameter $a$ (560; 12). Mathematically, the nonlinearities $\mathcal{F}_{j}(\mbox{\boldmath\(x\)})$ may take essentially arbitrary form; physiologically speaking, it is clearly wise to choose $\mathcal{F}_{j}(\mbox{\boldmath\(x\)})$ to reflect known facts about the anatomy and physiology of the system (e.g., $\mathcal{F}_{j}(\mbox{\boldmath\(x\)})$ might model inputs from a presynaptic layer whose responses are better-characterized than are those of the neuron of interest (450)).

Second, in many cases it is reasonable to include terms in $X$ that we may not be able to observe or calculate directly (e.g., intracellular noise, or the dynamical state of the network); fitting the model parameters in this case requires that we properly account for these ‘latent,’ unobserved variables in our likelihood. While inference in the presence of these hidden parameters is beyond the scope of this chapter, it is worth noting that this type of model may be fit tractably using generalizations of the methods described here, at the cost of increased computational complexity, but the benefit of enhanced model flexibility and realism (484; 563; 536).