Automatic Construction of Predictive Neuron Models through Large Scale Assimilation of Electrophysiological Data

We report on the construction of neuron models by assimilating electrophysiological data with large-scale constrained nonlinear optimization. The method implements interior point line parameter search to determine parameters from the responses to intracellular current injections of zebra finch HVC neurons. We incorporated these parameters into a nine ionic channel conductance model to obtain completed models which we then use to predict the state of the neuron under arbitrary current stimulation. Each model was validated by successfully predicting the dynamics of the membrane potential induced by 20–50 different current protocols. The dispersion of parameters extracted from different assimilation windows was studied. Differences in constraints from current protocols, stochastic variability in neuron output, and noise behave as a residual temperature which broadens the global minimum of the objective function to an ellipsoid domain whose principal axes follow an exponentially decaying distribution. The maximum likelihood expectation of extracted parameters was found to provide an excellent approximation of the global minimum and yields highly consistent kinetics for both neurons studied. Large scale assimilation absorbs the intrinsic variability of electrophysiological data over wide assimilation windows. It builds models in an automatic manner treating all data as equal quantities and requiring minimal additional insight.


II.
Comparing model activation/inactivation functions used in data assimilation and in Physiology

CaT: Gate variables
Model functions Empirical fitting functions

Recovery times Model functions in nonlinear optimization
Empirical fitting functions    Figure S2: A comparison of the model functions used in nonlinear optimization (Eqs.6-9 of the manuscript, full green lines) with those used to fit the activation (blue symbols) and inactivation (red symbols) of thalamocortical neurons 1 . Both functions and their parameters are given in the above table.

K2:
Gate variables Model functions Empirical fitting functions

Recovery times
Model functions Empirical fitting functions

III.
Frequency spectrum of the current protocols used to assimilate models of N1 and N2 Figure S3: Power spectrum of the current protocols used to extract the parameters of N1 and N2.
The low frequency limit (1Hz) corresponds to the width of the assimilation window (1s). The high frequency limit (50kHz) is the sampling frequency of the data. Both current protocols exhibit a low pass spectrum with a cut off frequency of 200Hz. Hence the predictive power of models constructed by assimilating these data is expected to decrease for currents protocols varying on very short timescales, smaller than 5ms -see Fig.3(f).  The width of each graph (190ms-1990ms) is the assimilation window.

IX. Covariance matrices of N1 and N2
Figure S10: Exponential decay of eigenvalue spectrum The best fit line (blue line) is a sum of three exponentials with different rates of decay.

Figure S11: The covariance matrices of neurons N1 and N2
Each covariance matrix was constructed from a statistical sample of 84 sets of parameters extracted over a range of assimilation windows. Parameters were ordered along the rows and columns of matrices in a slightly different order from Table II and as follows: 2. gNa,3. gNaP,4. ENa,5. gK1,6. gK2,7. gK3,8. EK,9. gL,10. EL,11. gout,12. gHCN,13,A,14. Vm (NaT) Finite off-diagonal elements (horizontal and vertical green lines in the covariance matrices) indicate the correlated parameters. Juxtaposing the covariance matrices of N1 and N2, suggests that correlated parameters do not occur at random but happen to be the same in N1 and N2 (horizontal arrows). These are the parameters that define the recovery time constant of ion channels: (NaT) m & m and to a lesser extent the conductance of the transient potassium channel: (K1) gK1.

Figure S12: Minimum finding in a noisy environment
Experimental noise has the characteristics of a residual temperature T which prevents reaching the global minimum of the objective function through direct parameter search. Instead, the tip of random vector P lies on the surface of a K-dimensional ellipsoid c=kT. Each set of parameters P extracted from one assimilation window gives one point on this surface (red dots). (a) Once a sufficiently large statistical sample of parameters (R>L) maps the ellipsoid, the maximum likelihood expectation <P> yields the global minimum p * . The parameter search has converged near the global minimum because the covariance matrices and eigenvalue spectra exhibit identical patterns in N1 and N2. If the parameter search had converged to local minima instead, (b) the covariance matrices of N1 and N2 would have different structures reflecting the different environments around these minima.
Insert in the cost function: The cross term cancels as it is proportional to the noise average which is zero. The second term on the RHS gives the variance of the membrane voltage which when driven by thermal fluctuations is given by Nyquist theorem: