Galves-Löcherbach Model in NESTML

In vivo, neurons display variability in their spiking behavior: continuous fluctuations of the membrane potential can be observed, and they exhibit spontaneous firing and different responses to the same repeated input stimulus.

To capture these phenomena, the Galves-Löcherbach (GL) model assumes that the firing of the neuron is a random event, whose probability of occurrence in any time step is a firing function \(Φ(V_\text{m})\) of membrane potential \(V_\text{m}\). By subsuming all sources of randomness into a single function, the GL neuron model simplifies the analysis and simulation of noisy spiking neural networks [1-3].

Preliminaries

[1]:

%matplotlib inline
from typing import Dict, Optional

import matplotlib as mpl

mpl.rcParams['axes.grid'] = True
mpl.rcParams['grid.color'] = 'k'
mpl.rcParams['grid.linestyle'] = ':'
mpl.rcParams['grid.linewidth'] = 0.5

import matplotlib.pyplot as plt

import nest
import numpy as np

from pynestml.codegeneration.nest_code_generator_utils import NESTCodeGeneratorUtils
from pynestml.codegeneration.nest_tools import NESTTools


import logging
logging.getLogger('matplotlib.font_manager').disabled = True


              -- N E S T --
  Copyright (C) 2004 The NEST Initiative

 Version: 3.9.0-post0.dev0
 Built: Oct 10 2025 14:26:32

 This program is provided AS IS and comes with
 NO WARRANTY. See the file LICENSE for details.

 Problems or suggestions?
   Visit https://www.nest-simulator.org

 Type 'nest.help()' to find out more about NEST.

GL model with biophysical units

We formulate the GL model as a linear integrate-and-fire model with membrane time constant \(\tau_\text{m}\), resting potential \(V_\text{r}\) and an instantaneous jump in the postsynaptic membrane potential upon receiving a spike (Dirac delta function postsynaptic kernel).

[2]:

# Neuron parameters
params = {
    "tau_m"   : 10.0,
    "t_ref"   : 2.0,
    "C_m"     : 250.0,
    "V_r"     : -65.0,
    "V_reset" : -65.0,
    "a"       : 1.2,
    "b"       : 27.0,
    "V_b"     : -51.3
}

[3]:

nestml_gl_exp_model = """
model gl_exp_neuron:
    state:
        refr_spikes_buffer mV = 0 mV
        refr_tick integer = 0    # Counts number of tick during the refractory period
        V_m mV = V_r     # Membrane potential

    equations:
        kernel G = delta(t)
        V_m' = -(V_m - V_r) / tau_m + (mV / ms) * convolve(G, spikes) + (I_e + I_stim) / C_m

    parameters:
        tau_m ms = 10 ms                  # Membrane time constant
        C_m pF = 250 pF                   # Capacitance of the membrane
        t_ref ms = 2 ms                   # Duration of refractory period
        tau_syn ms = 0.5 ms               # Time constant of synaptic current
        V_r mV = -65 mV                   # Resting membrane potential
        V_reset mV = -65 mV               # Reset potential of the membrane
        b real = 27                       # Parameter for the exponential curve
        a mV = 1.2 mV                     # Parameter for the exponential curve
        V_b mV = -51.3 mV                 # Membrane potential at which Phi(V)=1/b
        I_e pA = 0 pA                     # Constant external input current
        reset_after_spike boolean = true  # Whether to reset membrane potential after a spike is emitted

    internals:
        RefractoryCounts integer = steps(t_ref) # refractory time in steps

    input:
        spikes <- spike
        I_stim pA <- continuous

    output:
        spike

    function Phi(V_m mV) real:
        return (1 / b) * exp((V_m - V_b) / a)

    update:
        if refr_tick == 0:
            # neuron is not refractory
            integrate_odes()
        else:
            # neuron is absolute refractory
            refr_tick -= 1

        if random_uniform(0, 1) <= 1E-3 * resolution() * Phi(V_m):    # 1E-3 to convert ms to s
            # fire a spike!
            refr_tick = RefractoryCounts
            if reset_after_spike:
                V_m = V_reset

            emit_spike()
"""

[4]:

module_name, neuron_model_name = NESTCodeGeneratorUtils.generate_code_for(
    nestml_neuron_model=nestml_gl_exp_model,
    logging_level="ERROR"  # try "DEBUG" for more debug information
)

nest.Install(module_name)


              -- N E S T --
  Copyright (C) 2004 The NEST Initiative

 Version: 3.9.0-post0.dev0
 Built: Oct 10 2025 14:26:32

 This program is provided AS IS and comes with
 NO WARRANTY. See the file LICENSE for details.

 Problems or suggestions?
   Visit https://www.nest-simulator.org

 Type 'nest.help()' to find out more about NEST.


Feb 26 17:38:57 Install [Info]:
    loaded module nestml_f5ac9d5fccad41dea9f58ffd89f5fcb9_module

Firing rate

This should correspond to the \(\Phi(V_\text{m})\) function in the neuron (see plot below for the theoretical curve). Please note that in the literature, \(\Phi\) is sometimes defined as a probability of firing function with values between 0 and 1, but we here define it as giving the firing rate (in units of s\({}^{-1}\)).

[5]:

def measure_numerical_Phi_function(neuron_model_name, module_name, V_min=0., V_max=10., neuron_model_params=None, neuron_membrane_potential_name="V_m"):
    nest.ResetKernel()
    NESTTools.set_nest_verbosity("ERROR")

    nest.print_time = False
    nest.Install(module_name)
    nest.resolution = 1.   # check that results are independent of resolution...

    t_stop = 25000.

    V_range = np.linspace(V_min, V_max, 12)
    n_spikes = np.nan * np.ones_like(V_range)
    for i, V_m in enumerate(V_range):
        neuron = nest.Create(neuron_model_name)
        if neuron_model_params:
            neuron.set(neuron_model_params)

        neuron.set({neuron_membrane_potential_name: V_m})

        sr = nest.Create("spike_recorder")
        nest.Connect(neuron, sr)

        assert neuron.get(neuron_membrane_potential_name) == V_m
        nest.Simulate(t_stop)
        assert neuron.get(neuron_membrane_potential_name) == V_m

        n_spikes[i] = len(sr.events["times"])

    spike_rate = n_spikes / (t_stop / 1E3)

    return V_range, spike_rate

[6]:

# theoretical Phi vs V_m
V_range_theory = np.linspace(-60., -45., 100)
Phi_of_V_theory = (1 / params["b"]) * np.exp((V_range_theory - params["V_b"]) / params["a"])

# numerical Phi vs V_m
V_range_numeric, spike_rate_numeric = measure_numerical_Phi_function(neuron_model_name=neuron_model_name,
                                                                     module_name=module_name,
                                                                     V_min=-60.,
                                                                     V_max=-45.,
                                                                     neuron_model_params={"reset_after_spike": False,
                                                                                          "a": params["a"],
                                                                                          "b": params["b"],
                                                                                          "V_b": params["V_b"],
                                                                                          "tau_m": 1E99})

fig, ax = plt.subplots()
ax.plot(V_range_theory, Phi_of_V_theory, label="theory")
ax.plot(V_range_numeric, spike_rate_numeric, marker="o", label="numeric")
ax.legend()
ax.set_xlabel("$U$")
ax.set_ylabel("Firing rate [spikes/s]")

[6]:

Text(0, 0.5, 'Firing rate [spikes/s]')
../../_images/tutorials_gl_model_gl_model_tutorial_9_1.png

Computing reliability of the single neuron model

Here we plot the spike times for the neurons as a raster plot to study the stochastic behavior of the neuron model. First, we give the neurons a constant input current and then a frozen Poissonian input. Frozen Poissonian means that the spike times are random, but are exactly the same between trials.

We observe that with the constant input current, the spike times across trials are highly variable, and the firing rate at the bottom is “flat” indicating that the spike times are not reliably reproduced over trials. With the frozen Poissonian input, the corresponding firing rate panel shows several peaks indicating that the spike times tend to be more repeatable across trials.

[7]:

def evaluate_neuron(neuron_name, module_name, nneurons=1, neuron_parms=None, stimulus_type="constant", poisson_fr=0.0,
                    mu=500., sigma=0., t_sim=300., plot=False, rseed=1000, dt=0.1, input_freq=0.0):
    """
    Run a simulation in NEST for the specified neuron. Inject a stepwise
    current and plot the membrane potential dynamics and spikes generated.
    """
    nest.ResetKernel()
    NESTTools.set_nest_verbosity("ERROR")
    nest.print_time = False
    nest.Install(module_name)
    nest.SetKernelStatus({'rng_seed': rseed,
                          'resolution': dt})
    neuron = nest.Create(neuron_name, nneurons)
    if neuron_parms:
        for k, v in neuron_parms.items():
            nest.SetStatus(neuron, k, v)
    nest.SetStatus(neuron,{'V_m':np.random.uniform(-65.0, -50.0)})

    if stimulus_type == "noise":
        # Create a noise generator
        noise = nest.Create("noise_generator")
        # Set the parameters of the noise generator
        noise_params = {"mean": mu,
                        "std": sigma,
                        "dt": dt,
                        "frequency":input_freq}
        nest.SetStatus(noise, noise_params)
        nest.Connect(noise, neuron)

    elif stimulus_type == "poisson_spikes":
        # Create a Poisson generator device
        poisson = nest.Create("poisson_generator", params={"rate": poisson_fr})

        # Create a parrot neuron
        parrot = nest.Create("parrot_neuron")

        # Connect the Poisson generator to the parrot neuron
        nest.Connect(poisson, parrot)

        # Connect the parrot neuron to each neuron
        nest.Connect(parrot, neuron)
    else:
        assert stimulus_type == "constant"
        nest.SetStatus(neuron, "I_e", mu)

    multimeter = nest.Create("multimeter")
    multimeter.set({"record_from": ["V_m"],
                    "interval": dt})
    spike_recorder = nest.Create("spike_recorder")
    nest.Connect(multimeter, neuron)
    nest.Connect(neuron, spike_recorder)

    nest.Simulate(t_sim)

    dmm = nest.GetStatus(multimeter)[0]
    Voltages = dmm["events"]["V_m"]
    tv = dmm["events"]["times"]

    dSD = nest.GetStatus(spike_recorder, keys="events")[0]
    ns = dSD["senders"]
    ts = dSD["times"]

    _idx = [np.argmin((tv - spike_time)**2) - 1 for spike_time in ts]
    V_m_at_spike_times = Voltages[_idx]

    if plot:
        fig, ax = plt.subplots()
        ax.plot(tv, Voltages)
        ax.scatter(ts, V_m_at_spike_times)
        ax.set_xlabel("Time [ms]")
        ax.set_ylabel("V_m [mV]")
        ax.grid()

    return ts, ns

[8]:

dt = 0.1
nneurons = 50

ts_const_input, ns_const = evaluate_neuron(neuron_model_name,
                                     module_name,
                                     nneurons=nneurons,
                                     neuron_parms=params,
                                     stimulus_type="constant",
                                     mu=550.,
                                     t_sim=500.0,
                                     dt=dt,
                                     poisson_fr=0.0)

ts_poisson_input, ns_noise = evaluate_neuron(neuron_model_name,
                                     module_name,
                                     nneurons=nneurons,
                                     neuron_parms=params,
                                     stimulus_type="poisson_spikes",
                                     t_sim=500.0,
                                     poisson_fr=2000.0)

# rate averaging parameters
t_start = 0.0
t_stop = 500.0
t_step = 5.0
t_bins = np.arange(t_start, t_stop + t_step, t_step)

# Calculate the average rate
rate_const_input, _ = np.histogram(ts_const_input, bins=t_bins)
rate_const_input = rate_const_input / nneurons / (t_step * 1e-3)      # per trial, per second

rate_poisson_input, _ = np.histogram(ts_poisson_input, bins=t_bins)
rate_poisson_input = rate_poisson_input / nneurons / (t_step * 1e-3)  # per trial, per second

min_rate = min(np.amin(rate_const_input), np.amin(rate_poisson_input))
max_rate = max(np.amax(rate_const_input), np.amax(rate_poisson_input))

plt.figure(figsize=(10,4))

# Raster plots
plt.subplot2grid((3,2),(0,0), rowspan=2)
plt.plot(ts_const_input, ns_const, ".")
plt.xlim(100, 500)
plt.ylabel("#Trial")
plt.title("Constant input current")

plt.subplot2grid((3,2),(0,1), rowspan=2)
plt.plot(ts_poisson_input, ns_noise, ".")
plt.xlim(100, 500)
plt.ylabel("#Trial")
plt.title("Frozen Poisson input")

# Firing rate plots
plt.subplot2grid((3,2),(2,0), rowspan=1)
plt.plot(t_bins[:-1], rate_const_input)
plt.xlabel("Time [ms]")
plt.ylabel("Firing rate [spikes/s]")
plt.xlim(100, 500)
plt.ylim(min_rate, max_rate)

plt.subplot2grid((3,2),(2,1), rowspan=1)
plt.plot(t_bins[:-1], rate_poisson_input)
plt.xlabel("Time [ms]")
plt.ylabel("Firing rate [spikes/s]")
plt.xlim(100, 500)
plt.ylim(min_rate, max_rate)
plt.tight_layout()
../../_images/tutorials_gl_model_gl_model_tutorial_12_0.png

Acknowledgements

The authors would like to thank Renan Shimoura, Christophe Pouzat, Antonio Roque and Antonio Galves for their kind and helpful discussion and feedback.

References

[1] Brochini, L., de Andrade Costa, A., Abadi, M. et al. Phase transitions and self-organized criticality in networks of stochastic spiking neurons. Sci Rep 6, 35831 (2016). https://doi.org/10.1038/srep35831

[2] Lima, V., Pena, R. F. O., Shimoura, R. O., Kamiji, N. L., Ceballos, C. C., Borges, F. S., Higa, G. S. V., de Pasquale, R., & Roque, A. C. (2021). Modeling and characterizing stochastic neurons based on in vitro voltage-dependent spike probability functions. The European Physical Journal Special Topics 2021 230:14, 230(14), 2963–2972. https://doi.org/10.1140/EPJS/S11734-021-00160-7

[3] Christophe Pouzat (2020). Spike Train Analysis and Modeling. (Syllabus) Latin-American School on Computational Neuroscience (LASCON) 2020, São Paulo, Brazil.