Estimating Time-Varying Applied Current in the Hodgkin-Huxley Model

The classic Hodgkin-Huxley model is widely used for understanding the electrophysiological dynamics of a single neuron. While applying a constant current to the system results in a single voltage spike, it is possible to produce more interesting dynamics by applying time-varying currents, which may not be experimentally measurable. The aim of this work is to estimate time-varying applied currents of different deterministic forms given noisy voltage data. In particular, we utilize an augmented ensemble Kalman filter with parameter tracking to estimate four different deterministic applied currents, analyzing how the model dynamics change in each case. We test the efficiency of the parameter tracking algorithm in this setting by exploring the effects of changing the standard deviation of the parameter drift and the frequency of data available on the resulting time-varying applied current estimates and related uncertainty.


Introduction
The Hodgkin-Huxley model is a classical system of differential equations that is widely-used for understanding the electrophysiological dynamics of a single neuron [1]. The model is based on a simple circuit analogy, where each piece of the circuit corresponds to an electrophysiological component, representing the resistance of an electrically charged ion channel as a function of time and voltage [2,3]. While the Hodgkin-Huxley equations can be used to model the total current resulting from an applied voltage, the model can also be used to predict voltage given an externally applied current. The latter is particularly useful in experimental settings where voltage measurements are obtainable, making it possible to estimate the applied current based on observed voltage data.
In the standard model formulation, applying a constant current to the system results in a single voltage spike. However, by applying time-varying currents, it is possible to produce more interesting dynamics that include multiple voltage spikes. The work in this paper focuses on how the Hodgkin-Huxley system dynamics change when various deterministic, time-varying currents are applied. More specifically, the aim of this work is to estimate time-varying applied currents of different deterministic functional forms given some noisy observations of voltage. To tackle this inverse problem, we utilize an augmented ensemble Kalman filter (EnKF) with parameter tracking to estimate the Hodgkin-Huxley model states and time-varying applied current parameter.
Various methods have been used in the literature to estimate certain constant (or static) parameters in the Hodgkin-Huxley equations; see, e.g., [4][5][6][7]. However, the focus of this work is on estimating the time-varying applied current parameter, which we assume is unmeasurable with unknown dynamics. While many parameter estimation methods are available in the literature, the EnKF is particularly useful ! " # ℓ # %& # ' ((*) , ℓ , ' , %& Figure 1. The Hodgkin-Huxley model represented as a circuit. Here the capacitor (C M ) represents the charge storage capacity, the resistors (g , g K , and g Na ) act as the ion channels, and the batteries (V , V K , and V Na ) act as the electrochemical potentials.
for the problem at hand due to the sequential nature of the algorithm's updating scheme, which corrects the model prediction with the available data one point at a time [8,9]. If the time-varying parameter changes more slowly than the system dynamics, it is possible to track the change in the parameter over time using a random walk [10][11][12]. Further, since unknowns are treated as random variables in the Bayesian framework, there is a natural measure of uncertainty in the resulting parameter estimates, which lies in the estimated ensemble covariances of the underlying posterior probability distributions [13,14]. We analyze this problem using synthetic voltage data generated by applying four different deterministic functions as the applied current -a constant current, a step function with one long step, a step function with multiple shorter steps, and a sinusoidal function -in order to track and observe how the model dynamics change in each of these cases. We further test the efficiency of the EnKF with parameter tracking algorithm in these cases by performing numerical experiments first to establish baseline results, then to explore the effects of changing the standard deviation of the drift term in the parameter tracking algorithm as well as the frequency of data available on the resulting applied current parameter estimates.
The paper is organized as follows. Section 2 gives a brief review of the Hodgkin-Huxley model, summarizing the relevant equations. Section 3 describes the parameter estimation inverse problem and outlines the ensemble Kalman filtering algorithm, describing in particular time-varying parameter estimation using the EnKF with parameter tracking. Section 4 gives the numerical results, including the generation of synthetic data and the numerical experiments relating to estimating the time-varying applied current parameter, and Section 5 features a discussion of the results and future work.

Review of the Hodgkin-Huxley Model Equations
The Hodgkin-Huxley model provided the first quantitative description of electrical excitability in nerve cells [15], involving detailed mathematical equations to describe the voltage-dependent and time-dependent properties of the sodium and potassium conductances [1]. Each piece of the circuit shown in Figure 1 corresponds to a different electrophysiological component of the model. Capacitors represent the charge storage capacity of each gating variable; resistors represent the sodium, potassium, and leakage ion channels in the neuron; and batteries represent the electrochemical potentials that each gating variable has to let ions in and out of the charged cell. From this analogy, the Hodgkin-Huxley equation modeling total membrane current is given by where I = I(t) is the total membrane current in the axon (with a positive inward current), C M is the membrane capacity (assumed to be constant), and V = V(t) is the displacement of the membrane potential from its resting value (assumed to have a negative depolarization). Table 1 lists each model component, along with its corresponding units. For simplicity of terminology, we will refer interchangeably to V(t) as the voltage within this paper. Note that the voltage V is related to the membrane potential E via the relationship V = E − E r , where E r denotes the absolute value of the resting potential [1]. The ionic current density I ion is represented as the sum of the three currents where I Na , I K , and I model the currents relating to the sodium, potassium, and leakage channels occurring in the neuron, respectively. The form of the current for each ion channel follows from Ohm's law, where for i = Na, K, . Here g i (t, V) represents the gate for each channel generally as a function of time and voltage, and V − V i represents the difference between the overall voltage V of the system and the channel-specific voltages V i . We describe each of the three ionic currents in more detail as follows. Sodium current. The sodium current is given by where the sodium gate g Na (t, V) = m 3 hḡ Na is impacted by depolarization, which causes an increase in sodium conductance [1]. Hereḡ Na is a constant (conductance/cm 2 ), m = m(t) is the proportion of active sodium gates open (dimensionless variable which varies over time between 0 and 1), and h = h(t) is the proportion of inactive gates open (similarly dimensionless, varying between 0 and 1). The sodium voltage is given by V Na = E Na − E r , where E Na is an equilibrium potential for sodium. Table 2 lists the constant values of V Na andḡ Na .
The dynamics of the sodium gating variables m(t) and h(t) are governed by the following differential equations: where the voltage-dependent rate constants (msec -1 ) α m and α h represent the rate of flow of ions into the cell and β m and β h represent the flow out. The rate constants are modeled using the following equations, derived from Hodgkin and Huxley's experimental results [1]: Note that setting α m to its limit value of 1 at V = −25 mV avoids the discontinuity at that point. Potassium current. The potassium current is given by with potassium gate equation based on the assumption that potassium ions must have four similar particles in order to cross the membrane. Hereḡ K is a constant (conductance/cm 2 ) and n = n(t) is the proportion of potassium gates open (dimensionless, varying between 0 and 1). The potassium voltage is given by V K = E K − E r , where E K is an equilibrium potential for potassium, sensitive to the overall outside concentration of charged ions [16]. Table 2 lists the constant values of V K andḡ K . The dynamics of the gating variable n(t) are similarly modeled using the differential equation where the rate constants were also derived using experimental data [1]. Note the discontinuity in α n when V = −10 mV can be avoided by setting it equal to its limit value of 0.1 at that point. Leakage current. The leakage current is a small combined current, accounting mostly for chloride but also other ions. The leakage current is given by with constant conductanceḡ and leakage voltage V = E − E r . Here E is the potential at which the leak current is zero. The leakage voltage V is needed for any calculation for threshold, but it is unlikely to give any information about the nature of charged particles [16]. Table 2 lists the constant values of V andḡ . Model summary. In summary, the Hodgkin-Huxley model comprises the total membrane current equation (1), which depends on time, voltage, and the solutions to the transfer equations (6), (7), and (14). Note that when a constant voltage is applied, dV dt = 0 and (1) simplifies to I = I ion . In this case, equations (6), (7), and (14) can be solved independently to compute the total ionic current.
However, when voltage changes with time due to an applied current, dV dt = 0 and all four equations must be solved simultaneously. The complete system of coupled ordinary differential equations is given by Note that in this case, the applied current I = I(t) in (18) drives the system dynamics. We explore how different choices of deterministic, time-varying functions for I(t) affect the dynamics of the system in the numerical results.

Parameter Estimation and the Ensemble Kalman Filter
Given measurements of voltage, our aim is to estimate the time-varying applied current I(t) that best fits the available data. This is a parameter estimation inverse problem, where the parameter of interest is a time-varying deterministic function with assumed unknown dynamics. More specifically, we assume here that we cannot directly measure the time-varying applied current and that we do not have equations available to explain its dynamics.
The set-up for this inverse problem is similar to the standard set-up for estimating parameters in initial value problems of the form where x = x(t) denotes the model states and θ denotes the model parameters [17]. Given some discrete, noisy system measurements the inverse problem is to estimate the model states x(t) and parameters θ. Most classical approaches addressing this problem tend to focus on the case when the parameters are constants, i.e., when dθ dt = 0. In this case, however, θ = θ(t) and dθ dt is some unknown function. To estimate the time-varying applied current in this work, we use a version of the ensemble Kalman filter (EnKF) with parameter tracking [12,17]. The EnKF is an extension of the classical Kalman filter adapted to work with models that are not necessarily linear or Gaussian [8,9]. As a Bayesian statistical algorithm, the EnKF treats all unknowns as random variables with corresponding probability distributions. The filter uses a random sample to represent the current probability distribution of states and parameters, then utilizes ensemble statistics along with model predictions and observed data to update the sample at each discrete time point.
While the original EnKF was implemented for state estimation, the augmented EnKF allows for simultaneous state and parameter estimation [18]. The steps of the augmented EnKF are summarized as follows. At time j, the sample gives a discrete representation of the probability distribution, which is then updated using a two-step process to time j + 1. In the first step (the prediction step), we solve the system (22) to predict the state values at time j + 1. In this work, (22) is the Hodgkin-Huxley model given in (18)- (21). The state prediction ensemble is computed using the equation where F is the solution to (22) which are used to compute ensemble statistics. The prediction ensemble mean is computed using the formulaz and the prior covariance matrix by Note that while the parameters θ p j are not updated in the prediction step, their cross-correlation information with the predicted states is embedded in the prior covariance matrix and is used in the next step to update the posterior sample.
In the second step (the observation update), the predicted values are compared with the observed data y j+1 at time j + 1. The observation ensemble where w p j+1 ∼ N (0, D) represents the observation error, is compared to the observation model predictionŝ computed using the observation model G as in (23). In this work, G is a linear observation function measuring only the voltage in the Hodgkin-Huxley system. The combined posterior ensemble is then given by z for each p = 1, 2, . . . , N. The Kalman gain matrix K j+1 incorporates the cross-covariance of state and model predications, the forecast error covariance, and the observation noise covariance. For additional implementation details, see [17,19].
In the above formulation of the augmented EnKF, θ is assumed to be constant and is evolved artificially with time in order to obtain an estimate. When θ is time-varying, as in the case we are considering, the augmented EnKF as presented requires additional modification. If θ = θ(t) changes more slowly than the dynamics of the system, it is possible to track the changes in θ(t) by incorporating a random walk in the prediction step. To implement this, a new random variable ξ is introduced and added to current estimate of θ at each prediction step, allowing the previously fixed parameters to take a random walk of the form for each p = 1, 2, . . . , N. Parameter tracking of this type has been used in various data assimilation problems; see, e.g., [10][11][12]20]. Here we note that the choice of the standard deviation σ ξ in the drift term of the random walk is important in the accuracy and uncertainty of the resulting parameter estimate. We will explore this further in the numerical experiments.

Numerical Experiments
In this section we describe the numerical experiments performed using the augmented EnKF with parameter tracking to estimate the time-varying applied current in the Hodgkin-Huxley model. All experiments were performed using the MATLAB R programming language. In particular, we used the built-in solver ode15s to numerically solve (18)-(21) due to the potential stiffness of the system when applying different time-varying currents. We first describe how synthetic data was generated using four deterministic applied currents. We then provide the numerical results obtained using parameter tracking to estimate the applied current in each case, testing the effects of changing the standard deviation of the parameter drift in (32) and the amount of available data.

Synthetic Data Generation
To generate synthetic data, four different deterministic functions for the applied current I(t) were run through the Hodgkin-Huxley equations (18)- (21). The applied currents considered were: (a) a constant current, where I(t) = 2 mA/cm 2 (33) (b) a step function with one long step, such that (c) a "pulsing" step function with multiple short steps, such that  Figure 2 shows each data sets along with the corresponding applied current. Note that the gating variables n(t), m(t), and h(t) are unobserved states.

Estimating Time-Varying Applied Current via Parameter Tracking
To establish baseline results for parameter tracking in each of the four data sets described in Section 4.1, the augmented EnKF was employed using N = 100 ensemble members with the standard deviation of the parameter random walk in (32) set to σ ξ = 1. The results in Figures 3-6 show the parameter tracking estimates of I(t), along with the time series estimates of the Hodgkin-Huxley model states, for each of the four cases. Note that in each case, the filter mean is able to track the underlying true applied current, along with the unmeasurable system states, with uncertainty bounds represented by the ±2 estimated standard deviation curves.
As noted in Section 3, a carefully chosen standard deviation σ ξ for the random walk in (32) is crucial in maintaining the accuracy of the time-varying parameter estimate and avoiding filter divergence; see [19] and references therein. To demonstrate the sensitivity of the algorithm to this choice, we tested the effects of changing σ ξ by letting σ ξ = 10, 2, 1, 0.5, 0.25, and 0.1 in estimating the applied current, compared with σ ξ = 1 used in the baseline results above. For sake of demonstration, we focused on the data generated using the sinusoidal current defined in (36); similar results hold in the other cases.
The results in Figure 7 show that when σ ξ = 10, although the EnKF mean estimate was able to well-track the true solution, the estimated ±2 standard deviation curves around the mean are very large around the mean, reflecting a lack of confidence in the estimate. On the other hand, when σ ξ = 0.1, the filter is unable to well-track the true parameter and eventually diverges. In this case, while the filter is unable to track the true parameter, the estimated ±2 standard deviation curves are very tight around the mean, implying a high confidence in an incorrect estimate. For this example, the choice of σ ξ = 0.5 in the parameter random walk visually captures the true underlying parameter the best out of the values considered, as it keeps the EnKF estimated mean closely with a small standard deviation around the mean. In addition to the choice of σ ξ , the frequency of time-series data available also has a significant effect on the resulting parameter tracking estimates. To analyze this on the problem at hand, we subsampled the  . This means that data was considered in the observation step in the EnKF every first, second, or fifth millisecond (as compared to every 0.1 msec using the full data). Note that the filter still performed the prediction step every 0.1 msec as before, but the observation step was only preformed if data was available at that time point. The results in Figures 8 and 9 show that as data becomes more and more sparse, the parameter tracking estimate of the applied current parameter loses more and more characteristics of the underlying deterministic function. In particular, in Figure 8, as less data is available, the parameter tracking estimate of the pulsing step function begins lagging in estimating the steps and is unable to well maintain the shape. Similar results are seen in Figure 9 for the sinusoidal function, where the parameter tracking estimate loses its periodicity as less data becomes available. In both figures, it is clear that the sparser the data set, the more challenging it becomes for the filter to track the true underlying applied current function.

Discussion
The aim of this work was to utilize the EnKF with parameter tracking in estimating the time-varying applied current parameter in the Hodgkin-Huxley model. In particular, the parameter tracking algorithm was analyzed in estimating four different deterministic applied currents using synthetically-generated voltage data. We first verified that the algorithm was able to successfully track the underlying applied current, along with the unobserved model states, for each of the four test cases. In addition to tracking the applied current I(t), baseline results show that the filter is able to accurately estimate the unobserved states n(t), m(t), and h(t) from the generated voltage data. Further numerical experiments were conducted to analyze how the parameter tracking estimates of the applied current parameter are affected under different implementation conditions, namely, when changing the standard deviation of the parameter drift term in (32) and when the algorithm is provided increasingly less amounts of data.
In general, using the augmented EnKF with parameter tracking as described in Section 3, we were able to well track the underlying applied current functions in each of the four cases considered, establishing the baseline results shown in Figures 3-6. From these results, we were able to further explore the effects of two different aspects of the parameter tracking implementation: the choice of the standard deviation of the parameter drift term in (32), and the availability of time-series data.
We found that using different values for σ ξ resulted in drastically different levels of accuracy and confidence in the resulting time-varying parameter estimates for the applied current. As shown in Figure  7, when using σ ξ = 10 here for the parameter drift, the resulting mean estimate tracked well, however the resulting confidence in the estimate was low, as reflected in the wide range between the estimated ±2 standard deviation curves. Oppositely, when using σ ξ = 0.1, the parameter tracking estimate diverged, returning an estimate of the parameter which was far off from the true solution despite having high confidence in the estimate. Based on the results in Figure 7, it was determined that for the data considered, σ ξ = 0.5 was the best choice out of the tested values, since the algorithm was able to accurately track the mean with less uncertainty reflected in the resulting ±2 standard deviation curves around the mean.
We also explored how the accuracy of the parameter tracking estimate changed when the generated data was subsampled, resulting in fewer data points being used as updating information in the filter. The results in Figures 8 and 9 show that as less data was made available, the parameter tracking algorithm had increasing difficulty in estimating the true underlying applied current function, losing structural features such as the step onset of the pulsing step current and the periodicity of the sinusoidal current.
While the focus of this work was on estimating four different deterministic forms of the time-varying applied current parameter, in future work we aim to estimate stochastic forms of the current, as well as estimating applied currents relating to networks of neurons. In addition to the applied current, we are also interested in applying parameter tracking methodology to estimate the α(V) and β(V) rate functions as additional unknown parameters that vary with voltage, which would be particularly useful in applications of the Hodgkin-Huxley model where these rate functions do not necessarily share the same parameterized forms as in (8)-(11) and (15)- (16).
Additional future work includes comparing our results for the Hodgkin-Huxley model with results using other single neuron models, such as the FitzHugh-Nagumo and Hindmarsh-Rose models [21][22][23]. We also aim to apply our results to further biomedical applications utilizing Hodgkin-Huxley dynamics to model various neurodegenerative diseases affecting the function of neurons and ionic channels, such as Alzheimer's disease and amyotrophic lateral sclerosis [24][25][26].