eglif_multirec - Conductance based extended generalized leaky integrate and fire neuron model

Description

eglif_multirec is the generalized leaky integrate and fire neuron according to De Grazia et al. (2026) [1], with specific post-synaptic conductances for AMPA, NMDA and GABA receptors.

The membrane potential is given by the following differential equation:

\[\begin{split}\begin{cases} C_m \dfrac{dV_m}{dt} = \dfrac{C_m}{\tau_m}(V_m - E_L) - I_{adap} + I_{dep} + I_e + I_{stim} + I_{syn}\\ \dfrac{dI_{adap}}{dt} = k_{adap}(V_m - E_L) - k_2 \cdot I_{adap}\\ \dfrac{dI_{dep}}{dt} = -k_1 \cdot I_{dep}\\ \end{cases}\end{split}\]

where the synaptic current \(I_{syn}\) integrates the input from postsynaptic receptors:

\[I_{syn} = \sum_{i=1} ^{4} g_i (V_m - E_{rev,i})\]

Here, the synapse i is excitatory or inhibitory depending on the value of \(E_{rev,i}\).

The receptor-specific synaptic currents are defined as:

\[ \begin{align}\begin{aligned}\begin{cases}\\I_{AMPA} = -\left[ g_{AMPA,1}(t)(V_m - E_{AMPA}) + g_{AMPA,2}(t)(V_m - E_{AMPA}) \right]\\\begin{split}\\\end{split}\\I_{NMDA} = -\left[ g_{NMDA,1}(t)(V_m - E_{NMDA})B_{Mg}(V_m) + g_{NMDA,2}(t)(V_m - E_{NMDA})B_{Mg}(V_m) \right]\\\begin{split}\\\end{split}\\I_{GABA} = -\left[ g_{GABA,1}(t)(V_m - E_{GABA}) + g_{GABA,2}(t)(V_m - E_{GABA}) \right]\\\begin{split}\\\end{split}\\I_{AMPA2} = -\left[ g_{AMPA2,1}(t)(V_m - E_{AMPA2}) + g_{AMPA2,2}(t)(V_m - E_{AMPA2}) \right]\\\begin{split}\\\end{split}\\I_{AMPA3} = -\left[ g_{AMPA3,1}(t)(V_m - E_{AMPA3}) + g_{AMPA3,2}(t)(V_m - E_{AMPA3}) \right]\\\end{cases}\end{aligned}\end{align} \]

The total synaptic current is:

\[I_{syn} = I_{AMPA} + I_{NMDA} + I_{GABA} + I_{AMPA2} + I_{AMPA3}\]

The NMDA current is modulated by the voltage-dependent magnesium block:

\[B_{Mg}(V_m) = \frac{1}{1 + \exp\left(\frac{v_{0,block} - V_{m,2}}{k_{block}}\right)}\]

Neuron produces spikes stochastically according to a point process with the firing intensity:

\[\lambda = \lambda_0 \cdot e^{\dfrac{V_m - V_{th}}{\tau_V}}\]

In case of spike emission, the spike-triggered adaptation currents \(I_{adap}\) and \(I_{dep}\) are respectively increased and set by their respective constants (which can be positive or negative):

\[\begin{split}V_m &= V_{reset}\\ I_{dep} &= A1\\ I_{adap} &= I_{adap} + A2\end{split}\]

Warning

The model described here is not matching the other LIF based models because of the sign in the membrane potential equations: the leak current should drive the membrane potential towards the resting state and not the opposite.

Neuron parameters

Spike Stochasticity

The generation of action potentials can be configured to behave either deterministically or stochastically through the boolean parameter stochastic_spiking (default: false).

  • Deterministic Mode (stochastic_spiking = false): A spike is generated subitaneously whenever the membrane potential crosses the fixed threshold value:

    \[V_m \ge V_{th}\]

  • Stochastic Mode (stochastic_spiking = true): The neuron produces spikes stochastically according to a point process. The instantaneous firing intensity (escape rate) \(\lambda\) depends exponentially on the distance between the current membrane potential and the threshold:

    \[\lambda = \lambda_0 \cdot e^{\dfrac{V_m - V_{th}}{\tau_V}}\]

    Where \(\lambda_0\) is the base escape rate and \(\tau_V\) is the voltage scale parameter. The probability of emitting a spike within a simulation time step \(\Delta t\) is given by:

    \[P(\text{spike}) = 1 - e^{-\lambda \cdot \Delta t}\]

In-vitro state

The automatic tuning procedure was applied to the parameters of the E-GLIF model (\(A_1\), \(A_2\), \(k_1\), \(k_2\), \(I_e\), \(k_{adap}\)), yielding the cell-specific values. The rest of the biophysical parameters were defined in accordance with the electrophysiological studies for each neuronal population [2].

For each cell type, the goal was to preserve the main electrophysiological features of the corresponding multicompartmental reference model within the reduced point-neuron description. Reference responses were generated through depolarizing and hyperpolarizing current-injection protocols in NEURON. The same protocols were then simulated with the E-GLIF model, and electrophysiological descriptors were extracted from both model responses.

These descriptors were used to define cell-specific objective (fitness) functions, yielding a multi-objective minimization problem that captured the \(f-I\) relationship together with relevant dynamical features such as:

  • Rheobase

  • Tonic firing rate

  • Coefficient of Variation (CV)

  • \(f-I\) response

  • \(f-I\) curve

  • Post-inhibitory pauses (depending on the cell type)

  • Rebound responses

See Table 1 for details in [1].

The parameter search was implemented in a custom Python workflow based on the DEAP evolutionary framework, using the NSGA-II algorithm as a multi-objective optimizer. This approach was chosen to handle the simultaneous calibration of multiple electrophysiological features over a nonlinear parameter space, a setting in which evolutionary methods have been widely used for neuronal model tuning.

Since NSGA-II returns a Pareto set of non-dominated solutions, the final parameter set for each cell type was selected using an achievement scalarizing function, identifying a balanced trade-off among the different objectives [1].

The resulting parameters for all optimized cells are summarized in Table 2 and fitness errors are reported in Table 3 of the reference paper [1].

Synaptic parameters

Static Synapse

In-vitro state

For the static synapse configuration, synaptic transmission integrates receptor-specific kinetics for AMPA, NMDA, and GABA receptors. To capture multi-timescale synaptic profiles, the total conductance \(g_{X}(t)\) for each receptor is modeled as the sum of a fast and a slow exponential contribution:

\[g_{X}(t) = g_{X,\text{fast}}(t) + g_{X,\text{slow}}(t)\]

where \(X \in \{\text{AMPA}, \text{NMDA}, \text{GABA}\}\).

The underlying kernel parameters regulating these kinetic profiles, i.e., specifically the initial conductance \(g_{\text{init}}\), the rise scaling factor \(A_r\), the decay scaling factors (\(A_{d1}\), \(A_{d2}\)), and the respective time constants (\(\tau_r\), \(\tau_{d1}\), \(\tau_{d2}\)) were calibrated using a multi-objective genetic algorithm (NSGA-II).

For the mathematical details underlying AMPA, NMDA, GABA synaptic kernels, see Equations 4–7 in [1].

The optimization was targeted against reference synaptic conductance traces extracted from multicompartmental NEURON simulations under voltage-clamp protocols. The resulting time constants and scaling factors for each connection and postsynaptic cell type are illustrated in Figure 5 of the reference paper [1].

Furthermore, to compensate for the loss of dendritic attenuation inherent to the reduction from spatially extended structures to point-neuron representations, a network-level weight calibration was performed. The final effective static synaptic weights optimized for each pathway are summarized in Table 4 of [1].

Tsodyks Markram Synapse

Short-term plasticity (STP) is incorporated into the network using the Tsodyks-Markram formalism, with parameters assigned according to each specific connection type. Circuits based on the De Grazia et al. model [1] leverage the tsodyks_synapse version of the model.

In-vitro state

The connection-specific parameters regulating these dynamic synapses are defined in accordance with S1 Table of the reference paper [1].

To preserve consistency between static and plastic network simulations, an automated synaptic weight rescaling procedure was applied. A rescaling factor relative to the corresponding static weight was computed for each connection and receptor type, ensuring that the area under the synaptic conductance curve matches that of the static model under low-frequency stimulation (5 Hz), where STP effects are assumed to be negligible [1].

The final calibrated effective STP synaptic weights resulting from this connection-specific rescaling are reported in Table 4 of the reference paper [1].

References