Spike Timing-dependent Plasticity Curve

This tutorial uses a single cell and reproduces this Brian2 example. We aim to reproduce a spike timing-dependent plastivity curve which arises from stochastic calcium-based synapse dynamics described in Graupner and Brunel [1].

The synapse is modeled as synaptic efficacy variable, \(\rho\), which is a function of the calcium concentration, \(c(t)\). There are two stable states at \(\rho=0\) (DOWN) and \(\rho=1\) (UP), while \(\rho=\rho^\ast=0.5\) represents a third unstable state between the two stable states. The calcium concentration dynamics are represented by a simplified model which uses a linear sum of individual calcium transients elicited by trains of pre- and postsynaptic action potentials:

\[\begin{split}\begin{align*} c^\prime (t) &= - \frac{1}{\tau_{Ca}}c + C_{pre} \sum_i \delta(t-t_i-D) + C_{post} \sum_j \delta(t-t_j), \\ \rho^\prime(t) &= - \frac{1}{\tau}\left [ \rho (1 - \rho) (\rho^\ast - \rho) -\gamma_p (1-\rho) H\left(c - \theta_p \right) + \gamma_d \rho H\left(c - \theta_d \right) \right ] + N, \\ N &= \frac{\sigma}{\sqrt{\tau}} \sqrt{H\left( c - \theta_p \right) + H\left( c - \theta_d \right)} W. \end{align*}\end{split}\]

Here, the sums over \(i\) and \(j\) represent the contributions from all pre and postsynaptic spikes, respectively, with \(C_{pre}\) and \(C_{pre}\) denoting the jumps in concentration after a spike. The jump after the presynaptic spike is delayed by \(D\). The calcium decay time is assumed to be much faster than the synaptic time scale, \(\tau_{Ca} \ll \tau\). The subscripts \(p\) and \(d\) represent potentiation (increase in synaptic efficacy) and depression (decrease in synaptic efficacy), respectively, with \(\gamma\) and \(\theta\) being the corresponding rates and thresholds. \(H(x)\) is the right-continuous heaviside step function (\(H(0)=1\)).

This mechanism is stochastic, \(W\) represents a white noise process, and therefore our simulation needs to

  • use a stochastic synapse mechanism,

  • accumulate statistics over a large enough ensemble of initial states.

Implementation of a Stochastic Mechanism

Implementing a stochastic mechanism which is given by a stochastic differential equation (SDE) as above is straightforward to implement in Arbor’s NMODL dialect. Let’s examine the mechanism code in the Arbor repository.

The main difference compared to a deterministic (ODE) description is the additional WHITE_NOISE block,


which declares the white noise process \(W\), and the specification of the stochastic solver method,

    SOLVE state METHOD stochastic

This is sufficient to inform Arbor about the stochasticity of the mechanism. For more information about Arbor’s strategy to solve SDEs, please consult this overview, while details about the numerical solver can be found in the developers guide.

The Model

In this tutorial, the neuron model itself is simple with only passive (leaky) membrane dynamics, and it receives regular synaptic current input in one arbitrary chosen control volume (CV) to trigger regular spikes.

First we import some required modules:

import arbor
import random
import multiprocessing
import numpy  # You may have to pip install these.
import pandas  # You may have to pip install these.
import seaborn  # You may have to pip install these.

Next we set the simulation parameters in order to reproduce the plasticity curve:

# (1) Set simulation paramters

# Spike response delay (ms)
D = 13.7
# Spike frequency in Hertz
f = 1.0
# Number of spike pairs
num_spikes = 30
# time lag resolution
stdp_dt_step = 20.0
# Maximum time lag
stdp_max_dt = 100.0
# Ensemble size per initial value
ensemble_per_rho_0 = 100
# Simulation time step
dt = 0.1
# List of initial values for 2 states
rho_0 = [0] * ensemble_per_rho_0 + [1] * ensemble_per_rho_0
# We need a synapse for each sample path
num_synapses = len(rho_0)
# Time lags between spike pairs (post-pre: < 0, pre-post: > 0)
stdp_dt = numpy.arange(-stdp_max_dt, stdp_max_dt + stdp_dt_step, stdp_dt_step)

The time lag resolution, together with the maximum time lag, determine the number of cases we want to simulate. For each such case, however, we need to run many simulations in order to get a statistically meaningful result. The number of simulations per case is given by the ensemble size and the initial conditions. In our case, we have two inital states, \(\rho(0)=0\) and \(\rho(0)=1\), and for each initial state we want to run \(100\) simulations. We note, that the stochastic synapse mechanism does not alter the state of the cell, but couples one-way only by reacting to spikes. Therefore, we are allowed to simply place \(100\) synapses per initial state onto the cell without worrying about interference. Moreover, this has the benefit of exposing parallelism that Arbor can take advantage of.

Thus, we create a simple cell with a midpoint at which we place our mechanisms:

# (2) Make the cell

# Create a morphology with a single (cylindrical) segment of length=diameter=6 μm
tree = arbor.segment_tree()
tree.append(arbor.mnpos, arbor.mpoint(-3, 0, 0, 3), arbor.mpoint(3, 0, 0, 3), tag=1)

# Define the soma and its midpoint
labels = arbor.label_dict({"soma": "(tag 1)", "midpoint": "(location 0 0.5)"})

# Create and set up a decor object
decor = (
    .paint('"soma"', arbor.density("pas"))
    .place('"midpoint"', arbor.synapse("expsyn"), "driving_synapse")
    .place('"midpoint"', arbor.threshold_detector(-10), "detector")
for i in range(num_synapses):
    mech = arbor.mechanism("calcium_based_synapse")
    mech.set("rho_0", rho_0[i])
    decor.place('"midpoint"', arbor.synapse(mech), f"calcium_synapse_{i}")

# Create cell
cell = arbor.cable_cell(tree, decor, labels)

Since our stochastic mechanism calcium_based_synapse is not within Arbor’s default set of mechanism, we need to extend the mechanism catalogue within the cable cell properties:

# (3) Create extended catalogue including stochastic mechanisms

cable_properties = arbor.neuron_cable_properties()
cable_properties.catalogue = arbor.default_catalogue()
cable_properties.catalogue.extend(arbor.stochastic_catalogue(), "")

Our cell and cell properties can then later be used to create a simple recipe:

# (4) Recipe

class stdp_recipe(arbor.recipe):
    def __init__(self, cell, props, gens):
        self.the_cell = cell
        self.the_props = props
        self.the_gens = gens

    def num_cells(self):
        return 1

    def cell_kind(self, gid):
        return arbor.cell_kind.cable

    def cell_description(self, gid):
        return self.the_cell

    def global_properties(self, kind):
        return self.the_props

    def probes(self, gid):
        return [
            arbor.cable_probe_point_state_cell("calcium_based_synapse", "rho", "rho")

Note, that the recipe takes a cell, cell properties and a list of event generators as constructor arguments and returns them with its corresponding methods. Furthermore, the recipe also returns a list of probes which contains only one item: A query for our mechanism’s state variable \(\rho\). Since we placed a number of these mechanisms on our cell, we will receive a vector of values when probing.

Next we set up the simulation logic:

# (5) run simulation for a given time lag

def run(time_lag):
    # Time between stimuli
    T = 1000.0 / f

    # Simulation duration
    t1 = num_spikes * T

    # Time difference between post and pre spike including delay
    d = -time_lag + D

    # Stimulus and sample times
    t0_post = 0.0 if d >= 0 else -d
    t0_pre = d if d >= 0 else 0.0
    stimulus_times_post = numpy.arange(t0_post, t1, T)
    stimulus_times_pre = numpy.arange(t0_pre, t1, T)
    sched_post = arbor.explicit_schedule(stimulus_times_post)
    sched_pre = arbor.explicit_schedule(stimulus_times_pre)

    # Create strong enough driving stimulus
    generators = [arbor.event_generator("driving_synapse", 1.0, sched_post)]

    # Stimulus for calcium synapses
    for i in range(num_synapses):
        # Zero weight -> just modify synaptic weight via stdp

The pre- and postsynaptic events are generated at explicit schedules, where the presynaptic event is shifted in time by \(D -\text{time lag}\) with respect to the presynaptic event, which in turn is generated regularly with the frequency \(f\). The postsynaptic events are driven by the deterministic synapse with weight 1.0, while the presynaptic events are generated at the stochastic calcium synapses. The postsynaptic weight can be set arbitrarily as long as it is large enough to trigger the spikes.

Thus, we have all ingredients to create the recipe

        generators.append(arbor.event_generator(f"calcium_synapse_{i}", 0.0, sched_pre))

Now, we need to initialize the simulation, register a probe and run the simulation:

    recipe = stdp_recipe(cell, cable_properties, generators)

    # Select one thread and no GPU
    alloc = arbor.proc_allocation(threads=1, gpu_id=None)
    context = arbor.context(alloc, mpi=None)
    domains = arbor.partition_load_balance(recipe, context)

    # Get random seed
    random_seed = random.getrandbits(64)

    # Create simulation
    sim = arbor.simulation(recipe, context, domains, random_seed)

    # Register prope to read out stdp curve
    handle = sim.sample((0, "rho"), arbor.explicit_schedule([t1 - dt]))

Since we are interested in the long-term average value, we only query the probe at the end of the simulation.

After the simulation is finished, we calculate the change in synaptic strength by evaluating the transition probabilies from initial DOWN state to final UP state and vice versa.

    sim.run(t1, dt)

    # Process sampled data
    data, meta = sim.samples(handle)[0]
    data_down = data[-1, 1 : ensemble_per_rho_0 + 1]
    data_up = data[-1, ensemble_per_rho_0 + 1 :]
    # Initial fraction of synapses in DOWN state
    beta = 0.5
    # Synaptic strength ratio UP to DOWN (w1/w0)
    b = 5
    # Transition indicator form DOWN to UP
    P_UA = (data_down > 0.5).astype(float)
    # Transition indicator from UP to DOWN
    P_DA = (data_up < 0.5).astype(float)
    # Return change in synaptic strength
    ds_A = (
        (1 - P_UA) * beta
        + P_DA * (1 - beta)
        + b * (P_UA * beta + (1 - P_DA) * (1 - beta))

Since we need to run our simulation for each time lag case anew, we spawn a bunch of threads to carry out the work in parallel:

The collected results can then be plotted:


Comparison of this simulation with reference simulation [1]; for a simulation duration of 60 spikes at 1 Hertz, ensemble size of 2000 per initial state and time step dt=0.01 ms. The shaded region indicates the 95% confidence interval.

The full code

You can find the full code of the example at python/examples/calcium_stdp.py.