Brette et al. 2007 Benchmark Framework

---

Run this example as a Jupyter notebook:

JupyterHub service

See our guide for more information and troubleshooting.

---

This module provides a common framework for running the Brette et al. 2007 simulator review benchmarks. The benchmarks create sparsely coupled networks of excitatory and inhibitory neurons which exhibit self-sustained activity after an initial stimulus. The model is based on the Vogels & Abbott network model [1].

This framework is used by the individual benchmark scripts: - coba.py (Benchmark 1: Conductance-based synapses) - cuba.py (Benchmark 2: Current-based synapses) - hh_coba.py (Benchmark 3: Hodgkin-Huxley neurons) - cuba_ps.py (Benchmark 4: Precise spiking)

References

[1]

Vogels TP, Abbott LF. 2005. Signal propagation and logic gating in networks of integrate-and-fire neurons. Journal of Neuroscience. 25(46):10786-10795. https://doi.org/10.1523/JNEUROSCI.3508-05.2005

import nest


def compute_rate(spike_recorder, n_rec, simtime, num_processes=1):
    """
    Compute average firing rate per neuron from spike recorder.

    Parameters
    ----------
    spike_recorder : NodeCollection
        Spike recorder device
    n_rec : int
        Number of neurons recorded from
    simtime : float
        Simulation time [ms]
    num_processes : int
        Number of MPI processes (default: 1)

    Returns
    -------
    float
        Average firing rate per neuron [spikes/s]
    """
    # Assume evenly distributed neurons across processes, which is valid here but not in general.
    n_spikes = spike_recorder.n_events
    n_neurons = n_rec / num_processes
    rate = n_spikes / (n_neurons * simtime) * 1000.0  # convert to spikes/s
    return rate


def build_network(params):
    """
    Build the benchmark network.

    Parameters
    ----------
    params : dict
        Dictionary containing all benchmark parameters

    Returns
    -------
    tuple
        (E_neurons, I_neurons, E_stimulus, E_recorder, I_recorder, build_time)
    """

    # Set kernel parameters
    nest.SetKernelStatus(
        {
            "resolution": params["dt"],
            "total_num_virtual_procs": params["virtual_processes"],
            "overwrite_files": True,
        }
    )

    # Set default neuron parameters
    nest.SetDefaults(params["model"], params["model_params"])

    # Create populations
    print("Creating excitatory population ...")
    E_neurons = nest.Create(params["model"], params["NE"])

    print("Creating inhibitory population ...")
    I_neurons = nest.Create(params["model"], params["NI"])

    # Create stimulus generator
    print("Creating excitatory stimulus generator ...")
    E_stimulus = nest.Create(params["stimulus"], params["stimulus_params"])

    # Create spike recorders
    print("Creating spike recorders ...")
    E_recorder = nest.Create(params["recorder"], params["recorder_params"])
    I_recorder = nest.Create(params["recorder"], params["recorder_params"])

    # Calculate connection counts
    CE = int(params["NE"] * params["epsilon"])  # excitatory connections per neuron
    CI = int(params["NI"] * params["epsilon"])  # inhibitory connections per neuron

    # Create custom synapse models
    nest.CopyModel("static_synapse", "syn_ex", params["E_synapse_params"])
    nest.CopyModel("static_synapse", "syn_in", params["I_synapse_params"])
    nest.SetDefaults("syn_ex", {"delay": params["delay"]})
    nest.SetDefaults("syn_in", {"delay": params["delay"]})

    # Connect nodes
    print("Connecting excitatory population ...")
    nest.Connect(
        E_neurons,
        E_neurons,
        conn_spec={"rule": "fixed_indegree", "indegree": CE},
        syn_spec="syn_ex",
    )

    nest.Connect(
        I_neurons,
        E_neurons,
        conn_spec={"rule": "fixed_indegree", "indegree": CI},
        syn_spec="syn_in",
    )

    print("Connecting inhibitory population ...")
    nest.Connect(
        E_neurons,
        I_neurons,
        conn_spec={"rule": "fixed_indegree", "indegree": CE},
        syn_spec="syn_ex",
    )

    nest.Connect(
        I_neurons,
        I_neurons,
        conn_spec={"rule": "fixed_indegree", "indegree": CI},
        syn_spec="syn_in",
    )

    print("Connecting Poisson stimulus ...")
    nest.Connect(
        E_stimulus,
        E_neurons[: params["Nstim"]],
        conn_spec="all_to_all",
        syn_spec="syn_ex",
    )

    print("Connecting spike recorders ...")
    nest.Connect(E_neurons[: params["Nrec"]], E_recorder)
    nest.Connect(I_neurons[: params["Nrec"]], I_recorder)

    return E_neurons, I_neurons, E_recorder, I_recorder


def run_simulation(params):
    """
    Build and run the benchmark simulation.

    Parameters
    ----------
    params : dict
        Dictionary containing all benchmark parameters

    Returns
    -------
    dict
        Dictionary containing simulation results and timing information
    """
    nest.ResetKernel()

    # Build network
    E_neurons, I_neurons, E_recorder, I_recorder = build_network(params)

    # Run simulation
    print("Simulating ...")
    nest.Simulate(params["simtime"])

    # Get timing from kernel status
    kernel_status = nest.GetKernelStatus()
    build_time = kernel_status["network_build_time"]
    sim_time = kernel_status["simulation_time"]

    # Calculate number of synapses
    N = len(E_neurons) + len(I_neurons)
    CE = int(params["NE"] * params["epsilon"])
    CI = int(params["NI"] * params["epsilon"])
    Nsyn = (CE + CI) * N + params["Nrec"] * 2 + params["Nstim"]

    # Compute rates
    num_processes = nest.num_processes
    E_rate = compute_rate(E_recorder, params["Nrec"], params["simtime"], num_processes)
    I_rate = compute_rate(I_recorder, params["Nrec"], params["simtime"], num_processes)

    # Print summary
    print("\n" + "=" * 60)
    print("Simulation summary")
    print("=" * 60)
    print(f"Number of Neurons : {N}")
    print(f"Number of Synapses: {Nsyn}")
    print(f"Excitatory rate   : {E_rate:.2f} spikes/s")
    print(f"Inhibitory rate   : {I_rate:.2f} spikes/s")
    print(f"Building time     : {build_time:.2f} s")
    print(f"Simulation time   : {sim_time:.2f} s")
    print("=" * 60 + "\n")

    return {
        "N": N,
        "Nsyn": Nsyn,
        "E_rate": E_rate,
        "I_rate": I_rate,
        "build_time": build_time,
        "sim_time": sim_time,
    }