Algorithmic Parameter Estimation and Uncertainty Quantification for Hodgkin-Huxley Neuron Models

Experimental data-based parameter search for Hodgkin–Huxley-style (HH) neuron models is a major challenge for neuroscientists and neuroengineers. Current search strategies are often computationally expensive, are slow to converge, have difficulty handling nonlinearities or multimodalities in the objective function, or require good initial parameter guesses. Most important, many existing approaches lack quantification of uncertainties in parameter estimates even though such uncertainties are of immense biological significance. We propose a novel method for parameter inference and uncertainty quantification in a Bayesian framework using the Markov chain Monte Carlo (MCMC) approach. This approach incorporates prior knowledge about model parameters (as probability distributions) and aims to map the prior to a posterior distribution of parameters informed by both the model and the data. Furthermore, using the adaptive parallel tempering strategy for MCMC, we tackle the highly nonlinear, noisy, and multimodal loss function, which depends on the HH neuron model. We tested the robustness of our approach using the voltage trace data generated from a 9-parameter HH model using five levels of injected currents (0.0, 0.1, 0.2, 0.3, and 0.4 nA). Each test consisted of running the ground truth with its respective currents to estimate the model parameters. To simulate the condition for fitting a frequency-current (F-I) curve, we also introduced an aggregate objective that runs MCMC against all five levels simultaneously. We found that MCMC was able to produce many solutions with acceptable loss values (e.g., for 0.0 nA, 889 solutions were within 0.5% of the best solution and 1,595 solutions within 1% of the best solution). Thus, an adaptive parallel tempering MCMC search provides a “landscape” of the possible parameter sets with acceptable loss values in a tractable manner. Our approach is able to obtain an intelligently sampled global view of the solution distributions within a search range in a single computation. Additionally, the advantage of uncertainty quantification allows for exploration of further solution spaces, which can serve to better inform future experiments.

Is manual tuning of model parameters still practicable? 69 Fitting experimental data to neuron models has been a major challenge for neuroscientists. In fact, 70 the parameters of the first biologically realistic quantitative description of the neuronal membrane 71 were hand fitted by Hodgkin and Huxley (1952). This approach has remained popular; for exam-72 ple, the peak conductances of ion channels in a model of an elemental leech heartbeat oscillator 73 were hand tuned by Nadim et al. (1995). Likewise, both the single compartment and network pa-  Computational inference for neural models 85 The body of research on parameter estimation for models of neural dynamics for single cells or 86 circuits spans across various scientific communities, approaches, and neuron models. A gap ex- 87 ists, however, between researchers with rich physiological knowledge, on the one hand, and re-  ., 2011). The latter include simulated annealing, differential evolution, and genetic algorithms.

113
These approaches have the disadvantage of slowly converging to an optimal set of parameters 114 and thus being computationally expensive. Using gradient-based optimization to accelerate con-     Furthermore, the key advantage of parallel tempering MCMC is its ability to expose multimodality    In this study we address the issue of multiple realizability of models in neuroscience. Specifically, 173 we target the computational inference for HH-based neuron models with a state-of-the-art MCMC 174 algorithm. Our goal is to estimate the parameters in such models and quantify their uncertainties. 175 We design the parameter estimation problem in a Bayesian framework: (i) parameter sets in a 176 model are treated as probability distributions rather than one single optimal set; and (ii) the distri-177 butions of parameter sets can exhibit dependencies between individual values of the set, rather 178 than assuming each parameter to be independent from the others. This Bayesian framework for 179 parameter estimation does not assume that a single independent optimum can be reached. The 180 likelihood is given implicitly in the form of a loss function between observational data and model 181 output, hence requiring the numerical solution of the model for each new set of parameters. As 182 a solution algorithm we utilize an extension MCMC, the parallel tempering MCMC method, which 183 is critical for recovering distributions of parameter sets that are multimodal (i.e., multiple sets of 184 parameters of a neuron model are recovered for a single set of observational data). Further, we 185 visualize solution maps of the multimodal posterior that enable us to quantify which parameters 186 can reproduce the observational data more closely.

188
The results of our study show that our design of the inverse problem combined with parallel tem-189 pering MCMC for solving the inverse problem is able to successfully overcome the inference chal- solution points (i.e., a single set of "optimal" parameters) obtained from optimization algorithms. 200 We illustrate the results of the parallel tempering MCMC method in a simple inverse problem setup, 201 which nevertheless exhibits the difficulties of more complex inference setups in neural dynamics. 202 We use a relatively simple and well-known HH model (

208
The first goal of the inverse problem is to recover the true sodium and potassium conductances.

209
To demonstrate the challenges of the problem, we visualize the loss function that would have to 210 be minimized to find the optimum in Figure 1A, where the colors indicate the (positive) loss value.

211
The true parameter set is located in a "valley" in the loss "landscape," whereas the loss is large 212 for parameters that produce large discrepancies between model output and data such as shown 213 in graph S1 corresponding to point S1 in Figure 1A. However, multiple sets of parameters exist 214 that are local minima in this landscape, and these multiple local minima pose major challenges 215 because optimization algorithms may present one of them as the "optimal" solution. These local 216 minima will not give a sufficiently good fit of model outputs vs. data, as is illustrated in graph S2 of 217 Figure 1 that corresponds to point S2 in the landscape ( Figure 1A). Such local minima are known 218 to be problematic for numerical optimization algorithms (Nocedal and Wright, 2006). 219 The second goal of the inverse problem is to, additionally, quantify the uncertainty with respect 220 to these parameters and understand the sensitivities of the parameters. For instance, in Figure 1, 221 the graphs S3 and S4 illustrate how different parameters produce voltage traces similar to the data, 222 while at the same time this uncertainty in the parameters is visible in the landscape ( Figure 1A) as 223 a valley containing points S3 and S4. In order to recover these uncertainties, the formulation of 224 the inverse problem in a Bayesian framework becomes critical.

225
The numerical solution of the Bayesian inverse problem with parallel tempering MCMC success-226 fully provides a probability density spanned by the two-dimensional parameter space. The density 227 is large, where the loss function in Figure 1A has its major valley. We show the progress of the 228 MCMC algorithm in Figure 1B as the number of collected samples increases from 2,000 to 8,000.  Figure 1. Loss function corresponding to the inverse problem with a simple HH model and how MCMC is able to successfully recover optimal parameters and quantify uncertainties and parameter trade-offs. A: Loss function landscape (colors) over the parameters, sodium ( ) and ( ) potassium conductances. The "true" parameter set represents a global minimizer of the loss; the valley of the loss along points S3 and S4 (in dark blue color of loss) translates to parameter uncertainties (or trade-offs). Graphs S1-S4 depict voltage traces at corresponding points in A. B: Densities of MCMC samples, where the dark purple color of highest density overlaps with the true parameter set. As the number of MCMC samples increases, the method recovers the valley of the loss landscape A and hence quantifies uncertainties in the parameters.
As the sample count (i.e., iteration count of MCMC) grows, the algorithm generates a longer tail These results serve as a proof of concept that the parallel tempering MCMC algorithm can suc-236 cessfully tackle multimodal losses. Next we transition to a more complex model for neural dynam-237 ics, which will be the main result of this paper.

311
In the next section we formalize these intersections using the Wasserstein distance metric.   The approximate Wasserstein distance provides a quantitative metric of the total cost required to 318 transform one probability distribution to another probability distribution (Givens and Shortt, 1984). 319 The results are displayed in Table 2    Single-Solution Approach: As detailed in the Results, a parameter set exists that will produce the 344 lowest loss value for each of the individual and aggregate constraints (see Table 1). Many sets of 345 solutions also will provide extremely low loss values. For instance, we found, for the solution space 346 for the injected current of 0.0 nA, a total of 889 solutions within 0.5% of the best solution and a 347 total of 1,595 solutions with 1% of the best solution (see Table 3). 348 Even though the loss values are small for many parameter sets, the voltage traces produced by 349 these sets are different. To illustrate these differences, in Figure 6 we plotted the first five voltage  Table 3. Number of solutions given the injected current as a function of percentage difference from the best solution. For instance, for a 0.0 nA injected current, 889 solutions are 0.5% away from the best 0.0 nA solution. MCMC is able to find these multiple acceptable solutions. Figure 6 visualizes an example set of voltage traces generated from these solution parameter sets.
Multimodal Approach: Many more nearly-optimally fitting solutions exist, however. One can which is out of the scope of this work.

392
Physiological meaning 393 Extracting physiological information from the solution maps is the ultimate goal of this exercise.

394
For instance, one can quickly rank the importance of a parameter by looking at the 2D marginals.

395
Looking at the bottom row of Figure 2, we observe that is highly correlated to while  Alonso and Marder showed that multiple solutions exist for solving for the same neural out-401 come. MCMC complements this finding by showing that a much larger solution set can be found. 402 We believe that this multimodal solution set is the norm rather than the exception. As models in- stiff biological ODE models. LSODA has been shown to perform better in biological systems than 426 many modern ODE solvers (Städter et al., 2021). The voltage trace generated by the ground truth 427 parameters is shown in Figure 5 in orange color.

428
Loss function for measuring fit between data and model outputs 429 The main loss function used in the present study involving the AM model is the loss in Equation be dominant, associated with , the burst frequency mismatch, and with , the duty cycle mismatch.

440
Our objective is to test the multimodality of the inverse problem with parallel tempering MCMC 441 methods.

442
Loss function for individual and aggregate constraint 443 We use the loss function (detailed in Loss function for measuring fit between data and model out- In this work we choose the Bayesian formulation (also called Bayesian framework), which is cen-456 tered around Bayes' theorem: where denotes the parameters of the ODE model and is the data. On the left-hand side of (1),

458
( | ) is the unknown posterior (i.e., conditional probability for the model parameters given