Biophysically interpretable inference of single neuron dynamics

Abstract

Identification of key ionic channel contributors to the overall dynamics of a neuron is an important problem in experimental neuroscience. Such a problem is challenging since even in the best cases, identification relies on noisy recordings of membrane potential only, and strict inversion to the constituent channel dynamics is mathematically ill-posed. In this work, we develop a biophysically interpretable, learning-based strategy for data-driven inference of neuronal dynamics. In particular, we propose two optimization frameworks to learn and approximate neural dynamics from an observed voltage trajectory. In both the proposed strategies, the membrane potential dynamics are approximated as a weighted sum of ionic currents. In the first strategy, the ionic currents are represented using voltage dependent channel conductances and membrane potential in a parametric form, while in the second strategy, the currents are represented as a linear combination of generic basis functions. A library of channel activation/inactivation and time-constant curves describing prototypical channel kinetics are used to provide estimates of the channel variables to approximate the ionic currents. Finally, a linear optimization problem is solved to infer the weights/scaling variables in the membrane-potential dynamics. In the first strategy, the weights can be used to recover the channel conductances, and the reversal potentials while in the second strategy, using the estimated weights, active channels can be inferred and the trajectory of the gating variables are recovered, allowing for biophysically salient inference. Our results suggest that the complex nonlinear behavior of the neural dynamics over a range of temporal scales can be efficiently inferred in a data-driven manner from noisy membrane potential recordings.

Appendix

1.1 A.1 Computational procedure for strategy 1

We will describe strategy 1 using Example 1 without loss of generality. This procedure is also summarized as an algorithm (see Algorithm 1).

Library selection

We first constructed a library of the activation curves and the time-constant curves following Eq. (4) where the functional forms of $\alpha _{n_{j}}(V)$ and $\beta _{n_{j}}(V)$ were selected as given in Izhikevich (2007) with unknown parameters $d,e,f,\bar {d},\bar {e},\bar {f}$. These functions are:

$$ \begin{array}{@{}rcl@{}} \alpha_{n_{j}}(V)&=&\frac{d(V-e)}{(exp((-V+e)/f)-1)},\\ \beta_{n_{j}}(V)&=&\bar{d} exp\left( \frac{-(V-\bar{e})}{\bar{f}}\right), \end{array} $$

(14)

for the activation channels and

$$ \begin{array}{@{}rcl@{}} \alpha_{n_{j}}(V)&=&\bar{d} exp\left( \frac{-(V-\bar{e})}{\bar{f}}\right), \\ \beta_{n_{j}}(V)&=&\frac{d}{(1+exp((-V+e)/f))}, \end{array} $$

(15)

for the inactivation channels. The interval [V_b,V_th] is empirically observed to be [− 75, 30]mV with resting potential close to − 65mV (Fig. 13). Therefore, we constrained our design parameters

$$ e\in\{-35,-40, -45,-50, -55,-60,-65,-70\}, $$

$\bar {e} = e$, d ∈ (0, 5], f ∈ [1, 30], $\bar {d}= d$, and $\bar {f} = f$.

Fig. 13

Schematic of Strategy I. a The voltage trace (red color) and its time derivative (blue color) are used to determine the function domain for the activation/inactivation and time-constant curves. b The voltage dependent functions $\alpha _{n_{j}}(V)$ and $\beta _{n_{j}}(V)$ are drawn from the literature and used to generate the library of voltage dependent activation curves and time-constant curves. Here, blue solid lines are the curves associated with activating channels while the red solid lines are the curves associated with inactivating channels. c The current Eq. (2) are forward simulated to produce candidate gating variable trajectories, i.e., m(t) and h(t). Here, blue solid lines are the curves associated with sampled gating variables obtained from the library while the red solid lines are the gating variables corresponding to the ground truth model of the 4th order HH model. d We sample a set of candidate gating variable trajectories. Regression is then performed to produce a candidate model

Channel pruning

The design parameters were uniformly randomly sampled from these intervals to form a set of $\alpha _{n_{j}}(V)$, $\beta _{n_{j}}(V)$. Because these functions are the channel opening and closing rates, we expect them to attain non-zero values over voltage ranges in which the observed voltage trace displays a non-zero derivative. In particular, we prune those functions that meet the criteria

$$ \max_{V \in [-75,-50]} \left|\alpha_{n_{j}}(V)\right| < 10^{-5} $$

(16)

and

$$ \max_{V \in [-5,5]} \left|\beta_{n_{j}}(V)\right| < 10^{-5} $$

(17)

After pruning, we are left with a set of nontrivial opening and closing rate functions. From these, we construct the final library of activation and time-constant via (4).

Gating variables

Using the constructed library, the channel dynamical Eq. (3) were forward simulated to obtain a set of corresponding gating variable trajectories (see Fig. 13). From this, we built a set of regression functions using these gating variable trajectories in the current Eq. (2). Here, the values of a,b in the current equations were constrained to be in the set {0, 1, 2, 3, 4}, and

$$ \begin{array}{@{}rcl@{}} {\Phi}(V,\hat{I}_{x},I)(t)&=& [{m_{1}^{a}}(t){h_{1}^{b}}(t)V(t),-{m_{1}^{a}}(t){h_{1}^{b}}(t),\ldots, \\ &&{m_{l}^{a}}(t){h_{l}^{b}}(t)V(t),-{m_{l}^{a}}(t){h_{l}^{b}}(t) ]. \end{array} $$

(18)

Monte Carlo sampling and regression

To perform the regression and estimate the unknown parameters, we performed a Monte Carlo analysis. For Example 1, we used 150ms of (simulated) data with a sampling period of 0.01ms. Next, we picked in each iteration $\phi (V,\hat {I}_{x},I) \in \mathbb {R}^{15000 \times 11}$ (l = 11), a subset of Φ(t), to make the regression problem sufficiently over-determined. These functions were picked randomly with uniform probability, at each iteration. Then we solved the linear regression problem as described in Eq. (11) using a standard SVD algorithm. During each iteration, the estimated parameters were recorded along with their residual/regression error (R_n) and the condition number (C_n) of the regression matrix (see Fig. 14). We performed this iterative procedure for 1000 iterations.

Fig. 14

The procedure from Fig. 13 is repeated multiple times. a The Monte Carlo procedure produces many candidate models, each associated with a residual error R_n and condition number C_n. We select a model that minimizes error to within a standard deviation of the mean condition number, here indicated by the red arrow. b In the case of Example 1, we depict the final set of channels (red: passive; blue: active) among all candidate channels over Monte Carlo iterations

Model selection

The final model is selected as the one that minimizes the residual error to within one standard deviation of the mean condition number over all iterations.

1.2 A.2 Computational procedure for strategy 2

In strategy 2, the library selection step and the channel pruning steps are the same as described in Appendix A. We will describe the rest of the procedure next.

Gating variables

Using the library of activation/inactivation curves and time-constant curves, the channel dynamical Eq. (3) were forward simulated to obtain a set of corresponding gating variable trajectories. However, instead of using functional forms for the current equations as in strategy 1, we here use an (artificial) neural network architecture (Lewis et al. 1998) to transform the gating variable trajectories into currents.

Neural network (NN) design

We use a random vector functional link network (Lewis et al. 1998) with random weights between the input and hidden layers. The weights from the hidden layer to the output layer W_x are free and to be regressed. The inputs to the NN were the voltage trajectory, the extrinsic current and the gating variables obtained from the previous step. We constructed the NN inputs in order to distinguish transient, persistent and leak currents. Specifically, we used three blocks of inputs : (i) {V (t),m(t),h(t)} (for approximating transient currents); (ii) {V (t),m(t)} (for approximating persistent currents); and (iii) V (t) (to approximate leak current) (see Fig. 10). The exogenous current I(t) was applied as the input to the output neuron scaled by a bias (B) which was free to be tuned.

Monte Carlo sampling and weight tuning

At each iteration, we picked six gating variables (five activation channels and one inactivation channel). We updated the weights of the hidden to output layer via linear regression as in strategy 1. We used 30 hidden layer neurons for each input block. We performed this analysis over 1500 iterations and retained the model that minimized the residual error.

1.3 A.3 Online learning rule and convergence analysis

Let the following equations describe the approximated neuronal dynamics as

$$ \begin{array}{@{}rcl@{}} \dot{\hat{V}}(t)&=&\hat{W}^{T}{\phi}(t)+\tilde{A}\tilde{V}(t)\\ \dot{\hat{n}}_{j}(t)&=&\frac{n_{j\infty}(\hat V)-\hat n_{j}}{\tau_{n_{j}}(\hat V)} \end{array} $$

(19)

where $\hat {W}$ consists of the estimated parameters and ϕ is calculated using the estimated gating variables and measured voltage, $\tilde {V}$ is the voltage estimation error, and $\tilde {A}> 0$ is a design parameter. To bring the estimated parameters close to the true parameters, update the estimated parameters using the adaptation rule as follows: First define the performance measure as

$$ \begin{array}{@{}rcl@{}} E_{o} &= &\frac{1}{2}\beta(\dot{V}(t)-(\dot{\hat{V}}(t)-\tilde{A}\tilde{V}(t)))^{2}\\ && + \frac{1}{2}\alpha(V(t)-\hat{V}(t))^{2}+\frac{1}{2}\gamma\hat{W}^{T}(t)\hat{W}(t). \end{array} $$

(20)

Using the performance measure (20), define the learning rule as

$$ \begin{array}{@{}rcl@{}} \dot{\hat{W}}(t)&=& \alpha \phi(t)\tilde V(t)+ \beta \phi(t)[V(t+h)\\ & &-V(t)-h \hat{W}^{T} \phi]-\gamma\hat W(t) \end{array} $$

(21)

with $\alpha , \beta , h>0, \hat W(0)=\hat W_{0}$ and $\tilde {V}(t)=V(t)-\hat {V}(t)$. Define the error in the estimated parameters as $\tilde W(t) = W(t)-\hat {W}(t)$. To ensure persistency of excitation, the data (history) can be stacked to get the regression vector ϕ, voltage estimation error, $\tilde {V}$ in the update equation (by using the historical data, V (0),…,V (t) and $\hat {V}(0),\ldots , \hat {V}(t)$ ). Here, γ is the sigma modification term (Lewis et al. 1998), which performs the role of a regularizer, while avoiding parameter drift .

To analyse the convergence properties, consider a quadratic energy function,

$$ L(t)=\frac{\alpha}{2}\tilde{V}^{2}(t)+\frac{1}{2}\tilde{W}^{T}(t)\tilde{W}(t). $$

Taking the time derivative, we get

$$ \frac{d}{dt}L(t) = \alpha\tilde{V}(t)\dot{\tilde{V}}(t)+\tilde{W}^{T}(t)\dot{\tilde{W}}(t). $$

Using the adaptation rule (21) and the approximated dynamcis (19), we have

$$ \begin{array}{@{}rcl@{}} \frac{d}{dt}L(t) &=& \alpha\tilde{V}(t)(W^{T}(t)\phi(t)-\hat{W}^{T}(t){\phi}(t)-\tilde{A}\tilde{V}(t)) \\ &&+\tilde{W}^{T}(t)(-\alpha \phi(t)\tilde V(t) -\beta \phi(t)[V(t+h)\\ &&-V(t)-h \hat{W}^{T} \phi]\\ &&+\gamma\hat W). \end{array} $$

(22)

This is simplified as

$$ \begin{array}{@{}rcl@{}} \frac{d}{dt}L(t) &= &-\tilde{A}\alpha\tilde{V}^{2}(t)-h\beta \tilde{W}^{T}(t)\phi(t)\tilde{W}^{T}\phi\\ &&+\tilde{W}^{T}(t)\gamma W-\gamma\tilde{W}^{T}(t)\tilde{W}(t). \end{array} $$

The slope of the energy function is thus bounded by

$$ \begin{array}{@{}rcl@{}} \frac{d}{dt}L(t) &\le & -\tilde{A}\alpha\tilde{V}^{2}(t) -h\beta \tilde{W}^{T}(t)\phi(t)\phi^{T}(t)\tilde{W}\\ &&+\frac{\gamma^{2} \epsilon}{2} W^{T}W-\left( \gamma-\frac{1}{2\epsilon}\right)\tilde{W}^{T}(t)\tilde{W}. \end{array} $$

This implies that the voltage estimation error and the parameter estimation errors converge to a bound defined by the term $\frac {\gamma ^{2}\epsilon }{2} W^{T}W$.

If we consider the voltage estimation error dynamics, $\dot {\tilde {V}}(t)=\dot {V}(t)-\dot {\hat { V}}(t)$. Since we assumed that the exact channel dynamics are known in the analysis, the error dynamics is simplified as $\dot {\tilde {V}}(t)=\tilde {W}^{T}(t)\phi (t)-\tilde {A}\tilde V(t)$. The linear term $\tilde {A}\tilde V(t)$ is required to stabilize the model during the learning period. Furthermore, since the channel dynamics are not known, in general, we will have $\dot {\tilde {V}}(t)=\tilde {W}^{T}(t)\phi (t)+\hat W(t)(\phi (t)-\hat {\phi }(t))-\tilde {A}\tilde V(t)$. Here the regression vector, $\hat {\phi }(t)$, is calculated using a set of gating variables characterized by their activation curves, time-constants. As the errors in the gating variables converge, the regression vector, $\hat {\phi }(t) \to \phi (t)$.

By representing the dynamics in parametric form, it can also be noticed that the current, I(t), is required to be a non-constant stimulus inorder for the matrix ϕ(t)ϕ^T(t) to be nonsingular for the existence of a solution to the linear regression problem (due to the presence of leak channel). This is a potential source for numerical errors in the linear regression problem as the regression matrix can be poorly conditioned.

In the second strategy, due to the approximation, there will be a reconstruction/approximation error, o(V,n_j), i.e.,

$$ \dot{V}(t) = {\sum\limits_{x \in A} W_{x} {\Phi}_{x}(V(t),n_{j}(t))}+I(t) + o(V,n_{j}). $$

Therefore, the voltage estimation error dynamics can be obtained in this case as, $\dot {\tilde {V}}(t)=\tilde {W}^{T}(t)\phi (t)+o(V,n_{j})-\tilde {A}\tilde V(t)$. This approximation error will increase the error bound obtained earlier. By increasing the number of basis functions in the approximation, the approximation error can be reduced, however, at the expense of additional computations (Lewis et al. 1998). Since a random vector function link network is used, the regression problem is still linear.

It should also be noted that in the analysis above, the target parameters are assumed to be constants. The design parameters used in the learning algorithm (strategy 1) are α = 500,β = 750,γ = 0.0001, $\tilde {A} = 0.01$ and α = 250,γ = 0.001,β = 350,A = 1 for strategy 2, yielded results that are similar to that presented in the Examples (but not shown).