diff --git a/arbor/include/arbor/mechcat.hpp b/arbor/include/arbor/mechcat.hpp index bfe5113a57..45b30f83d3 100644 --- a/arbor/include/arbor/mechcat.hpp +++ b/arbor/include/arbor/mechcat.hpp @@ -112,6 +112,7 @@ class ARB_ARBOR_API mechanism_catalogue { ARB_ARBOR_API const mechanism_catalogue& global_default_catalogue(); ARB_ARBOR_API const mechanism_catalogue& global_allen_catalogue(); ARB_ARBOR_API const mechanism_catalogue& global_bbp_catalogue(); +ARB_ARBOR_API const mechanism_catalogue& global_stochastic_catalogue(); // Load catalogue from disk. ARB_ARBOR_API const mechanism_catalogue load_catalogue(const std::string&); diff --git a/doc/concepts/mechanisms.rst b/doc/concepts/mechanisms.rst index 42518b71ce..02d05ab4b4 100644 --- a/doc/concepts/mechanisms.rst +++ b/doc/concepts/mechanisms.rst @@ -77,6 +77,14 @@ Two catalogues are provided that collect mechanisms associated with specific pro * ``bbp_catalogue`` For models published by the Blue Brain Project (BBP). * ``allen_catalogue`` For models published on the Allen Brain Atlas Database. +A fourth catalogue ``stochastic_catalogue`` provides mechanisms expressed as stochastic differential +equations: + +* *ou_input* Synapse mechanism that can stochastically account for a population of *ou_input* + synapses. The mechanism is similar to *expsyn_curr* but with the exponential decay being subject + to noise due to a Ornstein-Uhlenbeck process. + + .. _mechanisms_dynamic: Adding Catalogues to Arbor @@ -316,6 +324,50 @@ Euler-Maruyama method. For specifics about the notation to define stochastic processes, please consult the :ref:`Arbor-specific NMODL extension `. +.. note:: + + While the units of :math:`\textbf{f}(t, \textbf{X}(t))` represent the deterministic rate of + change (per millisecond), + + .. math:: + + \left[\textbf{f}(t, \textbf{X}(t))\right] = \frac{\left[\textbf{X}(t)\right]}{ms}, + + the stochastic terms scale with :math:`t^{-1/2}`, + + .. math:: + + \left[\textbf{l}_i(t, \textbf{X}(t))\right] = \frac{\left[\textbf{X}(t)\right]}{\sqrt{ms}}. + + +**Example:** The popular Ornstein-Uhlenbeck process is described by a scalar linear mean-reverting SDE +and can be written as + +.. math:: + + X^\prime = -\frac{1}{\tau} (X - \mu) + \sqrt{\frac{2}{τ}} \sigma W, + +with white noise :math:`W`, and constant model parameters :math:`\tau`, :math:`\mu` and +:math:`\sigma`. The relaxation time :math:`\tau` determines how fast the process reverts back to its +mean value :math:`\mu`, and :math:`\sigma` controls the volatility (:math:`\mu` and :math:`\sigma` +have the same units as :math:`X`). The expected value and variance can be computed analytically and +yield + +.. math:: + + \begin{align*} + E[X] &= \mu - \left( \mu - X_0\right) e^{-t/\tau}, \\ + Var[X] &= \sigma^2 \left( 1 - e^{-2 t/\tau} \right), + \end{align*} + +which in the limit :math:`t \rightarrow \infty` converge to + +.. math:: + + \begin{align*} + E[X] &= \mu, \\ + Var[X] &= \sigma^2. + \end{align*} API --- diff --git a/example/CMakeLists.txt b/example/CMakeLists.txt index f139326839..b341b25736 100644 --- a/example/CMakeLists.txt +++ b/example/CMakeLists.txt @@ -14,3 +14,4 @@ add_subdirectory(probe-demo) add_subdirectory(plasticity) add_subdirectory(lfp) add_subdirectory(diffusion) +add_subdirectory(ornstein_uhlenbeck) diff --git a/example/ornstein_uhlenbeck/CMakeLists.txt b/example/ornstein_uhlenbeck/CMakeLists.txt new file mode 100644 index 0000000000..aa4d7ef98f --- /dev/null +++ b/example/ornstein_uhlenbeck/CMakeLists.txt @@ -0,0 +1,23 @@ +make_catalogue_standalone( + NAME ou + SOURCES "${CMAKE_CURRENT_SOURCE_DIR}" + MOD ornstein_uhlenbeck + CXX + CXX_FLAGS_TARGET ${ARB_CXX_FLAGS_TARGET_FULL} + VERBOSE ON) + +add_executable(ou EXCLUDE_FROM_ALL ou.cpp) +add_dependencies(ou ou-catalogue) +target_compile_options(ou PRIVATE ${ARB_CXX_FLAGS_TARGET_FULL}) +target_include_directories(ou PRIVATE "${CMAKE_CURRENT_BINARY_DIR}/generated/ou") + +target_link_libraries(ou PRIVATE arbor arborio ou-catalogue) +if (ARB_USE_BUNDLED_FMT) + target_include_directories(ou PRIVATE "${CMAKE_CURRENT_SOURCE_DIR}/../../ext/fmt/include") + target_compile_definitions(ou PRIVATE FMT_HEADER_ONLY) +else() + find_package(fmt REQUIRED) + target_link_libraries(ou PRIVATE fmt::fmt-header-only) +endif() + +add_dependencies(examples ou) diff --git a/example/ornstein_uhlenbeck/ornstein_uhlenbeck.mod b/example/ornstein_uhlenbeck/ornstein_uhlenbeck.mod new file mode 100644 index 0000000000..60f9e3ad56 --- /dev/null +++ b/example/ornstein_uhlenbeck/ornstein_uhlenbeck.mod @@ -0,0 +1,38 @@ +: Ornstein-Uhlenbeck: linear mean reverting process +: ================================================= +: dS(t) = kappa [mu - S(t)] dt + sigma dW(t) +: E[S(t)] = mu - (mu-S_0)*e^(-kappa t) +: Var[S(t)] = sigma^2/(2 kappa) [1 - e^(-2 kappa t)] +: in the limit: +: lim t->∞ E[S(t)] = mu +: lim t->∞ Var[S(t)] = sigma^2/(2 kappa) + +NEURON { + SUFFIX ornstein_uhlenbeck + GLOBAL kappa, mu, sigma +} + +PARAMETER { + kappa = 0.1 + mu = 1 + sigma = 0.1 +} + +STATE { + S +} + +INITIAL { + S=2 +} + +BREAKPOINT { + SOLVE state METHOD stochastic +} + +WHITE_NOISE { W } + +DERIVATIVE state { + S' = kappa*(mu - S) + sigma*W +} + diff --git a/example/ornstein_uhlenbeck/ou.cpp b/example/ornstein_uhlenbeck/ou.cpp new file mode 100644 index 0000000000..ad8f9c47b6 --- /dev/null +++ b/example/ornstein_uhlenbeck/ou.cpp @@ -0,0 +1,181 @@ +#include + +#include + +#include +#include +#include +#include +#include + +// forward declaration +arb::mechanism_catalogue build_catalogue(); + +// a single-cell recipe with probes +class recipe: public arb::recipe { +public: + recipe(unsigned ncvs) { + using namespace arb; + using namespace arborio::literals; + + // build catalogue with stochastic mechanism + cell_gprop_.catalogue = global_default_catalogue(); + cell_gprop_.catalogue.import(build_catalogue(), ""); + cell_gprop_.default_parameters = neuron_parameter_defaults; + + // paint the process on the whole cell + decor dec; + double const cv_size = 1.0; + dec.set_default(cv_policy_max_extent(cv_size)); + dec.paint("(all)"_reg , density("hh")); + dec.paint("(all)"_reg , density("ornstein_uhlenbeck")); + + // single-cell tree with ncvs control volumes + segment_tree tree; + tree.append(mnpos, {0, 0, 0.0, 4.0}, {0, 0, ncvs*cv_size, 4.0}, 1); + cell_ = cable_cell(morphology(tree), dec); + } + + arb::cell_size_type num_cells() const override { return 1; } + + arb::util::unique_any get_cell_description(arb::cell_gid_type) const override { return cell_; } + + arb::cell_kind get_cell_kind(arb::cell_gid_type) const override { return arb::cell_kind::cable; } + + std::vector get_probes(arb::cell_gid_type) const override { return probes_; } + + std::any get_global_properties(arb::cell_kind) const override { return cell_gprop_; } + + void add_probe(arb::probe_tag tag, std::any address) { probes_.emplace_back(std::move(address), tag); } + +protected: + std::vector probes_; + arb::cable_cell_global_properties cell_gprop_; + arb::cable_cell cell_; +}; + +// sampler for vector probes +struct sampler { + sampler(std::vector& data, std::size_t n_cvs, std::size_t n_steps): + n_cvs_{n_cvs} , n_steps_{n_steps}, data_{data} { + data_.resize(n_cvs*n_steps); + } + + void operator()(arb::probe_metadata pm, std::size_t n, const arb::sample_record* samples) { + const auto* m = arb::util::any_cast(pm.meta); + arb_assert(n_cvs_ == m->size()); + arb_assert(n_steps_ == n); + + for (std::size_t i=0; i(samples[i].data); + auto [lo, hi] = *data; + arb_assert(static_cast(hi-lo) == n_cvs_); + for (std::size_t j=0; j& data_; +}; + +// compute mean and variance online +// uses Welford-Knuth algorithm for the variance +struct accumulator { + std::size_t n_ = 0; + double mean_ = 0; + double var_ = 0; + + accumulator& operator()(double sample) { + double const delta = sample - mean_; + mean_ += delta / (++n_); + var_ += delta * (sample - mean_); + return *this; + } + + std::size_t n() const noexcept { return n_; } + double mean() const noexcept { return mean_; } + double variance() const noexcept { return n_ > 1 ? var_/(n_-1) : 0; } +}; + +int main(int argc, char** argv) { + + unsigned ncvs = 5000; // number of control volumes + double const dt = 1.0/1024; // time step + unsigned nsteps = 500; // number of time steps + + // create recipe and add probes + recipe rec(ncvs); + rec.add_probe(1, arb::cable_probe_density_state_cell{"ornstein_uhlenbeck", "S"}); + + // make context and simulation objects + auto context = arb::make_context({1, -1}); + arb::simulation sim = arb::simulation::create(rec) + .set_context(context) + .set_seed(137); + + // setup sampler and add it to the simulation with regular schedule + std::vector data; + sampler s{data, ncvs, nsteps}; + auto all_probes = [](arb::cell_member_type) { return true; }; + sim.add_sampler(all_probes, arb::regular_schedule(dt), s, arb::sampling_policy::lax); + + // run the simulation + sim.run(nsteps*dt, dt); + + // evaluate the mean for each time step across the ensembe of realizations + // (each control volume is a independent realization of the Ornstein-Uhlenbeck process) + std::vector acc(nsteps); + for (std::size_t t=0; t std::pair { + const double mu = 1.0; + const double S_0 = 2.0; + const double kappa = 0.1; + const double sigma = 0.1; + return { + mu - (mu-S_0)*std::exp(-kappa*t), + (sigma*sigma/(2*kappa))*(1.0 - std::exp(-2*kappa*t)) + }; + }; + + // print mean and expectation + for (std::size_t t=0; t(ptr); + for (int i=0; i(ty, *ic)); + if (ig) cat.register_implementation(nm, std::make_unique(ty, *ig)); + } + return cat; +} diff --git a/example/ornstein_uhlenbeck/readme.md b/example/ornstein_uhlenbeck/readme.md new file mode 100644 index 0000000000..a92066879d --- /dev/null +++ b/example/ornstein_uhlenbeck/readme.md @@ -0,0 +1,7 @@ +## Miniapp for demonstrating stochastic processes + +The simulation consists of a single cell with a stochastic process (Ornstein-Uhlenbeck) painted on +its control volumes. The stochastic process is described by a linear mean-reverting stochastic +differential equation which is specified in the accompanying NMODL file. All processes start from +the same initial condition and are averaged over the control volumes at each time step to generate an +ensemble statistic. These results are then compared to the analytical solution. diff --git a/mechanisms/CMakeLists.txt b/mechanisms/CMakeLists.txt index 9f1c14f5dc..56c0e20836 100644 --- a/mechanisms/CMakeLists.txt +++ b/mechanisms/CMakeLists.txt @@ -21,6 +21,12 @@ make_catalogue( VERBOSE ${ARB_CAT_VERBOSE} ADD_DEPS ON) +make_catalogue( + NAME stochastic + MOD ou_input + VERBOSE ${ARB_CAT_VERBOSE} + ADD_DEPS ON) + # This re-exports set(arbor-builtin-mechanisms ${arbor-builtin-mechanisms} PARENT_SCOPE) diff --git a/mechanisms/stochastic/ou_input.mod b/mechanisms/stochastic/ou_input.mod new file mode 100644 index 0000000000..2c63d75a7c --- /dev/null +++ b/mechanisms/stochastic/ou_input.mod @@ -0,0 +1,86 @@ +: Input current modeled by an Ornstein-Uhlenbeck process +: Sets non-specific current given mean, volatility and relaxation time +: +: Adapted from Jannik Luboeinski's example here: https://github.com/jlubo/arbor_ou_lif_example +: +: I_ou input current (nA) +: μ input current mean (nA) +: σ input current volatility (nA) +: τ relaxation time (ms) +: +: dI_ou/dt = -(1/τ) * (I_ou - μ) + √(2/τ) * σ * W +: +: statistic properties: +: E[I_ou] = μ - (μ - I_ou_0) * e^(-t/τ) +: Var[I_ou] = σ^2 * (1 - e^(-2 * t/τ)) +: +: in the limit t->∞: +: E[I_ou] = μ +: Var[I_ou] = σ^2 + +NEURON { + POINT_PROCESS ou_input + RANGE mu, sigma, tau + NONSPECIFIC_CURRENT I +} + +UNITS { + (ms) = (milliseconds) + (nA) = (nanoampere) +} + +PARAMETER { + mu = 1 (nA) : mean of the stochastic process + sigma = 1 (nA) : volatility of the stochastic process + tau = 1 (ms) : relaxation time +} + +STATE { + I_ou (nA) : instantaneous state + active : indicates if the process is currently enabled +} + +ASSIGNED { + alpha + beta +} + +WHITE_NOISE { + W +} + +INITIAL { + I_ou = 0 + active = -1 + alpha = 1.0/tau + beta = sigma * (2.0/tau)^0.5 +} + +BREAKPOINT { + SOLVE state METHOD stochastic + I = -I_ou +} + +DERIVATIVE state { + I_ou' = heaviside(active) * (alpha * (mu - I_ou) + beta * W) +} + +NET_RECEIVE(weight) { + if (weight >= 0) { : indicates that stimulation begins + I_ou = mu : initialize the process at the mean value + active = 1 + } + else { : indicates that stimulation ends + I_ou = 0 : switch off the process + active = -1 + } +} + +FUNCTION heaviside(x) { + if (x >= 0) { + heaviside = 1 + } + else { + heaviside = 0 + } +} diff --git a/python/example/ou_lif/Readme.md b/python/example/ou_lif/Readme.md new file mode 100644 index 0000000000..dedc039c29 --- /dev/null +++ b/python/example/ou_lif/Readme.md @@ -0,0 +1,37 @@ +# Ornstein-Uhlenbeck input to LIF neuron + +This example is taken from [Jannik Luboeinski's repo](https://github.com/jlubo/arbor_ou_lif_example) +and adapted for inclusion into Arbor. + +The simulation consists of a single compartment cable cell with simple leaky integrate-and-fire (LIF) +dynamics with Ornstein-Uhlenbeck (OU) input current. The OU input current can approximate the +synaptic input from a population of neurons. This is mathematically based on the diffusion +approximation, which analytically describes the input current that arises from presynaptic Poisson +spiking if the spiking activity is sufficiently high and the evoked postsynaptic potentials are +sufficiently small compared to the threshold potential [1,2]. Moreover, the exponentially decaying +autocorrelation function of the OU process resembles that of currents evoked by the spiking activity +of many neocortical neurons [3]. + +Here, we use the Ornstein-Uhlenbeck process in the following form, `dI/dt = -(1/τ) * (I - μ) + +√(2/τ) * σ * W`. Assuming homogeneous Poisson spiking at frequency *f* for a putative population of +*N* neurons, the mean current is `μ = N * f * w_out`, and the volatility (standard +deviation) is `σ = √(10^3*N*f/(2*τ))*w_out` where *f* is the firing rate in Hertz, *w_out* is the +synaptic weight in nC, and *τ* is the synaptic time constant in ms (cf. [4, 5]). + +![Resulting traces: voltage, spikes, currents](traces.svg) + +## References: + +[1] Ricciardi, LM. 1977. _Diffusion Processes and Related Topics in Biology_. +Springer, Berlin, Germany. + +[2] Moreno-Bote, R; Parga, N. 2010. Response of integrate-and-fire neurons to noisy inputs filtered +by synapses with arbitrary timescales: Firing rate and correlations. _Neural Comput._ 22. + +[3] Destexhe, A; Rudolph, M; Paré, D. 2003. The high-conductance state of neocortical neurons in +vivo. _Nat. Rev. Neurosci._ 4. + +[4] Dayan, P; Abbott, LF. 2001. _Theoretical Neuroscience_. MIT Press, Cambridge/MA, USA. + +[5] Luboeinski, J; Tetzlaff, C. 2021. Memory consolidation and improvement by synaptic tagging and +capture in recurrent neural networks. _Commun. Biol._ 275. diff --git a/python/example/ou_lif/lif.mod b/python/example/ou_lif/lif.mod new file mode 100644 index 0000000000..9271d6dcd2 --- /dev/null +++ b/python/example/ou_lif/lif.mod @@ -0,0 +1,72 @@ +: Leaky-integrate and fire neuron + +: If the threshold potential is crossed, the membrane resistance is set to a very low value +: to pull the membrane potential towards the reset potential for the duration of the refractory period. +: After that, the membrane resistance is set to the comparably higher leak resistance again. +: The threshold potential of Arbor's threshold detector should match the threshold potential here. + +: authors: Sebastian Schmitt, Jannik Luboeinski +: based on https://github.com/tetzlab/FIPPA/blob/main/STDP/mechanisms/lif.mod + +NEURON { + SUFFIX lif + NONSPECIFIC_CURRENT i + RANGE R_reset, R_leak, V_reset, V_rev, V_th, I_0, i_factor, t_ref +} + +UNITS { + (ms) = (milliseconds) + (mV) = (millivolt) + (MOhm) = (megaohm) +} + +STATE { + refractory_counter +} + +INITIAL { + refractory_counter = t_ref + 1 : start not refractory + + v = V_rev : membrane potential +} + +PARAMETER { + R_leak = 10 (MOhm) : membrane resistance during leakage (standard case) + R_reset = 1e-10 (MOhm) : membrane resistance during refractory period (very small to yield fast reset) + I_0 = 0 : constant input current (in nA) + i_factor = 1 : conversion factor from nA to mA/cm^2; for point neurons + + V_reset = -70 (mV) : reset potential + V_rev = -65 (mV) : reversal potential + V_th = -55 (mV) : threshold potential to be crossed for spiking + t_ref = 2 (ms) : refractory period +} + +BREAKPOINT { + SOLVE state METHOD cnexp + + LOCAL R_mem + LOCAL E + + : threshold crossed -> start refractory counter + if (v > V_th) { + refractory_counter = 0 + } + + : choose between leak and reset potential + if (refractory_counter <= t_ref) { : in refractory period - strong drive of v towards V_reset + R_mem = R_reset + E = V_reset + } else { : outside refractory period + R_mem = R_leak + E = V_rev + } + + i = (((v - E) / R_mem) - I_0) * i_factor : current density in units of mA/cm^2 +} + +DERIVATIVE state { + if (refractory_counter <= t_ref) { + refractory_counter' = 1 + } +} diff --git a/python/example/ou_lif/ou_lif.py b/python/example/ou_lif/ou_lif.py new file mode 100644 index 0000000000..c308ec18f1 --- /dev/null +++ b/python/example/ou_lif/ou_lif.py @@ -0,0 +1,342 @@ +#!/bin/python3 + +# adapted from Jannik Luboeinski's example here: https://github.com/jlubo/arbor_ou_lif_example + +import random +import subprocess +import arbor as arb +import numpy as np +import matplotlib.pyplot as plt + + +def make_catalogue(): + # build a new catalogue + out = subprocess.getoutput("arbor-build-catalogue ou_lif . --cpu True") + print(out) + # load the new catalogue and extend it with builtin stochastic catalogue + cat = arb.load_catalogue("./ou_lif-catalogue.so") + cat.extend(arb.stochastic_catalogue(), "") + return cat + + +def make_cell(): + + # cell morphology + # =============== + + tree = arb.segment_tree() + radius = 1e-10 # radius of cylinder (in µm) + height = 2 * radius # height of cylinder (in µm) + tree.append( + arb.mnpos, + arb.mpoint(-height / 2, 0, 0, radius), + arb.mpoint(height / 2, 0, 0, radius), + tag=1, + ) + labels = arb.label_dict({"center": "(location 0 0.5)"}) + + # LIF density mechanism + # ===================== + + # surface area of the cylinder in m^2 (excluding the circle-shaped ends, since Arbor does + # not consider current flux there) + area_m2 = 2 * np.pi * (radius * 1e-6) * (height * 1e-6) + area_cm2 = area_m2 * 1e4 + # conversion factor from nA to mA/cm^2; for point neurons + i_factor = (1e-9 / 1e-3) / area_cm2 + # neuronal capacitance in F + C_mem = 1e-9 + # specific capacitance in F/m^2, computed from absolute capacitance of point neuron + c_mem = C_mem / area_m2 + # reversal potential in mV + V_rev = -65.0 + # spiking threshold in mV + V_th = -55.0 + # initial synaptic weight in nC + h_0 = 4.20075 + # leak resistance in MOhm + R_leak = 10.0 + # membrane time constant in ms + tau_mem = 2.0 + + lif = arb.mechanism("lif") + lif.set("R_leak", R_leak) + lif.set("R_reset", 1e-10) + # set to initial value to zero (background input is applied via stochastic ou_bg_mech) + lif.set("I_0", 0) + lif.set("i_factor", i_factor) + lif.set("V_rev", V_rev) + lif.set("V_reset", -70.0) + lif.set("V_th", V_th) + # refractory time in ms + lif.set("t_ref", tau_mem) + + # stochastic background input current (point mechanism) + # ===================================================== + + # synaptic time constant in ms + tau_syn = 5.0 + # mean in nA + mu_bg = 0.15 + # volatility in nA + sigma_bg = 0.5 + # instantiate mechanism + ou_bg = arb.mechanism("ou_input") + ou_bg.set("mu", mu_bg) + ou_bg.set("sigma", sigma_bg) + ou_bg.set("tau", tau_syn) + + # stochastic stimulus input current (point mechanism) + # =================================================== + + # number of neurons in the putative population + N = 25 + # firing rate of the putative neurons in Hz + f = 100 + # synaptic weight in nC + w_out = h_0 / R_leak + # mean in nA + mu_stim = N * f * w_out + # volatility in nA + sigma_stim = np.sqrt((1000.0 * N * f) / (2 * tau_syn)) * w_out + # instantiate mechanism + ou_stim = arb.mechanism("ou_input") + ou_stim.set("mu", mu_stim) + ou_stim.set("sigma", sigma_stim) + ou_stim.set("tau", tau_syn) + + # paint and place mechanisms + # ========================== + + decor = arb.decor() + decor.set_property(Vm=V_rev, cm=c_mem) + decor.paint("(all)", arb.density(lif)) + decor.place('"center"', arb.synapse(ou_stim), "ou_stim") + decor.place('"center"', arb.synapse(ou_bg), "ou_bg") + decor.place('"center"', arb.threshold_detector(V_th), "spike_detector") + + return arb.cable_cell(tree, decor, labels) + + +class ou_recipe(arb.recipe): + def __init__(self, cell, cat): + arb.recipe.__init__(self) + + # simulation runtime parameters in ms + self.runtime = 20000 + self.dt = 0.2 + + # initialize catalogue and cell properties + self.the_props = arb.neuron_cable_properties() + self.the_props.catalogue = cat + self.the_cell = cell + + # protocol for stimulation + # "time_start": starting time, "scheme": protocol type + self.bg_prot = {"time_start": 0, "scheme": "FULL", "label": "ou_bg"} + self.stim1_prot = {"time_start": 10000, "scheme": "TRIPLET", "label": "ou_stim"} + self.stim2_prot = { + "time_start": 15000, + "scheme": "ONEPULSE", + "label": "ou_stim", + } + + def cell_kind(self, gid): + return arb.cell_kind.cable + + def cell_description(self, gid): + return self.the_cell + + def global_properties(self, kind): + return self.the_props + + def num_cells(self): + return 1 + + def probes(self, gid): + # probe membrane potential, total current, and external input currents + return [ + arb.cable_probe_membrane_voltage('"center"'), + arb.cable_probe_total_ion_current_cell(), + arb.cable_probe_point_state_cell("ou_input", "I_ou"), + ] + + def event_generators(self, gid): + gens = [] + # OU stimulus input #1 + gens.extend(self.get_generators(self.stim1_prot)) + # OU stimulus input #2 + gens.extend(self.get_generators(self.stim2_prot)) + # OU background input + gens.extend(self.get_generators(self.bg_prot)) + return gens + + # Returns arb.event_generator instances that describe the specifics of a given + # input/stimulation protocol for a mechanism implementing an Ornstein-Uhlenbeck process. + # Here, if the value of the 'weight' parameter is 1, stimulation is switched on, + # whereas if it is -1, stimulation is switched off. + def get_generators(self, protocol): + + prot_name = protocol["scheme"] # name of the protocol (defining its structure) + start_time = protocol["time_start"] # time at which the stimulus starts in s + label = protocol["label"] # target synapse (mechanism label) + + if prot_name == "ONEPULSE": + # create regular schedules to implement a stimulation pulse that lasts for 0.1 s + stim_on = arb.event_generator( + label, + 1, + arb.regular_schedule(start_time, self.dt, start_time + self.dt), + ) + stim_off = arb.event_generator( + label, + -1, + arb.regular_schedule( + start_time + 100, self.dt, start_time + 100 + self.dt + ), + ) + return [stim_on, stim_off] + + elif prot_name == "TRIPLET": + # create regular schedules to implement pulses that last for 0.1 s each + stim1_on = arb.event_generator( + label, + 1, + arb.regular_schedule(start_time, self.dt, start_time + self.dt), + ) + stim1_off = arb.event_generator( + label, + -1, + arb.regular_schedule( + start_time + 100, self.dt, start_time + 100 + self.dt + ), + ) + stim2_on = arb.event_generator( + label, + 1, + arb.regular_schedule( + start_time + 500, self.dt, start_time + 500 + self.dt + ), + ) + stim2_off = arb.event_generator( + label, + -1, + arb.regular_schedule( + start_time + 600, self.dt, start_time + 600 + self.dt + ), + ) + stim3_on = arb.event_generator( + label, + 1, + arb.regular_schedule( + start_time + 1000, self.dt, start_time + 1000 + self.dt + ), + ) + stim3_off = arb.event_generator( + label, + -1, + arb.regular_schedule( + start_time + 1100, self.dt, start_time + 1100 + self.dt + ), + ) + return [stim1_on, stim1_off, stim2_on, stim2_off, stim3_on, stim3_off] + + elif prot_name == "FULL": + # create a regular schedule that lasts for the full runtime + stim_on = arb.event_generator( + label, + 1, + arb.regular_schedule(start_time, self.dt, start_time + self.dt), + ) + return [stim_on] + + else: + return [] + + +if __name__ == "__main__": + + # set up and run simulation + # ========================= + + # create recipe + cell = make_cell() + cat = make_catalogue() + recipe = ou_recipe(cell, cat) + + # get random seed + random_seed = random.getrandbits(64) + print("random_seed = " + str(random_seed)) + + # select one thread and no GPU + alloc = arb.proc_allocation(threads=1, gpu_id=None) + context = arb.context(alloc, mpi=None) + domains = arb.partition_load_balance(recipe, context) + + # create simulation + sim = arb.simulation(recipe, context, domains, seed=random_seed) + + # create schedule for recording + reg_sched = arb.regular_schedule(0, recipe.dt, recipe.runtime) + + # set handles to probe membrane potential and currents + gid = 0 + handle_mem = sim.sample((gid, 0), reg_sched) # membrane potential + handle_tot_curr = sim.sample((gid, 1), reg_sched) # total current + handle_curr = sim.sample((gid, 2), reg_sched) # input current + + sim.record(arb.spike_recording.all) + sim.run(tfinal=recipe.runtime, dt=recipe.dt) + + # get traces and spikes from simulator + # ==================================== + + # there is only one location (single control volume) + loc = 0 + data_mem, data_curr = [], [] + + data, _ = sim.samples(handle_mem)[loc] + times = data[:, 0] + + # read out neuron data + data_mem.append(data[:, 1]) + + # get total membrane current + data_tot_curr, _ = sim.samples(handle_tot_curr)[loc] + data_curr.append(np.negative(data_tot_curr[:, 1])) + + # get Ornstein-Uhlenbeck currents + data_input_curr, _ = sim.samples(handle_curr)[loc] + data_curr.append(data_input_curr[:, 1]) + data_curr.append(data_input_curr[:, 2]) + + # get spikes + spike_times = sim.spikes()["time"] + + # assemble, store, and plot data + # ============================== + + data_header = "Time, V, I_stim, I_bg" + data_stacked = np.column_stack( + [times, np.transpose(data_mem), np.transpose(data_curr)] + ) + + np.savetxt("traces.txt", data_stacked, fmt="%.4f", header=data_header) + np.savetxt("spikes.txt", spike_times, fmt="%.4f") + + fig, axes = plt.subplots(nrows=4, ncols=1, sharex=False, figsize=(10, 10)) + axes[0].set_title("random_seed = " + str(random_seed)) + axes[0].plot(data_stacked[:, 0], data_stacked[:, 1], label="V", color="C3") + axes[0].plot( + spike_times, -55 * np.ones(len(spike_times)), ".", color="blue", markersize=1 + ) + axes[1].plot(data_stacked[:, 0], data_stacked[:, 2], label="I_tot", color="C1") + axes[2].plot(data_stacked[:, 0], data_stacked[:, 3], label="I_stim", color="C2") + axes[3].plot(data_stacked[:, 0], data_stacked[:, 4], label="I_bg", color="C0") + axes[3].set_xlabel("Time (ms)") + axes[0].set_ylabel("V (mV)") + axes[1].set_ylabel("I_tot (nA)") + axes[2].set_ylabel("I_stim (nA)") + axes[3].set_ylabel("I_bg (nA)") + fig.tight_layout() + fig.savefig("traces.svg", bbox_inches="tight") diff --git a/python/example/ou_lif/traces.svg b/python/example/ou_lif/traces.svg new file mode 100644 index 0000000000..13a3c49c53 --- /dev/null +++ b/python/example/ou_lif/traces.svg @@ -0,0 +1,22364 @@ + + + + + + + + 2022-11-09T15:50:52.860599 + image/svg+xml + + + Matplotlib v3.6.1, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/python/mechanism.cpp b/python/mechanism.cpp index a5ee801fc4..471d241609 100644 --- a/python/mechanism.cpp +++ b/python/mechanism.cpp @@ -200,6 +200,7 @@ void register_mechanisms(pybind11::module& m) { m.def("default_catalogue", [](){return arb::global_default_catalogue();}); m.def("allen_catalogue", [](){return arb::global_allen_catalogue();}); m.def("bbp_catalogue", [](){return arb::global_bbp_catalogue();}); + m.def("stochastic_catalogue", [](){return arb::global_stochastic_catalogue();}); m.def("load_catalogue", [](pybind11::object fn) { return arb::load_catalogue(util::to_string(fn)); }); // arb::mechanism_desc diff --git a/scripts/run_cpp_examples.sh b/scripts/run_cpp_examples.sh index efbcf8157c..5af6591ce7 100755 --- a/scripts/run_cpp_examples.sh +++ b/scripts/run_cpp_examples.sh @@ -26,7 +26,7 @@ check () { fi } -for ex in bench brunel gap_junctions generators lfp ring single-cell "probe-demo v" plasticity +for ex in bench brunel gap_junctions generators lfp ring single-cell "probe-demo v" plasticity ou do echo " - $ex" dir=`echo $ex | tr ' ' '_'` diff --git a/test/unit/CMakeLists.txt b/test/unit/CMakeLists.txt index 2a634d44cb..36892930db 100644 --- a/test/unit/CMakeLists.txt +++ b/test/unit/CMakeLists.txt @@ -117,11 +117,13 @@ set(unit_sources test_piecewise.cpp test_pp_util.cpp test_probe.cpp + test_rand.cpp test_range.cpp test_recipe.cpp test_ratelem.cpp test_schedule.cpp test_scope_exit.cpp + test_sde.cpp test_segment_tree.cpp test_simd.cpp test_simulation.cpp @@ -143,7 +145,6 @@ set(unit_sources test_unique_any.cpp test_vector.cpp test_version.cpp - test_sde.cpp # unit test driver test.cpp @@ -229,4 +230,4 @@ target_compile_definitions(unit PRIVATE "-DDATADIR=\"${CMAKE_CURRENT_SOURCE_DIR} target_compile_definitions(unit PRIVATE "-DLIBDIR=\"${PROJECT_BINARY_DIR}/lib\"") target_include_directories(unit PRIVATE "${CMAKE_CURRENT_BINARY_DIR}") target_include_directories(unit PRIVATE "${CMAKE_CURRENT_BINARY_DIR}/generated/testing") -target_link_libraries(unit PRIVATE ext-gtest arbor arborenv arborio arborio-private-headers arbor-private-headers arbor-sup) +target_link_libraries(unit PRIVATE ext-gtest ext-random123 arbor arborenv arborio arborio-private-headers arbor-private-headers arbor-sup) diff --git a/test/unit/test_rand.cpp b/test/unit/test_rand.cpp new file mode 100644 index 0000000000..8cee82b9d5 --- /dev/null +++ b/test/unit/test_rand.cpp @@ -0,0 +1,54 @@ +#include +#include +#include "backends/rand_impl.hpp" + +bool verify_uniform(arb::cbprng::array_type counter, arb::cbprng::array_type key, + std::array expected) { + const auto r = arb::cbprng::generator{}(counter,key); + return + (r[0] == expected[0]) && + (r[1] == expected[1]) && + (r[2] == expected[2]) && + (r[3] == expected[3]); +} + +bool verify_normal(std::array uniform, std::array expected) { + const auto [a0, a1] = r123::boxmuller(uniform[0], uniform[1]); + const auto [a2, a3] = r123::boxmuller(uniform[2], uniform[3]); + return + ((long long)(a0*1000000) == expected[0]) && + ((long long)(a1*1000000) == expected[1]) && + ((long long)(a2*1000000) == expected[2]) && + ((long long)(a3*1000000) == expected[3]); +} + +TEST(cbprng, uniform) { + EXPECT_TRUE(verify_uniform( + {0ull, 0ull, 0ull, 0ull}, + {0ull, 0ull, 0ull, 0ull}, + {29492327419918145ull, 4614276061115832004ull, 16925429801461668750ull, 5660986226915721659ull})); + + EXPECT_TRUE(verify_uniform( + {1ull, 2ull, 3ull, 4ull}, + {5ull, 6ull, 7ull, 8ull}, + {1161205111990270403ull, 7490910075796229492ull, 3916163298891572586ull, 624723054006169054ull})); + + EXPECT_TRUE(verify_uniform( + {6666ull, 77ull, 500ull, 8363839ull}, + {999ull, 137ull, 0xdeadf00dull, 0xdeadbeefull}, + {3717439728375325370ull, 14259638735392226729ull, 6108569366204981687ull, 8675995625794245694ull})); +} + +TEST(cbprng, normal) { + EXPECT_TRUE(verify_normal( + {29492327419918145ull, 4614276061115832004ull, 16925429801461668750ull, 5660986226915721659ull}, + {16723ll, 1664687ll, -761307ll, 1335286ll})); + + EXPECT_TRUE(verify_normal( + {1161205111990270403ull, 7490910075796229492ull, 3916163298891572586ull, 624723054006169054ull}, + {517262ll, 1238884ll, 2529374ll, 610685ll})); + + EXPECT_TRUE(verify_normal( + {3717439728375325370ull, 14259638735392226729ull, 6108569366204981687ull, 8675995625794245694ull}, + {684541ll, 215202ll, 1072054ll, -599461ll})); +} diff --git a/test/unit/test_sde.cpp b/test/unit/test_sde.cpp index da5274990a..6de99c3920 100644 --- a/test/unit/test_sde.cpp +++ b/test/unit/test_sde.cpp @@ -205,9 +205,9 @@ void test_statistics(const std::vector& ordered_samples, double mu, double si EXPECT_TRUE(chi_2_test_variance(sigma*sigma, acc.mean(), acc.variance(), acc.n())); } -// =========================== -// Recipe and global variables -// =========================== +// ================================= +// Declarations and global variables +// ================================= using namespace arb; using namespace arborio::literals; @@ -225,6 +225,10 @@ void advance_mean_reverting_stochastic_density_process2(arb_mechanism_ppack* pp) void advance_mean_reverting_stochastic_process(arb_mechanism_ppack* pp); void advance_mean_reverting_stochastic_process2(arb_mechanism_ppack* pp); +// =========================================== +// Recipe with multiple mechanism replacements +// =========================================== + // helper macro for replacing a mechanism's advance method #define REPLACE_IMPLEMENTATION(MECH) \ { \