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.
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 conductance1 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
The voltage equation 2.1 is the classical Kirchhoff relationship of a biophysical model: the current is the total ionic current intrinsic to the cellular membrane composition. Its voltage dependence is modeled through the voltage variable and lagged variables (, , ) that model the voltage filtered in distinct timescales (slow, ultraslow, ). 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 .
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 to 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 from the response of the biophysical model to voltage step inputs.
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.
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 . In this section, we determine 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.
The response of the total ionic current in the biophysical model to a voltage-clamp step is then obtained by direct substitution of and (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 .
The details of this simple idea are provided in the next sections.
3.1 Two-Timescale Models
The interpretation of this expression is that the value of the ion current function of the multiscale model evaluated at and can be read at from a voltage-clamp step experiment from to . Given that and can be chosen arbitrarily, this gives a simple method to evaluate on the biophysical model and use its value in the function of the multiscale model.
3.2 Three-Timescale Models
We apply the same idea to three-timescale models to find 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 . This can be resolved by devising a slightly more complex voltage clamp step experiment (see Figure 3): starting at , stepping to and, finally, after .
Following the same reasoning as before, after the fast gating variables will be approximately at their steady-state value . Similarly the slow gating variables will approximately be at the steady-state value , while the ultraslow gating variables will still be close to . The value of the current at time will therefore be a good approximation of of the multiscale model, assuming .
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 (or for three-timescale models). This is illustrated in Figure 4.
3.4 Calcium Dynamics
4 Iterative Optimization of the Structural Parameters
The previous section showed how to identify the function 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.
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 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 , the only parameters left to identify in model 2.4 to 2.6 are the membrane capacitance ; the reset parameters , , and ; and the cutoff voltage . Those parameters are easily estimated from current clamp simulations of the biophysical model.
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 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 is always taken to be equal to . The reset values of the slower voltages are parameters to be optimized but are physiologically constrained. A reasonable initialization is to set sufficiently (e.g., 20 mV) above . The parameter is constrained to be positive, since 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 is initially set to 1 F cm, as is often the case in a conductance-based model.
4.3 Local Parameter Optimization
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 -nullcline and the -nullcline . In other words, the level curves of the identified ion current 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 as a bifurcation parameter and by studying the family of phase portraits parameterized by . 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.
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 . 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 and . They reflect the timescale separation 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 displayed in Figure 7 (left). For a specific value of the applied current , the -nullcline intersects the -nullcline at the transcritical singularity (Franci, Drion, & Sepulchre, 2012). As the current increases from below to above the value , 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 . For , 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 and the structural parameters using the new value for . 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 than the original model, but both f-I curves remain close to each other afterward.
On the other hand, increasing above its nominal value of 47.7 mS cm to 200 mS cm results in a bistable phase portrait (see Figure 8, bottom right): for a specific range of , a stable hybrid limit cycle on the upper branch coexists with a stable fixed point on the lower branch of the -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 , 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 , , and . 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 , 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 , we did not find a similar improvement for (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.
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).
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 - space are shown for different values of . 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 , 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, 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.
When the current is lowered to , 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 -nullcline (see Figure 12, bottom left). The subsequent spiking, on the limit cycle on the upper branch of the -nullcline, causes to increase. This increase moves the two branches of the -nullcline closer, resulting in bistability between a limit cycle and a stable fixed point (see Figure 12, bottom center). As 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 and are a single variable in the reduced model.
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 is simply replaced by the associated voltage of the steady-state function: . 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
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 at which the DIC is calculated. The function 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.
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 is then a sum of functions of , , and . The I-V curves are defined as the sum of all functions acting on their specific timescale and those faster—for example, for the slow I-V curve.
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.
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: , , , , , , , , . The parameters for the model of Drion et al. (2018) are the means of the parameters used in their network simulations: , , , , , , , , , , , , , .
|7 and 8 (center)||0.58||0.022||6.7||−40||−25|
|7 and 8 (center)||0.58||0.022||6.7||−40||−25|
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 , the time constants, and the reset voltages are, respectively, F cm, ms, and mV.
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.
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.