Inferring the dynamics of ionic currents from recursive piecewise data assimilation of approximate neuron models

We construct neuron models from data by transferring information from an observed time series to the state variables and parameters of Hodgkin-Huxley models. When the learning period completes, the model will predict additional observations and its parameters uniquely characterise the complement of ion channels. However, the assimilation of biological data, as opposed to model data, is complicated by the lack of knowledge of the true neuron equations. Reliance on guessed conductance models is plagued with multi-valued parameter solutions. Here, we report on the distributions of parameters and currents predicted with intentionally erroneous models, over-specified models, and an approximate model fitting hippocampal neuron data. We introduce a recursive piecewise data assimilation (RPDA) algorithm that converges with near-perfect reliability when the model is known. When the model is unknown, we show model error introduces correlations between certain parameters. The ionic currents reconstructed from these parameters are excellent predictors of true currents and carry a higher degree of confidence,>95.5%, than underlying parameters,>53%. Unexpressed ionic currents are correctly filtered out even in the presence of mild model error. When the model is unknown, the covariance eigenvalues of parameter estimates are found to be a good gauge of model error. Our results suggest that biological information may be retrieved from data by focussing on current estimates rather than parameters.

Here, we show that ionic current predictions provide more reliable information on the biology than parameters and we propose that the covariance matrix of parameter estimates constitutes a suitable metric of model error when the biological model is unknown.We introduce Recursive Piecewise Data Assimilation (RPDA) as a novel parameter search algorithm.RPDA biases the state vector towards the true solution by re-injecting membrane voltage data in the state vector at regular intervals across the assimilation window.When the model is known, RPDA obtains the true solution with near perfect reliability and minimal error (< 0.1%), independently of the choice of initial state vectors and stimulation protocols.The uniqueness of solutions in turn establishes that observability and identifiability criteria are fulfilled.We then introduce model error by detuning a gate exponent and study the dispersion of parameters estimated by two variants (ErrM 1, ErrM 2) of the model used to generate the data.Two further model variants adding a supplementary ion channel to the original model (ErrM 3) and to its ErrM 1 variant (ErrM 4) gave the parameter and current dispersions of over-specified models.We find that model error introduces correlations between blocks of parameters defining the same ion channel type.The calculation of ionic currents integrates these correlations with the effect that mean current predictions deviate by less than 8% from their true values compared to over 100% for some parameters.Uncertainty on current estimates is less than 4.5% whereas it is 47% on parameters.Over-specifying conductance models by adding extra ion channels to them is not found to induce significant compensation between parameters.Unexpressed ion channels are correctly filtered out while those present are as-signed the correct current values, even by the mildly erroneous model ErrM 4. We finally compare the uncertainty on currents estimated with intentionally erroneous models (ErrM 1-ErrM 4) to those of a guessed CA1 model trained on CA1 hippocampal neuron data.We find that the uncertainty on the reconstructed CA1 currents is 14% compared to 4.5% on ErrM 2 currents.As the uncertainty on ionic currents increases with model error, we estimate that the error in the CA1 model is about 3 times greater than the error introduced by detuning a single gate exponent in ErrM 2. Our findings suggest that the predictive accuracy of RPDA is sufficient to detect ion channel dysfunction [51][52][53][54][55][56][57][58][59] and to infer the action of ion channel antagonists [60].

II. RECURSIVE PIECEWISE DATA ASSIMILATION
We use data assimilation [1,2] to optimize the state variables and parameters of Hodgkin-Huxley type models [61] by minimizing a least-square cost function [9,10,19].This function measures the misfit between a state variable representing the membrane voltage, x 1 , and the experimental membrane voltage, V mem , recorded at discrete times t i , i ∈ [0, N ] spanning the assimilation window of duration T : (1) The state of the neuron is described by a vector ⃗ x with L + K + 2 components.Vector components are x 1 , the membrane voltage; x 2 , . . ., x L , the gate variables of ion channels; x L+1 and x L+2 , the Tikhonov regularization variable [62][63][64] and its time derivative; and x L+2 , . . ., x L+K+2 , the model parameters.The model parameters do not depend on time.The cost function is minimized subject to both equality constraints: specified by the neuron model and the zero time derivatives of model parameters expressed as follows: and inequality constraints: specifying the range of variation of membrane voltage, gate variables, regularization term and parameters, , p max l l=L+3,...,L+K+2.
(5) J = J N a +J K +• • •−J inj is the current per unit area of the neuron membrane.This includes all ionic currents (Na, K ...) and the current injected to drive neuron oscillations, J inj .C is the membrane capacitance, τ l the recovery time of ionic gate l, and x l∞ is the steady-state value of gate variable x l .The user sets the parameter search range, [p min l , p max l ], to the widest biologically plausible range for each parameter.Data assimilation outputs the optimal parameters and the state variables at t = 0 as ⃗ x(0).
The convergence of data assimilation is compromised by badly conditioned problems.When the model is known, convergence may fail when the observability and identifiability criteria are unmet or, if the conductance model is large (L > 12), when the parameter search risks getting stuck in a local minimum of the cost function.These issues have been partially remedied by either increasing the embedding space [8, 34-36, 42, 65], improving the design of the injected current waveform, or using noise regularization [16,66].In addition, when the model is unknown, as in real neurons, the global minimum may not coincide with the true solution.Parameter searches starting from different initial state vectors may get stuck at different minima, giving multi-valued solutions.
The Recursive Piecewise Data Assimilation (RPDA) algorithm, we present in Fig. 1a, resolves these issues by biassing the parameter search towards the local minimum closest to the true solution in the ill-posed case, and towards the global minimum in the well-posed case.The parameter search is biased by re-injecting membrane voltage data into the model every M data-points.This is done by substituting x 1 with V mem at time points t 0 , t M −1 , t 2M −1 . . . in the linearized expression of the equality constraints [66]: l=1,...,L+K+2, i=0,...,N −2. ( where ∆t = T /N .At other times the state vector is propagated normally from t i to t i+2 by Eqs.6.The substitution of V mem also impacts the first and second derivatives of the objective function with respect to x 1 .Where the cost function ceases to depend on x 1 , its derivatives with respect to x 1 vanish.We explicitly set these derivatives to zero at times t 0 , t M −1 , t 2M −1 . . . .This produces a discontinuity in x 1 (t) at the beginning of each M -block (Fig. 1b) which is the tradeoff for imposing the bias to the solution.The RPDA method has similarities with the multiple shooting approach proposed by Zimmer and Sahle [3,67], which also fits data in a piecewise manner.However our RPDA method performs piecewise assimilation of blocks of data through redefined constraints, whereas the Zimmer and Sahle's approach redefines the cost function.
The RPDA algorithm recursively improves the accuracy on parameter estimates by relaxing the bias and reinjecting data over longer M -intervals, while restarting the parameter search from the previous estimate.The initial block size is M 0 = 2.The estimated state vector is then used as the new initial state in the next parameter search when M is incremented from M 0 to M 0 + 2. The search terminates when M becomes greater than N (N = 50, 001).The iterations ultimately restore the continuity of x 1 (t) across the assimilation window.When the initial bias is too strong the parameter search may fail.The algorithm then restarts the parameter search from the next larger block size (M 0 = 4).This process is depicted in Fig. 1c where the parameter search is restarted twice after convergence failed starting from M 0 = 2 then 4 to succeed with M 0 = 6.
Previous work [16,66] has explored the use of additive noise as a regularization method, showing that the stochastic perturbation of the fitting landscape by noise can disrupt local minima and improve convergence to the global minimum.In the RPDA method, each segment of length M similarly perturbs the fitting landscape.RPDA has the further advantage that the solution will remain a good solution for any value of M thanks to the bias from re-injected data.

III. RESULTS
We then examined the accuracy of RPDA in three cases of model error in the assimilating model: (i) the model is known and is correct; (ii) the model is erroneous but error is known; (iii) the model is erroneous and the error is unknown.In the first case RPDA converges 100% of the time towards to the true solution.In the second case, we assimilate the same data with 4 variants of the original model incorporating either an erroneous gate exponent, redundant ion channels, or a combination of both.In the third case, we use a guessed conductance model to assimilate membrane voltage recordings of hippocampal neurons.

A. Well-posed problems
Our reference model is a single compartment model of the rostral ventrolateral medulla (RVLM).This model is exemplar of the models used to predict neuron oscillations [9,10].It comprises five common types of ion channels: transient sodium (NaT), delayed-rectifier potas-sium (K), low threshold calcium (CaT), hyperpolarized cyclic nucleotide (HCN), and leak (Appendix A).The model has L = 7 state variables and K = 40 adjustable parameters whose reference values, p k (k = l − (L + 2)), are listed in Table .I.The RVLM model configured with these parameters was used to forward-integrate the current protocol I inj (t) in Fig. 2 and generate voltage oscillations V mem (t) over 1000ms.The sequence of depolarising and hyperpolarising steps of varying amplitudes and durations was designed to elicit a response from all ionic currents in the neuron.

Uniqueness of solutions
The convergence of RPDA was probed by initializing the state vector ⃗ x at 28 different locations: 18 initial parameter values were chosen at random in the parameter search range, the other 10 had parameters set every 1/10th of the search interval.In this way, the initial state vectors mapped the entire parameter space allowing the parameter search to approach the solution from different directions.All assimilation runs used the same 200ms long data set labelled w 1 in Fig. 2.
RPDA successfully converged to the true solution from all 28 initial conditions (100% convergence rate).23 of these assimilation runs converged from M 0 = 2. 3 more required one restart from M 0 = 4.The last 2 required two consecutive restarts from M 0 = 4 and 6, as shown in Fig. 1c.The mean block size at which parameter search arrived within 0.1% of true values was ⟨M ⟩ = 22.This shows convergence is rapid once the algorithm has selected a suitable starting block size M 0 .
The mean values of parameter estimates, µ k , and their standard deviations, σ k , were calculated over the 28 sets (Table.I).A majority of µ k (34/40) are within < 0.1% of the true parameter values.Outliers such as t 3 , t 5 , t 7 , ϵ 7 relate to gate recovery times.They nevertheless remain within 1% of the true values.In general the ultra narrow dispersions (σ k /µ k < 0.01%) confirm that parameters are fully constrained and the system is observable.
With a 100% convergence rate, the RPDA method improves over DA and DA with noise regularization which respectively achieve 67% and 94% convergence rates [66].We next demonstrate that the current protocol fulfills the identifiability criterion by sampling different data sets.

Identifiability
Identifiability was investigated by probing the dispersion of parameters estimated from 11 different sections of the stimulation protocol.These are the time windows labelled w 1 , . . ., w 11 in Fig. 2. Each window is 200ms long (10,001 data points) starting at 0ms, 40ms, 80ms, . . ., 400ms.The vector components of the initial state vector were set at the midpoint of the search range in all assimilation runs.k from true values of less than 0.1%.Unsurprisingly, parameters related to gate kinetics show the largest standard deviations.Yet the deviations of their mean values remain small suggesting parameters retain good predictive power.For example the standard deviations on t 6 and t 5 are 8.1% and 17% respectively whereas the deviations of mean predicted values from true values are 2.2% and 4.5%.Closer examination of the 11 parameter sets shows that the greater deviations occur in time windows which happen to include fewer action potentials, such as the 180-300ms interval in Fig. 2. Excluding these windows, the standard deviations fall in the normal range reported in the previous section.This is a reminder of the importance of including a minimum of 5-6 action potentials in the assimilation widow for the identifiability condition to be fulfilled.

Size of data set
We have carried out assimilations using RPDA over larger windows of 20,001, 30,001, 40,001, and 49,999 points and found that the well-posed model converges in all of these cases.The accuracy on parameter estimates is very similar to the accuracy reported in Table .I.This suggests that larger windows do not produce further gains in accuracy in well-posed problems once they include a minimum of ≈ 5 action potentials.

B. Ill-posed problems
We now study the dispersion of parameters inferred by four variants of the RVLM model which keep parametrization the same.The ErrM 1 variant had an erroneous gate exponent in the equation of the sodium current J N aT (Eq.8)where the gate probability x  II from ill-posed problems.These windows are also 200ms long but offset by 2ms so that all encompass the same action potential labelled (*) at 88ms.over-specified models was useful to determine the degree to which RPDA was able to filter out an unexpressed ion channel.A hypothesis which we will be testing here is whether RPDA assigns a finite conductance to the Achannel to compensate for model error in the Na gate probability, or whether RPDA succeeds in disentangling the contributions of each ion channel in spite of model error.
For each erroneous model, we estimated 41 sets of parameters from 41 assimilation windows, ω 1 , . . ., ω 41 in Fig. 2. We then calculated the mean values and standard deviations of the ErrM 1 and ErrM 2 parameters which we compared to the true parameters values of the RVLM model (Table II).The covariance matrices of the RVLM, ErrM 1 and ErrM 2 parameters are plotted in Fig. 3.These show the emergence of parameter correlations as soon as model error is switched on in ErrM 1 and ErrM 2. We then completed the ErrM 1...ErrM 4 models with the parameters estimated by RPDA and forward-integrated these models to predict the Na, K, HCN and A-type current waveforms (Fig. 4).These are compared to the true waveforms of the original RVLM model (dashed lines) to determine the error in the currents predicted by erroneous models.The degree of confidence on predictions was determined by computing the min and max current waveforms (Fig. 4) and the standard  The p k are the true parameters used to construct the data being assimilated (Fig. 2).The µ k and σ k are the mean values and standard deviations of parameters estimated from different initial guesses of the state vector and assimilation window w1.The µ ′ k and σ ′ k are the mean values and standard deviations of parameters estimated from assimilation widows w1, . . ., w11 starting from the same state vector.Deviations from mean greater than 1% are highlighted in bold.We set k = l − L − 2.
deviations on the integral charge under these curves (Table III).The Na, K and A-type current waveforms were reconstructed at the site of one and the same action potential (labelled * in Fig. 2) which was chosen for being common to all 41 assimilation windows used to generate the statistical sample of parameters.least well constrained parameters such as t 5 ,ϵ 5 or t 6 , ϵ 6 are the most sensitive to model error.The (t 6 ,ϵ 6 ) pair is an extreme example where standard deviations jump from (0.5%,0.05%) (RVLM), to (35%,32%) (ErrM 1), and (47%,42%) ErrM 2.
Parameter correlations are also evident in the finite off-diagonal components of the ErrM 1 and ErrM 2 co-variance matrices (Fig. 3) compared to the RVLM covariance matrix which has none -except for t 5 (parameter 28).The eigenvalues of the RVLM covariance matrix (Fig. 3a) are vanishingly small.In contrast, the ErrM 1 and ErrM 2 covariances have 6 non-zero eigenvalues.These indicate the lengths of the principal semiaxes along which parameter compensation occurs.The largest correlations are within the parameters of the CaT and HCN currents (k > 25) and they increase with increasing model error (Fig. 3b-d).An important point to note is the clustering of off-diagonal terms in blocks linking the parameters of individual ionic currents.This raises the prospect that parameter correlations may compensate each other in the calculation of ionic currents to make more accurate predictions than by relying on parameters alone.
Fig. 4 shows the range of current waveforms reconstructed from our 41 sets of 40 parameters.For clarity we only plot the min and max waveforms of this range.The min and max current waveforms predicted by the RVLM model are virtually identical to the original waveforms (Fig. 4a).Quite remarkably, the ErrM 1 model still predicts nearly identical min and max waveforms (Fig. 4b) despite the order of magnitude larger uncertainty in the underlying parameters (Table II).Only in the ErrM 2 model does the dispersion in predicted currents become noticeable (Fig. 4c).Note also the excellent agreement between the currents predicted by ErrM 1 and ErrM 2 and the true RVLM current waveform (dashed lines).These results demonstrate that approximate models can still predict the true currents with a high degree of confidence despite large uncertainty in underlying parameters.
We then calculated the total Na and K charge transferred under the predicted current waveforms.The mean values Q and standard deviations Σ are given in Ta- ble III.The standard deviations Σ 1 / Q1 (ErrM 1) and Σ 2 / Q2 (ErrM 2) remain under 5% for all currents.This is considerably lower than the uncertainty on underlying parameters which can be as high as 50% (Table III).If we examine how close the average estimates are to true values we find that the Q1 is within 5% (Na), 4% (K) of the original model prediction and Q2 is within 9% (Na), 7% (K).The ionic charges predicted by erroneous models ErrM 1 and ErrM 2 are therefore an excellent predictor of the true ionic charges.In our example, decreasing the Na gate exponent effectively decreases the slope of the activation curve.RPDA compensates for the subsequent widening of the activation curve by reducing the width parameter δV 2 from 10 (RVLM) to 8.1 (ErrM 1) and 5.2 (ErrM 2) in Table II.The large deviation (50%) of such parameter from true value is a second reason why ionic currents constitute more stable predictors of true values when the model is unknown.
2. Over-specified models: ErrM 3 and ErrM 4 It is sometimes argued that model over-specification causes multi-valued solutions.
We investigate this hypothesis by constructing two over-specified models, adding a super-numerary A-type current to the RVLM model (ErrM 3), and to the ErrM 1 model (ErrM 4).The A-type current density is J A = g A x 14 x 15 (x 1 − E K ), where x 14 and x 15 are the activation and inactivation variables respectively, and E K is the potassium reversal potential.
We find that the ErrM 3 problem has a single-valued solution that yields the true RVLM parameters.In the process, RPDA correctly filters out the A-type current by assigning a negligible value to g A .The standard deviations of the ErrM 3 parameters are also small (< 0.01%), and comparable to those of the RVLM model (Table II).The narrow dispersion of parameter estimates explains the narrow dispersion of the predicted current waveforms and their similarity to the true current waveforms (dashed lines, Fig. 4d).The mean ion volume discharged per action potential, Q3 , is within 1.5% (Na) and 3% (K) of true values hence is also an excellent predictor of true ion discharge.The standard deviations, Σ 3 , are similar to those of the RVLM model (Table III).These results show that large models still obtain the correct parameters and currents despite over-specification.Now turning to the ErrM 4 model, one hypothesis is that the A-channel parameters might compensate for the erroneous Na gate exponent.Examination of parameter distributions show nothing of the sort occurs.RPDA assigns a vanishing conductance, g A , to the supernumerary 'A-channel'.This can be seen in Fig. 4e where the A-current is zero within a standard deviation.The predicted Na and K current waveforms remains in excellent agreement with true waveforms (dashed lines) and carry a high degree of confidence in spite of model error.Turning to the mean charge transferred per action potential, Q4 , we find ErrM 4 predicts the true ionic discharge within 3.5% (Na) and 3% (K) of true values.The volume of A-type charge transfer, Q4 =82nC.cm−2 remains tiny compared to 1166nC.cm−2 for Na and 1131nC.cm−2 for K indicating the A-channel does not compensate for gate exponent error.The standard deviations on discharged ions, Σ 4 , are small and comparable to those of the RVLM model (Table III).
These remarkable results indicate that RPDA successfully disentangles the contributions of all ion current types even with approximate models that may include redundant ion channels.

Data error: assimilating noisy RVLM data
For completeness we examine the case where the V mem time series in Fig. 2a is corrupted by additive Gaussian noise.The noise r.m.s.amplitude, 0.1 mV, is comparable to the noise level in patch-clamp recordings.Fig. 4f shows the Na and K currents predicted by the RVLM model in this case.The current waveforms have a very narrow dispersion and are virtually identical to the true waveforms (dashed lines).Hence the impact of data error on predicted quantities will often be small compared to model error.TABLE IV.Na and K ion discharge per action potential of CA1 neuron Mean ion volumes discharged, Q, and standard deviation, Σ, at action potential (*) in Fig. 5.

Unknown model: hippocampal CA1 neuron
We complete our study of model error by investigating the distribution of parameters and ionic currents estimated using a guessed neuron model.Model error is now unknown which means that the mean parameter and current estimates can no longer be compared to the biological values which we are seeking.We used the model to assimilate the current-clamp recordings of a hippocampal CA1 neuron (Fig. 5a).We refer the reader to Abu-Hassan et al. [68] for details on the experimental protocol.The guessed CA1 model had 9 ion channels and 70 parameters (see Appendix B).It includes the NaT, NaP, K, A-type, Ca, BK, SK, HCN and leak ion channels believed to be present in the CA1 soma [69][70][71].We obtained 71 sets of 70 parameters from sliding time windows.Each window was 200ms long.Consecutive windows were offset by 2ms, same as in Fig. 2. The CA1 models completed with any of the 71 sets of parameters successfully predicted the observed membrane voltage (red trace, Fig. 5a).The good predictive power of the completed models therefore suggests that the error in the model equations is mild.
From these 71 sets of parameters, we reconstructed the ionic current waveforms at the site of the action potential labelled (*) in Fig. 5a.The predicted sodium currents (NaP, NaT) and potassium currents (K, A, SK, BK) are combined into Na and K waveforms in Fig. 5b.While the waveform shapes are similar to those of our earlier models (Fig. 4) the predictions cover a wider range and the potassium current has a slower rate of decay (Fig. 5b).The standard deviations on the ionic charge transferred per action potential are 14.4% (Na) and 14.6% (K) (Table IV).These are 7 times larger than in the ErrM 2 model suggesting that the discrepancy between the guessed CA1 model and the actual CA1 neuron is approximately 7 times greater than the discrepancy between the ErrM 2 model and the RVLM model.Another measure of model error are the eigenvalues of the parameter covariance matrix (Fig. 5c).Although the CA1, ErrM 1 and ErrM 2 eigenvalues fall in the same range (Figs.3b,c;5c) the CA1 distribution has a fatter tail of eigenvalues which suggests the CA1 model has more parameters with small functional overlap.

IV. DISCUSSION
In well-posed problems, the RPDA method achieves a remarkable 100% convergence rate from 28 different initial conditions (Table I) and with 41 current waveforms (Table II).Improvement over the 64% convergence rate reported by Taylor et al. [66] is achieved by re-injecting data in piecewise intervals and by recursively restarting parameter search from a larger piecewise interval if convergence fails.The well-posed case was used to validate the fulfilment of the identifiability and observability criteria when the true parameter solution was recovered starting from different initial conditions and current waveforms.With the caveat that a solution to such NP-hard inference problem is not known to exist in general, our empirical simulations suggest that in the special case of neuron-based conductance models the RPDA method obtains a solution for systems of up to at least 14 differential equations which describe most neuron types.The RVLM model is exemplar in this regard because any ionic current which may be added will have a similar mathematical form including sigmoidal activation and first order gate dynamics.
By intentionally introducing model error in well-posed data assimilation (ErrM 1, ErrM 2), multi-valued solutions were found to occur due to correlations between parameters (Fig. 3).Averaging the parameters estimated from sliding assimilation windows gave mean parameter estimates within 15% of true values except for gate activation times which deviated by a factor of two or more and had large uncertainty (up to 46%).In order to achieve a greater consistency in predictions, we calculated the ionic currents which integrate parameter correlations (Fig. 4).The mean ionic charge transferred per action potential deviated by less than 5% (ErrM 1) and 8% (ErrM 2) from true values and the coefficients of variation were very low at 0.26% (ErrM 1) and 4.5% (ErrM 2).It is therefore anticipated that ionic currents could be predicted with sufficient accuracy to resolve changes in individual ionic currents induced by inhibitory drugs or ion channel dysfunction.
Overspecifying conductance models (ErrM 3, ErrM 4) is no impediment to recovering the true parameters.The mean current waveforms estimated by ErrM 3 are identical to those of the RVLM model and their coefficient of variation are identical (< 1.5%).The ErrM 3 waveforms remain close to the true waveforms (dashed lines, Fig. 4) and their coefficient of variation is also small (< 1.2%).The charge transferred under the ErrM 4 waveforms was within 4% of true values and the uncertainty on predictions was less than 1.3% despite model error and overspecification.This demonstrates the ability of the RPDA method to disentangle the correct contribution of each ion channel to the membrane voltage without the overspecified current compensating for model error.The implication is that a universal conductance model could infer accurate information on ion channels without requiring any prior assumption on which ion channel might be expressed.This also suggests that reductionist models that minimize the number of parameters may not be the best route to improve accuracy on parameter estimates.Instead, focus should be on mitigating model error.
We found that standard deviations on parameters (Tables II) and current estimates (Tables III, IV) increase with increasing model error.Therefore standard deviations provide a metric to quantify how close a guessed model is to the unknown biological model.
A possible strategy for correcting model error has been proposed by Ye et al. [72,73].This consists in adding a model error term to the cost function and balancing data error and model error in such a way as to eliminate the bias of model error on the parameter solution.When model error is properly weighted, parameters are assigned true values without having to compensate for model error.The challenge is in the determination of the covariance matrix representing this weight [72].We have attempted to include a model error term in Eq.3 and weighted it with a single hyperparameter.However this approach has proved unpractical due to the difficulty of finding a finite hyper-parameter value that resolves all constraint violations.The optimal value of such hyperparameter was found to drift to zero or infinity during parameter search.Further progress will require evaluating the covariance matrix weighting model error.

V. CONCLUSIONS
We have introduced Recursive Piecewise Data Assimilation as an novel method for optimizing neuron-based conductance models.RPDA improves over earlier parameter estimation methods by biasing the parameter search with data and by iteratively improving the solution as the bias is gradually released.When the model is known, RPDA achieves a 100% convergence rate towards the true solution for a wide range of initial conditions and training data sets.When the model is unknown, model error introduces correlations between parameters estimates.The largest correlations occur between parameters that define each ionic current.By reconstructing the ionic currents, parameter correlations cancel out and the current waveforms are found to approximate very well the true current waveforms.The ionic charge transferred under these curves is also predicted with a high degree of confidence of 95.5% when model error affects a single ion channel and 85% when the true biological model is guessed (CA1 neuron).The increasing uncertainty of predictions with increasing model error suggests that the covariance matrix of parameters is a good metric of model error.We also found that model over-specification is not responsible for parameter sloppiness and that RPDA correctly disentangles all ion channel contributions to the membrane voltage data even when the model has a wrong gate exponent.Our work shows that combining variational inference with statistical analysis offers good prospects for extracting sensible biological information from data.
The nominal RVLM parameters are listed in the "True" value column of Table I.Eqs.7-10 were configured with these parameters to generate the membrane voltage data V mem in Fig. 2. The injected current waveform I inj (t) had 50,000 points with a 0.02 ms time step.The RVLM model was forward-integrated with the odeint function of the scipy package in Python 3. The resulting V mem data were then assimilated by the RVLM system and its Errm1-ErrM 4 variants.

VII. APPENDIX B: HIPPOCAMPAL NEURON MODEL
This hippocampal neuron model includes the ion channels in the soma of CA1 neurons [69-71, 74, 75]: persistent sodium current (NaP), a A-type potassium current (A) and calcium activated potassium currents (SK and BK) in addition to the NaT, K, Calcium, HCN and leak currents of the previous RVLM model.The model now has 9 ionic currents, 14 state variables (L = 14) and 69 adjustable parameters (K = 69).The membrane voltage dynamics is given by: The ionic current densities are: J N aT = g N aT x 2 x 3 3 (x 1 − E N a ), J N aP = g N aP x 4 (x 1 − E N a ), J K = g K x 4  5 (x 1 − E K ), J A = g A x 6 x 7 (x 1 − E K ), J Ca = g Ca x 2 8 x 9 (x 1 − E Ca ), J BK = g BK x 2 10 x 11 (x 1 − E K ), J SK = g SK x 12 (x 1 − E K ), J HCN = g H x 13 (x 1 − E HCN ), J L = g L (x 1 − E L ). (13) The BK current depends on both the membrane voltage x 1 and the internal calcium concentration [Ca 2+ ] i ≡ x 14 whereas the SK current depends on [Ca 2+ ] i only.
The BK current has fast activation which we take to be instantaneous in line with Warman et al. [69]: where the calcium equilibrium across the membrane is given by the following rate equation: driven by the calcium ion current J Ca .w is the thickness of the surface area across which Ca 2+ fluxes are calculated (w = 1µm).We modelled the dynamics of calcium gate variables x 8 and x 9 with Eqs.9 and 10 with the appropriate calcium gate parameters.The inactivation dynamics of the BK current has a long time constant τ 11 .It also follows rate equation: The activation dynamics of the SK current is fast (x 12 ≡ x 12,∞ ), and similar to the BK channel, the activation of the SK channel is described by: Gate variables x 2 (NaP channel) and x 6 , x 7 (A channel) follow the first order response in Eqs.9 and 10 with the appropriate NaP and A channel parameters.

FIG. 1 .
FIG. 1. RPDA algorithm (a) Algorithm flowchart; (b) The RPDA algorithm synchronizes the membrane voltage variable x1 to the time series data, Vmem, in blocks of M discrete time points.The algorithm re-injects Vmem in the model at the beginning of each block (red circle); (c) Example of recursive convergence starting from M0 = 6 and increasing the block size to M = N .This follows two "false starts" at M0 = 2 and M0 = 4.

3 2 x 3 FIG. 2 .
FIG. 2. Neuron oscillations and assimilation rangesMembrane voltage (black trace) generated by forward integration of the injected current protocol (blue trace) with the original RVLM model.w1, . . ., w11 are the 11 assimilation windows used to generate the data in Table I from well-posed problems.Each window is 200ms wide (yellow bands) and offset from the next by 40ms.Top panel : ω1, . . ., ω41 are the 41 assimilation windows used the generate the data in Table II from ill-posed problems.These windows are also 200ms long but offset by 2ms so that all encompass the same action potential labelled (*) at 88ms.

FIG. 3 .
FIG. 3. Covariance maps of parameters inferred with the RVLM, ErrM1 and ErrM2 models (a) Eigenvalues of the covariance matrix of parameters estimated with erroneous models ErrM1 (blue trace) and ErrM2 (red trace).The eigenvalues of the well-posed RVLM model (black trace) are shown for reference.Covariance maps of the (b) RVLM, (c) ErrM1, and (d) ErrM2 parameters.

FIG. 4 .%
FIG. 4. A comparison of ionic currents predicted by erroneous and exact models (a) Current waveforms predicted by the RV LM model.The min and max in the Na and K traces give the range of variation of currents reconstructed from the 41 assimilation windows.For reference, the true Na and K current waveforms predicted by the original RV LM model are shown as the black dashed lines in (b-f).Current waveforms predicted by erroneous models: (b) ErrM 1, (c) ErrM 2; an over-specified RV LM model: (d) ErrM 3; an over-specified erroneous model: (e) ErrM 4; and from noisy data, (f) N oisy.

FIG. 5 .
FIG. 5. Uncertainty on currents and parameters estimated from a hippocampal CA1 neuron (a) Membrane voltage oscillations (black trace) driven by injecting a calibrated current waveform (blue trace) into a hippocampal neuron (Wistar rat, CA1 neuron).Voltage oscillations predicted by the hippocampal neuron model completed with one set of parameters.The (*) labels the action potential whose ionic current waveforms we predict in (b).(b) Mini and maxi current waveforms (Na, K) predicted by 71 CA1 models completed with the parameters from 71 assimilation windows.(c) The eigenvalues of the 70 × 70 covariance matrix of estimated parameters.

TABLE I .
Dispersion of parameters estimated from different initial guesses and different assimilation windows

1 .
Erroneous model variants: ErrM 1 and ErrM 2Detuning the gate exponent from x 3 2 x 3 (RVLM) to x 2 2 x 3 (ErrM 1) and x 2 x 3 (ErrM 2) introduces correlations between parameters.Their standard deviations jump by an order of magnitude from RVLM to ErrM 1 (TableII) and by a smaller increment from ErrM 1 to ErrM 2. Gate recovery times previously identified as the

TABLE II .
Dispersion of parameters inferred from erroneous models ErrM1 and ErrM2The µ k and σ k are the mean values and standard deviations of parameters estimated from windows ω1, . . ., ω41 using the RVLM model (reference estimates).The µ 1k , σ 1k and µ 2k , σ 2k are the means and standard deviations of parameters estimated with the erroneous models, ErrM1 and ErrM2 respectively.Figures in bold (bold+italics) single out deviations from mean greater than 1% (10% respectively).

TABLE III .
Na, K and HCN ion discharge per action potential QNa, QK , and QHCN are the ionic volumes discharged during action potential (*).Q1 and Σ1 are the mean values and standard deviations derived from a statistical sample of 41 ion current waveforms.These waveforms are predicted by 41 ErrM 1 models completed with the 41 parameter sets from assimilation windows ω1, . . ., ω41.The same was done for RVLM, ErrM 2, ErrM 3 and ErrrM 4. QN and ΣN are mean charge and standard deviations obtained by assimilating noisy data with the RVLM model.