Arbor, "noise" sources", stochastic differential equations: desire for benchmarking and validation - help required

[mgiugliano]

Dear All, dear Developers,

together with my colleagues, I committed (by a grant with the EU/ERA-NET, called Arbor-IO) to investigate and validate SDEs and "noise sources" of Arbor.

While I am waiting for the next Arbor release, where the "noise sources" will be available, I need some input and help.

May I ask you to clarify whether, upon a change in dt, the variance of the Wiener process (i.e. the term premultiplying the gauss-distributed pseudorandom number on the right hand side of each SDE, scales with the square root of dt ?

I would like to publish a short and concise document, to be anyway submitted to EIC-EU Project Officer and Reviewers, including a

• comparison between the statistical properties an Ornstein-Uhlenbeck process implemented as a single-compartment "without HH mechanism" is implemented in Arbor and it is checked against a Python naive implementation and, similarly,
• a comparison between the statistics of spike trains of a single-compartmental HH model (e.g. as in this) and those implemented in NEURON, upon external stimulation by a current-clamp fluctuating random input waveform.

The piece of information on the scaling with 'dt' and 'sqrt(dt)' are key to us at this point. If it does not scale like 'sqrt(dt)' then we won't be able to match control simulations.

Thank you in advance
Michele GIUGLIANO

[thorstenhater]

Dear Michele,

SDE support is still unreleased, but the scaling is indeed with the square root of the
time step width. The documentation will explain that fact as well. The feature likely will
be merged into master today and appear in version 0.8 next month. If you so desire,
cloning the repository and executing pip install . at its top-level will install arbor with
SDE support (after the merge).

However, we do not have built-in support for random current sources, so a workaround
will be needed: Either pre-generating a random waveform for consumption in a arbor.iclamp's
envelope or writing an NMODL based current source/point process.

Looking forward to hearing about your progress.

Best regards

[mgiugliano]

Thank you! Your comment is precious to me.

May I still ask you to clarify a couple of remaining doubts?

• Wouldn't a simple "RC" (ordinary first-order differential equation) with a Wiener process - appearing on the right hand side of its charge balance equation - qualifies to be an SDE?

• Indeed, upon reading here, I assumed that reimplementing this NEURON's stochastic mechanism would have been a piece of cake.

I had thought to define a new custom mechanism in ARBOR, incorporating a WHITE_NOISE block, as indicated in the example:

WHITE_NOISE {
       a 
   }

and then simply allowing a reimplementation of the (e.g.) Orsntein-Uhlenbeck realisation like

DERIVATIVE state {
       x' = (m - x)/tau +  s/tau * a
   }

and of course adding

   BREAKPOINT {
       SOLVE state METHOD stochastic
   }

Would this strategy work? If it does, then it answers my original question with a passive and with an active compartment.

Generating offline a waveform and then "re-injecting it" as an arbor.iclamp input would be truly cumbersome (e.g. for large networks, as well as for mimicking channel noise with a diffusion approximation (see here in distinct compartments of multiple cells).

Any hints would be gold for us.

[boeschf]

@mgiugliano To answer your question: yes this would work, and it is now merged in master. Note, that the noise is generated per synapse or per control volume, and crucially, is independently distributed. Thus, at this point it is not straight forward to introduce spatially correlated white noise. For documentation see the following links: mechanism description, nmodl reference and implementation details.
Furthermore, the seed value for the random number generator can be set using the simulation class (python API, c++ API - see also the new builder pattern for the simulation class).

[thorstenhater]

Hi Michele,

from my reading of the code and documentation, I assume an implmentation like below would be a feasible option

NEURON {
  POINT_PROCESS o_u_current
  NONSPECIFIC_CURRENT i
  RANGE mu, theta, sigma
}

PARAMETER {
  mu
  theta
  sigma
}

STATE { x }

WHITE_NOISE { W }

DERIVATIVE dX {
  x' = (mu - x)*theta +  sigma*W
}

BREAKPOINT {
  SOLVE dX METHOD stochastic
  i = x
}

and give you a random-walk-like current. It could also be possible to use W
directly in the BREAKPOINT block, if needed.

See also #1987.

[boeschf]

just to clarify, you cannot use W directly in the BREAKPOINT block, since the multiplication with $\sqrt{\Delta t}$ happens inside the code generated by solving the DERIVATIVE block. Plain W is just a normally distributed random number with mean $\mu = 0$ and variance $\sigma^2 = 1$

[llandsmeer]

Made together with @boeschf and @thorstenhater, here is a multicore demo that simulates OU current noise and outputs numpy arrays with the voltage values.

There current version runs on boeschf SDE branch as it's not merged into arbor master yet. To build from source:

git clone --depth 1 --recursive --branch=SDE 'https://github.com/boeschf/arbor'
cd arbor ; mkdir build ; cd build
cmake -DCMAKE_INSTALL_PREFIX=$(realpath prefix) -DARB_USE_BUNDLED_LIBS=ON -DARB_WITH_PYTHON=ON ..
make -j 12 && make install

Then update the install path in the ARBOR_REPOSITORY repository.
If testing after the merge has happened (and not building from source), you can comment the lines following #### IF INSTALLED FROM SOURCE:

import sys
import importlib.util
import os.path

##### IF INSTALLED USING PIP:
#ARBOR_LOCATION = 'default'
ARBOR_BUILD_CATALOGUE = 'arbor-build-catalogue'

#### IF INSTALLED FROM SOURCE:
ARBOR_REPOSITORY = '/home/llandsmeer/tmp/a-e/arbor'
ARBOR_LOCATION = os.path.join(ARBOR_REPOSITORY,
        'build/prefix/lib/python3.10/site-packages/arbor/__init__.py')
ARBOR_BUILD_CATALOGUE = os.path.join(ARBOR_REPOSITORY,
        'build/prefix/bin/arbor-build-catalogue')

if ARBOR_LOCATION != 'default':
    spec = importlib.util.spec_from_file_location('arbor', ARBOR_LOCATION)
    module = importlib.util.module_from_spec(spec)
    sys.modules[spec.name] = module
    spec.loader.exec_module(module)
# END INSTALLATION OPTIONS

import numpy as np
import re
import subprocess
import tempfile
import arbor
import matplotlib.pyplot as plt

mechanism_source = '''
NEURON {
  SUFFIX noise
  NONSPECIFIC_CURRENT i
  RANGE mu, theta, sigma
}
PARAMETER {
  mu = 0
  theta = 1
  sigma = 1
}
STATE { x }
INITIAL { x=0 }
WHITE_NOISE { W }
DERIVATIVE dX {
  x' = (mu - x)*theta +  sigma*W
}
BREAKPOINT {
  SOLVE dX METHOD stochastic
  i = x
}
'''

def build_noise_catalogue():
    with tempfile.TemporaryDirectory() as tmpdir:
        with open(os.path.join(tmpdir, 'noise.mod'), 'w') as f:
            print(mechanism_source, file=f)
        out = subprocess.getoutput(f'{ARBOR_BUILD_CATALOGUE} noise {tmpdir}')
        m = re.search('[^ ]*\.so', out)
        if m is None:
            print(out)
            exit(1)
        return arbor.load_catalogue(m.group(0))


def make_cable_cell(gid):
    tree = arbor.segment_tree()
    soma = tree.append(arbor.mnpos, arbor.mpoint(-12, 0, 0, 6), arbor.mpoint(0, 0, 0, 6), tag=1)
    labels = arbor.label_dict(dict(soma="(tag 1)", root= "(root)"))
    decor = arbor.decor()
    decor.paint('"soma"', arbor.density("noise"))
    return arbor.cable_cell(tree, labels, decor)

class Recipe(arbor.recipe):
    def __init__(self, ncells=100):
        arbor.recipe.__init__(self)
        self.ncells = ncells
        self.props = arbor.neuron_cable_properties()
        self.props.catalogue.extend(build_noise_catalogue(), '')
    def num_cells(self): return self.ncells
    def cell_description(self, gid): return make_cable_cell(gid)
    def cell_kind(self, gid): return arbor.cell_kind.cable
    def connections_on(self, gid): return []
    def event_generators(self, gid): return []
    def probes(self, gid): return [arbor.cable_probe_membrane_voltage('"root"')]

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

# (11) Instantiate recipe
recipe = Recipe(ncells=10)
ctx = arbor.context('avail_threads')
sim = arbor.simulation(recipe, ctx)
handles = [sim.sample((gid, 0), arbor.regular_schedule(1)) for gid in range(recipe.num_cells())]

sim.run(100)
times = np.array(sim.samples(handles[0])[0][0][:,0])
results = np.array([sim.samples(handles[gid])[0][0][:, 1] for gid in range(recipe.num_cells())])

plt.plot(times,results.T)
plt.show()

[mgiugliano]

@llandsmeer @boeschf and @thorstenhater, the sample code unfortunately never worked.

I get the following error on a MacBook Pro, 16" 2019, 8 cores, Ventura macOs 13.3 (22E252), home-brew Python 3.11.2 and after having pip3-installed Arbor (version 0.8.1).

`
Traceback (most recent call last):
File "/.../testme.py", line 91, in
sim = arbor.simulation(recipe, ctx)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File ".../testme.py", line 80, in cell_description
def cell_description(self, gid): return make_cable_cell(gid)
^^^^^^^^^^^^^^^^^^^^
File ".../testme.py", line 71, in make_cable_cell
return arbor.cable_cell(tree, labels, decor)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
TypeError: init(): incompatible constructor arguments. The following argument types are supported:
1. arbor._arbor.cable_cell(morphology: arb::morphology, decor: arbor._arbor.decor, labels: Optional[arbor._arbor.label_dict] = Labels)
2. arbor._arbor.cable_cell(segment_tree: arb::segment_tree, decor: arbor._arbor.decor, labels: Optional[arbor._arbor.label_dict] = Labels)

Invoked with: <arbor._arbor.segment_tree object at 0x1122b1eb0>, (label_dict (region "soma" (tag 1)) (locset "root" (root))), <arbor._arbor.decor object at 0x112111130>`

print(arbor.__config__) returns
{'mpi': False, 'mpi4py': False, 'gpu': None, 'vectorize': False, 'profiling': False, 'neuroml': True, 'bundled': True, 'version': '0.8.1', 'source': '2022-12-22T15:48:03+01:00 28129e032a9dfef45568aea6697f995b4821cbbe', 'arch': 'none', 'prefix': '/Users/runner/work/arbor/arbor/_skbuild/macosx-10.15-universal2-3.11/cmake-install', 'binary_path': 'bin', 'lib_path': 'lib', 'data_path': 'share', 'CXX': '/Applications/Xcode_13.2.1.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin/c++', 'pybind-version': '2.10.1'}

[thorstenhater]

Hi @mgiugliano,

we had a breaking change to our APIs at some point between September and now. Just swap
the arguments to the cable cell constructor as mentioned in the error:

def make_cable_cell(gid):
    tree = arbor.segment_tree()
    soma = tree.append(arbor.mnpos, arbor.mpoint(-12, 0, 0, 6), arbor.mpoint(0, 0, 0, 6), tag=1)
    labels = arbor.label_dict(dict(soma="(tag 1)", root= "(root)"))
    decor = arbor.decor()
    decor.paint('"soma"', arbor.density("noise"))
    return arbor.cable_cell(tree, decor, labels)

[mgiugliano]

Thank you! Now I get the final plot.

[mgiugliano]

For posterity, here is Google Colab notebook that runs the demo too. https://colab.research.google.com/drive/1MQpfdrpxzTqGvRp-tc8QCpk2BKNEOR-b?usp=sharing