Inevitably, the data we collect is noisy, and equally the model we select to describe the production of those data has errors. This means we must, at the outset, address a conditional probability distribution P(X|Y) as our goal in the data assimilation transfer from Y to the model. In Abarbanel [3] we describe how to use the Markov nature of the model x(n) → x(n+1) =f(x(n), q) and the definition of conditional probabilities to derive the recursion relation:

where we have identified CMI(a,b|c)=log[P(a,b|c)P(a|c)P(b|c)]. This is Shannon's conditional mutual information [10] telling us how many bits (for log₂) we know about a when observing b conditioned on c. For us a = {y(n + 1)}, b = {x(n + 1), X(n)}, c = {Y(n)}. We can simplify this further with the assumption that an observation at any time depends only on the state of the system.

Using this recursion relation to move backwards through the observation window from t_F = t₀ + NΔt through the measurements at times τ_k to the start of the window at t₀, we may write, up to factors independent of X

If we now write P(X|Y) ∝ exp[−A(X)] where A(X), the negative of the log likelihood, we call the action, then conditional expected values for functions along the path X are defined by

dX=∏n=0NdDx(n), and all factors in the action independent of X cancel out here. The action takes the convenient expression

which is the sum of the terms which modify the conditional probability distribution when an observation is made at t = τ_k and the sum of the stochastic version of x(n) → x(n + 1) − f(x(n), q) and finally the distribution when the observation window opens at t₀.

What quantities G(X) are of interest? One natural one is the path G(X) = X_μ; μ = {a, n} itself; another is the covariance around that mean 〈Xμ〉=X̄μ=〈Xμ〉:〈(Xμ-X̄μ)(Xν-X̄ν)〉. Other moments are of interest, of course. If one has an anticipated form for the distribution at large X, then G(X) may be chosen as a parametrized version of that form and those parameters determined near the maximum of P(X|Y).

The action simplifies to what we call the “standard model” of data assimilation when (1) observations y are related to their model counterparts via Gaussian noise with zero mean and diagonal precision matrix R_m, and (2) model errors are associated with Gaussian errors of mean zero and diagonal precision matrix R_f:

If we have knowledge of the distribution P(x(0)) at t₀ we may add it to this action. If we have no knowledge of P(x(0)), we may take its distribution to be uniform over the dynamic range of the model variables, then it, as here, is absent, canceling numerator and denominator in Equation (6).

Our challenge is to perform integrals such as Equation (6). One should anticipate that the dominant contribution to the expected value comes from the maxima of P(X|Y) or, equivalently the minima of A(X). At such minima, the two contributions to the action, the measurement error and the model error, balance each other to accomplish the explicit transfer of information from the former to the latter.

We note, as before, that when f(x(n), q) is nonlinear in X, as it always is in interesting examples, the expected value integral Equation (6) is not Gaussian. So, some thinking is in order to approximate this high dimensional integral. We turn to that now. After consideration of methods to do the integral, we will return to a variety of examples taken from neuroscience.

The two generally useful methods available for evaluating this kind of high dimensional integral are Laplace's method [5] and Monte Carlo techniques [6, 11, 12]. We address them in order. We also add our own new and useful versions of the methods.

To locate the minima of the action A(X) = −log[P(X|Y)] we must seek paths X^(j); j = 0, 1, … such that ∂A(X)/∂X|X(j)=0, and then check that the second derivative at X^(j), the Hessian, is a positive definite matrix in path coordinates. The vanishing of the derivative is a necessary condition.

Laplace's method [5] expands the action around the X^(j) seeking the path X⁽⁰⁾ with the smallest value of A(X). The contribution of X⁽⁰⁾ to the integral Equation (6) is approximately exp[A(X⁽¹⁾) − A(X⁽⁰⁾)] bigger than that of the path with the next smallest action.

This sounds more or less straightforward; however, finding the global minimum of a nonlinear function such as A(X) is an NP-complete problem [13]. In a practical sense one cannot expect to succeed with such problems. However there is an attractive feature of the form of A(X) that permits us to accomplish more.

We now discuss two algorithmic approaches to implementing Laplace's method.

Looking at Equation (8) we see that if the precision of the model is zero, R_f = 0, the action is quadratic in the L measured variables x_l(n) and independent of the remaining states. The global minimum of such an action comes with x_l(τ_k) = y_l(τ_k) and any choice for the remaining states and parameters. Choose the path with these values of x(τ_k) and values from a uniform distribution of the other state variables and parameters covering the expected dynamic range of those, and call it path X_init. In practice, we recognize that the global minimum of A(X) is degenerate at R_f = 0, so we select many initial paths. We choose N_I of them, and initialize whatever numerical optimization program we have selected, to run on each of them. We continue to call the collection of N_I paths X_init.

We call this method precision annealing (PA) [14–17]. It slowly turns up the precision of the model collecting paths at each R_f that emerged from the degenerate global minimum at R_f = 0. In practice it is able to track N_I possible minima of A(X) at each R_f. When not enough information is presented to the model, that is L is too small, there are many local minima at all R_f. This is a manifestation of the NP-completeness of the minimization of A(X) problem. None of the minima may dominate the expected value integral of interest.

As L increases, and enough information is transmitted to the model, for large R_f one minimum appears to stand out as the global minimum, and the paths associated with that smallest minimum yields good predictions. We note that there are always paths, not just a single path, as we have a distribution of paths, N_I of which are sampled in the PA procedure, within a variation of size 1/Rm. A clear example of this is seen in Shirman [18] in a small, illustrative model.

In the even that the chosen model is inconsistent with the data, or there is too much noise in the model error term, a single minimum of the action will not appear for large R_f. As in the case of too few measurements, there will be multiple local minima. An example of this can be seen in Ye et al. [14].

In meteorology one approach to data assimilation is to add a term to the deterministic dynamics which move the state of a model toward the observations [19]

where u(n) > 0 and vanishes except where a measurement is available. This is referred to as “nudging.” It appears in an ad hoc, but quite useful, manner.

Within the structure we have developed, one may see that the “nudging term” arises through the balance between the measurement error term and the model error term in the action. This is easy to see when we look at the continuous time version of the data assimilation standard model

The extremum of this action is given by the Euler-Lagrange equations for the variational problem [20]

in which the right hand side is the “nudging” term appearing in a natural manner. Approximating the operator δabddt+∂Fb(x(t)∂xa(t) we can rewrite this Euler-Lagrange equation in “nudging” form

We will use both the full variation of the action, in discrete time, as well as its nudging form in our examples below.

Monte Carlo methods [6, 11, 17, 21] are well covered in the literature. We have not used them in the examples in this paper. However, the development of a precision annealing version of Monte Carlo techniques promises to address the difficulties with large matrices for the Jacobian and Hessians required in variational principles (Wong et al., unpublished). When one comes to network problems, about which we comment later, this method may be essential.

A twin experiment is a synthetic numerical experiment meant to mirror the conditions of a laboratory experiment. We use our mathematical model with some informed parameter choices in order to generate numerical data. Noise is added to observable variables in the model, here V(t). These data are then put through our SDA procedure to estimate parameters and unobserved states of the model. The neuron model is now calibrated or completed.

Using the parameters and the full state of the model at the end t_F of an observation window [t₀, t_F], we take a current waveform I_injected(t ≥ t_F) to drive the model neuron and predict the time course of all dynamical variables in the prediction window [t_F, …]. This validation of the model is the critical test of our SDA procedure, here PA. In a laboratory experiment we have no specific knowledge of the parameters in the model and, by definition, cannot observe the unobserved state variables; here we can do that. So, “fitting” the observed data within the observation window [t₀, t_F] is not enough, we must reproduce all states for t ≥ t_F to test our SDA methods.

We use the model laid out in the Supplementary Data Sheet 1. We assume that the neuron has a resting potential of −70 mV and set the initial values for the voltage and each gating variable accordingly. We assume that the internal calcium concentration of the cell is C_in = 0.1 μM. We use an integration time step of 0.02 ms and integrate forward in time using an adaptive Runge-Kutta method [6]. Noise is added to the voltage time course by sampling from a Gaussian distribution N(0,2) in units of mV.

The waveform of the injected current was chosen to have three key attributes: (1) It is strong enough to cause spiking in the neuron, (2) it dwells a long time in a hyperpolarizing region, and (3) its overall frequency content is low enough to not be filtered out by the neuron's RC time constant. On this last point, a neuron behaves like an RC circuit, it has a cut off frequency limited by the time constant of the system. Any input current which has a frequency higher than that of the cut off frequency won't be “seen” by the neuron. The time constant is taken to be the time it takes to spike and return back to 37% above its resting voltage. We chose a current where the majority of the power spectral density exists below 50 Hz.

The first two seconds of our chosen current waveform is a varying hyperpolarizing current. In order to characterize I_h(t) and I_CaT(t), it is necessary to thoroughly explore the region where the current is active. I_h(t) is only activated when the neuron is hyperpolarized. The activation of I_h(t) deinactivates I_CaT(t), thereby allowing its dynamics to be explored. In order to characterize I_Na(t) and I_K(t), it is necessary to cause spiking in the neuron. The depolarizing current must be strong enough to hit the threshold potential for spike activation.

The parameters used to generate the data used in the twin experiment are in Table 1, and the injection current data and the membrane voltage response may be seen in Figures 3A,B.

The numbers chosen for the data assimilation procedure in this paper are α = 1.4 and β ranging from 1 to 100. Rf,0,V=10-4 for voltage and R_{f, 0, j} = 1 for all gating variables. These numbers are chosen because the voltage range is 100 times large than the gating variable range. Choosing a single R_{f, 0} value would result in the gating variable equations being less enforced than the voltage equation by a factor of 10⁴. The α and β numbers are chosen because we seek to make RfRf0 sufficiently large. The α and β values chosen allow RfRf0 to reach 10¹⁵.

During estimation we instructed our methods to estimate the inverse capacitance and estimate the ratio g′=gCm instead of g and C_m independently. That separation can present a challenge to numerical procedures. We also estimated the reversal potentials as a check on the SDA method.

Within our computational capability we can reasonably perform estimates on 50,000 data points. This captures a second of data when Δt = 0.02ms. However, there are time constants in the model neuron which are on order 1 second. In order for us to estimate the behavior of these parameters accurately, we need to see multiple instances of the full response. We need a window on the order of 2–3 s. We can obtain this by downsampling the data. We know from previous results that downsampling can lead to better estimations [38]. We take every i^th data point, depending on the level of downsampling we want to do. In this data assimilation run, we downsampled by a factor of 4 to incorporate 4 s of data in the estimation window.

We look at a plot of the action as a function of β; that is, log[R_f/R_f0]. We expect to see a leveling off of the action [16] as a function of R_f. If the action becomes independent of R_f, we can then explore how well our parameter estimations perform when integrating them forward as predictions of the calibrated model. Looking at the action plot in Figure 4, we can see there is a region in which the action appears to level off, around β = 40. It is in this region where we look for our parameter estimates.

We examine all solutions around this region of β and utilize their parameter estimates in our predictions. We compare our numerical prediction to the “real" data from our synthetic experiment. We evaluate good predictions by finding the correlation coefficient between these two curves. This metric is chosen instead of a simple root mean square error because slight variations in spike timings yield a high amount of error even if the general spiking pattern is correct. The prediction plot and parameters for the best prediction can be seen in Figure 5 and Table 2. The voltage trace in red is the estimated voltage after data assimilation is completed. It is overlayed on the synthetic input data in black. The blue time course is a prediction, starting at the last time point of the red estimated V(t) trace and using the parameter estimates for t ≤ 4, 000ms.

The red curve matches the computed voltage trace quite well. There is no significant deviation in the frequency of spikes, spike amplitudes, or the hyperpolarized region of the cell. Looking at the prediction window, we can see that there is some deviation in the hyperpolarized voltage trace after 7,000 ms. Our our predicted voltage does not become nearly as hyperpolarized as the real data. This is an indication that our parameter estimates for currents activated in this region are not entirely correct. Comparing parameters in Table 2, we can see that E_h is estimated as lower than its actual value. Despite that, we still are able to reproduce neuron behavior fairly well.

Our next use of SDA employs the “nudging” method described in Eq. (9). In this section we used some of the data [Daou and Margoliash (in review)] taken in experiments on multiple HVC_X neurons from different zebra finches. The questions we asked was whether we could, using SDA, identify differences in biophysical characteristics of the birds. This question is motivated by prior biological observations [Daou and Margoliash (in review)].

Using the same HVC_X model as before, we estimated the maximal conductances {g_Na, g_K, g_CaT, g_SK, g_h} holding fixed other kinetic parameters and the Nernst/Reversal potentials. The baseline characteristics of an ion channel are set by the properties of the cell membrane and the complex proteins penetrating the membrane forming the physical channel. Differences among birds would then come from the density or numbers of various channels as characterized by the maximal conductances. If such differences were identified, it would promote further investigation of the biologically exciting proposition that these differences arise in relation to some aspect of the song learning experience of the birds [Daou and Margoliash (in review)].

In Figures 6A,B we display the stimulating current and membrane voltage response from one of 9 neurons in our large sample. The analysis using SDA was of four neurons from one bird and seven neurons from another. The results for {g_CaT, g_Na, g_SK} is displayed in Figure 7. The maximal conductances from one bird are shown in blue and from the other bird, in red. There is a striking difference between the distributions of maximal conductances.

We do not propose here to delve into the biological implications of these results. Nevertheless, we note that the neurons from each bird occupy a small but distinct region of the parameter space (Figure 7). This result and its implications for birdsong learning, and more broadly for neuroscience, are described in Daou and Margoliash (in review). Here, however, we display this result as an example of the power of SDA to address a biologically important question in a systematic, principled manner beyond what is normally achieved in analyses of such data.

To fully embrace the utility of SDA for these experiments, however, further work is needed. A limitation of the present result is that the SDA estimates for g_SK for a subset of the neurons/observations for Bird One reach the bounds of the observation window (Figure 7). Addressing such issues would be prelude to the exciting possibility of estimating more parameters than just the principle ion currents in the Hodgkin-Huxley equations. This could use SDA numerical techniques to calculate, over hours or days, estimates of parameters that could require months or years of work to measure with traditional biological and biophysical approaches, in some cases requiring specialized equipment beyond that available for most in vitro recording set ups. In contrast, applying SDA to such data sets requires only a computer.