# 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,

WHITE_NOISE {
W
}


which declares the white noise process $$W$$, and the specification of the stochastic solver method,

BREAKPOINT {
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 = (
arbor.decor()
.set_property(Vm=-40)
.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):
arbor.recipe.__init__(self)
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")]

def event_generators(self, gid):
return self.the_gens


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
generators.append(arbor.event_generator(f"calcium_synapse_{i}", 0.0, sched_pre))


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

    # Create recipe
recipe = stdp_recipe(cell, cable_properties, generators)


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

    # Select one thread and no GPU
context = arbor.context(alloc, mpi=None)

# 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, 0), arbor.explicit_schedule([t1 - dt]))

# Run simulation
sim.run(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.

    # 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))
) / (beta + (1 - beta) * b)
return pandas.DataFrame({"ds": ds_A, "ms": time_lag, "type": "Arbor"})


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:

with multiprocessing.Pool() as p:
results = p.map(run, stdp_dt)


The collected results can then be plotted:

## The full code¶

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

## References¶

1(1,2)

Graupner and Brunel, PNAS 109 (10): 3991-3996 (2012); https://doi.org/10.1073/pnas.1109359109, https://www.pnas.org/doi/10.1073/pnas.1220044110.