Sde Tutorial (#2044)

Implementation of stochastic calcium plasticity curves, originally
authored by Sebastian Schmitt @schmitts, and adapted here for Arbor
inclusion.
- python example plasticity_stdp.py
- stochastic mechanism calcium_based_synapse.mod
- tutorial

Closes #1987

doc/dev/index.rst

@@ -27,4 +27,4 @@ Here we document internal components of Arbor. These pages can be useful if you'
    mechanism_abi
    util
    version
-.. numerics
+   numerics

doc/dev/numerics.rst

@@ -18,3 +18,5 @@ Mechanisms
 Exponential Euler `cnexp`.
 
 Implicit Euler `sparse`.
+
+Euler-Maruyama (explicit Euler) `stochastic`

doc/tutorial/calcium_stdp_curve.rst

+.. _tutorial_calcium_stpd_curve:
+
+Spike Timing-dependent Plasticity Curve
+=======================================
+
+This tutorial uses a single cell and reproduces `this Brian2 example
+<https://brian2.readthedocs.io/en/latest/examples/frompapers.Graupner_Brunel_2012.html>`_.  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, :math:`\rho`, which is a function of the
+calcium concentration, :math:`c(t)`. There are two stable states at :math:`\rho=0` (DOWN) and
+:math:`\rho=1` (UP), while :math:`\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:
+
+.. math::
+
+   \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*}
+
+Here, the sums over :math:`i` and :math:`j` represent the contributions from all pre and
+postsynaptic spikes, respectively, with :math:`C_{pre}` and :math:`C_{pre}` denoting the jumps in
+concentration after a spike. The jump after the presynaptic spike is delayed by :math:`D`.  The
+calcium decay time is assumed to be much faster than the synaptic time scale,
+:math:`\tau_{Ca} \ll \tau`. The subscripts :math:`p` and :math:`d` represent potentiation (increase
+in synaptic efficacy) and depression (decrease in synaptic efficacy), respectively, with
+:math:`\gamma` and :math:`\theta` being the corresponding rates and thresholds. :math:`H(x)` is the
+right-continuous heaviside step function (:math:`H(0)=1`).
+
+This mechanism is stochastic, :math:`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 :ref:`Arbor's NMODL dialect <format-sde>`. Let's examine
+the mechanism code in the `Arbor repository
+<https://github.com/arbor-sim/arbor/mechanisms/stochastic/calcium_based_synapse.mod>`_.
+
+The main difference compared to a deterministic (ODE) description is the additional `WHITE_NOISE`
+block,
+
+.. code:: none
+
+   WHITE_NOISE {
+       W
+   }
+
+which declares the white noise process :math:`W`, and the specification of the `stochastic` solver
+method,
+
+.. code:: none
+
+   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 :ref:`this overview <mechanisms-sde>`, while
+details about the numerical solver can be found in the :ref:`developers guide <sde>`.
+
+
+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:
+
+.. literalinclude:: ../../python/example/calcium_stdp.py
+   :language: python
+   :lines: 13-18
+
+Next we set the simulation parameters in order to reproduce the plasticity curve:
+
+.. literalinclude:: ../../python/example/calcium_stdp.py
+   :language: python
+   :lines: 20-41
+
+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, :math:`\rho(0)=0` and
+:math:`\rho(0)=1`, and for each initial state we want to run :math:`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 :math:`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:
+
+.. literalinclude:: ../../python/example/calcium_stdp.py
+   :language: python
+   :lines: 44-67
+
+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:
+
+.. literalinclude:: ../../python/example/calcium_stdp.py
+   :language: python
+   :lines: 70-74
+
+Our cell and cell properties can then later be used to create a simple recipe:
+
+.. literalinclude:: ../../python/example/calcium_stdp.py
+   :language: python
+   :lines: 77-103
+
+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
+:math:`\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:
+
+.. literalinclude:: ../../python/example/calcium_stdp.py
+   :language: python
+   :lines: 106-134
+
+The pre- and postsynaptic events are generated at explicit schedules, where the presynaptic event
+is shifted in time by :math:`D -\text{time lag}` with respect to the presynaptic event, which in
+turn is generated regularly with the frequency :math:`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
+
+.. literalinclude:: ../../python/example/calcium_stdp.py
+   :language: python
+   :lines: 136-137
+
+Now, we need to initialize the simulation, register a probe and run the simulation:
+
+.. literalinclude:: ../../python/example/calcium_stdp.py
+   :language: python
+   :lines: 139-154
+
+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.
+
+.. literalinclude:: ../../python/example/calcium_stdp.py
+   :language: python
+   :lines: 156-174
+
+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:
+
+.. literalinclude:: ../../python/example/calcium_stdp.py
+   :language: python
+   :lines: 177-178
+
+The collected results can then be plotted:
+
+.. figure:: calcium_stdp.svg
+    :width: 1600
+    :align: center
+
+    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``.
+
+References
+----------
+.. [1] 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>`_.

doc/tutorial/index.rst

@@ -52,6 +52,14 @@ Probes
 
    probe_lfpykit
 
+Stochastic Mechanisms
+---------------------
+
+.. toctree::
+   :maxdepth: 1
+
+   calcium_stdp_curve
+
 Hardware
 --------
 

mechanisms/CMakeLists.txt

@@ -23,7 +23,7 @@ make_catalogue(
 
 make_catalogue(
   NAME stochastic
-  MOD ou_input
+  MOD ou_input calcium_based_synapse
   VERBOSE ${ARB_CAT_VERBOSE}
   ADD_DEPS ON)
 

mechanisms/stochastic/calcium_based_synapse.mod

+: Calcium-based plasticity model
+: Based on the work of 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
+:
+: Author: Sebastian Schmitt
+:
+: The synapse is modeled as synaptic efficacy variable, ρ, which is a function of the calcium
+: concentration, c(t). The synapse model features two stable states at ρ=0 (DOWN) and ρ=1 (UP),
+: while ρ=ρ_star=0.5 represents a third unstable state between the two stable states.
+: The calcium concentration dynamics are represented by a simplified model which ueses a linear sum
+: of individual calcium transients elicited by trains of pre- and postsynaptic action potentials.
+:
+: drho/dt = -(1/τ)ρ(1-ρ)(ρ_star-ρ) + (γ_p/τ)(1-ρ) H[c(t)-θ_p] - (γ_d/τ)ρ H[c(t)-θ_d] + N
+:       N = (σ/√τ) √(H[c(t)-θ_p] + H[c(t)-θ_d]) W
+:
+:   dc/dt = -(1/τ_Ca)c + C_pre Σ_i δ(t-t_i-D) + C_post Σ_j δ(t-t_j)
+:
+: rho      synaptic efficacy variable (unit-less)
+: rho_star second root of cubic polynomial (unit-less), rho_star=0.5
+: rho_0    initial value (unit-less)
+: tau      synaptic time constant (ms), order of seconds to minutes
+: gamma_p  rate of synaptic increase (unit-less)
+: theta_p  potentiaton threshold (concentration)
+: gamma_d  rate of synaptic decrease (unit-less)
+: theta_d  depression threshold (concentration)
+: sigma    noise amplitude
+: W        white noise
+: c        calcium concentration (concentration)
+: C_pre    concentration jump after pre-synaptic spike (concentration)
+: C_post   concentration jump after post-synaptic spike (concentration)
+: tau_Ca   Calcium decay time constant (ms), order of milliseconds
+: H        right-continuous heaviside step function ( H[x]=1 for x>=0; H[x]=0 otherwise )
+: t_i      presynaptic spike times
+: t_j      postsynaptic spike times
+: D        time delay
+
+NEURON {
+    POINT_PROCESS calcium_based_synapse
+    RANGE rho_0, tau, theta_p, gamma_p, theta_d, gamma_d, C_pre, C_post, tau_Ca, sigma
+}
+
+STATE {
+    c
+    rho
+}
+
+PARAMETER {
+    rho_star = 0.5
+    rho_0 = 1
+    tau = 150000 (ms)
+    gamma_p = 321.808
+    theta_p = 1.3
+    gamma_d = 200
+    theta_d = 1
+    sigma = 2.8248
+    C_pre = 1
+    C_post = 2
+    tau_Ca = 20 (ms)
+}
+
+ASSIGNED {
+    one_over_tau
+    one_over_tau_Ca
+    sigma_over_sqrt_tau
+}
+
+INITIAL {
+    c = 0
+    rho = rho_0
+
+    one_over_tau = 1/tau
+    one_over_tau_Ca = 1/tau_Ca
+    sigma_over_sqrt_tau = sigma/(tau^0.5)
+}
+
+BREAKPOINT {
+    SOLVE state METHOD stochastic
+}
+
+WHITE_NOISE {
+    W
+}
+
+DERIVATIVE state {
+    LOCAL hsp
+    LOCAL hsd
+    hsp = step_right(c - theta_p)
+    hsd = step_right(c - theta_d)
+    rho' = (-rho*(1-rho)*(rho_star-rho) + gamma_p*(1-rho)*hsp - gamma_d*rho*hsd)*one_over_tau + (hsp + hsd)^0.5*sigma_over_sqrt_tau*W
+    c' = -c*one_over_tau_Ca
+}
+
+NET_RECEIVE(weight) {
+    c = c + C_pre
+}
+
+POST_EVENT(time) {
+    c = c + C_post
+}

python/example/calcium_stdp.py

+#!/usr/bin/env python3
+# This script is included in documentation. Adapt line numbers if touched.
+#
+# Authors: Sebastian Schmitt
+#          Fabian Bösch
+#
+# Single-cell simulation: Calcium-based synapse which models synaptic efficacy as proposed by
+# 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.
+# The synapse dynamics is affected by additive white noise. The results reproduce the spike
+# timing-dependent plasticity curve for the DP case described in Table S1 (supplemental material).
+
+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.
+
+# (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)
+
+
+# (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)
+
+
+# (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(), "")
+
+
+# (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
+
+
+# (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))
+
+    # Create recipe
+    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, 0), arbor.explicit_schedule([t1 - dt]))
+
+    # Run simulation
+    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))
+    ) / (beta + (1 - beta) * b)
+    return pandas.DataFrame({"ds": ds_A, "ms": time_lag, "type": "Arbor"})
+
+
+with multiprocessing.Pool() as p:
+    results = p.map(run, stdp_dt)
+
+ref = numpy.array(
+    [
+        [-100, 0.9793814432989691],
+        [-95, 0.981715028725338],
+        [-90, 0.9932274542583821],
+        [-85, 0.982392230227282],
+        [-80, 0.9620761851689686],
+        [-75, 0.9688482001884063],
+        [-70, 0.9512409611378684],
+        [-65, 0.940405737106768],
+        [-60, 0.9329565205853866],
+        [-55, 0.9146720800329048],
+        [-50, 0.8896156244609853],
+        [-45, 0.9024824529979171],
+        [-40, 0.8252814817763271],
+        [-35, 0.8171550637530018],
+        [-30, 0.7656877496052755],
+        [-25, 0.7176064429672677],
+        [-20, 0.7582385330838939],
+        [-15, 0.7981934216985763],
+        [-10, 0.8835208109434913],
+        [-5, 0.9390513341028807],
+        [0, 0.9927519271849183],
+        [5, 1.2354639175257733],
+        [10, 1.2255075694250952],
+        [15, 1.1760718597832],
+        [20, 1.1862298823123565],
+        [25, 1.1510154042112806],
+        [30, 1.125958948639361],
+        [35, 1.1205413366238108],
+        [40, 1.0812636495110723],
+        [45, 1.0717828284838595],
+        [50, 1.0379227533866708],
+        [55, 1.0392771563905585],
+        [60, 1.023024320343908],
+        [65, 1.046049171409996],
+        [70, 1.040631559394446],
+        [75, 1.0257331263516831],
+        [80, 1.0013538722817072],
+        [85, 1.0121890963128077],
+        [90, 1.0013538722817072],
+        [95, 1.0094802903050326],
+        [100, 0.9918730512544945],
+    ]
+)
+df_ref = pandas.DataFrame({"ds": ref[:, 1], "ms": ref[:, 0], "type": "Reference"})
+
+seaborn.set_theme()
+df = pandas.concat(results)
+df = pandas.concat([df, df_ref])
+plt = seaborn.relplot(kind="line", data=df, x="ms", y="ds", hue="type")
+plt.set_xlabels("lag time difference (ms)")
+plt.set_ylabels("change in synaptic strenght (after/before)")
+plt._legend.set_title("")
+plt.savefig("calcium_stdp.svg")

scripts/run_python_examples.sh

@@ -32,3 +32,4 @@ $PREFIX python python/example/network_two_cells_gap_junctions.py
 $PREFIX python python/example/diffusion.py
 $PREFIX python python/example/plasticity.py
 $PREFIX python python/example/v-clamp.py
+$PREFIX python python/example/calcium_stdp.py