Continuation methods as a tool for parameter inference in electrophysiology modeling

TL;DR

This work tackles the computational bottleneck of calibrating cardiac electrophysiology models to stable periodic data by employing continuation methods to track limit cycles as parameters vary. The authors implement a Noble-model–based framework extended with ion concentrations and a convergence-rate parameter, and demonstrate that continuation can reduce the cost of obtaining converged action potentials by approximately $70\%$ during MCMC-based parameter inference. They compare Standard, Tracking, and Continuation strategies, showing that continuation delivers significant speed-ups while preserving posterior accuracy, and discuss practical limitations such as multistability and possible extinction of limit cycles. The approach broadens the toolkit for fast, robust parameter inference in fast-slow dynamical systems and can facilitate more extensive Bayesian analyses and experimental design in electrophysiology.

Abstract

Parameterizing mathematical models of biological systems often requires fitting to stable periodic data. In cardiac electrophysiology this typically requires converging to a stable action potential through long simulations. We explore this problem through the theory of dynamical systems, bifurcation analysis and continuation methods; under which a converged action potential is a stable limit cycle. Various attempts have been made to improve the efficiency of identifying these limit cycles, with limited success. We demonstrate that continuation methods can more efficiently infer the converged action potential as proposed model parameter sets change during optimization or inference routines. In an example electrophysiology model this reduces parameter inference computation time by 70%. We also discuss theoretical considerations and limitations of continuation method use in place of time-consuming model simulations. The application of continuation methods allows more robust optimization by making extra runs from multiple starting locations computationally tractable, and facilitates the application of inference methods such as Markov Chain Monte Carlo to gain more information on the plausible parameter space.

Figures (10)

• Figure 1: Convergence to a limit cycle or 'stable action potential'. Simulations used the O'Hara-Rudy CiPA model of duttaOptimizationSilicoCardiac2017. A: Two action potentials at different stages of convergence. An unconverged AP in pink, and a more converged AP (after 30) in black. B: Phase diagram of intracellular calcium concentration against membrane voltage, showing the convergence from a perturbed state (pink) to a stable action potential (black).
• Figure 2: Benchmarking computation times for differently sized parameter perturbations. Box plots show the quartiles across repeated simulations with whiskers are the maximum and minimum. The 0th, 25th, 50th, and 75th percentiles overlap. Computation times reported in seconds for: 'ODE - Tracking' ODE simulation from previous limit cycle (Table \ref{['tab:ic']} - Converged), 'ODE - Standard' ODE simulation from standard non-converged initial conditions, 'Cont - Shooting' Shooting method continuation to follow a previously converged limit cycle (Table \ref{['tab:ic']} --- Converged). A: Small perturbation $(\theta_1=1, \theta_2=1, \theta_3=1)\rightarrow(\theta_1=1.1,\theta_2=1,\theta_3=1)$. B: Large perturbation $(\theta_1=1, \theta_2=1, \theta_3=1)\rightarrow(\theta_1=1.5,\theta_2=1.2,\theta_3=0.8)$
• Figure 3: APs for the Standard and Tracking approaches, aligned on the first upcrossing of $V=\qty{-20}{\mV}$. 'Start - Tracking' --- the converged AP for $\theta_1=1,\theta_2=1,\theta_3=1$ which is used to initiate the Tracking approach (Table \ref{['tab:ic']}). 'Start - Standard' --- the non-converged initial AP for the Standard approach (Table \ref{['tab:ic']}). 'End - Small Perturbation' --- the converged AP for the small perturbation parameters, $\theta_1=1.1,\theta_2=1,\theta_3=1$. 'End - Large Perturbation' --- the converged AP for the large perturbation parameters, $\theta_1=1.5,\theta_2=1.2,\theta_3=0.8$.
• Figure 4: MCMC estimated posterior. The Noble concentration model, Eq. \ref{['equ:nobleConc']}, is fitted to a single noisy converged AP from the data generating parameters $\theta_1=1, \theta_2=1,\theta_3=1$. The method for generating the data and calculating the likelihood are summarized in Section \ref{['sec:mcmcMethods']} with further details given in Appendix \ref{['sup:mcmc']}. The lower triangle shows scatter plots (blue: strong correlation, yellow: weak correlation). The upper triangle shows heatmaps of the bivariate densities. The diagonal contains histograms of the marginal densities. The data generating parameters are given in red. Results shown are from the Continuation approach.
• Figure 5: Examples of possible problems with using continuation methods when multiple limit cycles are possible. A: Different methods can converge to different limit cycles. If the Continuation approach initially converges to the limit cycle at $p_1$ (starting from the pink initial condition line), and then follows this limit cycle to $p_2$ (blue trajectory), it identifies a different limit cycle to the Standard approach which is run from the initial conditions at $p_2$ (red trajectory). B: The limit cycle at $p_1$ may not exist at $p_2$. The Continuation approach (blue trajectory) finds a limit cycle at $p_1$ which it attempts to follow towards $p_2$, but the limit cycle doesn't exist at $p_2$. As such, the Continuation approach instead follows the unstable branch, while this cannot happen to the Standard approach (red trajectory). LC: Limit cycle, IC: Initial condition.