NESTML Ornstein-Uhlenbeck noise tutorial

In this tutorial, we will formulate the Ornstein-Uhlenbeck (O-U) noise process in NESTML and simulate it in NEST Simulator.

[1]:
%matplotlib inline
import matplotlib.pyplot as plt
import nest
import numpy as np
import os
from pynestml.frontend.pynestml_frontend import generate_nest_target

NEST_SIMULATOR_INSTALL_LOCATION = nest.ll_api.sli_func("statusdict/prefix ::")

              -- N E S T --
  Copyright (C) 2004 The NEST Initiative

 Version: master@1510a5790
 Built: Mar 21 2023 14:28:29

 This program is provided AS IS and comes with
 NO WARRANTY. See the file LICENSE for details.

 Problems or suggestions?
   Visit https://www.nest-simulator.org

 Type 'nest.help()' to find out more about NEST.

The Ornstein-Uhlenbeck process

The Ornstein-Uhlenbeck process is often used as a source of noise because it is well understood and has convenient properties (it is a Gaussian process, has the Markov property, and is stationary). Let the O-U process, denoted \(U(t)\) (with \(t\geq 0\)) , be defined as the solution of the following stochastic differential equation:

\begin{align} \frac{dU}{dt} = \frac{\mu - U}{\tau} + \sigma\sqrt{\frac 2 \tau} \frac{dB(t)}{dt} \end{align}

The first right-hand side term is a “drift” term which is deterministic and slowly reverts \(U_t\) to the mean \(\mu\), with time constant \(\tau\). The second term is stochastic as \(B(t)\) is the Brownian motion (also “Wiener”) process, and \(\sigma>0\) is the standard deviation of the noise.

It turns out that the infinitesimal step in Brownian motion is white noise, that is, an independent and identically distributed sequence of Gaussian \(\mathcal{N}(0, 1)\) random variables. The noise \(dB(t)/dt\) can be sampled at time \(t\) by drawing a sample from that Gaussian distribution, so if the process is sampled at discrete intervals of length \(h\), we can write (equation 2.47 from [1]):

\begin{align} U(t + h) = (U(t) - \mu)\exp(-h/\tau) + \sigma\sqrt{(1 - \exp(-2h / \tau ))} \cdot\mathcal{N}(0, 1) \end{align}

Simulating the Ornstein-Uhlenbeck process with NESTML

Formulating the model in NESTML

The O-U process is now defined as a stepwise update equation, that is carried out at each timestep. All statements contained in the NESTML update block are run at every timestep, so this is a natural place to place the O-U updates.

Updating the O-U process state requires knowledge about how much time has elapsed. As we will update the process once every simulation timestep, we pass this value by invoking the resolution() function.

[2]:
nestml_ou_model = '''
neuron ornstein_uhlenbeck_noise:

    parameters:
        mean_noise real = 500    # mean of the noise
        sigma_noise real = 50    # std. dev. of the noise
        tau_noise ms = 20 ms     # time constant of the noise

    internals:
        A_noise real = sigma_noise * ((1 - exp(-2 * resolution() / tau_noise)))**.5

    state:
        U real = mean_noise   # set the initial condition

    update:
        U = mean_noise \
          + (U - mean_noise) * exp(-resolution() / tau_noise) \
          + A_noise * random_normal(0, 1)
'''

Save to a temporary file and make the model available to instantiate in NEST:

[3]:
with open("ornstein_uhlenbeck_noise.nestml", "w") as nestml_model_file:
    print(nestml_ou_model, file=nestml_model_file)

generate_nest_target(input_path="ornstein_uhlenbeck_noise.nestml",
                     target_path="/tmp/nestml-ou-noise-target",
                     module_name="nestml_ou_module",
                     suffix="_nestml",
                     logging_level="ERROR",
                     codegen_opts={"nest_path": NEST_SIMULATOR_INSTALL_LOCATION})
nest.Install("nestml_ou_module")
Warning: PyGSL is not available. The stiffness test will be skipped.
Warning: No module named 'pygsl'



              -- N E S T --
  Copyright (C) 2004 The NEST Initiative

 Version: master@1510a5790
 Built: Mar 21 2023 14:28:29

 This program is provided AS IS and comes with
 NO WARRANTY. See the file LICENSE for details.

 Problems or suggestions?
   Visit https://www.nest-simulator.org

 Type 'nest.help()' to find out more about NEST.

CMake Warning:
  Ignoring empty string ("") provided on the command line.


-- The CXX compiler identification is AppleClang 14.0.0.14000029
-- Detecting CXX compiler ABI info
-- Detecting CXX compiler ABI info - done
-- Check for working CXX compiler: /Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin/c++ - skipped
-- Detecting CXX compile features
-- Detecting CXX compile features - done

-------------------------------------------------------
nestml_ou_module Configuration Summary
-------------------------------------------------------

C++ compiler         : /Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin/c++
Build static libs    : OFF
C++ compiler flags   :
NEST compiler flags  : -Wextra -Wno-unknown-pragmas -D_GLIBCXX_ASSERTIONS -std=c++11 -Wall  -O2 -g
NEST include dirs    :  -I/Users/pooja/nest/include/nest -I/usr/local/include -I/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX13.1.sdk/usr/include -I/usr/local/Cellar/gsl/2.7/include -I/Users/pooja/miniconda3/envs/nestml_3.10/include -I/Users/pooja/miniconda3/envs/nestml_3.10/include
NEST libraries flags : -L/Users/pooja/nest/lib/nest -lnestutil -lsli -lnestkernel  /usr/local/lib/libltdl.dylib /Users/pooja/miniconda3/envs/nestml_3.10/lib/libreadline.dylib /Users/pooja/miniconda3/envs/nestml_3.10/lib/libncurses.dylib /usr/local/Cellar/gsl/2.7/lib/libgsl.dylib /usr/local/Cellar/gsl/2.7/lib/libgslcblas.dylib   /Users/pooja/miniconda3/envs/nestml_3.10/lib/libmpi.dylib

-------------------------------------------------------

You can now build and install 'nestml_ou_module' using
  make
  make install

The library file libnestml_ou_module.so will be installed to
  /Users/pooja/nest/lib/nest
The module can be loaded into NEST using
  (nestml_ou_module) Install       (in SLI)
  nest.Install(nestml_ou_module)   (in PyNEST)

CMake Warning (dev) in CMakeLists.txt:
  No cmake_minimum_required command is present.  A line of code such as

    cmake_minimum_required(VERSION 3.23)

  should be added at the top of the file.  The version specified may be lower
  if you wish to support older CMake versions for this project.  For more
  information run "cmake --help-policy CMP0000".
This warning is for project developers.  Use -Wno-dev to suppress it.

-- Configuring done
-- Generating done
-- Build files have been written to: /tmp/nestml-ou-noise-target
Consolidate compiler generated dependencies of target nestml_ou_module_module
[ 33%] Building CXX object CMakeFiles/nestml_ou_module_module.dir/nestml_ou_module.o
[ 66%] Building CXX object CMakeFiles/nestml_ou_module_module.dir/ornstein_uhlenbeck_noise_nestml.o
In file included from In file included from /tmp/nestml-ou-noise-target/ornstein_uhlenbeck_noise_nestml.cpp:33:
In file included from /Users/pooja/nest/include/nest/kernel_manager.h:/tmp/nestml-ou-noise-target/nestml_ou_module.cpp:26:
In file included from 27/Users/pooja/nest/include/nest/connection_manager_impl.h::
26In file included from :
In file included from /Users/pooja/nest/include/nest/connection_manager.h:/Users/pooja/nest/include/nest/connection_manager.h36::
36In file included from :
/Users/pooja/nest/include/nest/connector_base.hIn file included from :/Users/pooja/nest/include/nest/connector_base.h:35:
In file included from /Users/pooja/nest/include/nest/sort.h:35:
37:
In file included from In file included from /usr/local/include/boost/sort/spreadsort/spreadsort.hpp/Users/pooja/nest/include/nest/sort.h:37:27:
In file included from :
/usr/local/include/boost/sort/spreadsort/integer_sort.hpp:26In file included from :
/usr/local/include/boost/sort/spreadsort/spreadsort.hpp:27:
In file included from /usr/local/include/boost/sort/spreadsort/detail/integer_sort.hpp/usr/local/include/boost/sort/spreadsort/integer_sort.hpp::43826:29:
: warning: unused parameter 'shift' [-Wunused-parameter]
/usr/local/include/boost/sort/spreadsort/detail/integer_sort.hpp:438:29: warning: unused parameter 'shift' [-Wunused-parameter]
                Right_shift shift, Compare comp)
                            ^
                Right_shift shift, Compare comp)
                            ^
/usr/local/include/boost/sort/spreadsort/detail/integer_sort.hpp/usr/local/include/boost/sort/spreadsort/detail/integer_sort.hpp::482482:29: warning: unused parameter 'shift' [-Wunused-parameter]:29
: warning: unused parameter 'shift' [-Wunused-parameter]
                Right_shift shift)
                            ^
                Right_shift shift)
                            ^
In file included from /tmp/nestml-ou-noise-target/nestml_ou_module.cppIn file included from /tmp/nestml-ou-noise-target/ornstein_uhlenbeck_noise_nestml.cpp:33:
In file included from /Users/pooja/nest/include/nest/kernel_manager.h::2726:
:
In file included from In file included from /Users/pooja/nest/include/nest/connection_manager_impl.h/Users/pooja/nest/include/nest/connection_manager.h:36:
In file included from /Users/pooja/nest/include/nest/connector_base.h::3526:
:
In file included from In file included from /Users/pooja/nest/include/nest/sort.h/Users/pooja/nest/include/nest/connection_manager.h::3736:
:
In file included from In file included from /Users/pooja/nest/include/nest/connector_base.h/usr/local/include/boost/sort/spreadsort/spreadsort.hpp:35:
In file included from /Users/pooja/nest/include/nest/sort.h:37:
In file included from /usr/local/include/boost/sort/spreadsort/spreadsort.hpp:28:
:In file included from 28/usr/local/include/boost/sort/spreadsort/float_sort.hpp:
:In file included from 25/usr/local/include/boost/sort/spreadsort/float_sort.hpp:25:
/usr/local/include/boost/sort/spreadsort/detail/float_sort.hpp:778:28: warning: unused parameter 'rshift' [-Wunused-parameter]
:
/usr/local/include/boost/sort/spreadsort/detail/float_sort.hpp:778:28: warning: unused parameter 'rshift' [-Wunused-parameter]
               Right_shift rshift)
                           ^
               Right_shift rshift)
                           ^
/usr/local/include/boost/sort/spreadsort/detail/float_sort.hpp:820:28: warning: unused parameter 'rshift' [-Wunused-parameter]
               Right_shift rshift, Compare comp)
                           ^
/usr/local/include/boost/sort/spreadsort/detail/float_sort.hpp:820:28: warning: unused parameter 'rshift' [-Wunused-parameter]
               Right_shift rshift, Compare comp)
                           ^
In file included from /tmp/nestml-ou-noise-target/nestml_ou_module.cpp:26:
In file included from /Users/pooja/nest/include/nest/connection_manager_impl.h:26:
In file included from In file included from /tmp/nestml-ou-noise-target/ornstein_uhlenbeck_noise_nestml.cpp:33:
In file included from /Users/pooja/nest/include/nest/connection_manager.h/Users/pooja/nest/include/nest/kernel_manager.h::3627:
:
In file included from In file included from /Users/pooja/nest/include/nest/connection_manager.h/Users/pooja/nest/include/nest/connector_base.h::3635:
:
In file included from In file included from /Users/pooja/nest/include/nest/connector_base.h/Users/pooja/nest/include/nest/sort.h::3537:
:
In file included from In file included from /Users/pooja/nest/include/nest/sort.h/usr/local/include/boost/sort/spreadsort/spreadsort.hpp::3729:
:
In file included from In file included from /usr/local/include/boost/sort/spreadsort/spreadsort.hpp/usr/local/include/boost/sort/spreadsort/string_sort.hpp::2923:
:
In file included from /usr/local/include/boost/sort/spreadsort/string_sort.hpp:23/usr/local/include/boost/sort/spreadsort/detail/string_sort.hpp:
:/usr/local/include/boost/sort/spreadsort/detail/string_sort.hpp:752752::2626::  warning: warning: unused parameter 'get_character' [-Wunused-parameter]unused parameter 'get_character' [-Wunused-parameter]

                Get_char get_character, Get_length length, Unsigned_char_type)
                Get_char get_character, Get_length length, Unsigned_char_type)
                         ^
                         ^/usr/local/include/boost/sort/spreadsort/detail/string_sort.hpp:752:52: warning: unused parameter 'length' [-Wunused-parameter]
                Get_char get_character, Get_length length, Unsigned_char_type)
                                                   ^

/usr/local/include/boost/sort/spreadsort/detail/string_sort.hpp:752:52: warning: unused parameter 'length' [-Wunused-parameter]
                Get_char get_character, Get_length length, Unsigned_char_type)
                                                   ^
/usr/local/include/boost/sort/spreadsort/detail/string_sort.hpp:780:18: warning: unused parameter 'get_character' [-Wunused-parameter]
        Get_char get_character, Get_length length, Compare comp, Unsigned_char_type)
                 ^
/usr/local/include/boost/sort/spreadsort/detail/string_sort.hpp:780:44: warning: unused parameter 'length' [-Wunused-parameter]
/usr/local/include/boost/sort/spreadsort/detail/string_sort.hpp:780:18:        Get_char get_character, Get_length length, Compare comp, Unsigned_char_type)
warning:                                            ^unused parameter 'get_character' [-Wunused-parameter]

        Get_char get_character, Get_length length, Compare comp, Unsigned_char_type)
                 ^
/usr/local/include/boost/sort/spreadsort/detail/string_sort.hpp:780:44: warning: unused parameter 'length' [-Wunused-parameter]
        Get_char get_character, Get_length length, Compare comp, Unsigned_char_type)
                                           ^
/usr/local/include/boost/sort/spreadsort/detail/string_sort.hpp:807:18: warning: unused parameter 'get_character' [-Wunused-parameter]
        Get_char get_character, Get_length length, Compare comp, Unsigned_char_type)
                 ^
/usr/local/include/boost/sort/spreadsort/detail/string_sort.hpp:807:44: warning: unused parameter 'length' [-Wunused-parameter]
        Get_char get_character, Get_length length, Compare comp, Unsigned_char_type)
                                           ^
/usr/local/include/boost/sort/spreadsort/detail/string_sort.hpp:807:18: warning: unused parameter 'get_character' [-Wunused-parameter]
        Get_char get_character, Get_length length, Compare comp, Unsigned_char_type)
                 ^
/usr/local/include/boost/sort/spreadsort/detail/string_sort.hpp:807:44: warning: unused parameter 'length' [-Wunused-parameter]
        Get_char get_character, Get_length length, Compare comp, Unsigned_char_type)
                                           ^
/tmp/nestml-ou-noise-target/nestml_ou_module.cpp:100:41: warning: unused parameter 'i' [-Wunused-parameter]
nestml_ou_module::init( SLIInterpreter* i )
                                        ^
/tmp/nestml-ou-noise-target/ornstein_uhlenbeck_noise_nestml.cpp:99:16: warning: unused variable '__resolution' [-Wunused-variable]
  const double __resolution = nest::Time::get_resolution().get_ms();  // do not remove, this is necessary for the resolution() function
               ^
11 warnings generated.
11 warnings generated.
[100%] Linking CXX shared module nestml_ou_module.so
[100%] Built target nestml_ou_module_module
Consolidate compiler generated dependencies of target nestml_ou_module_module
[100%] Built target nestml_ou_module_module
Install the project...
-- Install configuration: ""
-- Installing: /Users/pooja/nest/lib/nest/nestml_ou_module.so

Mar 23 12:27:41 Install [Info]:
    loaded module nestml_ou_module

Running the simulation in NEST

Let’s define a function that will instantiate and set parameters for the O-U model, run a simulation, and plot and return the results.

[4]:
def evaluate_ou_process(h: float=.1, t_sim:float=100., neuron_parms=None, title=None, plot=True):
    """
    h : float
        timestep in ms
    t_sim : float
        total simulation time in ms
    """
    nest.ResetKernel()
    nest.SetKernelStatus({"resolution": h})
    neuron = nest.Create("ornstein_uhlenbeck_noise_nestml")

    if neuron_parms:
        for k, v in neuron_parms.items():
            nest.SetStatus(neuron, k, v)

    multimeter = nest.Create("multimeter")
    multimeter.set({"record_from": ["U"],
                    "interval": h})
    nest.Connect(multimeter, neuron)

    nest.Simulate(t_sim)

    dmm = nest.GetStatus(multimeter)[0]
    U = dmm["events"]["U"]
    timevec = dmm["events"]["times"]

    if plot:
        fig, ax = plt.subplots(figsize=(12., 5.))
        if title is not None:
            fig.suptitle(title)
        ax.plot(timevec, U)
        ax.set_xlabel("Time [ms]")
        ax.set_ylabel("U")
        ax.set_xlim(0., t_sim)
        ax.grid()

    return timevec, U

Then, we can run the process. Hint: play with the parameters a bit here and see the effects it has on the returned timeseries.

[5]:
timevec, U = evaluate_ou_process(h=1.,
                                 t_sim=1000.,
                                 neuron_parms={"U" : -2500.,
                                               "mean_noise": -3333.,
                                               "tau_noise": 20.,
                                               "sigma_noise": 100.},
                                title=r"Ornstein-Uhlenbeck process")

Mar 23 12:27:47 correlation_detector [Info]:
    Default for delta_tau changed from 0.5 to 5 ms

Mar 23 12:27:47 correlomatrix_detector [Info]:
    Default for delta_tau changed from 0.5 to 5 ms

Mar 23 12:27:47 correlospinmatrix_detector [Info]:
    Default for delta_tau changed from 0.1 to 1 ms

Mar 23 12:27:47 SimulationManager::set_status [Info]:
    Temporal resolution changed from 0.1 to 1 ms.

Mar 23 12:27:47 NodeManager::prepare_nodes [Info]:
    Preparing 2 nodes for simulation.

Mar 23 12:27:47 SimulationManager::start_updating_ [Info]:
    Number of local nodes: 2
    Simulation time (ms): 1000
    Not using OpenMP
    Number of MPI processes: 1

Mar 23 12:27:47 SimulationManager::run [Info]:
    Simulation finished.
../../_images/tutorials_ornstein_uhlenbeck_noise_nestml_ou_noise_tutorial_10_1.png

Correctness test based on predicting variance

Assuming that the initial value of the process is picked equal to the process mean, we can predict the variance of the timeseries as ([1], eq. 2.26):

\begin{align} \text{Var}(U) = \sigma^2 \end{align}

We now run a consistency check across parameters, to make sure that the prediction matches the result derived from the timeseries generated by sampling the process.

[6]:
h = [.01, .1, 1.]
tau_noise = [10., 100., 1000.]
sigma_noise = [0., 10., 100., 1000.]

max_rel_error = .25   # yes, this is pretty terrible!
                      # Need a very long integration time (esp. for large tau)
                      # to get a tighter bound.

for _h in h:
    for _tau_noise in tau_noise:
        for _sigma_noise in sigma_noise:
            print("For h = " + str(_h) + ", tau_noise = " + str(_tau_noise) + ", sigma_noise = " + str(_sigma_noise))
            c = (_sigma_noise * np.sqrt(2 / _tau_noise))**2
            timevec, U = evaluate_ou_process(h=_h,
                                             t_sim=25000.,
                                             neuron_parms={"U" : 0.,
                                                           "mean_noise": 0.,
                                                           "tau_noise": _tau_noise,
                                                           "sigma_noise": _sigma_noise},
                                            title=r"Ornstein-Uhlenbeck process ($\tau$=" + str(_tau_noise) + ", $\sigma$=" + str(_sigma_noise) + ")",
                                            plot=False)
            var_actual = np.var(U)
            print("Actual variance: " + str(var_actual))
            var_expected = _sigma_noise**2
            print("Expected variance: " + str(var_expected))
            if var_actual < 1E-15:
                assert var_expected < 1E-15
            else:
                assert np.abs(var_expected - var_actual) / (var_expected + var_actual) < max_rel_error
For h = 0.01, tau_noise = 10.0, sigma_noise = 0.0

Mar 23 12:27:50 SimulationManager::set_status [Info]:
    Temporal resolution changed from 0.1 to 0.01 ms.

Mar 23 12:27:50 NodeManager::prepare_nodes [Info]:
    Preparing 2 nodes for simulation.

Mar 23 12:27:50 SimulationManager::start_updating_ [Info]:
    Number of local nodes: 2
    Simulation time (ms): 25000
    Not using OpenMP
    Number of MPI processes: 1

Mar 23 12:27:51 SimulationManager::run [Info]:
    Simulation finished.
Actual variance: 0.0
Expected variance: 0.0
For h = 0.01, tau_noise = 10.0, sigma_noise = 10.0

Mar 23 12:27:51 SimulationManager::set_status [Info]:
    Temporal resolution changed from 0.1 to 0.01 ms.

Mar 23 12:27:51 NodeManager::prepare_nodes [Info]:
    Preparing 2 nodes for simulation.

Mar 23 12:27:51 SimulationManager::start_updating_ [Info]:
    Number of local nodes: 2
    Simulation time (ms): 25000
    Not using OpenMP
    Number of MPI processes: 1

Mar 23 12:27:51 SimulationManager::run [Info]:
    Simulation finished.
Actual variance: 102.21030336065829
Expected variance: 100.0
For h = 0.01, tau_noise = 10.0, sigma_noise = 100.0

Mar 23 12:27:51 SimulationManager::set_status [Info]:
    Temporal resolution changed from 0.1 to 0.01 ms.

Mar 23 12:27:51 NodeManager::prepare_nodes [Info]:
    Preparing 2 nodes for simulation.

Mar 23 12:27:51 SimulationManager::start_updating_ [Info]:
    Number of local nodes: 2
    Simulation time (ms): 25000
    Not using OpenMP
    Number of MPI processes: 1

Mar 23 12:27:51 SimulationManager::run [Info]:
    Simulation finished.
Actual variance: 10221.030336065825
Expected variance: 10000.0
For h = 0.01, tau_noise = 10.0, sigma_noise = 1000.0

Mar 23 12:27:51 SimulationManager::set_status [Info]:
    Temporal resolution changed from 0.1 to 0.01 ms.

Mar 23 12:27:51 NodeManager::prepare_nodes [Info]:
    Preparing 2 nodes for simulation.

Mar 23 12:27:51 SimulationManager::start_updating_ [Info]:
    Number of local nodes: 2
    Simulation time (ms): 25000
    Not using OpenMP
    Number of MPI processes: 1

Mar 23 12:27:52 SimulationManager::run [Info]:
    Simulation finished.
Actual variance: 1022103.0336065828
Expected variance: 1000000.0
For h = 0.01, tau_noise = 100.0, sigma_noise = 0.0

Mar 23 12:27:52 SimulationManager::set_status [Info]:
    Temporal resolution changed from 0.1 to 0.01 ms.

Mar 23 12:27:52 NodeManager::prepare_nodes [Info]:
    Preparing 2 nodes for simulation.

Mar 23 12:27:52 SimulationManager::start_updating_ [Info]:
    Number of local nodes: 2
    Simulation time (ms): 25000
    Not using OpenMP
    Number of MPI processes: 1

Mar 23 12:27:52 SimulationManager::run [Info]:
    Simulation finished.
Actual variance: 0.0
Expected variance: 0.0
For h = 0.01, tau_noise = 100.0, sigma_noise = 10.0

Mar 23 12:27:52 SimulationManager::set_status [Info]:
    Temporal resolution changed from 0.1 to 0.01 ms.

Mar 23 12:27:52 NodeManager::prepare_nodes [Info]:
    Preparing 2 nodes for simulation.

Mar 23 12:27:52 SimulationManager::start_updating_ [Info]:
    Number of local nodes: 2
    Simulation time (ms): 25000
    Not using OpenMP
    Number of MPI processes: 1

Mar 23 12:27:53 SimulationManager::run [Info]:
    Simulation finished.
Actual variance: 114.93624532212093
Expected variance: 100.0
For h = 0.01, tau_noise = 100.0, sigma_noise = 100.0

Mar 23 12:27:53 SimulationManager::set_status [Info]:
    Temporal resolution changed from 0.1 to 0.01 ms.

Mar 23 12:27:53 NodeManager::prepare_nodes [Info]:
    Preparing 2 nodes for simulation.

Mar 23 12:27:53 SimulationManager::start_updating_ [Info]:
    Number of local nodes: 2
    Simulation time (ms): 25000
    Not using OpenMP
    Number of MPI processes: 1

Mar 23 12:27:53 SimulationManager::run [Info]:
    Simulation finished.
Actual variance: 11493.62453221208
Expected variance: 10000.0
For h = 0.01, tau_noise = 100.0, sigma_noise = 1000.0

Mar 23 12:27:53 SimulationManager::set_status [Info]:
    Temporal resolution changed from 0.1 to 0.01 ms.

Mar 23 12:27:53 NodeManager::prepare_nodes [Info]:
    Preparing 2 nodes for simulation.

Mar 23 12:27:53 SimulationManager::start_updating_ [Info]:
    Number of local nodes: 2
    Simulation time (ms): 25000
    Not using OpenMP
    Number of MPI processes: 1

Mar 23 12:27:53 SimulationManager::run [Info]:
    Simulation finished.
Actual variance: 1149362.453221211
Expected variance: 1000000.0
For h = 0.01, tau_noise = 1000.0, sigma_noise = 0.0

Mar 23 12:27:54 SimulationManager::set_status [Info]:
    Temporal resolution changed from 0.1 to 0.01 ms.

Mar 23 12:27:54 NodeManager::prepare_nodes [Info]:
    Preparing 2 nodes for simulation.

Mar 23 12:27:54 SimulationManager::start_updating_ [Info]:
    Number of local nodes: 2
    Simulation time (ms): 25000
    Not using OpenMP
    Number of MPI processes: 1

Mar 23 12:27:54 SimulationManager::run [Info]:
    Simulation finished.
Actual variance: 0.0
Expected variance: 0.0
For h = 0.01, tau_noise = 1000.0, sigma_noise = 10.0

Mar 23 12:27:54 SimulationManager::set_status [Info]:
    Temporal resolution changed from 0.1 to 0.01 ms.

Mar 23 12:27:54 NodeManager::prepare_nodes [Info]:
    Preparing 2 nodes for simulation.

Mar 23 12:27:54 SimulationManager::start_updating_ [Info]:
    Number of local nodes: 2
    Simulation time (ms): 25000
    Not using OpenMP
    Number of MPI processes: 1

Mar 23 12:27:54 SimulationManager::run [Info]:
    Simulation finished.
Actual variance: 94.22834590929104
Expected variance: 100.0
For h = 0.01, tau_noise = 1000.0, sigma_noise = 100.0

Mar 23 12:27:54 SimulationManager::set_status [Info]:
    Temporal resolution changed from 0.1 to 0.01 ms.

Mar 23 12:27:54 NodeManager::prepare_nodes [Info]:
    Preparing 2 nodes for simulation.

Mar 23 12:27:54 SimulationManager::start_updating_ [Info]:
    Number of local nodes: 2
    Simulation time (ms): 25000
    Not using OpenMP
    Number of MPI processes: 1

Mar 23 12:27:55 SimulationManager::run [Info]:
    Simulation finished.
Actual variance: 9422.834590928938
Expected variance: 10000.0
For h = 0.01, tau_noise = 1000.0, sigma_noise = 1000.0

Mar 23 12:27:55 SimulationManager::set_status [Info]:
    Temporal resolution changed from 0.1 to 0.01 ms.

Mar 23 12:27:55 NodeManager::prepare_nodes [Info]:
    Preparing 2 nodes for simulation.

Mar 23 12:27:55 SimulationManager::start_updating_ [Info]:
    Number of local nodes: 2
    Simulation time (ms): 25000
    Not using OpenMP
    Number of MPI processes: 1

Mar 23 12:27:55 SimulationManager::run [Info]:
    Simulation finished.
Actual variance: 942283.4590929038
Expected variance: 1000000.0
For h = 0.1, tau_noise = 10.0, sigma_noise = 0.0

Mar 23 12:27:55 SimulationManager::set_status [Info]:
    Temporal resolution changed from 0.1 to 0.1 ms.

Mar 23 12:27:55 NodeManager::prepare_nodes [Info]:
    Preparing 2 nodes for simulation.

Mar 23 12:27:55 SimulationManager::start_updating_ [Info]:
    Number of local nodes: 2
    Simulation time (ms): 25000
    Not using OpenMP
    Number of MPI processes: 1

Mar 23 12:27:55 SimulationManager::run [Info]:
    Simulation finished.
Actual variance: 0.0
Expected variance: 0.0
For h = 0.1, tau_noise = 10.0, sigma_noise = 10.0

Mar 23 12:27:55 SimulationManager::set_status [Info]:
    Temporal resolution changed from 0.1 to 0.1 ms.

Mar 23 12:27:55 NodeManager::prepare_nodes [Info]:
    Preparing 2 nodes for simulation.

Mar 23 12:27:55 SimulationManager::start_updating_ [Info]:
    Number of local nodes: 2
    Simulation time (ms): 25000
    Not using OpenMP
    Number of MPI processes: 1

Mar 23 12:27:55 SimulationManager::run [Info]:
    Simulation finished.
Actual variance: 99.44588939205474
Expected variance: 100.0
For h = 0.1, tau_noise = 10.0, sigma_noise = 100.0

Mar 23 12:27:55 SimulationManager::set_status [Info]:
    Temporal resolution changed from 0.1 to 0.1 ms.

Mar 23 12:27:55 NodeManager::prepare_nodes [Info]:
    Preparing 2 nodes for simulation.

Mar 23 12:27:55 SimulationManager::start_updating_ [Info]:
    Number of local nodes: 2
    Simulation time (ms): 25000
    Not using OpenMP
    Number of MPI processes: 1

Mar 23 12:27:55 SimulationManager::run [Info]:
    Simulation finished.
Actual variance: 9944.588939205476
Expected variance: 10000.0
For h = 0.1, tau_noise = 10.0, sigma_noise = 1000.0

Mar 23 12:27:55 SimulationManager::set_status [Info]:
    Temporal resolution changed from 0.1 to 0.1 ms.

Mar 23 12:27:55 NodeManager::prepare_nodes [Info]:
    Preparing 2 nodes for simulation.

Mar 23 12:27:55 SimulationManager::start_updating_ [Info]:
    Number of local nodes: 2
    Simulation time (ms): 25000
    Not using OpenMP
    Number of MPI processes: 1

Mar 23 12:27:55 SimulationManager::run [Info]:
    Simulation finished.
Actual variance: 994458.8939205472
Expected variance: 1000000.0
For h = 0.1, tau_noise = 100.0, sigma_noise = 0.0

Mar 23 12:27:55 SimulationManager::set_status [Info]:
    Temporal resolution changed from 0.1 to 0.1 ms.

Mar 23 12:27:55 NodeManager::prepare_nodes [Info]:
    Preparing 2 nodes for simulation.

Mar 23 12:27:55 SimulationManager::start_updating_ [Info]:
    Number of local nodes: 2
    Simulation time (ms): 25000
    Not using OpenMP
    Number of MPI processes: 1

Mar 23 12:27:56 SimulationManager::run [Info]:
    Simulation finished.
Actual variance: 0.0
Expected variance: 0.0
For h = 0.1, tau_noise = 100.0, sigma_noise = 10.0

Mar 23 12:27:56 SimulationManager::set_status [Info]:
    Temporal resolution changed from 0.1 to 0.1 ms.

Mar 23 12:27:56 NodeManager::prepare_nodes [Info]:
    Preparing 2 nodes for simulation.

Mar 23 12:27:56 SimulationManager::start_updating_ [Info]:
    Number of local nodes: 2
    Simulation time (ms): 25000
    Not using OpenMP
    Number of MPI processes: 1

Mar 23 12:27:56 SimulationManager::run [Info]:
    Simulation finished.
Actual variance: 96.98034018669549
Expected variance: 100.0
For h = 0.1, tau_noise = 100.0, sigma_noise = 100.0

Mar 23 12:27:56 SimulationManager::set_status [Info]:
    Temporal resolution changed from 0.1 to 0.1 ms.

Mar 23 12:27:56 NodeManager::prepare_nodes [Info]:
    Preparing 2 nodes for simulation.

Mar 23 12:27:56 SimulationManager::start_updating_ [Info]:
    Number of local nodes: 2
    Simulation time (ms): 25000
    Not using OpenMP
    Number of MPI processes: 1

Mar 23 12:27:56 SimulationManager::run [Info]:
    Simulation finished.
Actual variance: 9698.034018669543
Expected variance: 10000.0
For h = 0.1, tau_noise = 100.0, sigma_noise = 1000.0

Mar 23 12:27:56 SimulationManager::set_status [Info]:
    Temporal resolution changed from 0.1 to 0.1 ms.

Mar 23 12:27:56 NodeManager::prepare_nodes [Info]:
    Preparing 2 nodes for simulation.

Mar 23 12:27:56 SimulationManager::start_updating_ [Info]:
    Number of local nodes: 2
    Simulation time (ms): 25000
    Not using OpenMP
    Number of MPI processes: 1

Mar 23 12:27:56 SimulationManager::run [Info]:
    Simulation finished.
Actual variance: 969803.4018669543
Expected variance: 1000000.0
For h = 0.1, tau_noise = 1000.0, sigma_noise = 0.0

Mar 23 12:27:56 SimulationManager::set_status [Info]:
    Temporal resolution changed from 0.1 to 0.1 ms.

Mar 23 12:27:56 NodeManager::prepare_nodes [Info]:
    Preparing 2 nodes for simulation.

Mar 23 12:27:56 SimulationManager::start_updating_ [Info]:
    Number of local nodes: 2
    Simulation time (ms): 25000
    Not using OpenMP
    Number of MPI processes: 1

Mar 23 12:27:56 SimulationManager::run [Info]:
    Simulation finished.
Actual variance: 0.0
Expected variance: 0.0
For h = 0.1, tau_noise = 1000.0, sigma_noise = 10.0

Mar 23 12:27:56 SimulationManager::set_status [Info]:
    Temporal resolution changed from 0.1 to 0.1 ms.

Mar 23 12:27:56 NodeManager::prepare_nodes [Info]:
    Preparing 2 nodes for simulation.

Mar 23 12:27:56 SimulationManager::start_updating_ [Info]:
    Number of local nodes: 2
    Simulation time (ms): 25000
    Not using OpenMP
    Number of MPI processes: 1

Mar 23 12:27:56 SimulationManager::run [Info]:
    Simulation finished.
Actual variance: 137.92623396537726
Expected variance: 100.0
For h = 0.1, tau_noise = 1000.0, sigma_noise = 100.0

Mar 23 12:27:56 SimulationManager::set_status [Info]:
    Temporal resolution changed from 0.1 to 0.1 ms.

Mar 23 12:27:56 NodeManager::prepare_nodes [Info]:
    Preparing 2 nodes for simulation.

Mar 23 12:27:56 SimulationManager::start_updating_ [Info]:
    Number of local nodes: 2
    Simulation time (ms): 25000
    Not using OpenMP
    Number of MPI processes: 1

Mar 23 12:27:56 SimulationManager::run [Info]:
    Simulation finished.
Actual variance: 13792.623396537845
Expected variance: 10000.0
For h = 0.1, tau_noise = 1000.0, sigma_noise = 1000.0

Mar 23 12:27:56 SimulationManager::set_status [Info]:
    Temporal resolution changed from 0.1 to 0.1 ms.

Mar 23 12:27:56 NodeManager::prepare_nodes [Info]:
    Preparing 2 nodes for simulation.

Mar 23 12:27:56 SimulationManager::start_updating_ [Info]:
    Number of local nodes: 2
    Simulation time (ms): 25000
    Not using OpenMP
    Number of MPI processes: 1

Mar 23 12:27:56 SimulationManager::run [Info]:
    Simulation finished.
Actual variance: 1379262.33965378
Expected variance: 1000000.0
For h = 1.0, tau_noise = 10.0, sigma_noise = 0.0

Mar 23 12:27:56 correlation_detector [Info]:
    Default for delta_tau changed from 0.5 to 5 ms

Mar 23 12:27:56 correlomatrix_detector [Info]:
    Default for delta_tau changed from 0.5 to 5 ms

Mar 23 12:27:56 correlospinmatrix_detector [Info]:
    Default for delta_tau changed from 0.1 to 1 ms

Mar 23 12:27:56 SimulationManager::set_status [Info]:
    Temporal resolution changed from 0.1 to 1 ms.

Mar 23 12:27:56 NodeManager::prepare_nodes [Info]:
    Preparing 2 nodes for simulation.

Mar 23 12:27:56 SimulationManager::start_updating_ [Info]:
    Number of local nodes: 2
    Simulation time (ms): 25000
    Not using OpenMP
    Number of MPI processes: 1

Mar 23 12:27:56 SimulationManager::run [Info]:
    Simulation finished.
Actual variance: 0.0
Expected variance: 0.0
For h = 1.0, tau_noise = 10.0, sigma_noise = 10.0

Mar 23 12:27:56 correlation_detector [Info]:
    Default for delta_tau changed from 0.5 to 5 ms

Mar 23 12:27:56 correlomatrix_detector [Info]:
    Default for delta_tau changed from 0.5 to 5 ms

Mar 23 12:27:56 correlospinmatrix_detector [Info]:
    Default for delta_tau changed from 0.1 to 1 ms

Mar 23 12:27:56 SimulationManager::set_status [Info]:
    Temporal resolution changed from 0.1 to 1 ms.

Mar 23 12:27:56 NodeManager::prepare_nodes [Info]:
    Preparing 2 nodes for simulation.

Mar 23 12:27:56 SimulationManager::start_updating_ [Info]:
    Number of local nodes: 2
    Simulation time (ms): 25000
    Not using OpenMP
    Number of MPI processes: 1

Mar 23 12:27:56 SimulationManager::run [Info]:
    Simulation finished.
Actual variance: 102.32554029453819
Expected variance: 100.0
For h = 1.0, tau_noise = 10.0, sigma_noise = 100.0

Mar 23 12:27:56 correlation_detector [Info]:
    Default for delta_tau changed from 0.5 to 5 ms

Mar 23 12:27:56 correlomatrix_detector [Info]:
    Default for delta_tau changed from 0.5 to 5 ms

Mar 23 12:27:56 correlospinmatrix_detector [Info]:
    Default for delta_tau changed from 0.1 to 1 ms

Mar 23 12:27:56 SimulationManager::set_status [Info]:
    Temporal resolution changed from 0.1 to 1 ms.

Mar 23 12:27:56 NodeManager::prepare_nodes [Info]:
    Preparing 2 nodes for simulation.

Mar 23 12:27:56 SimulationManager::start_updating_ [Info]:
    Number of local nodes: 2
    Simulation time (ms): 25000
    Not using OpenMP
    Number of MPI processes: 1

Mar 23 12:27:56 SimulationManager::run [Info]:
    Simulation finished.
Actual variance: 10232.554029453817
Expected variance: 10000.0
For h = 1.0, tau_noise = 10.0, sigma_noise = 1000.0

Mar 23 12:27:56 correlation_detector [Info]:
    Default for delta_tau changed from 0.5 to 5 ms

Mar 23 12:27:56 correlomatrix_detector [Info]:
    Default for delta_tau changed from 0.5 to 5 ms

Mar 23 12:27:56 correlospinmatrix_detector [Info]:
    Default for delta_tau changed from 0.1 to 1 ms
Actual variance: 1023255.4029453819
Expected variance: 1000000.0
For h = 1.0, tau_noise = 100.0, sigma_noise = 0.0

Mar 23 12:27:56 SimulationManager::set_status [Info]:
    Temporal resolution changed from 0.1 to 1 ms.

Mar 23 12:27:56 NodeManager::prepare_nodes [Info]:
    Preparing 2 nodes for simulation.

Mar 23 12:27:56 SimulationManager::start_updating_ [Info]:
    Number of local nodes: 2
    Simulation time (ms): 25000
    Not using OpenMP
    Number of MPI processes: 1

Mar 23 12:27:56 SimulationManager::run [Info]:
    Simulation finished.

Mar 23 12:27:56 correlation_detector [Info]:
    Default for delta_tau changed from 0.5 to 5 ms

Mar 23 12:27:56 correlomatrix_detector [Info]:
    Default for delta_tau changed from 0.5 to 5 ms

Mar 23 12:27:56 correlospinmatrix_detector [Info]:
    Default for delta_tau changed from 0.1 to 1 ms

Mar 23 12:27:56 SimulationManager::set_status [Info]:
    Temporal resolution changed from 0.1 to 1 ms.

Mar 23 12:27:56 NodeManager::prepare_nodes [Info]:
    Preparing 2 nodes for simulation.

Mar 23 12:27:56 SimulationManager::start_updating_ [Info]:
    Number of local nodes: 2
    Simulation time (ms): 25000
    Not using OpenMP
    Number of MPI processes: 1

Mar 23 12:27:56 SimulationManager::run [Info]:
    Simulation finished.
Actual variance: 0.0
Expected variance: 0.0
For h = 1.0, tau_noise = 100.0, sigma_noise = 10.0

Mar 23 12:27:56 correlation_detector [Info]:
    Default for delta_tau changed from 0.5 to 5 ms

Mar 23 12:27:56 correlomatrix_detector [Info]:
    Default for delta_tau changed from 0.5 to 5 ms

Mar 23 12:27:56 correlospinmatrix_detector [Info]:
    Default for delta_tau changed from 0.1 to 1 ms

Mar 23 12:27:56 SimulationManager::set_status [Info]:
    Temporal resolution changed from 0.1 to 1 ms.

Mar 23 12:27:56 NodeManager::prepare_nodes [Info]:
    Preparing 2 nodes for simulation.

Mar 23 12:27:56 SimulationManager::start_updating_ [Info]:
    Number of local nodes: 2
    Simulation time (ms): 25000
    Not using OpenMP
    Number of MPI processes: 1

Mar 23 12:27:56 SimulationManager::run [Info]:
    Simulation finished.
Actual variance: 102.5676581333639
Expected variance: 100.0
For h = 1.0, tau_noise = 100.0, sigma_noise = 100.0

Mar 23 12:27:56 correlation_detector [Info]:
    Default for delta_tau changed from 0.5 to 5 ms

Mar 23 12:27:56 correlomatrix_detector [Info]:
    Default for delta_tau changed from 0.5 to 5 ms

Mar 23 12:27:56 correlospinmatrix_detector [Info]:
    Default for delta_tau changed from 0.1 to 1 ms

Mar 23 12:27:56 SimulationManager::set_status [Info]:
    Temporal resolution changed from 0.1 to 1 ms.

Mar 23 12:27:56 NodeManager::prepare_nodes [Info]:
    Preparing 2 nodes for simulation.

Mar 23 12:27:56 SimulationManager::start_updating_ [Info]:
    Number of local nodes: 2
    Simulation time (ms): 25000
    Not using OpenMP
    Number of MPI processes: 1

Mar 23 12:27:56 SimulationManager::run [Info]:
    Simulation finished.
Actual variance: 10256.76581333639
Expected variance: 10000.0
For h = 1.0, tau_noise = 100.0, sigma_noise = 1000.0

Mar 23 12:27:56 correlation_detector [Info]:
    Default for delta_tau changed from 0.5 to 5 ms

Mar 23 12:27:56 correlomatrix_detector [Info]:
    Default for delta_tau changed from 0.5 to 5 ms

Mar 23 12:27:56 correlospinmatrix_detector [Info]:
    Default for delta_tau changed from 0.1 to 1 ms

Mar 23 12:27:56 SimulationManager::set_status [Info]:
    Temporal resolution changed from 0.1 to 1 ms.

Mar 23 12:27:56 NodeManager::prepare_nodes [Info]:
    Preparing 2 nodes for simulation.

Mar 23 12:27:56 SimulationManager::start_updating_ [Info]:
    Number of local nodes: 2
    Simulation time (ms): 25000
    Not using OpenMP
    Number of MPI processes: 1

Mar 23 12:27:56 SimulationManager::run [Info]:
    Simulation finished.
Actual variance: 1025676.581333639
Expected variance: 1000000.0
For h = 1.0, tau_noise = 1000.0, sigma_noise = 0.0

Mar 23 12:27:56 correlation_detector [Info]:
    Default for delta_tau changed from 0.5 to 5 ms

Mar 23 12:27:56 correlomatrix_detector [Info]:
    Default for delta_tau changed from 0.5 to 5 ms

Mar 23 12:27:56 correlospinmatrix_detector [Info]:
    Default for delta_tau changed from 0.1 to 1 ms

Mar 23 12:27:56 SimulationManager::set_status [Info]:
    Temporal resolution changed from 0.1 to 1 ms.

Mar 23 12:27:56 NodeManager::prepare_nodes [Info]:
    Preparing 2 nodes for simulation.

Mar 23 12:27:56 SimulationManager::start_updating_ [Info]:
    Number of local nodes: 2
    Simulation time (ms): 25000
    Not using OpenMP
    Number of MPI processes: 1

Mar 23 12:27:56 SimulationManager::run [Info]:
    Simulation finished.
Actual variance: 0.0
Expected variance: 0.0
For h = 1.0, tau_noise = 1000.0, sigma_noise = 10.0

Mar 23 12:27:56 correlation_detector [Info]:
    Default for delta_tau changed from 0.5 to 5 ms

Mar 23 12:27:56 correlomatrix_detector [Info]:
    Default for delta_tau changed from 0.5 to 5 ms

Mar 23 12:27:56 correlospinmatrix_detector [Info]:
    Default for delta_tau changed from 0.1 to 1 ms

Mar 23 12:27:56 SimulationManager::set_status [Info]:
    Temporal resolution changed from 0.1 to 1 ms.

Mar 23 12:27:56 NodeManager::prepare_nodes [Info]:
    Preparing 2 nodes for simulation.

Mar 23 12:27:56 SimulationManager::start_updating_ [Info]:
    Number of local nodes: 2
    Simulation time (ms): 25000
    Not using OpenMP
    Number of MPI processes: 1

Mar 23 12:27:56 SimulationManager::run [Info]:
    Simulation finished.
Actual variance: 65.70587968208686
Expected variance: 100.0
For h = 1.0, tau_noise = 1000.0, sigma_noise = 100.0

Mar 23 12:27:56 correlation_detector [Info]:
    Default for delta_tau changed from 0.5 to 5 ms

Mar 23 12:27:56 correlomatrix_detector [Info]:
    Default for delta_tau changed from 0.5 to 5 ms

Mar 23 12:27:56 correlospinmatrix_detector [Info]:
    Default for delta_tau changed from 0.1 to 1 ms

Mar 23 12:27:56 SimulationManager::set_status [Info]:
    Temporal resolution changed from 0.1 to 1 ms.

Mar 23 12:27:56 NodeManager::prepare_nodes [Info]:
    Preparing 2 nodes for simulation.

Mar 23 12:27:56 SimulationManager::start_updating_ [Info]:
    Number of local nodes: 2
    Simulation time (ms): 25000
    Not using OpenMP
    Number of MPI processes: 1

Mar 23 12:27:56 SimulationManager::run [Info]:
    Simulation finished.
Actual variance: 6570.587968208678
Expected variance: 10000.0
For h = 1.0, tau_noise = 1000.0, sigma_noise = 1000.0

Mar 23 12:27:56 correlation_detector [Info]:
    Default for delta_tau changed from 0.5 to 5 ms

Mar 23 12:27:56 correlomatrix_detector [Info]:
    Default for delta_tau changed from 0.5 to 5 ms

Mar 23 12:27:56 correlospinmatrix_detector [Info]:
    Default for delta_tau changed from 0.1 to 1 ms

Mar 23 12:27:56 SimulationManager::set_status [Info]:
    Temporal resolution changed from 0.1 to 1 ms.

Mar 23 12:27:56 NodeManager::prepare_nodes [Info]:
    Preparing 2 nodes for simulation.

Mar 23 12:27:56 SimulationManager::start_updating_ [Info]:
    Number of local nodes: 2
    Simulation time (ms): 25000
    Not using OpenMP
    Number of MPI processes: 1

Mar 23 12:27:56 SimulationManager::run [Info]:
    Simulation finished.
Actual variance: 657058.7968208689
Expected variance: 1000000.0

Autocorrelogram

If we plot the autocorrelogram, we should find high correlations only for lags lower than \(\tau\).

Let’s look at a process with \(\tau=1000\):

[7]:
timevec, U = evaluate_ou_process(h=_h,
                                 t_sim=10000.,
                                 neuron_parms={"U" : 0.,
                                               "mean_noise": 0.,
                                               "tau_noise": 1000.,
                                               "sigma_noise": 1000.},
                                 title=r"Ornstein-Uhlenbeck process ($\tau$=1000, $\sigma$=1000)")

Mar 23 12:28:02 correlation_detector [Info]:
    Default for delta_tau changed from 0.5 to 5 ms

Mar 23 12:28:02 correlomatrix_detector [Info]:
    Default for delta_tau changed from 0.5 to 5 ms

Mar 23 12:28:02 correlospinmatrix_detector [Info]:
    Default for delta_tau changed from 0.1 to 1 ms

Mar 23 12:28:02 SimulationManager::set_status [Info]:
    Temporal resolution changed from 0.1 to 1 ms.

Mar 23 12:28:02 NodeManager::prepare_nodes [Info]:
    Preparing 2 nodes for simulation.

Mar 23 12:28:02 SimulationManager::start_updating_ [Info]:
    Number of local nodes: 2
    Simulation time (ms): 10000
    Not using OpenMP
    Number of MPI processes: 1

Mar 23 12:28:02 SimulationManager::run [Info]:
    Simulation finished.
../../_images/tutorials_ornstein_uhlenbeck_noise_nestml_ou_noise_tutorial_14_1.png
[8]:
def autocorr(x):
    return np.correlate(x, x, mode='full')

fig, ax = plt.subplots(figsize=(12., 5.))
U_autocorr = autocorr(U)
lags = np.arange(len(U_autocorr)) - len(U_autocorr) / 2
ax.plot(lags, U_autocorr)
ax.set_xlabel("Lag [ms]")
ax.set_ylabel("Autocorr[U]")
ax.grid()
../../_images/tutorials_ornstein_uhlenbeck_noise_nestml_ou_noise_tutorial_15_0.png

Integrating Ornstein-Uhlenbeck noise into a NESTML integrate-and-fire neuron

Now, the O-U noise process is integrated into an existing NESTML model, as an extra somatic current \(I_{noise}\). In this example, we use the integrate-and-fire neuron with current based synapses and an exponentially shaped synaptic kernel, and no refractoriness mechanism.

This neuron can be written in just a few lines of NESTML:

[50]:
iaf_psc_exp = '''
neuron iaf_psc_exp:

    state:
        V_m mV = E_L
        I_noise pA = mean_noise

    equations:
        kernel psc_kernel = exp(-t / tau_syn)
        V_m' = -(V_m - E_L) / tau_m + (convolve(psc_kernel, spikes) + I_e + I_noise) / C_m

    parameters:
        E_L mV = -65 mV      # resting potential
        I_e pA = 0 pA        # constant external input current
        tau_m ms = 25 ms     # membrane time constant
        tau_syn ms = 5 ms    # synaptic time constant
        C_m pF = 250 pF        # membrane capacitance
        V_theta mV = -30 mV  # threshold potential
        mean_noise pA = 0.057 pA       # mean of the noise current
        sigma_noise pA = 0.003 pA      # standard deviation of the noise current
        tau_noise ms = 10 ms              # time constant of the noise process

    internals:
        A_noise real = sigma_noise * ((1 - exp(-2 * resolution() / tau_noise)))**.5

    input:
        spikes pA <- spike

    output:
        spike

    update:
        I_noise = mean_noise \
                  + (I_noise - mean_noise) * exp(-resolution() / tau_noise) \
                  + A_noise * random_normal(0, 1)

        integrate_odes()

        if V_m > V_theta:
            V_m = E_L
            emit_spike()
'''
[51]:
with open("iaf_psc_exp.nestml", "w") as nestml_model_file:
    print(iaf_psc_exp, file=nestml_model_file)

generate_nest_target(input_path="iaf_psc_exp.nestml",
                     target_path="/tmp/nestml-target",
                     module_name="nestml_iaf_module",
                     suffix="_nestml", logging_level="ERROR",dev=True,
                     codegen_opts={"nest_path": NEST_SIMULATOR_INSTALL_LOCATION})


              -- N E S T --
  Copyright (C) 2004 The NEST Initiative

 Version: master@1510a5790
 Built: Mar 21 2023 14:28:29

 This program is provided AS IS and comes with
 NO WARRANTY. See the file LICENSE for details.

 Problems or suggestions?
   Visit https://www.nest-simulator.org

 Type 'nest.help()' to find out more about NEST.

CMake Warning:
  Ignoring empty string ("") provided on the command line.


-- The CXX compiler identification is AppleClang 14.0.0.14000029
-- Detecting CXX compiler ABI info
-- Detecting CXX compiler ABI info - done
-- Check for working CXX compiler: /Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin/c++ - skipped
-- Detecting CXX compile features
-- Detecting CXX compile features - done

-------------------------------------------------------
nestml_iaf_module Configuration Summary
-------------------------------------------------------

C++ compiler         : /Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin/c++
Build static libs    : OFF
C++ compiler flags   :
NEST compiler flags  : -Wextra -Wno-unknown-pragmas -D_GLIBCXX_ASSERTIONS -std=c++11 -Wall  -O2 -g
NEST include dirs    :  -I/Users/pooja/nest/include/nest -I/usr/local/include -I/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX13.1.sdk/usr/include -I/usr/local/Cellar/gsl/2.7/include -I/Users/pooja/miniconda3/envs/nestml_3.10/include -I/Users/pooja/miniconda3/envs/nestml_3.10/include
NEST libraries flags : -L/Users/pooja/nest/lib/nest -lnestutil -lsli -lnestkernel  /usr/local/lib/libltdl.dylib /Users/pooja/miniconda3/envs/nestml_3.10/lib/libreadline.dylib /Users/pooja/miniconda3/envs/nestml_3.10/lib/libncurses.dylib /usr/local/Cellar/gsl/2.7/lib/libgsl.dylib /usr/local/Cellar/gsl/2.7/lib/libgslcblas.dylib   /Users/pooja/miniconda3/envs/nestml_3.10/lib/libmpi.dylib

-------------------------------------------------------

You can now build and install 'nestml_iaf_module' using
  make
  make install

The library file libnestml_iaf_module.so will be installed to
  /Users/pooja/nest/lib/nest
The module can be loaded into NEST using
  (nestml_iaf_module) Install       (in SLI)
  nest.Install(nestml_iaf_module)   (in PyNEST)

CMake Warning (dev) in CMakeLists.txt:
  No cmake_minimum_required command is present.  A line of code such as

    cmake_minimum_required(VERSION 3.23)

  should be added at the top of the file.  The version specified may be lower
  if you wish to support older CMake versions for this project.  For more
  information run "cmake --help-policy CMP0000".
This warning is for project developers.  Use -Wno-dev to suppress it.

-- Configuring done
-- Generating done
-- Build files have been written to: /tmp/nestml-target
Consolidate compiler generated dependencies of target nestml_iaf_module_module
[ 66%] Building CXX object CMakeFiles/nestml_iaf_module_module.dir/iaf_psc_exp_nestml.o
[ 66%] Building CXX object CMakeFiles/nestml_iaf_module_module.dir/nestml_iaf_module.o
In file included from /tmp/nestml-target/iaf_psc_exp_nestml.cpp:33:
In file included from /tmp/nestml-target/nestml_iaf_module.cpp:26In file included from :
In file included from /Users/pooja/nest/include/nest/connection_manager_impl.h:/Users/pooja/nest/include/nest/kernel_manager.h26::
27In file included from :
/Users/pooja/nest/include/nest/connection_manager.hIn file included from :/Users/pooja/nest/include/nest/connection_manager.h36::
36:
In file included from /Users/pooja/nest/include/nest/connector_base.hIn file included from :35:
In file included from /Users/pooja/nest/include/nest/sort.h:37:
In file included from /Users/pooja/nest/include/nest/connector_base.h/usr/local/include/boost/sort/spreadsort/spreadsort.hpp:27:
In file included from /usr/local/include/boost/sort/spreadsort/integer_sort.hpp:26:
:/usr/local/include/boost/sort/spreadsort/detail/integer_sort.hpp35::
438In file included from :/Users/pooja/nest/include/nest/sort.h29::37 :
warning: In file included from unused parameter 'shift' [-Wunused-parameter]/usr/local/include/boost/sort/spreadsort/spreadsort.hpp
:27:
In file included from /usr/local/include/boost/sort/spreadsort/integer_sort.hpp:26:
/usr/local/include/boost/sort/spreadsort/detail/integer_sort.hpp:438:29: warning: unused parameter 'shift' [-Wunused-parameter]
                Right_shift shift, Compare comp)
                            ^
                Right_shift shift, Compare comp)
                            ^
/usr/local/include/boost/sort/spreadsort/detail/integer_sort.hpp:482:29: warning: unused parameter 'shift' [-Wunused-parameter]
/usr/local/include/boost/sort/spreadsort/detail/integer_sort.hpp:482                Right_shift shift):
29                            ^:
 warning: unused parameter 'shift' [-Wunused-parameter]
                Right_shift shift)
                            ^
In file included from /tmp/nestml-target/iaf_psc_exp_nestml.cpp:33:
In file included from /Users/pooja/nest/include/nest/kernel_manager.h:27:
In file included from /Users/pooja/nest/include/nest/connection_manager.h:36:
In file included from /tmp/nestml-target/nestml_iaf_module.cpp:26:
In file included from /Users/pooja/nest/include/nest/connection_manager_impl.h:26:
In file included from /Users/pooja/nest/include/nest/connection_manager.h:In file included from /Users/pooja/nest/include/nest/connector_base.h:35:
36:
In file included from In file included from /Users/pooja/nest/include/nest/connector_base.h/Users/pooja/nest/include/nest/sort.h::3537:
:
In file included from In file included from /Users/pooja/nest/include/nest/sort.h/usr/local/include/boost/sort/spreadsort/spreadsort.hpp::3728:
:
In file included from In file included from /usr/local/include/boost/sort/spreadsort/spreadsort.hpp/usr/local/include/boost/sort/spreadsort/float_sort.hpp::2825:
:
In file included from /usr/local/include/boost/sort/spreadsort/float_sort.hpp:25/usr/local/include/boost/sort/spreadsort/detail/float_sort.hpp:778:28::
/usr/local/include/boost/sort/spreadsort/detail/float_sort.hpp:778 warning: unused parameter 'rshift' [-Wunused-parameter]
:28:                Right_shift rshift)
                           ^
warning: unused parameter 'rshift' [-Wunused-parameter]
               Right_shift rshift)
                           ^
/usr/local/include/boost/sort/spreadsort/detail/float_sort.hpp:820:28: warning: unused parameter 'rshift' [-Wunused-parameter]
               Right_shift rshift, Compare comp)
                           ^
/usr/local/include/boost/sort/spreadsort/detail/float_sort.hpp:820:28: warning: unused parameter 'rshift' [-Wunused-parameter]
               Right_shift rshift, Compare comp)
                           ^
In file included from /tmp/nestml-target/nestml_iaf_module.cpp:26:
In file included from /Users/pooja/nest/include/nest/connection_manager_impl.h:26:
In file included from /Users/pooja/nest/include/nest/connection_manager.h:36:
In file included from /Users/pooja/nest/include/nest/connector_base.h:35:
In file included from /Users/pooja/nest/include/nest/sort.h:37:
In file included from /usr/local/include/boost/sort/spreadsort/spreadsort.hpp:In file included from /tmp/nestml-target/iaf_psc_exp_nestml.cpp:33:
29In file included from /Users/pooja/nest/include/nest/kernel_manager.h:
:In file included from 27/usr/local/include/boost/sort/spreadsort/string_sort.hpp:
:In file included from 23/Users/pooja/nest/include/nest/connection_manager.h:
:36/usr/local/include/boost/sort/spreadsort/detail/string_sort.hpp:752:26::
 warning: unused parameter 'get_character' [-Wunused-parameter]
In file included from /Users/pooja/nest/include/nest/connector_base.h:35:
In file included from /Users/pooja/nest/include/nest/sort.h                Get_char get_character, Get_length length, Unsigned_char_type):
37                         ^:

In file included from /usr/local/include/boost/sort/spreadsort/spreadsort.hpp:29:
In file included from /usr/local/include/boost/sort/spreadsort/string_sort.hpp/usr/local/include/boost/sort/spreadsort/detail/string_sort.hpp:752:52:23:
:/usr/local/include/boost/sort/spreadsort/detail/string_sort.hpp warning: unused parameter 'length' [-Wunused-parameter]
:                Get_char get_character, Get_length length, Unsigned_char_type)
                                                   ^
752:26: warning: unused parameter 'get_character' [-Wunused-parameter]
                Get_char get_character, Get_length length, Unsigned_char_type)
                         ^
/usr/local/include/boost/sort/spreadsort/detail/string_sort.hpp:752:52: warning: unused parameter 'length' [-Wunused-parameter]
                Get_char get_character, Get_length length, Unsigned_char_type)
                                                   ^
/usr/local/include/boost/sort/spreadsort/detail/string_sort.hpp:780:18: warning: unused parameter 'get_character' [-Wunused-parameter]
        Get_char get_character, Get_length length, Compare comp, Unsigned_char_type)
                 ^
/usr/local/include/boost/sort/spreadsort/detail/string_sort.hpp:780:44: warning: unused parameter 'length' [-Wunused-parameter]
        Get_char get_character, Get_length length, Compare comp, Unsigned_char_type)
                                           ^
/usr/local/include/boost/sort/spreadsort/detail/string_sort.hpp:780:18: warning: unused parameter 'get_character' [-Wunused-parameter]
        Get_char get_character, Get_length length, Compare comp, Unsigned_char_type)
                 ^
/usr/local/include/boost/sort/spreadsort/detail/string_sort.hpp:780:44: warning: unused parameter 'length' [-Wunused-parameter]
        Get_char get_character, Get_length length, Compare comp, Unsigned_char_type)
                                           ^
/usr/local/include/boost/sort/spreadsort/detail/string_sort.hpp:807:18: warning: unused parameter 'get_character' [-Wunused-parameter]
        Get_char get_character, Get_length length, Compare comp, Unsigned_char_type)
                 ^
/usr/local/include/boost/sort/spreadsort/detail/string_sort.hpp:807:44: warning: unused parameter 'length' [-Wunused-parameter]
/usr/local/include/boost/sort/spreadsort/detail/string_sort.hpp:807:        Get_char get_character, Get_length length, Compare comp, Unsigned_char_type)18
:                                           ^
warning: unused parameter 'get_character' [-Wunused-parameter]
        Get_char get_character, Get_length length, Compare comp, Unsigned_char_type)
                 ^
/usr/local/include/boost/sort/spreadsort/detail/string_sort.hpp:807:44: warning: unused parameter 'length' [-Wunused-parameter]
        Get_char get_character, Get_length length, Compare comp, Unsigned_char_type)
                                           ^
In file included from /tmp/nestml-target/iaf_psc_exp_nestml.cpp:43:
/tmp/nestml-target/iaf_psc_exp_nestml.h:233:17: warning: 'iaf_psc_exp_nestml::get_C_m' hides overloaded virtual function [-Woverloaded-virtual]
  inline double get_C_m() const
                ^
/Users/pooja/nest/include/nest/node.h:692:18: note: hidden overloaded virtual function 'nest::Node::get_C_m' declared here: different number of parameters (1 vs 0)
  virtual double get_C_m( int comp );
                 ^
In file included from /tmp/nestml-target/nestml_iaf_module.cpp:47:
/tmp/nestml-target/iaf_psc_exp_nestml.h:233:17: warning: 'iaf_psc_exp_nestml::get_C_m' hides overloaded virtual function [-Woverloaded-virtual]
  inline double get_C_m() const
                ^
/Users/pooja/nest/include/nest/node.h:692:18: note: hidden overloaded virtual function 'nest::Node::get_C_m' declared here: different number of parameters (1 vs 0)
  virtual double get_C_m( int comp );
                 ^
/tmp/nestml-target/nestml_iaf_module.cpp:100:42: warning: unused parameter 'i' [-Wunused-parameter]
nestml_iaf_module::init( SLIInterpreter* i )
                                         ^
/tmp/nestml-target/iaf_psc_exp_nestml.cpp:105:16: warning: unused variable '__resolution' [-Wunused-variable]
  const double __resolution = nest::Time::get_resolution().get_ms();  // do not remove, this is necessary for the resolution() function
               ^
/tmp/nestml-target/iaf_psc_exp_nestml.cpp:258:24: warning: comparison of integers of different signs: 'long' and 'const size_t' (aka 'const unsigned long') [-Wsign-compare]
    for (long i = 0; i < NUM_SPIKE_RECEPTORS; ++i)
                     ~ ^ ~~~~~~~~~~~~~~~~~~~
12 warnings generated.
13 warnings generated.
[100%] Linking CXX shared module nestml_iaf_module.so
[100%] Built target nestml_iaf_module_module
Consolidate compiler generated dependencies of target nestml_iaf_module_module
[100%] Built target nestml_iaf_module_module
Install the project...
-- Install configuration: ""
-- Installing: /Users/pooja/nest/lib/nest/nestml_iaf_module.so

Now, the NESTML model is ready to be used in a simulation.

[64]:
def evaluate_neuron(neuron_name, neuron_parms=None, I_e=0.,
                    mu=0., sigma=0., t_sim=300., plot=True):
    """
    Run a simulation in NEST for the specified neuron. Inject a stepwise
    current and plot the membrane potential dynamics and spikes generated.
    """
    dt = .1   # [ms]

    nest.ResetKernel()
    try:
        nest.Install("nestml_iaf_module")
    except :
        pass
    neuron = nest.Create(neuron_name)
    if neuron_parms:
        for k, v in neuron_parms.items():
            nest.SetStatus(neuron, k, v)
    nest.SetStatus(neuron, "I_e", I_e)
    nest.SetStatus(neuron, "mean_noise", mu)
    nest.SetStatus(neuron, "sigma_noise", sigma)

    multimeter = nest.Create("multimeter")
    multimeter.set({"record_from": ["V_m"],
                    "interval": dt})
    sr = nest.Create("spike_recorder")
    nest.Connect(multimeter, neuron)
    nest.Connect(neuron, sr)

    nest.Simulate(t_sim)

    dmm = nest.GetStatus(multimeter)[0]
    Voltages = dmm["events"]["V_m"]
    tv = dmm["events"]["times"]
    dSD = nest.GetStatus(sr, keys='events')[0]
    spikes = dSD['senders']
    ts = dSD["times"]

    _idx = [np.argmin((tv - spike_time)**2) - 1 for spike_time in ts]
    V_m_at_spike_times = Voltages[_idx]

    if plot:
        fig, ax = plt.subplots()
        ax.plot(tv, Voltages)
        ax.scatter(ts, V_m_at_spike_times)
        ax.set_xlabel("Time [ms]")
        ax.set_ylabel("V_m [mV]")
        ax.grid()

    return ts

If we stimulate the neuron with a constant current of 300 pA alone, it does not spike:

[74]:
spike_times = evaluate_neuron("iaf_psc_exp_nestml", mu=300, sigma=0.)

Mar 24 12:22:41 NodeManager::prepare_nodes [Info]:
    Preparing 3 nodes for simulation.

Mar 24 12:22:41 SimulationManager::start_updating_ [Info]:
    Number of local nodes: 3
    Simulation time (ms): 300
    Not using OpenMP
    Number of MPI processes: 1

Mar 24 12:22:41 SimulationManager::run [Info]:
    Simulation finished.
../../_images/tutorials_ornstein_uhlenbeck_noise_nestml_ou_noise_tutorial_22_1.png

However, for the same \(\mu\)=300 pA, but setting \(\sigma\)=200 pA/√ms, the effect of the noise can be clearly seen in the membrane potential, and the occasional spiking of the neuron:

[68]:
spike_times = evaluate_neuron("iaf_psc_exp_nestml", mu=300, sigma=200)
assert spike_times.size > 0

Mar 23 15:14:20 NodeManager::prepare_nodes [Info]:
    Preparing 3 nodes for simulation.

Mar 23 15:14:20 SimulationManager::start_updating_ [Info]:
    Number of local nodes: 3
    Simulation time (ms): 300
    Not using OpenMP
    Number of MPI processes: 1

Mar 23 15:14:20 SimulationManager::run [Info]:
    Simulation finished.
../../_images/tutorials_ornstein_uhlenbeck_noise_nestml_ou_noise_tutorial_24_1.png

To further quantify the effects of the noise, we can plot a distribution of interspike intervals:

[73]:
spike_times = evaluate_neuron("iaf_psc_exp_nestml",
                              mu=300,
                              sigma=200,
                              t_sim=25000.,
                              plot=False)

ISI = np.diff(spike_times)
assert ISI.size > 0

ISI_mean = np.mean(ISI)
ISI_std = np.std(ISI)

print(str(len(spike_times)) + " spikes recorded")
print("Mean ISI: " + str(ISI_mean))
print("ISI std. dev.: " + str(ISI_std))

count, bin_edges = np.histogram(ISI)

fig, ax = plt.subplots()
ax.bar(bin_edges[:-1], count, width=.8 * (bin_edges[1] - bin_edges[0]))
ylim = ax.get_ylim()
ax.plot([ISI_mean, ISI_mean], ax.get_ylim(), c="#11CC22", linewidth=3)
ax.plot([ISI_mean - ISI_std, ISI_mean - ISI_std], ylim, c="#CCCC11", linewidth=3)
ax.plot([ISI_mean + ISI_std, ISI_mean + ISI_std], ylim, c="#CCCC11", linewidth=3)
ax.set_ylim(ylim)
ax.set_xlabel("ISI [ms]")
ax.set_ylabel("Count")
ax.grid()

Mar 23 15:15:19 NodeManager::prepare_nodes [Info]:
    Preparing 3 nodes for simulation.

Mar 23 15:15:19 SimulationManager::start_updating_ [Info]:
    Number of local nodes: 3
    Simulation time (ms): 25000
    Not using OpenMP
    Number of MPI processes: 1

Mar 23 15:15:19 SimulationManager::run [Info]:
    Simulation finished.
265 spikes recorded
Mean ISI: 94.38219696969698
ISI std. dev.: 88.85184438847593
../../_images/tutorials_ornstein_uhlenbeck_noise_nestml_ou_noise_tutorial_26_1.png

Further directions

  • Calculate and plot the time-varying expected variance of the OU process when its initial value is unequal to the process mean ([1], eq. 2.26; see Fig. 2 for a visual example)

  • Make an extension of the neuron model, that stimulates the cell with an inhibitory as well as excitatory noise current.

  • Instead of a current-based noise, make a neuron model that contains a conductance-based noise.

References

[1] D.T. Gillespie, “The mathematics of Brownian motion and Johnson noise”, Am. J. Phys. 64, 225 (1996); doi: 10.1119/1.18210

Acknowledgements

Thanks to Tobias Schulte to Brinke, Barna Zajzon, Renato Duarte, Claudia Bachmann and all participants of the CNS2020 tutorial on NEST Desktop & NESTML!

This software was developed in part or in whole in the Human Brain Project, funded from the European Union’s Horizon 2020 Framework Programme for Research and Innovation under Specific Grant Agreements No. 720270 and No. 785907 (Human Brain Project SGA1 and SGA2).

License

This notebook (and associated files) is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 2 of the License, or (at your option) any later version.

This notebook (and associated files) is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.