SDE

[boeschf]

Stochastic Differential Equations

Introduces sources of white noise to arbor (and nmodl) and addresses issues #1643 and #1655. Both point and density mechanisms may now use white noise as part of the state updates, turning the ODEs effectively into SDEs. The noise sources are provided per connection end point (point mech.) or control volume (density mech.) and are uncorrelated (spatially, temporally and across different mechanism instantiations)

Background

SDE's can be described in their differential form (informal engineering notation) as

$$ X^\prime = f\left(X(t),t\right) + \sum_i g_i\left( X(t), t\right) W_i(t)$$

with (standard) white noises $W_i$. More accurately, this can be expressed as integral equation,

$$ X(t+s) = X(t) + \int_{t}^{t+s} f\left( X(\tau), \tau\right) d\tau + \sum_i \int_t^{t+s} g_i\left( X(\tau), \tau \right)dB_i(\tau)$$

with (standard) Brownian motions $B_i$ (also referred to as Wiener processes). The white noise $W_i$ are stationary, Gaussian processes with $μ=0$ and $σ^2=1$. The stochastic integrals are interpreted in the Itô sense.

Systems of stochastic differential equations can be solved using a few different methods. We chose the Euler-Maruyama method here (based on Taylor-Itô expansion) for its simplicity and generality. The main drawback of this approach is accuracy (first order only). While there exist other methods (analytical and numerical), these are not general and not straight-forward to extend to systems of SDEs. Future work could be done to investigate if applying more accurate methods in corner-cases makes sense.

nmodl

The white noise sources can be defined in the model files using a WHITE_NOISE block:

   WHITE_NOISE {
       a b 
       c
   }

Arbitrary white noise variables can be declared afterwards (a, b, c in the example above).

If the state is updated by involving at least one of the declared white noise variables the system is considered to be stochastic:

   DERIVATIVE state {
       s' = f + g*a
   }

The solver method must then accordingly set to stochastic:

   BREAKPOINT {
       SOLVE state METHOD stochastic
   }

Simulation

The simulation class takes an optional seed argument at construction which sets the global seed for the random number generator. The default value is 0. In order to make construction of the simulation object more convenient, a builder class is added.

modcc

modcc gains a new solver method for Euler-Maruyama. The implementation required extension of the expressions and tokens to handle white noise and seed values. The code generation in the cpu and gpu printers is pretty straight-forward and akin to looking up state variables.

ABI

ABI changes by extending the struct arb_mechanism_ppack to accomodate for the random numbers. This requires a bump of the ABI version.

Arbor

The main changes in Arbor concern the backend implementations of the shared_state. Generally, we want to generate a batch of new random numbers before the advance_state procedure is invoked. As those random numbers need to be independent, we use the counter-based pseudo random number generators from Random123. We can prime the counters and keys with the following identifiers to guarantee independent random values:

• global seed value (per simulation)
• global cell id
• per-cell connection end point index (point mech.) or per-cell control volume index (density mech)
• mechanism id
• per-mechanism variable identifier
• per-mechanism simulation time (in the form of an integral counter)

The random123 library usually provides more than one independent random values per input (4 in our case as of now). In order to make use of these values, the backend caches the number and updates the numbers only every 4 steps. Cpu, cpu-vectorized and gpu variants should produce the same random values. In order to guarantee correct random numbers, synapse collapsing is disabled.

[thorstenhater]

So, the GPU pRNG is within $4\epsilon$ of the CPU version?

[boeschf]

yes (absolute tolerance)

[thorstenhater]

This is my only painpoint left:
Could we make it so that shared state doesn't need to know about this flag? Currently it's used in three
four different places, which might cause issues down the road. How about delegating these things to the
pRNG impl

• writing the values (here)
• setting gid/idx (ll 358)
• allocating (ll 396)/ filling (ll 477) the indices

via the use of an appropriate interface. In said interface we can make the decision once.

[thorstenhater]

🎉 👍 🎉
A massive effort and a very cool addition to Arbor! Thanks!