ARTICLE - Ricerca sull'intelligenza artificiale specializzata al MIT

ARTICLE

Communicated by Alain Destexhe

From Biophysical to Integrate-and-Fire Modeling

Tomas Van Pottelbergh
tmjv2@cam.ac.uk
Department of Engineering, University of Cambridge, Cambridge CB2 1PZ, U.K.

Guillaume Drion
gdrion@ulg.ac.be
Department of Electrical Engineering and Computer Science, University of Liège,
4000 Liège, Belgium

Rodolphe Sepulchre
r.sepulchre@eng.cam.ac.uk
Department of Engineering, University of Cambridge, Cambridge CB2 1PZ, U.K.

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 maxi-
mal 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 neuromodula-
zione: the transition between type I and type II excitability and the transi-
tion between spiking and bursting.

1 introduzione

Integrate-and-fire models have a long history, dating back to the begin-
ning of the twentieth century (Lapicque, 1907). Because of their simplicity,
they have long served as phenomenological models for action poten-
tial generation in neurons. Over the past few decades, the original leaky
integrate-and-fire model (Hill, 1936; Stein, 1965; Knight, 1972) has gradu-
ally 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 dy-
namics 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;

Calcolo neurale 33, 563–589 (2021)
https://doi.org/10.1162/neco_a_01353

D
o
w
N
o
UN
D
e
D

F
R
o
M
H

T
T

:
/
/

D
io
R
e
C
T
.

io
T
.

e
D
tu
N
e
C
o
UN
R
T
io
C
e
–
P
D

F
/

3
3
3
5
6
3
1
8
8
9
4
0
5
N
e
C
o
_
UN
_
0
1
3
5
3
P
D

B
sì
G
tu
e
S
T

o
N
0
7
S
e
P
e
M
B
e
R
2
0
2
3

564

T. Van Pottelbergh, G. Drion, and R. Sepulchre

Brette & Gerstner, 2005; Drion, Franci, Seutin, & Sepulchre, 2012). Basso-
dimensional integrate-and-fire models have also been shown to be suffi-
cient 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 bio-
physical model (Van Pottelbergh, Drion, & Sepulchre, 2018). This is a se-
vere obstacle when studying the role of cellular neuromodulation in a
large network (Drion et al., 2018). While intrinsic neuromodulation is stud-
ied 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 mo-
tivation 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 ex-
perimental spiking data. Invece, 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 re-
quires 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 nega-
tive conductance1 in each timescale. This model has been quite successful
at capturing important modulation properties in a robust and qualitative
maniera.

We introduce a multiscale integrate-and-fire model in section 2 and high-
light 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. Sezione 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

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

D
o
w
N
o
UN
D
e
D

F
R
o
M
H

T
T

:
/
/

D
io
R
e
C
T
.

io
T
.

e
D
tu
N
e
C
o
UN
R
T
io
C
e
–
P
D

F
/

3
3
3
5
6
3
1
8
8
9
4
0
5
N
e
C
o
_
UN
_
0
1
3
5
3
P
D

B
sì
G
tu
e
S
T

o
N
0
7
S
e
P
e
M
B
e
R
2
0
2
3

From Biophysical to Integrate-and-Fire Modeling

565

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

− Iion(V, Vs, Vus, . . .),

C ˙V = Iapp
τs ˙Vs = V − Vs,
τus ˙Vus = V − Vus.
. . .

(2.1)

(2.2)

(2.3)

The voltage equation 2.1 is the classical Kirchhoff relationship of a bio-
physical model: the current Iion is the total ionic current intrinsic to the cel-
lular membrane composition. Its voltage dependence is modeled through
the voltage variable V and lagged variables (Vs, Vus, . . .) that model the volt-
age filtered in distinct timescales (slow, ultraslow, . . .). Those variables dif-
fer from the gating variables of traditional conductance-based models, Ma
closely relate to the equivalent potentials originally defined in Kepler, Ab-
bott, 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 Iion.

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

When a reset is added to the multiscale model, as in equations 2.4 A
2.6, 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 Pottel-
bergh et al., 2018):

− Iion(V, Vs, Vus, . . .)

C ˙V = Iapp
τs ˙Vs = V − Vs
τus ˙Vus = V − Vus
. . .

if V ≥ Vmax :
V ← Vr
Vs ← Vs,R
Vus ← Vus + (cid:3)Vus

(2.4)

(2.5)

(2.6)

D
o
w
N
o
UN
D
e
D

F
R
o
M
H

T
T

:
/
/

D
io
R
e
C
T
.

io
T
.

e
D
tu
N
e
C
o
UN
R
T
io
C
e
–
P
D

F
/

3
3
3
5
6
3
1
8
8
9
4
0
5
N
e
C
o
_
UN
_
0
1
3
5
3
P
D

B
sì
G
tu
e
S
T

o
N
0
7
S
e
P
e
M
B
e
R
2
0
2
3

566

T. Van Pottelbergh, G. Drion, and R. Sepulchre

Figura 1: Block diagram representation of the multiscale model 2.1 A 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 Iion (section 3). Questo
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, Tuttavia, that there is a direct correspon-
dence between model 2.1 A 2.3 and model 2.4 A 2.6 provided that model 2.1
A 2.3 generates spikes. Each spike in the continuous-time model is replaced
by a reset mechanism in model 2.4 A 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. Invece, the nu-
merical 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 (other-
wise 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 differ-
ential equation (see Van Pottelbergh et al., 2018, for details). This constraint
is key to retaining a direct correspondence to the physiology of a biophysi-
cal model.

An important property of the integrate-and-fire model 2.4 A 2.6 is that
the nonlinear function Iion only needs to be identified in the subthreshold
voltage range V < Vmax. 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 l D o w n o a d e d f r o m h t t p : / / d i r e c t . m i t . / e d u n e c o a r t i c e - p d / l f / / / / 3 3 3 5 6 3 1 8 8 9 4 0 5 n e c o _ a _ 0 1 3 5 3 p d . / f b y g u e s t t o n 0 7 S e p e m b e r 2 0 2 3 From Biophysical to Integrate-and-Fire Modeling 567 structural parameters of the model (the time constants and the reset pa- rameters) and the single scalar nonlinear function Iion. In this section, we determine Iion for a given set of structural parameters. In the next section, we address the determination of the structural parameters and how the de- sign 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 pro- pose to identify the ionic current Iion 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 ex- periments in the first place. It is also computationally tractable because a voltage-clamp step response can be calculated in closed form in both a bio- physical model and in the integrate-and-fire model. In a biophysical model, each gating variable obeys a differential equation of the form τxi (V ) ˙xi = xi,∞(V ) − xi , (3.1) which becomes linear for a fixed value of V. The solution after a step from V0 to Vstep at t = 0 is described by xi(t) = xi,∞(V0) + [xi,∞(Vstep) − xi,∞(V0)] · (cid:2) 1 − e −t/τxi (Vstep ) (cid:3) , (3.2) assuming the membrane potential is at equilibrium at t = 0. The response of the total ionic current Iion in the biophysical model to a voltage-clamp step is then obtained by direct substitution of V (t) = Vstep and xi(t) (given by equation 3.2) in the expression for the total ionic cur- rent. The same calculation holds in the integrate-and-fire model 2.4 to 2.6 to obtain an expression of Iion(V (t), Vs(t), Vus(t), . . .). The details of this simple idea are provided in the next sections. 3.1 Two-Timescale Models. We start by describing the method for bio- physical models that can be described well by a fast and slow timescale, which are well separated (τs (cid:5) τ f ). We consider a simple voltage-clamp step experiment as in Figure 2, stepping from the initial voltage V0 to the final voltage Vstep. Assuming τ f is a good approximation of the time con- stants of the fast gating variables, they will have approximately reached their steady-state value xi,∞(Vstep) 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 xi,∞(V0). Similarly, in the multiscale integrate-and-fire model, the fast voltage f . V f will have approximately reached its steady-state value Vstep at t = 3τ l D o w n o a d e d f r o m h t t p : / / d i r e c t . m i t . / e d u n e c o a r t i c e - p d / l f / / / / 3 3 3 5 6 3 1 8 8 9 4 0 5 n e c o _ a _ 0 1 3 5 3 p d . / f b y g u e s t t o n 0 7 S e p e m b e r 2 0 2 3 568 T. Van Pottelbergh, G. Drion, and R. Sepulchre Figure 2: 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 V0 to Vstep. Assuming V = V f , the value of I , V0) of the multiscale model. at 3τ f will give a good approximation of Iion(Vstep Under the assumption of τs (cid:5) τ f , the slow voltage Vs will still be approxi- mately equal to V0. After merging the variables V and V f in the multiscale model, the following observation can be made for the voltage-clamp exper- iment in Figure 2: I(3τ f ) ≈ Iion(Vstep , V0). (3.3) The interpretation of this expression is that the value of the ion current function of the multiscale model evaluated at V = Vstep and Vs = V0 can be read at t = 3τ f from a voltage-clamp step experiment from V0 to Vstep. Given that V0 and Vstep can be chosen arbitrarily, this gives a simple method to evaluate Iion on the biophysical model and use its value in the function Iion(V, Vs) 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 3.2, resulting in the expression = xi,∞(V0) + [xi,∞(Vstep) − xi,∞(V0)] · [1 − e −3τ f /τxi (Vstep )] xi (3.4) for each gating variable. Substituting this expression in the ion current , V0) of the integrate- equation of the biophysical model then makes Iion(Vstep 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 separa- tion, the formula becomes a simple substitution of the fast and slow gating l D o w n o a d e d f r o m h t t p : / / d i r e c t . m i t . / e d u n e c o a r t i c e - p d / l f / / / / 3 3 3 5 6 3 1 8 8 9 4 0 5 n e c o _ a _ 0 1 3 5 3 p d . / f b y g u e s t t o n 0 7 S e p e m b e r 2 0 2 3 From Biophysical to Integrate-and-Fire Modeling 569 Figure 3: Voltage-clamp step experiment on a three-timescale model—the model of Drion et al. (2018). The voltage is stepped from V0 to Vstep,1 and then to Vstep,2 after 3τ f will give a good approximation of Iion(Vstep,2 s. Assuming V = V f , the value of I at 3τ , V0) of the multiscale model. , Vstep,1 + 3τ s variables by their steady-state values at Vstep and V0, respectively: xi lim τxi →0 →+∞ xi limτxi = xi,∞(Vstep), = xi,∞(V0). (3.5) (3.6) 3.2 Three-Timescale Models. We apply the same idea to three- timescale models to find Iion(V, Vs, Vus) 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 , V0). This can be resolved voltage, and therefore only evaluate Iion(Vstep by devising a slightly more complex voltage clamp step experiment (see Figure 3): starting at V0, stepping to Vstep,1 and, finally, Vstep,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 xi,∞(Vstep,2). Sim- ilarly the slow gating variables will approximately be at the steady-state value xi,∞(Vstep,1), while the ultraslow gating variables will still be close to xi,∞(V0). The value of the current at time 3τs + 3τ f will therefore be a good , V0) of the multiscale model, assuming approximation of Iion(Vstep,2 V = V f . , Vstep,1 , V0 The necessary substitution for the gating variables in the ion current equation is given by l D o w n o a d e d f r o m h t t p : / / d i r e c t . m i t . / e d u n e c o a r t i c e - p d / l f / / / / 3 3 3 5 6 3 1 8 8 9 4 0 5 n e c o _ a _ 0 1 3 5 3 p d . / f b y g u e s t t o n 0 7 S e p e m b e r 2 0 2 3 + [xi,∞(Vstep,2) − ¯xi] · [1 − e = ¯xi = xi,∞(V0) + [xi,∞(Vstep,1) − xi,∞(V0)] · [1 − e /τxi (Vstep,2 )], −3τ f −3τs/τxi (Vstep,1 )], xi ¯xi (3.7) (3.8) 570 T. Van Pottelbergh, G. Drion, and R. Sepulchre 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: xi lim τxi →0 →+∞ xi limτxi = xi,∞(Vstep,2), = xi,∞(V0), while xi ≈ xi,∞(Vstep,1) if τxi ≈ τs and τs (cid:5) τ f . (3.9) (3.10) 3.3 Precompensation of Voltages in Absence of Timescale Separa- tion. The method in the previous sections assumes that each gating vari- able 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 assump- tion 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 de- scribed 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τ + 3τs for three- timescale models). This is illustrated in Figure 4. f (or t = 3τ f In the case of two timescales, only the initial voltage of the voltage-clamp f can be computed given step needs to be changed. The value of Vs at t = 3τ the initial voltage V0 and the step voltage Vstep: Vs(3τ f ) = V0 + (Vstep − V0) (cid:4) 1 − e −3τ /τs f (cid:5) . (3.11) For this voltage to reach the desired value V ∗ sion can be solved for V0: s at t = 3τ f , the above expres- = V0 (cid:6) ∗ V s − Vstep (cid:4) 1 − e −3τ (cid:5)(cid:7) /τs f −3τ /e /τs . f (3.12) Making this additional substitution in 3.4 results in the substitution expres- , V ∗ sion to obtain the value for Iion(Vstep s ). The same idea can be applied to three-timescale models. Vstep,1 is calcu- lated in the same way as V0 in the two-timescale case: Vstep,1 = (cid:6) ∗ V s − Vstep,2 (cid:4) 1 − e −3τ /τs f (cid:5)(cid:7) −3τ /e /τs . f (3.13) l D o w n o a d e d f r o m h t t p : / / d i r e c t . m i t . / e d u n e c o a r t i c e - p d / l f / / / / 3 3 3 5 6 3 1 8 8 9 4 0 5 n e c o _ a _ 0 1 3 5 3 p d . / f b y g u e s t t o n 0 7 S e p e m b e r 2 0 2 3 From Biophysical to Integrate-and-Fire Modeling 571 l D o w n o a d e d f r o m h t t p : / / d i r e c t . m i t . Figure 4: Illustration of the precompensation of step voltages to reach the de- sired values of the slow voltages at t = 3τ s. (Top) The choice f = V ∗ of V0 makes Vs f for two-timescale models. (Bottom) The choice of V0, Vstep,1 and Vstep,2 leads to Vs s for three- timescale models. us at t = 3τ s at t = 3τ f or t = 3τ s and Vus = V ∗ = V ∗ + 3τ + 3τ f This can be done because Vs can be assumed to be at steady state at t = 3τs. V0 can then be computed using Vstep,1. For simplicity, it is assumed that τus (cid:5) τ f , giving = V0 (cid:6) V ∗ us − Vstep,1 (cid:4) 1 − e −3τs/τus (cid:5)(cid:7) −3τs/τus . /e (3.14) If τus (cid:5) τ cated: f does not hold, the expression becomes slightly more compli- = V0 ¯V = (cid:4) (cid:6) ¯V − Vstep,1 (cid:6) − Vstep,2 V ∗ us 1 − e (cid:4) −3τs/τus −3τs/τus , (cid:5)(cid:7) /e (cid:5)(cid:7) 1 − e −3τ /τus f −3τ /τus . f /e (3.15) (3.16) / e d u n e c o a r t i c e - p d / l f / / / / 3 3 3 5 6 3 1 8 8 9 4 0 5 n e c o _ a _ 0 1 3 5 3 p d . / f b y g u e s t t o n 0 7 S e p e m b e r 2 0 2 3 After making these substitutions in equations 3.7 and 3.8, the substitu- tion expressions can be used to obtain the value for Iion(Vstep,2 , V ∗ s , V ∗ us). 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 572 T. Van Pottelbergh, G. Drion, and R. Sepulchre equations of the type ˙ τ [Ca2+] τxCa,i (V ) ˙xCa,i [Ca2+] = [Ca2+]∞(V, xCa,i = xCa,i,∞(V ) − xCa,i . . . , . . .) − [Ca2+], (3.17) (3.18) where xCa,i are the gating variables of the calcium channels. This will gener- ally not result in a simple substitution expression. However, under the as- sumption 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. [Ca2+](Vs) and [Ca2+](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 Iion 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 con- stants, 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 Iapp α 1 ≈ C ˙V + α 0 ∂Iion ∂V , α 2 = 2V f + α 3Vs + α 4Vus + α = 1V + α ∂Iion ∂V f , . . . where the coefficients α i depend on the membrane potential around which the experiment is performed. The contribution of C ˙V is easily removed as it only produces a sharp pulse at the time of the step. Repeating this exper- iment at multiple subthreshold voltages gives a more robust estimation, as some of the coefficients α i might be small around a single voltage. l D o w n o a d e d f r o m h t t p : / / d i r e c t . m i t . / e d u n e c o a r t i c e - p d / l f / / / / 3 3 3 5 6 3 1 8 8 9 4 0 5 n e c o _ a _ 0 1 3 5 3 p d . / f b y g u e s t t o n 0 7 S e p e m b e r 2 0 2 3 From Biophysical to Integrate-and-Fire Modeling 573 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 τ , τs, τus will hence be the negative inverse of the identified poles. f 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) (En- qvist & Ljung, 2005; Schoukens & Rolain, 2012). While this approach can result in more accurate estimates, it requires more data and the gaussian in- put signals used are less suited to the voltage-clamp setting. Furthermore, the time constants only have to be accurate enough to make the optimiza- tion 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 Iion, the only parameters left to identify in model 2.4 to 2.6 are the membrane capacitance C; the reset parameters Vr, Vs,r, and (cid:3)Vus; and the cutoff voltage Vmax. Those parameters are easily estimated from current clamp simulations of the biophysical model. Vmax 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 Vmax lies suffi- ciently 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 Vr is always taken to be equal to Vmax. The reset values of the slower voltages are parameters to be optimized but are physiologi- cally constrained. A reasonable initialization is to set Vs,r sufficiently (e.g., 20 mV) above Vr. The parameter (cid:3)Vus is constrained to be positive, since Vus 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 represen- tative current clamp data (see Figure 5, top). 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; l D o w n o a d e d f r o m h t t p : / / d i r e c t . m i t . / e d u n e c o a r t i c e - p d / l f / / / / 3 3 3 5 6 3 1 8 8 9 4 0 5 n e c o _ a _ 0 1 3 5 3 p d . / f b y g u e s t t o n 0 7 S e p e m b e r 2 0 2 3 574 T. Van Pottelbergh, G. Drion, and R. Sepulchre = Figure 5: Test data generated using the Connor-Stevens model with ¯gA 0 mS cm−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 lin-
ear low-pass filter and reset to Vs,r after every spike, and Iion is evaluated using
these V and Vs.

Lepora, Overton, & Gurney, 2012). The residual current for our model is
defined as

Ires

= C ˙V − Iapp

+ Iion(V, Vs, Vus)

(4.1)

and is zero when the membrane current equation is satisfied. The test data
provide V and Iapp, while the derivative ˙V can be obtained by numeric dif-
ferentiation. The term Iion depends on the chosen time constants, and to
evaluate it, the voltage trace is filtered by the corresponding first-order lin-
ear low-pass filters. The action potentials themselves are eliminated by re-
moving the data for which V > Vmax, as the model is only optimized in the
subthreshold regime. Invece, the different voltages are reset after every
spike using their respective reset parameters, as shown in Figure 5 (bot-
tom). The voltage traces of Figure 5 (bottom) can then be used to evaluate
Iion(V, Vs, Vus) 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 param-
eters in the optimization are C and the reset values Vs,r and (cid:3)Vus. The ad-
vantage 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 contrasto, a cost

D
o
w
N
o
UN
D
e
D

F
R
o
M
H

T
T

:
/
/

D
io
R
e
C
T
.

io
T
.

e
D
tu
N
e
C
o
UN
R
T
io
C
e
–
P
D

F
/

3
3
3
5
6
3
1
8
8
9
4
0
5
N
e
C
o
_
UN
_
0
1
3
5
3
P
D

B
sì
G
tu
e
S
T

o
N
0
7
S
e
P
e
M
B
e
R
2
0
2
3

From Biophysical to Integrate-and-Fire Modeling

575

s for the Connor-Stevens model with ¯gA

Figura 6: Contour plots of two least squares cost functions as a function of τ
F
= 0 mS cm−2. (Left) Cost function
and τ
based on the residual current, as used in this article. (Center and right: Cost
= V − V ∗, where V is the voltage
function based on the residual voltage (Vres
of the test data and V ∗ is the voltage of the simulation of the integrate-and-fire
modello).

function based on the residual voltage has many local minima and is much
more irregular.

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 compu-
tational. We now show that it is also amenable to phase portrait analysis,
which provides mathematical insight into the initial biophysical model, Rif-
gardless 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 Iion(V, Vs) =
Iapp and the Vs-nullcline Vs = V. In other words, the level curves of the iden-
tified ion current Iion 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 Vus as a bifurcation pa-
rameter and by studying the family of phase portraits parameterized by Vus.
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

D
o
w
N
o
UN
D
e
D

F
R
o
M
H

T
T

:
/
/

D
io
R
e
C
T
.

io
T
.

e
D
tu
N
e
C
o
UN
R
T
io
C
e
–
P
D

F
/

3
3
3
5
6
3
1
8
8
9
4
0
5
N
e
C
o
_
UN
_
0
1
3
5
3
P
D

B
sì
G
tu
e
S
T

o
N
0
7
S
e
P
e
M
B
e
R
2
0
2
3

576

T. Van Pottelbergh, G. Drion, and R. Sepulchre

Figura 7: Phase portraits of the integrate-and-fire model obtained from the
Connor-Stevens model with standard parameters (Connor et al., 1977). The V-
and Vs-nullclines are drawn as full and dashed lines, rispettivamente. (Left) IL
V-nullclines are drawn for different values of Iapp: low (light blue), Iapp 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 Iapp), together
with the corresponding trace of the limit cycle oscillation in orange.

Connor-Stevens model, whose behavior can be easily modulated by a
change of maximal conductance, and the model by Drion et al. (2018), ex-
hibiting a switchable slow negative conductance.

5.1 Modulation of Excitability Type in the Connor-Stevens Model.
The Connor-Stevens model (Connor & Stevens, 1971; Connor et al., 1977) È
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 IA. One of its characteristic features is type I excitabil-
ità: the spiking frequency approaches zero when the applied current ap-
proaches the rheobase. This is in contrast to the type II excitability of the
Hodgkin-Huxley model, whose spiking frequency makes a jump as the ap-
plied 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 Ke-
pler 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 Ke-
pler 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 pro-
cedure are τ
= 0.022 ms and τs = 6.7 ms. They reflect the timescale sepa-
ration of the gating variables, which span a range of 0.03 ms to 3 ms.

The phase portrait of the two-timescale continuous-time model is dis-
app, IL

played in Figure 7 (left). For a specific value of the applied current I∗

D
o
w
N
o
UN
D
e
D

F
R
o
M
H

T
T

:
/
/

D
io
R
e
C
T
.

io
T
.

e
D
tu
N
e
C
o
UN
R
T
io
C
e
–
P
D

F
/

3
3
3
5
6
3
1
8
8
9
4
0
5
N
e
C
o
_
UN
_
0
1
3
5
3
P
D

B
sì
G
tu
e
S
T

o
N
0
7
S
e
P
e
M
B
e
R
2
0
2
3

From Biophysical to Integrate-and-Fire Modeling

577

D
o
w
N
o
UN
D
e
D

F
R
o
M
H

T
T

:
/
/

D
io
R
e
C
T
.

io
T
.

e
D
tu
N
e
C
o
UN
R
T
io
C
e
–
P
D

F
/

3
3
3
5
6
3
1
8
8
9
4
0
5
N
e
C
o
_
UN
_
0
1
3
5
3
P
D

B
sì
G
tu
e
S
T

o
N
0
7
S
e
P
e
M
B
e
R
2
0
2
3

Figura 8: 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 (O-
ange) for different values of ¯gA. The V- and Vs-nullclines are drawn as full and
dashed lines, rispettivamente, the trajectories after losing stability in orange. (Left)
¯gA

= 47.7 mS cm−2. (Right) ¯gA

= 0 mS cm−2. (Center) ¯gA

= 200 mS cm−2.

Vs-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) bifur-
catione. 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
A 2.6 in orange. Since the model has a reset mechanism to replace the ac-
tion potential generation, the periodic spiking occurs due to a hybrid limit
cycle.

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

= 0 mS cm

The excitability type of the model can be changed easily by modulating
−2, the model is similar to
the maximal conductance ¯gA. For ¯gA
the Hodgkin-Huxley model, which exhibits type II excitability. We again
obtain a two-timescale integrate-and-fire approximation of this model by
finding Iion and the structural parameters using the new value for ¯gA. 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

578

T. Van Pottelbergh, G. Drion, and R. Sepulchre

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 Iapp than the original model, but both f-I curves remain close to each other
afterward.

D'altra parte, increasing ¯gA above its nominal value of 47.7 mS cm

−2
−2 results in a bistable phase portrait (Guarda la figura 8, bottom
A 200 mS cm
right): for a specific range of Iapp, a stable hybrid limit cycle on the up-
per 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. Questo era
called type II* excitability in Drion, O’Leary et al. (2015), which is simi-
lar to type II excitability, but has a hysteretic f-I curve. Figura 8 (top right)
shows that the integrate-and-fire model captures this hysteresis. Although
the model starts spiking at a lower value of Iapp, 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 re-
duction of the Connor-Stevens model neuron model is useful to obtain a
qualitative approximation of the behavior of the original model. The mod-
ulation of excitability type can be predicted from the phase portraits, E
it is possible to construct an integrate-and-fire model that replicates this
modulation. Tuttavia, the resulting integrate-and-fire model lacks a quan-
titative 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. IL
three timescales are not strongly separated, Tuttavia, 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 τ
= 0.037 ms, τs = 1.7 ms, and τus = 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 vari-
able of the current IA, which are different from those of the other gating
variables. Figura 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 (Guarda la figura 9, left) reveals very similar voltage traces.

−2, we did not find a similar improvement for ¯gA

While the f-I curves match well for the three-timescale model with ¯gA
= 200 mS cm

47.7 mS cm
(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.

=
−2

D
o
w
N
o
UN
D
e
D

F
R
o
M
H

T
T

:
/
/

D
io
R
e
C
T
.

io
T
.

e
D
tu
N
e
C
o
UN
R
T
io
C
e
–
P
D

F
/

3
3
3
5
6
3
1
8
8
9
4
0
5
N
e
C
o
_
UN
_
0
1
3
5
3
P
D

B
sì
G
tu
e
S
T

o
N
0
7
S
e
P
e
M
B
e
R
2
0
2
3

From Biophysical to Integrate-and-Fire Modeling

579

Figura 9: (Left) Voltage traces of the original Connor-Stevens model (top) E
the three-timescale integrate-and-fire model (bottom) for Iapp linearly increasing
from 8 A 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 Switch-
able 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 chan-
nels endow the model with slow regenerativity by having an activation that
is slower than the sodium channel activation. Inoltre, 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 exter-
nal current: it is switched on only by hyperpolarization. This current was
hypothesized in Drion et al. (2018) as a critical mechanism for the neuro-
modulation of network states. Figura 10 (left) shows how the switchable
negative conductance can be observed from a voltage clamp step experi-
ment. 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).

We use the method of section 3.2 to construct an integrate-and-fire ap-
proximation of the biophysical model of Drion et al. (2018). The calcium

D
o
w
N
o
UN
D
e
D

F
R
o
M
H

T
T

:
/
/

D
io
R
e
C
T
.

io
T
.

e
D
tu
N
e
C
o
UN
R
T
io
C
e
–
P
D

F
/

3
3
3
5
6
3
1
8
8
9
4
0
5
N
e
C
o
_
UN
_
0
1
3
5
3
P
D

B
sì
G
tu
e
S
T

o
N
0
7
S
e
P
e
M
B
e
R
2
0
2
3

580

T. Van Pottelbergh, G. Drion, and R. Sepulchre

Figura 10: Comparison of the response to a voltage-clamp step from a hyper-
polarized (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.

=
Figura 11: (Left) Voltage trace of the model of Drion et al. (2018) for Iapp
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 Vs-
nullclines are drawn as full and dashed lines, rispettivamente, the stable (hybrid)
limit cycle in orange.

dynamics are treated as explained in section 3.4. Its behavior can be stud-
ied using three timescales. To visualize the obtained reduced model, phase
portraits in the V-Vs space are shown for different values of Vus. Assuming
the ultraslow timescale is much slower than the slow timescale, the behav-
ior of the model can be analyzed by looking at the fast-slow system in these
phase portraits.
For Iapp

−2, the original model shows regular spiking, con
a lower spiking frequency than during the bursting (Guarda la figura 11, left).
As this is a two-timescale behavior, Vus can be approximated by a con-
stant value. Figura 11 (right) shows the phase portrait corresponding to this

= 0 μA cm

D
o
w
N
o
UN
D
e
D

F
R
o
M
H

T
T

:
/
/

D
io
R
e
C
T
.

io
T
.

e
D
tu
N
e
C
o
UN
R
T
io
C
e
–
P
D

F
/

3
3
3
5
6
3
1
8
8
9
4
0
5
N
e
C
o
_
UN
_
0
1
3
5
3
P
D

B
sì
G
tu
e
S
T

o
N
0
7
S
e
P
e
M
B
e
R
2
0
2
3

From Biophysical to Integrate-and-Fire Modeling

581

Figura 12: (Top) Voltage trace of a burst in the model of Drion et al. (2018) con
= −1.6 μA cm−2. (Bottom) Portions of fast-slow phase portraits of the three-
Iapp
dimensional integrate-and-fire approximation for different values of Vus. IL
blue curve is the V-nullcline, the dashed curve is the Vs-nullcline, and each or-
ange curve is a trajectory of the phase portrait that captures a piece of the burst-
ing attractor in the three-dimensional model.

situation. There is no bistability between the resting and spiking state, E
the phase portrait is the standard phase portrait of a spiking model.

= −1.6 μA cm

When the current is lowered to Iapp

−2, the model bursts
(Guarda la figura 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
(Guarda la figura 12, bottom left). The subsequent spiking, on the limit cycle on the
upper branch of the V-nullcline, causes Vus 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 (Guarda la figura 12, bottom center). As Vus
increases even further, the limit cycle is lost in a saddle-homoclinic bifur-
catione, moving the trajectory to the stable fixed point and thus terminating
the burst (Guarda la figura 12, bottom right).

Figura 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 re-
sponse 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.

D
o
w
N
o
UN
D
e
D

F
R
o
M
H

T
T

:
/
/

D
io
R
e
C
T
.

io
T
.

e
D
tu
N
e
C
o
UN
R
T
io
C
e
–
P
D

F
/

3
3
3
5
6
3
1
8
8
9
4
0
5
N
e
C
o
_
UN
_
0
1
3
5
3
P
D

B
sì
G
tu
e
S
T

o
N
0
7
S
e
P
e
M
B
e
R
2
0
2
3

582

T. Van Pottelbergh, G. Drion, and R. Sepulchre

Figura 13: Voltage trace of the model of Drion et al. (2018) (top) and its three-
dimensional integrate-and-fire approximation (bottom) for Iapp stepping from 0
to −1.6 μA cm−2 at t = 500 ms.

Again, the obtained reduction can be used to construct an integrate-
and-fire model. Figura 13 shows the voltage trace before and after a
hyperpolarizing step current for the original model (top) and the reduced
modello (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.

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 biophys-
ical 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, Tuttavia, 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 be-
low), 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 vari-
ables in a biophysical model into potentials or voltages. Every variable xi

D
o
w
N
o
UN
D
e
D

F
R
o
M
H

T
T

:
/
/

D
io
R
e
C
T
.

io
T
.

e
D
tu
N
e
C
o
UN
R
T
io
C
e
–
P
D

F
/

3
3
3
5
6
3
1
8
8
9
4
0
5
N
e
C
o
_
UN
_
0
1
3
5
3
P
D

B
sì
G
tu
e
S
T

o
N
0
7
S
e
P
e
M
B
e
R
2
0
2
3

From Biophysical to Integrate-and-Fire Modeling

583

= x

is simply replaced by the associated voltage of the steady-state function:
−1
Vxi
io,∞(xi). In doing so, only the description of the system is changed,
not its input-output dynamics. The new form, Tuttavia, 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 equiv-
alent potentials with similar dynamics and replacing them by a weighted
average of their group. The weights are found by optimizing the local ap-
proximation of the model. For the method to work, the dynamics of the
equivalent potentials should fall into groups with similar dynamics. IL
method of this article provides a simpler alternative based on voltage clamp
esperimenti, 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 A 2.3, almeno
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 care-
fully validated afterward.

6.2 Local Approximation by Multi-Quadratic Integrate-and-Fire
Model. The multi-quadratic integrate-and-fire (MQIF) modello (Drion et al.,
2012; Van Pottelbergh et al., 2018) is a special case of the multiscale integrate-
and-fire model 2.4 A 2.6, in which Iion is a sum of quadratic functions in V,
Vs, and Vus:

Iion(V, Vs, Vus) = ¯g f (V − V 0)2 + ¯gs(Vs − V 0

S )2 + ¯gus(Vus − V 0

us)2.

(6.1)

The quadratics in V and Vs provide a normal form of the transcritical sin-
gularity in the fast subsystem, organizing the rest-spike bistability. IL
quadratic in Vus models the ultraslow feedback necessary for bursting. IL
relative positions of the parameters V 0 and V 0
s 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 pos-
sible 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
modello 2.4 A 2.6 by finding the point (V 0, V 0

S ) for which

∂Iion(V, Vs)
∂V

∂Iion(V, Vs)
∂Vs

= 0,

(6.2)

together with the condition that the determinant of the Hessian at (V 0, V 0
S )
should be negative (or zero, requiring further investigation). Iion(V 0, V 0
S ) È
then the current offset necessary to have the transcritical bifurcation occur

D
o
w
N
o
UN
D
e
D

F
R
o
M
H

T
T

:
/
/

D
io
R
e
C
T
.

io
T
.

e
D
tu
N
e
C
o
UN
R
T
io
C
e
–
P
D

F
/

3
3
3
5
6
3
1
8
8
9
4
0
5
N
e
C
o
_
UN
_
0
1
3
5
3
P
D

B
sì
G
tu
e
S
T

o
N
0
7
S
e
P
e
M
B
e
R
2
0
2
3

584

T. Van Pottelbergh, G. Drion, and R. Sepulchre

< V 0 (left), type I for V 0 s Figure 14: Modulation of the excitability type in the MQIF model. Chang- ing V 0 in the MQIF model results in different types of excitability: type II for s = V 0 (center) and type II* (Drion, O’Leary et al., V 0 s > V 0 (right). The phase portraits show the V-nullclines just before
2015) for V 0
S
(light blue) and after the bifurcation (dark blue). The Vs-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.

at the same value of Iapp. The values of ¯g f and ¯gs are simply the elements on
the diagonal of the Hessian at (V 0, V 0
S ).

6.3 Connection to Dynamic Input Conductances. The analysis method
of this article has some similarities with that of the dynamic input conduc-
tances (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. Invece, this ar-
ticle 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 contrasto, the voltage-clamp simulations
considered here could in principle be replaced by actual voltage-clamp
esperimenti.

The other main difference between both methods is that DICs are de-
fined around a steady state: all variables are at their steady-state value for
the voltage V at which the DIC is calculated. The function Iion(V, Vs, Vus) In
Questo articolo, Tuttavia, 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 Iion(V, Vs, Vus)
come segue:

D
o
w
N
o
UN
D
e
D

F
R
o
M
H

T
T

:
/
/

D
io
R
e
C
T
.

io
T
.

e
D
tu
N
e
C
o
UN
R
T
io
C
e
–
P
D

F
/

3
3
3
5
6
3
1
8
8
9
4
0
5
N
e
C
o
_
UN
_
0
1
3
5
3
P
D

B
sì
G
tu
e
S
T

o
N
0
7
S
e
P
e
M
B
e
R
2
0
2
3

From Biophysical to Integrate-and-Fire Modeling

g f (V ) =

gs(V ) =

gus(V ) =

∂Iion(V f

, Vs, Vus)

∂V f

∂Iion(V f

, Vs, Vus)
∂Vs
, Vs, Vus)

∂Vus

(cid:8)
(cid:8)
(cid:8)
(cid:8)
V f
(cid:8)
(cid:8)
(cid:8)
(cid:8)
V f
(cid:8)
(cid:8)
(cid:8)
(cid:8)
V f

=Vs=Vus=V

585

(6.3)

(6.4)

(6.5)

where g f , gs, and gus are the equivalent of the fast, slow, and ultraslow DICs,
rispettivamente. Since all voltages are taken as equal, it is clear that the DICs
only reveal a subset of the information provided by Iion(V, Vs, Vus).

6.4 Connection to I-V Curve Analysis. Ribar and Sepulchre (2019) pro-
pose 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 func-
tions in Ribar and Sepulchre (2019) are univariate. The current Iion is then a
sum of functions of V, Vs, and Vus. The I-V curves are defined as the sum of
all functions acting on their specific timescale and those faster—for exam-
ple, I f (V ) + Is(V ) for the slow I-V curve.

7 Conclusione

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 un-
der 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. Tutto
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 bio-
physical models used in this article take the original published values, ex-
cept where a change of parameter is indicated. For the Connor-Stevens
−2, ¯gL =
modello (Connor et al., 1977), these are the following: C = 1 μF cm
−2,
20 mS cm

−2, ¯gNa = 120 mS cm

−2, ¯gK = 20 mS cm

= 47.7 mS cm

−2, ¯gA

D
o
w
N
o
UN
D
e
D

F
R
o
M
H

T
T

:
/
/

D
io
R
e
C
T
.

io
T
.

e
D
tu
N
e
C
o
UN
R
T
io
C
e
–
P
D

F
/

3
3
3
5
6
3
1
8
8
9
4
0
5
N
e
C
o
_
UN
_
0
1
3
5
3
P
D

B
sì
G
tu
e
S
T

o
N
0
7
S
e
P
e
M
B
e
R
2
0
2
3

586

T. Van Pottelbergh, G. Drion, and R. Sepulchre

Tavolo 1: Parameters of the Integrate-and-Fire Models.

Figura

7 E 8 (center)
8 (left)
8 (right)
9
10–13

0.58
0.86
0.3
1.2
0.82

0.022
0.03
0.027
0.037
0.89

τs

6.7
3
23
1.7
4.3

τus

Vr
Vs,R
−40 −25
−45 −24
−35 −24
−40 −20
7.5

2.8
278 −45

(cid:3)Vus

0
1.7

VL = −17 mV, VNa = 55 mV, VK = −72 mV, VA
= −75 mV. The parameters
for the model of Drion et al. (2018) are the means of the parameters used
−2, ¯gNa =
in their network simulations: C = 1 μF cm
−2,
−2, ¯gK,D = 40 mS cm
170 mS cm
=
¯gH = 0.01 mS cm
120 mV, KD = 170, k1

= 4 mS cm
= 0.55 mS cm
−2, VL = −55 mV, VNa = 50 mV, VK = −85 mV, VCa

−2, ¯gL = 0.055 mS cm

= 0.1, k2

−2, ¯gK,Ca

−2, ¯gCa,T

= 0.01.

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
testo. The parameters of the integrate-and-fire models were found using the
method described in section 4, with SciPy’s Trust Region Reflective algo-
rithm (Branch, Coleman, & Li, 1999) for the least-squares optimization. IL
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, rispettivamente,
μF cm−2, ms, and mV.

Ringraziamenti

We thank Ilario Cirillo, Thiago Burghi, Luka Ribar, and Christian Grussler
for helpful discussions. T.V.P received a fees scholarship from the Engineer-
ing and Physical Sciences Research Council (https://www.epsrc.ac.uk) un-
der 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.

Riferimenti

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.

D
o
w
N
o
UN
D
e
D

F
R
o
M
H

T
T

:
/
/

D
io
R
e
C
T
.

io
T
.

e
D
tu
N
e
C
o
UN
R
T
io
C
e
–
P
D

F
/

3
3
3
5
6
3
1
8
8
9
4
0
5
N
e
C
o
_
UN
_
0
1
3
5
3
P
D

B
sì
G
tu
e
S
T

o
N
0
7
S
e
P
e
M
B
e
R
2
0
2
3

From Biophysical to Integrate-and-Fire Modeling

587

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.

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.

Brunel, N., Hakim, V., & Richardson, M. J. (2014). Single neuron dynamics and com-

putation. Current Opinion in Neurobiology, 25, 149–155.

Burghi, T. B., Schoukens, M., & Sepulchre, R. (2019). Feedback for nonlinear sys-
tem identification. In Proceedings of the 18th European Control Conference (pag. 1344–
1349). Piscataway, NJ: IEEE.

Ceballos, C. C., Roque, UN. C., & Leão, R. M. (2017). The role of negative conductances
in neuronal subthreshold properties and synaptic integration. Biophysical Reviews,
9(5), 827–834.

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.

Connor, J. A., Walter, D., & McKown, R. (1977). Neural repetitive firing: Modifica-
tions of the Hodgkin-Huxley axon suggested by experimental results from crus-
tacean axons. Biophysical Journal, 18(1):81–102.

Drion, G., Dethier, J., Franci, A., & Sepulchre, R. (2018). Switchable slow cellular con-
ductances determine robustness and tunability of network states. PLOS Compu-
tational Biology, 14(4), e1006125.

Drion, G., Franci, A., Dethier, J., & Sepulchre, R. (2015). Dynamic input conductances

shape neuronal spiking. eNeuro, 2(1).

Drion, G., Franci, A., Seutin, V., & Sepulchre, R. (2012). A novel phase portrait for

neuronal excitability. PLOS One, 7(8), e41806.

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.

Enqvist, M., & Ljung, l. (2005). Linear approximations of nonlinear FIR systems for

separable input processes. Automatica, 41(3), 459–473.

Ermentrout, B. (1996). Type I membranes, phase resetting curves, and synchrony.

Calcolo neurale, 8(5), 979–1001.

FitzHugh, R. (1961). Impulses and physiological states in theoretical models of nerve

membrane. Biophysical Journal, 1(6), 445–466.

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.

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.

Franci, A., Drion, G., & Sepulchre, R. (2017). Robust and tunable bursting requires

slow positive feedback. Journal of Neurophysiology, 119, 1222–1234.

Hill, UN. V. (1936). Excitation and accommodation in nerve. In Proceedings of the Royal

Society of London B: Biological Sciences, 119(814), 305–355.

D
o
w
N
o
UN
D
e
D

F
R
o
M
H

T
T

:
/
/

D
io
R
e
C
T
.

io
T
.

e
D
tu
N
e
C
o
UN
R
T
io
C
e
–
P
D

F
/

3
3
3
5
6
3
1
8
8
9
4
0
5
N
e
C
o
_
UN
_
0
1
3
5
3
P
D

B
sì
G
tu
e
S
T

o
N
0
7
S
e
P
e
M
B
e
R
2
0
2
3

588

T. Van Pottelbergh, G. Drion, and R. Sepulchre

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.

Hodgkin, UN. L., & Huxley, UN. F. (1952). A quantitative description of membrane cur-
rent and its application to conduction and excitation in nerve. Journal of Physiol-
ogy, 117(4), 500–544.

Huys, Q. J. M., Ahrens, M. B., & Paninski, l. (2006). Efficient estimation of detailed

single-neuron models. Journal of Neurophysiology, 96(2), 872–890.

Izhikevich, E. M. (2003). Simple model of spiking neurons. IEEE Transactions on Neu-

ral Networks, 14(6), 1569–1572.

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.

Jolivet, R., Rauch, A., Lüscher, H.-R., & Gerstner, W. (2006). Predicting spike timing
of neocortical pyramidal neurons by simple threshold models. Journal of Compu-
tational Neuroscience, 21(1), 35–49.

Kepler, T. B., Abbott, l. F., & Marder, E. (1992). Reduction of conductance-based

neuron models. Biological Cybernetics, 66(5), 381–387.

Knight, B. W. (1972). Dynamics of encoding in a population of neurons. Journal of

General Physiology, 59(6), 734–766.

Kung, S. Y. (1978). A new identification and model reduction algorithm via singu-
lar value decomposition. In Proceedings of the 12th Asilomar Conference on Circuits,
Sistemi, and Computers (pag. 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.

Latham, P. E., Richmond, B. J., Nelson, P. G., & Nirenberg, S. (2000). Intrinsic dynam-
ics in neuronal networks. IO. Theory. Journal of Neurophysiology, 83(2), 808–827.
Lepora, N. F., Overton, P. G., & Gurney, K. (2012). Efficient fitting of conductance-
based model neurons from somatic current clamp. Journal of Computational Neu-
roscience, 32(1), 1–24.

Mauro, A., Conti, F., Dodge, F., & Schor, R. (1970). Subthreshold behavior and phe-
nomenological impedance of the squid giant axon. Journal of General Physiology,
55(4), 497–523.

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.

Mihalas, S., & Niebur, E. (2008). A generalized linear integrate-and-fire neural model

produces diverse spiking behaviors. Calcolo neurale, 21(3), 704–718.

Mugnaio, D. N., Hulett, J., McLaughlin, J., & de Callafon, R. UN. (2013). Thermal dynam-
ical identification of light-emitting diodes by step-based realization and convex
optimization. IEEE Transactions on Components, Packaging and Manufacturing Tech-
nology, 3(6), 997–1007.

Morse, T. M., Davison, UN. P., & Hines, M. (2001). Parameter space reduction in neu-
ron model optimization through minimization of residual voltage clamp current.
Society for Neuroscience Abstracts, 27.

D
o
w
N
o
UN
D
e
D

F
R
o
M
H

T
T

:
/
/

D
io
R
e
C
T
.

io
T
.

e
D
tu
N
e
C
o
UN
R
T
io
C
e
–
P
D

F
/

3
3
3
5
6
3
1
8
8
9
4
0
5
N
e
C
o
_
UN
_
0
1
3
5
3
P
D

B
sì
G
tu
e
S
T

o
N
0
7
S
e
P
e
M
B
e
R
2
0
2
3

From Biophysical to Integrate-and-Fire Modeling

589

Nagumo, J., Arimoto, S., & Yoshizawa, S. (1962). An active pulse transmission line

simulating nerve axon. Proceedings of the IRE, 50(10), 2061–2070.

Pozzorini, C., Naud, R., Mensi, S., & Gerstner, W. (2013). Temporal whitening by
power-law adaptation in neocortical neurons. Nature Neuroscience, 16(7), 942–948.
Ribar, L., & Sepulchre, R. (2019). Neuromodulation of neuromorphic circuits. IEEE

Transactions on Circuits and Systems I: Regular Papers, 66(8), 3028–3040.

Schoukens, M., & Rolain, Y. (2012). Parametric identification of parallel Wiener sys-
tems. IEEE Transactions on Instrumentation and Measurement, 61(10), 2825–2832.
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.

Stein, R. B. (1965). A theoretical analysis of neuronal variability. Biophysical Journal,

5(2), 173–194.

Van Pottelbergh, T., Drion, G., & Sepulchre, R. (2018). Robust modulation of

integrate-and-fire models. Calcolo neurale, 30(4), 987–1011.

Wehmeier, U., Dong, D., Koch, C., & Van Essen, D. (1989). Modeling the mammalian
visual system. In C. Koch & IO. Segev (Eds.), Methods in neuronal modeling (pag. 335–
359). Cambridge, MA: CON Premere.

Zeiger, H., & McEwen, UN. (1974). Approximate linear realizations of given dimen-
sion via Ho’s algorithm. IEEE Transactions on Automatic Control, 19(2), 153–153.

Received June 16, 2020; accepted September 24, 2020.

D
o
w
N
o
UN
D
e
D

F
R
o
M
H

T
T

:
/
/

D
io
R
e
C
T
.

io
T
.

e
D
tu
N
e
C
o
UN
R
T
io
C
e
–
P
D

F
/

3
3
3
5
6
3
1
8
8
9
4
0
5
N
e
C
o
_
UN
_
0
1
3
5
3
P
D

B
sì
G
tu
e
S
T

o
N
0
7
S
e
P
e
M
B
e
R
2
0
2
3
Scarica il pdf