From Biophysical to Integrate-and-Fire Modeling

Van Pottelbergh, Tomas; Drion, Guillaume; Sepulchre, Rodolphe

doi:10.1162/neco_a_01353

Abstract

This article proposes a methodology to extract a low-dimensional integrate-and-fire model from an arbitrarily detailed single-compartment biophysical model. The method aims at relating the modulation of maximal conductance parameters in the biophysical model to the modulation of parameters in the proposed integrate-and-fire model. The approach is illustrated on two well-documented examples of cellular neuromodulation: the transition between type I and type II excitability and the transition between spiking and bursting.

1 Introduction

Integrate-and-fire models have a long history, dating back to the beginning of the twentieth century (Lapicque, 1907). Because of their simplicity, they have long served as phenomenological models for action potential generation in neurons. Over the past few decades, the original leaky integrate-and-fire model (Hill, 1936; Stein, 1965; Knight, 1972) has gradually been extended and modified to explain more neurological data and phenomena. Important examples include the replacement of the linear function by a quadratic or exponential nonlinearity (Ermentrout, 1996; Latham, Richmond, Nelson, & Nirenberg, 2000; Fourcaud-Trocmé, Hansel, van Vreeswijk, & Brunel, 2003), the linearization of the subthreshold dynamics in biophysical models (Mauro, Conti, Dodge, & Schor, 1970; Brunel, Hakim, & Richardson, 2014), and adding variables modeling the effects of refractoriness and adaptation (Wehmeier, Dong, Koch, & Van Essen, 1989; Mihalas & Niebur, 2008; Jolivet, Rauch, Lüscher, & Gerstner, 2006; Mensi et al., 2012; Pozzorini, Naud, Mensi, & Gerstner, 2013; Izhikevich, 2003; Brette & Gerstner, 2005; Drion, Franci, Seutin, & Sepulchre, 2012). Low-dimensional integrate-and-fire models have also been shown to be sufficient to reproduce experimental spike trains with great accuracy (Jolivet, Lewis, & Gerstner, 2004; Badel et al., 2008).

The simplicity of integrate-and-fire models nevertheless comes with an important limitation: they lack the physiological interpretation of a biophysical model (Van Pottelbergh, Drion, & Sepulchre, 2018). This is a severe obstacle when studying the role of cellular neuromodulation in a large network (Drion et al., 2018). While intrinsic neuromodulation is studied via changes of maximal conductances in biophysical models, mapping this effect in an integrate-and-fire model with no biophysical connection is challenging.

The objective of this article is to provide a quantitative computational bridge between a detailed single-compartment biophysical model and its compact representation by a multiscale integrate-and-fire model. Our motivation is to allow for a systematic mapping between the physiological parameters of a biophysical model and the abstract parameters of its integrate-and-fire approximation. Our models are not optimized to fit experimental spiking data. Instead, we aim at simplified models amenable to qualitative phase-portrait analysis and able to capture the qualitative changes of neuronal excitability induced by neuromodulators, which requires nonlinear subthreshold dynamics.

We exploit the architecture of the multi-quadratic integrate-and-fire (MQIF) model introduced in Drion et al. (2012) and further studied in Van Pottelbergh et al. (2018). The state variables of the model have the interpretation of the membrane potential filtered in different timescales. The model parameters relate to a local balance between positive and negative conductance^¹ in each timescale. This model has been quite successful at capturing important modulation properties in a robust and qualitative manner.

We introduce a multiscale integrate-and-fire model in section 2 and highlight a key feature of the proposed approach: fitting the parameters from voltage-clamp data rather than from current-clamp data. The parameters of this integrate-and-fire model to be identified from the biophysical model can be split into two groups. Section 3 discusses how to identify the ion current function of the integrate-and-fire model. The other parameters are treated in section 4, which also discusses local optimization to improve the result. We then apply this method to model modulation in two biophysical models from the literature in section 5. We end by discussing the limitations of the method and its connection to other methods for the analysis of neural behavior.

2 A Multiscale Integrate-and-Fire Model Structure

We consider a general multiscale neuronal model of the form

C \dot{V} = I_{app} - I_{ion} (V, V_{s}, V_{u s}, \dots),

(2.1)

τ_{s} \dot{V_{s}} = V - V_{s},

(2.2)

τ_{u s} {\dot{V}}_{u s} = V - V_{u s} . \dots

(2.3)

The voltage equation 2.1 is the classical Kirchhoff relationship of a biophysical model: the current $I_{ion}$ is the total ionic current intrinsic to the cellular membrane composition. Its voltage dependence is modeled through the voltage variable $V$ and lagged variables (⁠ $V_{s}$ ⁠, $V_{u s}$ ⁠, $\dots$ ⁠) that model the voltage filtered in distinct timescales (slow, ultraslow, $\dots$ ⁠). Those variables differ from the gating variables of traditional conductance-based models, but closely relate to the equivalent potentials originally defined in Kepler, Abbott, and Marder (1992). Note that the dynamics of those filters is chosen to be linear. The nonlinearity of the model is entirely concentrated in the scalar function $I_{ion}$ ⁠.

While the model 2.1 to 2.3 can exhibit highly nonlinear behavior in current clamp, this is not the case in voltage clamp. The voltage-clamped relation from $V$ to $I_{app}$ is shown in Figure 1. It has the classic structure of a parallel Wiener system. The identification of this nonlinear fading memory model is much simpler than the direct identification of the current clamp dynamics (Burghi, Schoukens, & Sepulchre, 2019). In particular, our methodology identifies the nonlinear function $I_{ion}$ from the response of the biophysical model to voltage step inputs.

When a reset is added to the multiscale model, as in equations to , the model is converted into an integrate-and-fire model. This simple integrate-and-fire model structure includes many models in the literature: the one-dimensional leaky, quadratic, and exponential integrate-and-fire models, but also the multidimensional Izhikevich (Izhikevich, 2003), AdEx (Brette & Gerstner, 2005), and MQIF models (Drion et al., 2012; Van Pottelbergh et al., 2018):

if V \geq V_{max} : C \dot{V} = I_{app} - I_{ion} (V, V_{s}, V_{u s}, \dots) V \leftarrow V_{r}

(2.4)

τ_{s} {\dot{V}}_{s} = V - V_{s} V_{s} \leftarrow V_{s, r}

(2.5)

τ_{u s} {\dot{V}}_{u s} = V - V_{u s} V_{u s} \leftarrow V_{u s} + Δ V_{u s} \dots

(2.6)

Figure 1:

View large Download slide

Block diagram representation of the multiscale model 2.1 to 2.3, in voltage clamp. Note the additional voltage $V_{f}$ ⁠, accounting for the dynamics of the fast ion channels. While this voltage will be merged with $V$ in the final multiscale model, $τ_{f}$ is necessary for the identification of $I_{ion}$ (section 3). This is a parallel Wiener system with the linear systems on the left feeding into the nonlinearity on the right.

In the remainder of this article, we concentrate on the integrate-and-fire version of the model. Be aware, however, that there is a direct correspondence between model 2.1 to 2.3 and model 2.4 to 2.6 provided that model 2.1 to 2.3 generates spikes. Each spike in the continuous-time model is replaced by a reset mechanism in model 2.4 to 2.6.

The motivation for this substitution is computational: the reset avoids the stiff integration of a spike, which can be a significant computational gain in the numerical integration of a large-scale spiking model. Instead, the numerical integration is stopped at threshold and reinitialized at a novel initial condition specified by a discrete rule. We insist that the reset operation in our integrate-and-fire model has no other role than short-cutting the (otherwise tedious) numerical integration of the spike in the differential equation. This is in contrast with other integrate-and-fire models that design the reset rule to generate solutions that do not correspond to solutions of the differential equation (see Van Pottelbergh et al., 2018, for details). This constraint is key to retaining a direct correspondence to the physiology of a biophysical model.

An important property of the integrate-and-fire model 2.4 to 2.6 is that the nonlinear function $I_{ion}$ only needs to be identified in the subthreshold voltage range $V < V_{max}$ ⁠.

3 Identification of the Ion Current Function from a Biophysical Model

Our methodology approximates a given biophysical model by an integrate-and-fire model of the type 2.4 to 2.6. We make a distinction between the structural parameters of the model (the time constants and the reset parameters) and the single scalar nonlinear function $I_{ion}$ ⁠. In this section, we determine $I_{ion}$ for a given set of structural parameters. In the next section, we address the determination of the structural parameters and how the design can be optimized by iterating between the identification of the total ionic current and the identification of the structural parameters.

Given a biophysical model and choice of structural parameters, we propose to identify the ionic current

I_{ion}

so as to optimize the matching between voltage-clamp experiments on the biophysical model and on the integrate-and-fire approximation. We choose this criterion because it combines physiological relevance and computational tractability. It is physiologically relevant because biophysical models are developed from voltage-clamp experiments in the first place. It is also computationally tractable because a voltage-clamp step response can be calculated in closed form in both a biophysical model and in the integrate-and-fire model. In a biophysical model, each gating variable obeys a differential equation of the form

τ_{x_{i}} (V) {\dot{x}}_{i} = x_{i, \infty} (V) - x_{i},

(3.1)

which becomes linear for a fixed value of

V

⁠. The solution after a step from

V_{0}

to

V_{step}

at

t = 0

is described by

x_{i} (t) = x_{i, \infty} (V_{0}) + [x_{i, \infty} (V_{step}) - x_{i, \infty} (V_{0})] \cdot [1 - e^{- t / τ_{x_{i}} (V_{step})}],

(3.2)

assuming the membrane potential is at equilibrium at

t = 0

⁠.

The response of the total ionic current $I_{ion}$ in the biophysical model to a voltage-clamp step is then obtained by direct substitution of $V (t) = V_{step}$ and $x_{i} (t)$ (given by equation 3.2) in the expression for the total ionic current. The same calculation holds in the integrate-and-fire model 2.4 to 2.6 to obtain an expression of $I_{ion} (V (t), V_{s} (t), V_{u s} (t), \dots)$ ⁠.

The details of this simple idea are provided in the next sections.

3.1 Two-Timescale Models

We start by describing the method for biophysical models that can be described well by a fast and slow timescale, which are well separated (⁠

τ_{s} ≫ τ_{f}

⁠). We consider a simple voltage-clamp step experiment as in Figure 2, stepping from the initial voltage

V_{0}

to the final voltage

V_{step}

⁠. Assuming

τ_{f}

is a good approximation of the time constants of the fast gating variables, they will have approximately reached their steady-state value

x_{i, \infty} (V_{step})

at

t = 3 τ_{f}

⁠. At the same time, because of timescale separation, the slow gating variables will not have significantly changed from their steady-state value

x_{i, \infty} (V_{0})

⁠.

Figure 2:

View large Download slide

Voltage clamp step experiment on a two-timescale model—the Connor-Stevens model (Connor & Stevens, 1971; Connor, Walter, & McKown, 1977). The voltage is stepped from $V_{0}$ to $V_{step}$ ⁠. Assuming $V = V_{f}$ ⁠, the value of $I$ at $3 τ_{f}$ will give a good approximation of $I_{ion} (V_{step}, V_{0})$ of the multiscale model.

Similarly, in the multiscale integrate-and-fire model, the fast voltage

V_{f}

will have approximately reached its steady-state value

V_{step}

at

t = 3 τ_{f}

⁠. Under the assumption of

τ_{s} ≫ τ_{f}

⁠, the slow voltage

V_{s}

will still be approximately equal to

V_{0}

⁠. After merging the variables

V

and

V_{f}

in the multiscale model, the following observation can be made for the voltage-clamp experiment in Figure 2:

I (3 τ_{f}) \approx I_{ion} (V_{step}, V_{0}) .

(3.3)

The interpretation of this expression is that the value of the ion current function of the multiscale model evaluated at $V = V_{step}$ and $V_{s} = V_{0}$ can be read at $t = 3 τ_{f}$ from a voltage-clamp step experiment from $V_{0}$ to $V_{step}$ ⁠. Given that $V_{0}$ and $V_{step}$ can be chosen arbitrarily, this gives a simple method to evaluate $I_{ion}$ on the biophysical model and use its value in the function $I_{ion} (V, V_{s})$ of the multiscale model.

Instead of doing this matching using actual voltage-clamp experiments, we can use the analytical solution for voltage-clamp steps on biophysical models , resulting in the expression

x_{i} = x_{i, \infty} (V_{0}) + [x_{i, \infty} (V_{step}) - x_{i, \infty} (V_{0})] \cdot [1 - e^{- 3 τ_{f} / τ_{x_{i}} (V_{step})}]

(3.4)

for each gating variable. Substituting this expression in the ion current equation of the biophysical model then makes

I_{ion} (V_{step}, V_{0})

of the integrate-and-fire model a function of the equations for the biophysical model and the time constant

τ_{f}

⁠. It is clear that in the limit of infinite timescale separation, the formula becomes a simple substitution of the fast and slow gating variables by their steady-state values at

V_{step}

and

V_{0}

⁠, respectively:

lim_{τ_{x_{i}} \to 0} x_{i} = x_{i, \infty} (V_{step}),

(3.5)

lim_{τ_{x_{i}} \to + \infty} x_{i} = x_{i, \infty} (V_{0}) .

(3.6)

Figure 3:

View large Download slide

Voltage-clamp step experiment on a three-timescale model—the model of Drion et al. (2018). The voltage is stepped from $V_{0}$ to $V_{step, 1}$ and then to $V_{step, 2}$ after $3 τ_{s}$ ⁠. Assuming $V = V_{f}$ ⁠, the value of $I$ at $3 τ_{s} + 3 τ_{f}$ will give a good approximation of $I_{ion} (V_{step, 2}, V_{step, 1}, V_{0})$ of the multiscale model.

3.2 Three-Timescale Models

We apply the same idea to three-timescale models to find $I_{ion} (V, V_{s}, V_{u s})$ as a function of the solution of a voltage clamp experiment. The simple procedure for two-timescale models cannot be used anymore, as it would always couple the slow and ultraslow voltage, and therefore only evaluate $I_{ion} (V_{step}, V_{0}, V_{0})$ ⁠. This can be resolved by devising a slightly more complex voltage clamp step experiment (see Figure 3): starting at $V_{0}$ ⁠, stepping to $V_{step, 1}$ and, finally, $V_{step, 2}$ after $3 τ_{s}$ ⁠.

Following the same reasoning as before, after $3 τ_{s} + 3 τ_{f}$ the fast gating variables will be approximately at their steady-state value $x_{i, \infty} (V_{step, 2})$ ⁠. Similarly the slow gating variables will approximately be at the steady-state value $x_{i, \infty} (V_{step, 1})$ ⁠, while the ultraslow gating variables will still be close to $x_{i, \infty} (V_{0})$ ⁠. The value of the current at time $3 τ_{s} + 3 τ_{f}$ will therefore be a good approximation of $I_{ion} (V_{step, 2}, V_{step, 1}, V_{0})$ of the multiscale model, assuming $V = V_{f}$ ⁠.

The necessary substitution for the gating variables in the ion current equation is given by

x_{i} = {\bar{x}}_{i} + [x_{i, \infty} (V_{step, 2}) - {\bar{x}}_{i}] \cdot [1 - e^{- 3 τ_{f} / τ_{x_{i}} (V_{step, 2})}],

(3.7)

{\bar{x}}_{i} = x_{i, \infty} (V_{0}) + [x_{i, \infty} (V_{step, 1}) - x_{i, \infty} (V_{0})] \cdot [1 - e^{- 3 τ_{s} / τ_{x_{i}} (V_{step, 1})}],

(3.8)

which is simply the solution of the gating variable equations at

t = 3 τ_{s} + 3 τ_{f}

for the voltage-clamp experiment shown in Figure 3. The limits for the time constants going to zero and infinity are similar to those in the case of two-timescales:

lim_{τ_{x_{i}} \to 0} x_{i} = x_{i, \infty} (V_{step, 2}),

(3.9)

lim_{τ_{x_{i}} \to + \infty} x_{i} = x_{i, \infty} (V_{0}),

(3.10)

while

x_{i} \approx x_{i, \infty} (V_{step, 1})

if

τ_{x_{i}} \approx τ_{s}

and

τ_{s} ≫ τ_{f}

⁠.

3.3 Precompensation of Voltages in Absence of Timescale Separation

The method in the previous sections assumes that each gating variable of the biophysical model can be grouped into one of three categories: fast, slow, or ultraslow. Furthermore, the approximations rely on those timescales to be well separated from each other.

When the model does not show clear timescale separation, the assumption that the slower voltages will not significantly change after a step no longer holds. As long as the dynamics of the different gating variables can still be grouped in two or three timescales, it is possible to modify the described methods to compensate for the dynamics of the slower variables during the step experiment. In essence, the voltages used in the voltage clamp step experiments are modified so that the slower voltage(s) reach the desired values at the specified times $t = 3 τ_{f}$ (or $t = 3 τ_{f} + 3 τ_{s}$ for three-timescale models). This is illustrated in Figure 4.

In the case of two timescales, only the initial voltage of the voltage-clamp step needs to be changed. The value of

V_{s}

at

t = 3 τ_{f}

can be computed given the initial voltage

V_{0}

and the step voltage

V_{step}

⁠:

V_{s} (3 τ_{f}) = V_{0} + (V_{step} - V_{0}) (1 - e^{- 3 τ_{f} / τ_{s}}) .

(3.11)

For this voltage to reach the desired value

V_{s}^{*}

at

t = 3 τ_{f}

⁠, the above expression can be solved for

V_{0}

⁠:

V_{0} = [V_{s}^{*} - V_{step} (1 - e^{- 3 τ_{f} / τ_{s}})] / e^{- 3 τ_{f} / τ_{s}} .

(3.12)

Making this additional substitution in results in the substitution expression to obtain the value for

I_{ion} (V_{step}, V_{s}^{*})

⁠.

The same idea can be applied to three-timescale models.

V_{step, 1}

is calculated in the same way as

V_{0}

in the two-timescale case:

V_{step, 1} = [V_{s}^{*} - V_{step, 2} (1 - e^{- 3 τ_{f} / τ_{s}})] / e^{- 3 τ_{f} / τ_{s}} .

(3.13)

This can be done because

V_{s}

can be assumed to be at steady state at

t = 3 τ_{s}

⁠.

V_{0}

can then be computed using

V_{step, 1}

⁠. For simplicity, it is assumed that

τ_{u s} ≫ τ_{f}

⁠, giving

V_{0} = [V_{u s}^{*} - V_{step, 1} (1 - e^{- 3 τ_{s} / τ_{u s}})] / e^{- 3 τ_{s} / τ_{u s}} .

(3.14)

If

τ_{u s} ≫ τ_{f}

does not hold, the expression becomes slightly more complicated:

V_{0} = [\bar{V} - V_{step, 1} (1 - e^{- 3 τ_{s} / τ_{u s}})] / e^{- 3 τ_{s} / τ_{u s}},

(3.15)

\bar{V} = [V_{u s}^{*} - V_{step, 2} (1 - e^{- 3 τ_{f} / τ_{u s}})] / e^{- 3 τ_{f} / τ_{u s}} .

(3.16)

After making these substitutions in equations 3.7 and 3.8, the substitution expressions can be used to obtain the value for $I_{ion} (V_{step, 2}, V_{s}^{*}, V_{u s}^{*})$ ⁠.

Figure 4:

View large Download slide

Illustration of the precompensation of step voltages to reach the desired values of the slow voltages at $t = 3 τ_{f}$ or $t = 3 τ_{f} + 3 τ_{s}$ ⁠. (Top) The choice of $V_{0}$ makes $V_{s} = V_{s}^{*}$ at $t = 3 τ_{f}$ for two-timescale models. (Bottom) The choice of $V_{0}$ ⁠, $V_{step, 1}$ and $V_{step, 2}$ leads to $V_{s} = V_{s}^{*}$ and $V_{u s} = V_{u s}^{*}$ at $t = 3 τ_{f} + 3 τ_{s}$ for three-timescale models.

3.4 Calcium Dynamics

Many biophysical models also contain calcium-gated ion channels apart from the more common voltage-gated ion channels. The conductance of these channels is usually described by a (nonlinear) function of the calcium concentration instead of a product of gating variables. The calcium dynamics obey a system of differential equations of the type

τ_{[{Ca}^{2+}]} \dot{[{Ca}^{2+}]} = {[{Ca}^{2+}]}_{\infty} (V, x_{Ca, i}, \dots) - [{Ca}^{2+}],

(3.17)

τ_{x_{Ca, i}} (V) {\dot{x}}_{Ca, i} = x_{Ca, i, \infty} (V) - x_{Ca, i} \dots

(3.18)

where

x_{Ca, i}

are the gating variables of the calcium channels. This will generally not result in a simple substitution expression. However, under the assumption that the calcium concentration changes much more slowly than the calcium gating variables, an approximate expression can be obtained. Assuming the calcium gating variables are at steady state whenever the voltage is constant during the voltage-clamp experiment, the expressions of the previous sections can be reused.

[{Ca}^{2+}] (V_{s})

and

[{Ca}^{2+}] (V)

in these expressions are obtained by substituting the values of the calcium gating variables at

t = 3 τ_{s}

and

t = 3 τ_{s} + 3 τ_{f}

⁠, respectively. Because of its simplicity, this approximation will be used in the rest of this article.

4 Iterative Optimization of the Structural Parameters

The previous section showed how to identify the function $I_{ion}$ for a given set of structural parameters. In this section, we discuss how to initialize those parameters and iteratively optimize them using current clamp data.

4.1 Initialization of the Time Constants

To estimate the time constants, we assume that the biophysical model can be described by the multiscale model 2.1 to 2.3 in the subthreshold regime. This is of course an approximation, but it allows finding a simple procedure for an initialization of the time constants.

The Wiener structure of parallel linear systems followed by a static nonlinearity (see Figure 1) can be exploited to estimate the time constants of the model. Applying a sufficiently small voltage-clamp signal will reveal the dynamics of the linear systems, since

\begin{matrix} I_{app} & \approx & C \dot{V} + α_{0} + α_{1} V + α_{2} V_{f} + α_{3} V_{s} + α_{4} V_{u s} \\ α_{1} & = & \frac{\partial I_{ion}}{\partial V}, α_{2} = \frac{\partial I_{ion}}{\partial V_{f}}, \dots \end{matrix}

where the coefficients

α_{i}

depend on the membrane potential around which the experiment is performed. The contribution of

C \dot{V}

is easily removed as it only produces a sharp pulse at the time of the step. Repeating this experiment at multiple subthreshold voltages gives a more robust estimation, as some of the coefficients

α_{i}

might be small around a single voltage.

Classical system identification techniques can then be used to obtain the parameters of the linear systems from this voltage-clamp data. We choose an extension by Miller, Hulett, McLaughlin, and de Callafon (2013) of the Ho-Kalman-Kung algorithm (Ho & Kalman, 1966; Zeiger & McEwen, 1974; Kung, 1978) for step responses. The measured-voltage clamp step responses are treated as the different outputs of the system to a single step input. The algorithm allows selecting a model order (number of voltage variables) based on the Hankel singular values. The poles can also be restricted to be real and stable ones, to fit the structure of the multiscale model. The sought time constants $τ_{f}, τ_{s}, τ_{u s}$ will hence be the negative inverse of the identified poles.

The described method to obtain initial estimates for the time constants is by no means optimal or the only one possible. A common approach used in Wiener system identification is the best linear approximation (BLA) (Enqvist & Ljung, 2005; Schoukens & Rolain, 2012). While this approach can result in more accurate estimates, it requires more data and the gaussian input signals used are less suited to the voltage-clamp setting. Furthermore, the time constants only have to be accurate enough to make the optimization procedure of section 4.3 converge to a reasonable set of parameters.

4.2 Initialization of the Other Structural Parameters

Apart from the time constants and $I_{ion}$ ⁠, the only parameters left to identify in model 2.4 to 2.6 are the membrane capacitance $C$ ⁠; the reset parameters $V_{r}$ ⁠, $V_{s, r}$ ⁠, and $Δ V_{u s}$ ⁠; and the cutoff voltage $V_{max}$ ⁠. Those parameters are easily estimated from current clamp simulations of the biophysical model.

$V_{max}$ is taken to be just after the onset of the spike, around the maximum of the first derivative of the membrane potential. As long as $V_{max}$ lies sufficiently above threshold, its precise value should not influence the behavior of the integrate-and-fire model. It is therefore fixed during the optimization step.

The voltage reset $V_{r}$ is always taken to be equal to $V_{max}$ ⁠. The reset values of the slower voltages are parameters to be optimized but are physiologically constrained. A reasonable initialization is to set $V_{s, r}$ sufficiently (e.g., 20 mV) above $V_{r}$ ⁠. The parameter $Δ V_{u s}$ is constrained to be positive, since $V_{u s}$ would only increase after a spike in model 2.1 to 2.3. Therefore, it can be initialized at 0 mV.

Finally, the membrane capacitance $C$ is initially set to 1 $μ$ F cm $^{- 2}$ ⁠, as is often the case in a conductance-based model.

4.3 Local Parameter Optimization

The structural parameters of the model can be iteratively optimized to improve the matching of representative current clamp data (see Figure 5, top).

Figure 5:

Test data generated using the Connor-Stevens model with g¯A=0mScm-2. (Top) Current clamp experiment with a decreasing ramp of the current Iapp. (Bottom) Data used to evaluate the cost function, where points for which V>Vmax are removed. Vs is obtained by filtering V with a first-order linear low-pass filter and reset to Vs,r after every spike, and Iion is evaluated using these V and Vs.

View large Download slide

Test data generated using the Connor-Stevens model with ${\bar{g}}_{A} = 0 mS {cm}^{- 2}$ ⁠. (Top) Current clamp experiment with a decreasing ramp of the current $I_{app}$ ⁠. (Bottom) Data used to evaluate the cost function, where points for which $V > V_{max}$ are removed. $V_{s}$ is obtained by filtering $V$ with a first-order linear low-pass filter and reset to $V_{s, r}$ after every spike, and $I_{ion}$ is evaluated using these $V$ and $V_{s}$ ⁠.

We choose to minimize a cost function based on the residual current, which has been used before for the parameter estimation of biophysical models (Morse, Davison, & Hines, 2001; Huys, Ahrens, & Paninski, 2006; Lepora, Overton, & Gurney, 2012). The residual current for our model is defined as

I_{res} = C \dot{V} - I_{app} + I_{ion} (V, V_{s}, V_{u s})

(4.1)

and is zero when the membrane current equation is satisfied. The test data provide

V

and

I_{app}

⁠, while the derivative

\dot{V}

can be obtained by numeric differentiation. The term

I_{ion}

depends on the chosen time constants, and to evaluate it, the voltage trace is filtered by the corresponding first-order linear low-pass filters. The action potentials themselves are eliminated by removing the data for which

V > V_{max}

⁠, as the model is only optimized in the subthreshold regime. Instead, the different voltages are reset after every spike using their respective reset parameters, as shown in Figure 5 (bottom). The voltage traces of Figure 5 (bottom) can then be used to evaluate

I_{ion} (V, V_{s}, V_{u s})

and thus the residual current.

It was found that simple least-squares minimization was sufficient for the purpose of this article. Apart from the time constants, the free parameters in the optimization are

C

and the reset values

V_{s, r}

and

Δ V_{u s}

⁠. The advantage of using the residual current over the residual voltage for the cost function is that the former does not suffer from the extreme sensitivity to tiny variations in the spike timing. This can be observed in Figure 6: the cost function based on the residual current is locally convex. In contrast, a cost function based on the residual voltage has many local minima and is much more irregular.

Figure 6:

Contour plots of two least squares cost functions as a function of τf and τs for the Connor-Stevens model with g¯A=0mScm-2. (Left) Cost function based on the residual current, as used in this article. (Center and right: Cost function based on the residual voltage (Vres=V-V*, where V is the voltage of the test data and V* is the voltage of the simulation of the integrate-and-fire model).

View large Download slide

Contour plots of two least squares cost functions as a function of $τ_{f}$ and $τ_{s}$ for the Connor-Stevens model with ${\bar{g}}_{A} = 0 mS {cm}^{- 2}$ ⁠. (Left) Cost function based on the residual current, as used in this article. (Center and right: Cost function based on the residual voltage (⁠ $V_{res} = V - V^{*}$ ⁠, where $V$ is the voltage of the test data and $V^{*}$ is the voltage of the simulation of the integrate-and-fire model).

An important property of the multiscale integrate-and-fire model is that its structural parameters have a clear interpretation. This makes it easy to initialize them at a reasonable value by inspection of a few voltage and current clamp experiments. The additional optimization procedure is a straightforward least-squares local optimization.

5 Integrate-and-Fire Modeling and Phase Portrait Analysis

The benefit of the proposed integrate-and-fire model is not purely computational. We now show that it is also amenable to phase portrait analysis, which provides mathematical insight into the initial biophysical model, regardless of its dimension.

For a two-timescale analysis, the phase portrait of the integrate-and-fire model is entirely characterized by two curves: the $V$ -nullcline $I_{ion} (V, V_{s}) = I_{app}$ and the $V_{s}$ -nullcline $V_{s} = V$ ⁠. In other words, the level curves of the identified ion current $I_{ion}$ determine the phase portrait of the model.

When the integrate-and-fire model is three-dimensional, we can describe its dynamics by considering the ultraslow variable $V_{u s}$ as a bifurcation parameter and by studying the family of phase portraits parameterized by $V_{u s}$ ⁠. In that sense, it can be said that the proposed integrate-and-fire model maps an arbitrary biophysical model to a family of phase portraits. This is very convenient for a qualitative understanding of the dynamical properties of the model.

The following sections illustrate the identification procedure and the phase portrait analysis of the integrate-and-fire model. The model is identified on two biophysical models from the literature; the classical Connor-Stevens model, whose behavior can be easily modulated by a change of maximal conductance, and the model by Drion et al. (2018), exhibiting a switchable slow negative conductance.

Figure 7:

View large Download slide

Phase portraits of the integrate-and-fire model obtained from the Connor-Stevens model with standard parameters (Connor et al., 1977). The $V$ - and $V_{s}$ -nullclines are drawn as full and dashed lines, respectively. (Left) The $V$ -nullclines are drawn for different values of $I_{app}$ ⁠: low (light blue), $I_{app}$ at the transcritical bifurcation in fast subsystem (medium blue), and high (dark blue). (Right) Detailed portion of the phase portrait (for the high value of $I_{app}$ ⁠), together with the corresponding trace of the limit cycle oscillation in orange.

5.1 Modulation of Excitability Type in the Connor-Stevens Model

The Connor-Stevens model (Connor & Stevens, 1971; Connor et al., 1977) is a six-dimensional conductance-based model for gastropod neuron somas. It has all the variables of the Hodgkin-Huxley model in addition to an extra potassium current $I_{A}$ ⁠. One of its characteristic features is type I excitability: the spiking frequency approaches zero when the applied current approaches the rheobase. This is in contrast to the type II excitability of the Hodgkin-Huxley model, whose spiking frequency makes a jump as the applied current is increased.

5.1.1 A Two-Timescale Integrate-and-Fire Approximation of the Connor-Stevens Model

Considering that this model was the example used by Kepler et al. (1992) for the equivalent potentials method, the method of this article is expected to work well on this model. Similar to the work of Kepler et al., we start by a two-dimensional reduction of the model using the method of section 3.1. The time constants found after the optimization procedure are $τ_{f} = 0.022 ms$ and $τ_{s} = 6.7 ms$ ⁠. They reflect the timescale separation of the gating variables, which span a range of 0.03 ms to 3 ms.

Figure 8:

View large Download slide

Phase portraits of the integrate-and-fire model (bottom) and f-I curves (top) of the Connor-Stevens model (blue) and the integrate-and-fire model (orange) for different values of ${\bar{g}}_{A}$ ⁠. The $V$ - and $V_{s}$ -nullclines are drawn as full and dashed lines, respectively, the trajectories after losing stability in orange. (Left) ${\bar{g}}_{A} = 0 mS {cm}^{- 2}$ ⁠. (Center) ${\bar{g}}_{A} = 47.7 mS {cm}^{- 2}$ ⁠. (Right) ${\bar{g}}_{A} = 200 mS {cm}^{- 2}$ ⁠.

The phase portrait of the two-timescale continuous-time model is displayed in Figure 7 (left). For a specific value of the applied current $I_{app}^{*}$ ⁠, the $V_{s}$ -nullcline intersects the $V$ -nullcline at the transcritical singularity (Franci, Drion, & Sepulchre, 2012). As the current increases from below to above the value $I_{app}^{*}$ ⁠, stability is lost in a saddle-node on invariant circle (SNIC) bifurcation. This was shown to be the mechanism of type I excitability in this model in Drion, O'Leary, and Marder (2015). The phase portrait of Figure 7 (right) shows the trajectory of spiking in the integrate-and-fire model 2.4 to 2.6 in orange. Since the model has a reset mechanism to replace the action potential generation, the periodic spiking occurs due to a hybrid limit cycle.

Figure 8 (top center) illustrates the frequency-current (f-I) curves of the model. The bifurcation voltage is well predicted by the integrate-and-fire model. Although the onset of the f-I curve is sharper than for the original model, both curves eventually converge.

The excitability type of the model can be changed easily by modulating the maximal conductance ${\bar{g}}_{A}$ ⁠. For ${\bar{g}}_{A} = 0 mS {cm}^{- 2}$ ⁠, the model is similar to the Hodgkin-Huxley model, which exhibits type II excitability. We again obtain a two-timescale integrate-and-fire approximation of this model by finding $I_{ion}$ and the structural parameters using the new value for ${\bar{g}}_{A}$ ⁠. Its phase portrait is shown in Figure 8 (bottom left) and is qualitatively the same as that of the FitzHugh-Nagumo model (FitzHugh, 1961; Nagumo, Arimoto, & Yoshizawa, 1962). There is a subcritical Hopf bifurcation resulting in type II excitability, confirmed by the f-I curve in Figure 8 (top left). The integrate-and-fire model loses stability at a slightly higher value of $I_{app}$ than the original model, but both f-I curves remain close to each other afterward.

On the other hand, increasing ${\bar{g}}_{A}$ above its nominal value of 47.7 mS cm $^{- 2}$ to 200 mS cm $^{- 2}$ results in a bistable phase portrait (see Figure 8, bottom right): for a specific range of $I_{app}$ ⁠, a stable hybrid limit cycle on the upper branch coexists with a stable fixed point on the lower branch of the $V$ -nullcline. The fixed point loses stability in a saddle-node bifurcation, while the limit cycle disappears in a fold limit cycle bifurcation. This was called type II* excitability in Drion, O'Leary et al. (2015), which is similar to type II excitability, but has a hysteretic f-I curve. Figure 8 (top right) shows that the integrate-and-fire model captures this hysteresis. Although the model starts spiking at a lower value of $I_{app}$ ⁠, the saddle-node bifurcation occurs at the same point. This can be explained by the fact that the fold limit cycle bifurcation is a global bifurcation, which is harder to capture than the local saddle-node bifurcation.

5.1.2 A Three-Timescale Integrate-and-Fire Approximation of the Connor-Stevens Model

The results in section 5.1.1 show that the two-timescale reduction of the Connor-Stevens model neuron model is useful to obtain a qualitative approximation of the behavior of the original model. The modulation of excitability type can be predicted from the phase portraits, and it is possible to construct an integrate-and-fire model that replicates this modulation. However, the resulting integrate-and-fire model lacks a quantitative approximation of the voltage trace (not shown) and f-I curve. Kepler et al. improved the quality of their reduction of the Connor-Stevens model by adding a third equivalent potential.

In an analogous effort, we approximate the original model by a three-timescale integrate-and-fire model using the method of section 3.2. The three timescales are not strongly separated, however, and the method of section 3.3 was used to adjust for this lack of separation. This resulted in a model with the time constants $τ_{f} = 0.037 ms$ ⁠, $τ_{s} = 1.7 ms$ ⁠, and $τ_{u s} = 2.8 ms$ ⁠. It is hypothesized that, as in the equivalent potentials method, this third timescale is necessary to account for the dynamics of the inactivation variable of the current $I_{A}$ ⁠, which are different from those of the other gating variables. Figure 9 (right) shows that the f-I curve for this three-timescale model almost perfectly matches that of the original model. A comparison of a current clamp simulation with a linearly increasing applied current for both models (see Figure 9, left) reveals very similar voltage traces.

While the f-I curves match well for the three-timescale model with ${\bar{g}}_{A} = 47.7 mS {cm}^{- 2}$ ⁠, we did not find a similar improvement for ${\bar{g}}_{A} = 200 mS {cm}^{- 2}$ (not shown). More work is necessary to determine whether the method is able to find a better approximation by using different current clamp data for the optimization of the structural parameters.

Figure 9:

View large Download slide

(Left) Voltage traces of the original Connor-Stevens model (top) and the three-timescale integrate-and-fire model (bottom) for $I_{app}$ linearly increasing from 8 to $12 μ A {cm}^{- 2}$ ⁠. (Right) f-I curves for the original Connor-Stevens model (blue) and the three-timescale integrate-and-fire model (orange).

5.2 Modulation between Spiking and Bursting Due to a Switchable Slow Negative Conductance

Our second illustration uses the eight-dimensional conductance-based model introduced in Drion et al. (2018). A critical physiological feature of this model is a T-type calcium channel with low-threshold activation in the slow timescale. T-type calcium channels endow the model with slow regenerativity by having an activation that is slower than the sodium channel activation. Furthermore, the channels are inactivated at a relatively low threshold (on an even slower timescale). This results in a slow, negative conductance that is switchable by an external current: it is switched on only by hyperpolarization. This current was hypothesized in Drion et al. (2018) as a critical mechanism for the neuromodulation of network states. Figure 10 (left) shows how the switchable negative conductance can be observed from a voltage clamp step experiment. The slope of the current response in the slow timescale determines the absence or presence of slow regenerativity (Franci, Drion, & Sepulchre, 2017).

A distinctive behavior of the model in current clamp is hyperpolari-zation-induced bursting (HIB). This means the model can be switched from slow spiking to bursting by sufficiently lowering the input current, as shown in Figure 10 (top).

Figure 10:

View large Download slide

Comparison of the response to a voltage-clamp step from a hyperpolarized (blue) and depolarized (orange) state for the original model of Drion et al. (2018) (left) and its reduction (right). The voltage is stepped from $-$ 90 mV and $-$ 60 mV to $-$ 42 mV. The part of the responses highlighted in green shows the presence of a slow, negative conductance in the hyperpolarized state and its absence in the depolarized state.

Figure 11:

View large Download slide

(Left) Voltage trace of the model of Drion et al. (2018) for $I_{app} = 0 μ A {cm}^{- 2}$ ⁠. (Right) Phase portrait of the three-dimensional integrate-and-fire approximation showing the absence of rest-spike bistability. The $V$ - and $V_{s}$ -nullclines are drawn as full and dashed lines, respectively, the stable (hybrid) limit cycle in orange.

We use the method of section 3.2 to construct an integrate-and-fire approximation of the biophysical model of Drion et al. (2018). The calcium dynamics are treated as explained in section 3.4. Its behavior can be studied using three timescales. To visualize the obtained reduced model, phase portraits in the $V$ - $V_{s}$ space are shown for different values of $V_{u s}$ ⁠. Assuming the ultraslow timescale is much slower than the slow timescale, the behavior of the model can be analyzed by looking at the fast-slow system in these phase portraits.

For $I_{app} = 0 μ A {cm}^{- 2}$ ⁠, the original model shows regular spiking, with a lower spiking frequency than during the bursting (see Figure 11, left). As this is a two-timescale behavior, $V_{u s}$ can be approximated by a constant value. Figure 11 (right) shows the phase portrait corresponding to this situation. There is no bistability between the resting and spiking state, and the phase portrait is the standard phase portrait of a spiking model.

Figure 12:

View large Download slide

(Top) Voltage trace of a burst in the model of Drion et al. (2018) with $I_{app} = - 1.6 μ A {cm}^{- 2}$ ⁠. (Bottom) Portions of fast-slow phase portraits of the three-dimensional integrate-and-fire approximation for different values of $V_{u s}$ ⁠. The blue curve is the $V$ -nullcline, the dashed curve is the $V_{s}$ -nullcline, and each orange curve is a trajectory of the phase portrait that captures a piece of the bursting attractor in the three-dimensional model.

When the current is lowered to $I_{app} = - 1.6 μ A {cm}^{- 2}$ ⁠, the model bursts (see Figure 12, top). The phase portraits during the different phases of the burst are shown in Figure 12 (bottom). The burst is initiated by the loss of stability in a saddle-node bifurcation on the lower branch of the $V$ -nullcline (see Figure 12, bottom left). The subsequent spiking, on the limit cycle on the upper branch of the $V$ -nullcline, causes $V_{u s}$ to increase. This increase moves the two branches of the $V$ -nullcline closer, resulting in bistability between a limit cycle and a stable fixed point (see Figure 12, bottom center). As $V_{u s}$ increases even further, the limit cycle is lost in a saddle-homoclinic bifurcation, moving the trajectory to the stable fixed point and thus terminating the burst (see Figure 12, bottom right).

Figure 10 (right) shows the voltage-clamp response of the reduced model for steps from two different voltages. Although different from the full biophysical model of Figure 10 (left), the integrate-and-fire model retains the switchable negative conductance responsible for hyperpolarization-induced bursting (HIB). Only for the hyperpolarized step, the slow response has a negative slope (highlighted in green). The fast response is instantaneous because $V_{f}$ and $V$ are a single variable in the reduced model.

Again, the obtained reduction can be used to construct an integrate-and-fire model. Figure 13 shows the voltage trace before and after a hyperpolarizing step current for the original model (top) and the reduced model (bottom). While the result is not a perfect quantitative match, it is quite accurate given the reduction of the number of variables from eight to three.

Figure 13:

View large Download slide

Voltage trace of the model of Drion et al. (2018) (top) and its three-dimensional integrate-and-fire approximation (bottom) for $I_{app}$ stepping from 0 to $- 1.6 μ A {cm}^{- 2}$ at $t = 500 ms$ ⁠.

6 Discussion

6.1 Validity of the Method

The method of this article attempts to match the voltage-clamp response of an integrate-and-fire and a biophysical model. The idea of modeling the voltage-clamp response to obtain a description of the neural dynamics is exactly what Hodgkin and Huxley (1952) did in their seminal work.

An important difference, however, is that we only match the response to a series of steps at a specific time. While this results in an analytical expression to reduce a biophysical model, it also means that the voltage-clamp response of both models will not necessarily match for other inputs and at other times. Nevertheless, under certain assumptions (explained below), the proposed method can be seen as a reduction method similar to the method of the equivalent potentials of Kepler et al. (1992).

The notion of equivalent potentials is a way to convert the gating variables in a biophysical model into potentials or voltages. Every variable $x_{i}$ is simply replaced by the associated voltage of the steady-state function: $V_{x_{i}} = x_{i, \infty}^{- 1} (x_{i})$ ⁠. In doing so, only the description of the system is changed, not its input-output dynamics. The new form, however, can make it easier to discover relationships among different variables, necessary to reduce the system.

Kepler et al. (1992) propose such a reduction method by grouping equivalent potentials with similar dynamics and replacing them by a weighted average of their group. The weights are found by optimizing the local approximation of the model. For the method to work, the dynamics of the equivalent potentials should fall into groups with similar dynamics. The method of this article provides a simpler alternative based on voltage clamp experiments, by making the additional assumption that the dynamics of each group of slower equivalent potentials can be described sufficiently well by a first-order linear low-pass filter as in equations 2.2 to 2.3, at least in the subthreshold regime.

This assumption might not always fully hold in practice, but the method can then still produce an acceptable reduction, although this should be carefully validated afterward.

6.2 Local Approximation by Multi-Quadratic Integrate-and-Fire Model

The multi-quadratic integrate-and-fire (MQIF) model (Drion et al., 2012; Van Pottelbergh et al., 2018) is a special case of the multiscale integrate-and-fire model to , in which

I_{ion}

is a sum of quadratic functions in

V

⁠,

V_{s}

⁠, and

V_{u s}

⁠:

I_{ion} (V, V_{s}, V_{u s}) = {\bar{g}}_{f} {(V - V^{0})}^{2} + {\bar{g}}_{s} {(V_{s} - V_{s}^{0})}^{2} + {\bar{g}}_{u s} {(V_{u s} - V_{u s}^{0})}^{2} .

(6.1)

The quadratics in

V

and

V_{s}

provide a normal form of the transcritical singularity in the fast subsystem, organizing the rest-spike bistability. The quadratic in

V_{u s}

models the ultraslow feedback necessary for bursting. The relative positions of the parameters

V^{0}

and

V_{s}^{0}

determine whether there is rest-spike bistability or not, but also influence the excitability type. This is illustrated in Figure 14, which sketches the phase portraits for the three possible regimes.

We can regard the MQIF model as a local approximation of the model in this article around a transcritical singularity. This singularity is identified in model to by finding the point

(V^{0}, V_{s}^{0})

for which

\frac{\partial I_{ion} (V, V_{s})}{\partial V} = \frac{\partial I_{ion} (V, V_{s})}{\partial V_{s}} = 0,

(6.2)

together with the condition that the determinant of the Hessian at

(V^{0}, V_{s}^{0})

should be negative (or zero, requiring further investigation).

I_{ion} (V^{0}, V_{s}^{0})

is then the current offset necessary to have the transcritical bifurcation occur at the same value of

I_{app}

⁠. The values of

{\bar{g}}_{f}

and

{\bar{g}}_{s}

are simply the elements on the diagonal of the Hessian at

(V^{0}, V_{s}^{0})

⁠.

Figure 14:

View large Download slide

Modulation of the excitability type in the MQIF model. Changing $V_{s}^{0}$ in the MQIF model results in different types of excitability: type II for $V_{s}^{0} < V^{0}$ (left), type I for $V_{s}^{0} = V^{0}$ (center) and type II* (Drion, O'Leary et al., 2015) for $V_{s}^{0} > V^{0}$ (right). The phase portraits show the $V$ -nullclines just before (light blue) and after the bifurcation (dark blue). The $V_{s}$ -nullclines are drawn as dashed lines. The stable fixed points are represented by filled circles, the saddle points by crosses. The possible trajectories are drawn in orange, with the reset points represented by squares.

6.3 Connection to Dynamic Input Conductances

The analysis method of this article has some similarities with that of the dynamic input conductances (DICs) introduced in Drion, Franci, Dethier, and Sepulchre (2015). Both methods group the contributions of different ion channels into a fast, slow, and ultraslow timescale. Nonetheless, a crucial difference between the methods is that the DICs are differential properties: they represent the change in current with an infinitesimal change in voltage. Instead, this article considers the ion current itself as a function of voltages in different timescales. This is important because infinitesimal voltage clamp steps are hard to perform experimentally. In contrast, the voltage-clamp simulations considered here could in principle be replaced by actual voltage-clamp experiments.

The other main difference between both methods is that DICs are defined around a steady state: all variables are at their steady-state value for the voltage $V$ at which the DIC is calculated. The function $I_{ion} (V, V_{s}, V_{u s})$ in this article, however, allows the voltages in each timescale to be different, thus capturing the transient behavior of the total ionic current. This is an important difference, because it allows our method to identify an integrate-and-fire neuron model, which is not possible from DICs.

While not identical, because of the way variables are divided over timescales, a property similar to DICs can be derived from

I_{ion} (V, V_{s}, V_{u s})

as follows:

g_{f} (V) = {\frac{\partial I_{ion} (V_{f}, V_{s}, V_{u s})}{\partial V_{f}}|}_{V_{f} = V_{s} = V_{u s} = V},

(6.3)

g_{s} (V) = {\frac{\partial I_{ion} (V_{f}, V_{s}, V_{u s})}{\partial V_{s}}|}_{V_{f} = V_{s} = V_{u s} = V},

(6.4)

g_{u s} (V) = {\frac{\partial I_{ion} (V_{f}, V_{s}, V_{u s})}{\partial V_{u s}}|}_{V_{f} = V_{s} = V_{u s} = V},

(6.5)

where

g_{f}

⁠,

g_{s}

⁠, and

g_{u s}

are the equivalent of the fast, slow, and ultraslow DICs, respectively. Since all voltages are taken as equal, it is clear that the DICs only reveal a subset of the information provided by

I_{ion} (V, V_{s}, V_{u s})

⁠.

6.4 Connection to I-V Curve Analysis

Ribar and Sepulchre (2019) propose a similar model of the total ionic current from I-V curves in different timescales. The main difference between the models is that the current functions in Ribar and Sepulchre (2019) are univariate. The current $I_{ion}$ is then a sum of functions of $V$ ⁠, $V_{s}$ ⁠, and $V_{u s}$ ⁠. The I-V curves are defined as the sum of all functions acting on their specific timescale and those faster—for example, $I_{f} (V) + I_{s} (V)$ for the slow I-V curve.

7 Conclusion

This article introduced a method to obtain an integrate-and-fire model from a conductance-based model by matching voltage-clamp and current-clamp responses.

The proposed multiscale integrate-and-fire model has a simple structure but retains a close connection to the physiology of the biophysical model.

The proposed method is applicable to arbitrary biophysical models under the key assumption that the kinetics of the gating variables can be grouped in a few distinct timescales.

Appendix: Methods

A.1 Software

Simulations were performed in Python with the SciPy library using the equations stated in the text. Differential equations were solved using SciPy's backward differentiation formula (BDF) method. All figures were drawn using the Python packages Matplotlib and Tikzplotlib, and/or the Latex packages PGF/TikZ and PGFPlots.

A.2 Parameters of the Biophysical Models

The parameters of the biophysical models used in this article take the original published values, except where a change of parameter is indicated. For the Connor-Stevens model (Connor et al., 1977), these are the following: $C = 1 μ F {cm}^{- 2}$ ⁠, ${\bar{g}}_{L} = 20 mS {cm}^{- 2}$ ⁠, ${\bar{g}}_{N a} = 120 mS {cm}^{- 2}$ ⁠, ${\bar{g}}_{K} = 20 mS {cm}^{- 2}$ ⁠, ${\bar{g}}_{A} = 47.7 mS {cm}^{- 2}$ ⁠, $V_{L} = - 17 mV$ ⁠, $V_{N a} = 55 mV$ ⁠, $V_{K} = - 72 mV$ ⁠, $V_{A} = - 75 mV$ ⁠. The parameters for the model of Drion et al. (2018) are the means of the parameters used in their network simulations: $C = 1 μ F {cm}^{- 2}$ ⁠, ${\bar{g}}_{L} = 0.055 mS {cm}^{- 2}$ ⁠, ${\bar{g}}_{N a} = 170 mS {cm}^{- 2}$ ⁠, ${\bar{g}}_{K, D} = 40 mS {cm}^{- 2}$ ⁠, ${\bar{g}}_{C a, T} = 0.55 mS {cm}^{- 2}$ ⁠, ${\bar{g}}_{K, C a} = 4 mS {cm}^{- 2}$ ⁠, ${\bar{g}}_{H} = 0.01 mS {cm}^{- 2}$ ⁠, $V_{L} = - 55 mV$ ⁠, $V_{N a} = 50 mV$ ⁠, $V_{K} = - 85 mV$ ⁠, $V_{C a} = 120 mV$ ⁠, $K_{D} = 170$ ⁠, $k_{1} = 0.1$ ⁠, $k_{2} = 0.01$ ⁠.

Table 1:

Parameters of the Integrate-and-Fire Models.

Figure	$C$	$τ_{f}$	$τ_{s}$	$τ_{u s}$	$V_{r}$	$V_{s, r}$	$Δ V_{u s}$
7 and 8 (center)	0.58	0.022	6.7		−40	−25
8 (left)	0.86	0.03	3		−45	−24
8 (right)	0.3	0.027	23		−35	−24
9	1.2	0.037	1.7	2.8	−40	−20	0
10–13	0.82	0.89	4.3	278	−45	7.5	1.7

Figure	$C$	$τ_{f}$	$τ_{s}$	$τ_{u s}$	$V_{r}$	$V_{s, r}$	$Δ V_{u s}$
7 and 8 (center)	0.58	0.022	6.7		−40	−25
8 (left)	0.86	0.03	3		−45	−24
8 (right)	0.3	0.027	23		−35	−24
9	1.2	0.037	1.7	2.8	−40	−20	0
10–13	0.82	0.89	4.3	278	−45	7.5	1.7

A.3 Parameters of the Integrate-and-Fire Models

All simulations and phase portraits of the integrate-and-fire models are based on the published biophysical model equations and parameters unless stated otherwise in the text. The parameters of the integrate-and-fire models were found using the method described in section 4, with SciPy's Trust Region Reflective algorithm (Branch, Coleman, & Li, 1999) for the least-squares optimization. The current clamp test data were generated using linearly decreasing currents for the Connor-Stevens model and a hyperpolarizing step current for the model of Drion et al. (2018). The data were sampled at 100 kHz and the transient removed. The obtained values are given in Table 1 for each figure. The units of $C$ ⁠, the time constants, and the reset voltages are, respectively, $μ$ F cm $^{- 2}$ ⁠, ms, and mV.

Note

¹

These conductances are not absolute but differential conductances, that is, slopes of the I-V curve. The negative conductance is sometimes referred to as “negative slope conductance” in the literature, for example, in the early observation of Stafstrom, Schwindt, and Crill (1982). See also Ceballos, Roque, and Leão (2017) for a recent review.

Acknowledgments

We thank Ilario Cirillo, Thiago Burghi, Luka Ribar, and Christian Grussler for helpful discussions. T.V.P received a fees scholarship from the Engineering and Physical Sciences Research Council (https://www.epsrc.ac.uk) under grant 1611337. Both T.V.P. and R.S. were supported by the European Research Council (https://erc.europa.eu) under the Advanced ERC Grant Agreement 670645. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

References

Badel

,

L.

,

Lefort

,

S.

,

Brette

,

R.

,

Petersen

,

C. C. H.

,

Gerstner

,

W.

, &

Richardson

,

M. J. E.

(

2008

).

Dynamic I-V curves are reliable predictors of naturalistic pyramidal-neuron voltage traces

.

Journal of Neurophysiology

,

99

(

2

),

656

–

666

.

Google Scholar

Branch

,

M.

,

Coleman

,

T.

, &

Li

,

Y.

(

1999

).

A subspace, interior, and conjugate gradient method for large-scale bound-constrained minimization problems

.

SIAM Journal on Scientific Computing

,

21

(

1

),

1

–

23

.

Google Scholar

Brette

,

R.

, &

Gerstner

,

W.

(

2005

).

Adaptive exponential integrate-and-fire model as an effective description of neuronal activity

.

Journal of Neurophysiology

,

94

(

5

),

3637

–

3642

.

Google Scholar

Brunel

,

N.

,

Hakim

,

V.

, &

Richardson

,

M. J.

(

2014

).

Single neuron dynamics and computation

.

Current Opinion in Neurobiology

,

25

,

149

–

155

.

Google Scholar

Burghi

,

T. B.

,

Schoukens

,

M.

, &

Sepulchre

,

R.

(

2019

).

Feedback for nonlinear system identification

. In

Proceedings of the 18th European Control Conference

(pp.

1344

–

1349

).

Piscataway, NJ

:

IEEE

.

Google Scholar

Ceballos

,

C. C.

,

Roque

,

A. C.

, &

Leão

,

R. M.

(

2017

).

The role of negative conductances in neuronal subthreshold properties and synaptic integration

.

Biophysical Reviews

,

9

(

5

),

827

–

834

.

Google Scholar

Connor

,

J. A.

, &

Stevens

,

C. F.

(

1971

).

Prediction of repetitive firing behaviour from voltage clamp data on an isolated neurone soma

.

Journal of Physiology

,

213

(

1

),

31

–

53

.

Google Scholar

Connor

,

J. A.

,

Walter

,

D.

, &

McKown

,

R.

(

1977

).

Neural repetitive firing: Modifications of the Hodgkin-Huxley axon suggested by experimental results from crustacean axons

.

Biophysical Journal

,

18

(

1

):

81

–

102

.

Google Scholar

Drion

,

G.

,

Dethier

,

J.

,

Franci

,

A.

, &

Sepulchre

,

R.

(

2018

).

Switchable slow cellular conductances determine robustness and tunability of network states

.

PLOS Computational Biology

,

14

(

4

), e1006125.

Google Scholar

Drion

,

G.

,

Franci

,

A.

,

Dethier

,

J.

, &

Sepulchre

,

R.

(

2015

).

Dynamic input conductances shape neuronal spiking

.

eNeuro

,

2

(

1

).

Google Scholar

Drion

,

G.

,

Franci

,

A.

,

Seutin

,

V.

, &

Sepulchre

,

R.

(

2012

).

A novel phase portrait for neuronal excitability

.

PLOS One

,

7

(

8

), e41806.

Google Scholar

Drion

,

G.

,

O'Leary

,

T.

, &

Marder

,

E.

(

2015

).

Ion channel degeneracy enables robust and tunable neuronal firing rates

. In

Proceedings of the National Academy of Sciences

,

112

(

38

),

E5361

–

E5370

.

Google Scholar

Enqvist

,

M.

, &

Ljung

,

L.

(

2005

).

Linear approximations of nonlinear FIR systems for separable input processes

.

Automatica

,

41

(

3

),

459

–

473

.

Google Scholar

Ermentrout

,

B.

(

1996

).

Type I membranes, phase resetting curves, and synchrony

.

Neural Computation

,

8

(

5

),

979

–

1001

.

Google Scholar

FitzHugh

,

R.

(

1961

).

Impulses and physiological states in theoretical models of nerve membrane

.

Biophysical Journal

,

1

(

6

),

445

–

466

.

Google Scholar

Fourcaud-Trocmé

,

N.

,

Hansel

,

D.

,

van Vreeswijk

,

C.

, &

Brunel

,

N.

(

2003

).

How spike generation mechanisms determine the neuronal response to fluctuating inputs

.

Journal of Neuroscience

,

23

(

37

),

11628

–

11640

.

Google Scholar

Franci

,

A.

,

Drion

,

G.

, &

Sepulchre

,

R.

(

2012

).

An organizing center in a planar model of neuronal excitability

.

SIAM Journal on Applied Dynamical Systems

,

11

(

4

),

1698

–

1722

.

Google Scholar

Franci

,

A.

,

Drion

,

G.

, &

Sepulchre

,

R.

(

2017

).

Robust and tunable bursting requires slow positive feedback

.

Journal of Neurophysiology

,

119

,

1222

–

1234

.

Google Scholar

Hill

,

A. V.

(

1936

).

Excitation and accommodation in nerve

. In

Proceedings of the Royal Society of London B: Biological Sciences

,

119

(

814

),

305

–

355

.

Google Scholar

Ho

,

B. L.

, &

Kalman

,

R. E.

(

1966

).

Editorial: Effective construction of linear state-variable models from input/output functions

.

Automatisierungstechnik

,

14

(

1–12

),

545

–

548

.

Google Scholar

Hodgkin

,

A. L.

, &

Huxley

,

A. F.

(

1952

).

A quantitative description of membrane current and its application to conduction and excitation in nerve

.

Journal of Physiology

,

117

(

4

),

500

–

544

.

Google Scholar

Huys

,

Q. J. M.

,

Ahrens

,

M. B.

, &

Paninski

,

L.

(

2006

).

Efficient estimation of detailed single-neuron models

.

Journal of Neurophysiology

,

96

(

2

),

872

–

890

.

Google Scholar

Izhikevich

,

E. M.

(

2003

).

Simple model of spiking neurons

.

IEEE Transactions on Neural Networks

,

14

(

6

),

1569

–

1572

.

Google Scholar

Jolivet

,

R.

,

Lewis

,

T. J.

, &

Gerstner

,

W.

(

2004

).

Generalized integrate-and-fire models of neuronal activity approximate spike trains of a detailed model to a high degree of accuracy

.

Journal of Neurophysiology

,

92

(

2

),

959

–

976

.

Google Scholar

Jolivet

,

R.

,

Rauch

,

A.

,

Lüscher

,

H.-R.

, &

Gerstner

,

W.

(

2006

).

Predicting spike timing of neocortical pyramidal neurons by simple threshold models

.

Journal of Computational Neuroscience

,

21

(

1

),

35

–

49

.

Google Scholar

Kepler

,

T. B.

,

Abbott

,

L. F.

, &

Marder

,

E.

(

1992

).

Reduction of conductance-based neuron models

.

Biological Cybernetics

,

66

(

5

),

381

–

387

.

Google Scholar

Knight

,

B. W.

(

1972

).

Dynamics of encoding in a population of neurons

.

Journal of General Physiology

,

59

(

6

),

734

–

766

.

Google Scholar

Kung

,

S. Y.

(

1978

). A new identification and model reduction algorithm via singular value decomposition. In

Proceedings of the 12th Asilomar Conference on Circuits, Systems, and Computers

(pp.

705

–

714

).

Piscataway, NJ

:

IEEE

.

Lapicque

,

L.

(

1907

).

Recherches quantitatives sur l'excitation électrique des nerfs traitée comme une polarisation

.

Journal de Physiologie et de Pathologie Générale

,

9

(

1

),

620

–

635

.

Google Scholar

Latham

,

P. E.

,

Richmond

,

B. J.

,

Nelson

,

P. G.

, &

Nirenberg

,

S.

(

2000

).

Intrinsic dynamics in neuronal networks. I. Theory

.

Journal of Neurophysiology

,

83

(

2

),

808

–

827

.

Google Scholar

Lepora

,

N. F.

,

Overton

,

P. G.

, &

Gurney

,

K.

(

2012

).

Efficient fitting of conductance-based model neurons from somatic current clamp

.

Journal of Computational Neuroscience

,

32

(

1

),

1

–

24

.

Google Scholar

Mauro

,

A.

,

Conti

,

F.

,

Dodge

,

F.

, &

Schor

,

R.

(

1970

).

Subthreshold behavior and phenomenological impedance of the squid giant axon

.

Journal of General Physiology

,

55

(

4

),

497

–

523

.

Google Scholar

Mensi

,

S.

,

Naud

,

R.

,

Pozzorini

,

C.

,

Avermann

,

M.

,

Petersen

,

C. C. H.

, &

Gerstner

,

W.

(

2012

).

Parameter extraction and classification of three cortical neuron types reveals two distinct adaptation mechanisms

.

Journal of Neurophysiology

,

107

(

6

),

1756

–

1775

.

Google Scholar

Mihalas

,

S.

, &

Niebur

,

E.

(

2008

).

A generalized linear integrate-and-fire neural model produces diverse spiking behaviors

.

Neural Computation

,

21

(

3

),

704

–

718

.

Google Scholar

Miller

,

D. N.

,

Hulett

,

J.

,

McLaughlin

,

J.

, &

de Callafon

,

R. A.

(

2013

).

Thermal dynamical identification of light-emitting diodes by step-based realization and convex optimization

.

IEEE Transactions on Components, Packaging and Manufacturing Technology

,

3

(

6

),

997

–

1007

.

Google Scholar

Morse

,

T. M.

,

Davison

,

A. P.

, &

Hines

,

M.

(

2001

).

Parameter space reduction in neuron model optimization through minimization of residual voltage clamp current

.

Society for Neuroscience Abstracts

,

27

.

Google Scholar

Nagumo

,

J.

,

Arimoto

,

S.

, &

Yoshizawa

,

S.

(

1962

).

An active pulse transmission line simulating nerve axon

.

Proceedings of the IRE

,

50

(

10

),

2061

–

2070

.

Google Scholar

Pozzorini

,

C.

,

Naud

,

R.

,

Mensi

,

S.

, &

Gerstner

,

W.

(

2013

).

Temporal whitening by power-law adaptation in neocortical neurons

.

Nature Neuroscience

,

16

(

7

),

942

–

948

.

Google Scholar

Ribar

,

L.

, &

Sepulchre

,

R.

(

2019

).

Neuromodulation of neuromorphic circuits

.

IEEE Transactions on Circuits and Systems I: Regular Papers

,

66

(

8

),

3028

–

3040

.

Google Scholar

Schoukens

,

M.

, &

Rolain

,

Y.

(

2012

).

Parametric identification of parallel Wiener systems

.

IEEE Transactions on Instrumentation and Measurement

,

61

(

10

),

2825

–

2832

.

Google Scholar

Stafstrom

,

C. E.

,

Schwindt

,

P. C.

, &

Crill

,

W. E.

(

1982

).

Negative slope conductance due to a persistent subthreshold sodium current in cat neocortical neurons in vitro

.

Brain Research

,

236

(

1

),

221

–

226

.

Google Scholar

Stein

,

R. B.

(

1965

).

A theoretical analysis of neuronal variability

.

Biophysical Journal

,

5

(

2

),

173

–

194

.

Google Scholar

Van Pottelbergh

,

T.

,

Drion

,

G.

, &

Sepulchre

,

R.

(

2018

).

Robust modulation of integrate-and-fire models

.

Neural Computation

,

30

(

4

),

987

–

1011

.

Google Scholar

Wehmeier

,

U.

,

Dong

,

D.

,

Koch

,

C.

, &

Van Essen

,

D.

(

1989

). Modeling the mammalian visual system. In

C.

Koch

&

I.

Segev

(Eds.),

Methods in neuronal modeling

(pp.

335

–

359

).

Cambridge, MA

:

MIT Press

.

Zeiger

,

H.

, &

McEwen

,

A.

(

1974

).

Approximate linear realizations of given dimension via Ho's algorithm

.

IEEE Transactions on Automatic Control

,

19

(

2

),

153

–

153

.

Google Scholar

2021

Massachusetts Institute of Technology

From Biophysical to Integrate-and-Fire Modeling

Abstract

1 Introduction

2 A Multiscale Integrate-and-Fire Model Structure

3 Identification of the Ion Current Function from a Biophysical Model

3.1 Two-Timescale Models

3.2 Three-Timescale Models

3.3 Precompensation of Voltages in Absence of Timescale Separation

3.4 Calcium Dynamics

4 Iterative Optimization of the Structural Parameters

4.1 Initialization of the Time Constants

4.2 Initialization of the Other Structural Parameters

4.3 Local Parameter Optimization

5 Integrate-and-Fire Modeling and Phase Portrait Analysis

5.1 Modulation of Excitability Type in the Connor-Stevens Model

5.1.1 A Two-Timescale Integrate-and-Fire Approximation of the Connor-Stevens Model

5.1.2 A Three-Timescale Integrate-and-Fire Approximation of the Connor-Stevens Model

5.2 Modulation between Spiking and Bursting Due to a Switchable Slow Negative Conductance

6 Discussion

6.1 Validity of the Method

6.2 Local Approximation by Multi-Quadratic Integrate-and-Fire Model

6.3 Connection to Dynamic Input Conductances

6.4 Connection to I-V Curve Analysis

7 Conclusion

Appendix: Methods

A.1 Software

A.2 Parameters of the Biophysical Models

A.3 Parameters of the Integrate-and-Fire Models

Note

Acknowledgments

References

Email alerts

Cited By

A product of The MIT Press

MIT Press Direct

Information

MIT Press

Contact Us

From Biophysical to Integrate-and-Fire Modeling

Abstract

1 Introduction

2 A Multiscale Integrate-and-Fire Model Structure

3 Identification of the Ion Current Function from a Biophysical Model

3.1 Two-Timescale Models

3.2 Three-Timescale Models

3.3 Precompensation of Voltages in Absence of Timescale Separation

3.4 Calcium Dynamics

4 Iterative Optimization of the Structural Parameters

4.1 Initialization of the Time Constants

4.2 Initialization of the Other Structural Parameters

4.3 Local Parameter Optimization

5 Integrate-and-Fire Modeling and Phase Portrait Analysis

5.1 Modulation of Excitability Type in the Connor-Stevens Model

5.1.1 A Two-Timescale Integrate-and-Fire Approximation of the Connor-Stevens Model

5.1.2 A Three-Timescale Integrate-and-Fire Approximation of the Connor-Stevens Model

5.2 Modulation between Spiking and Bursting Due to a Switchable Slow Negative Conductance

6 Discussion

6.1 Validity of the Method

6.2 Local Approximation by Multi-Quadratic Integrate-and-Fire Model

6.3 Connection to Dynamic Input Conductances

6.4 Connection to I-V Curve Analysis

7 Conclusion

Appendix: Methods

A.1 Software

A.2 Parameters of the Biophysical Models

A.3 Parameters of the Integrate-and-Fire Models

Note

Acknowledgments

References

Email alerts

Cited By

Related Articles

Related Book Chapters

A product of The MIT Press

MIT Press Direct

Information

MIT Press

Contact Us

This Feature Is Available To Subscribers Only