AEP: Stochastic differential equations in mechanisms

[noraabiakar]

This issue documents the discussion and ideas for supporting stochastic differential equations for independent point sources of noise. Spatially distributed and correlated noise sources are not discussed here.

Representation in nmodl:

Expose a new symbol ζ that represents Gaussian white noise.
The Stochastic differential equation for one-dimensional (time) processes, can typically be described in the form:

d X(t) = f(X,t) dt + g(X,t) d W(t)

where W is the Wiener process, but with some abuse of notation, this can be written like an ODE:

X’ = f(X,t) + g(X,t)ζ(t)

where ζ represents (n-dimensional) Gaussian white noise. The NMODL representation could then piggyback on the ODE representation with a new special symbol representing ζ.

Solving the SDE:

Possibilities:

1. The Euler-Maruyama method, described here.
2. A second-order solution such as the Runge-Kutta method for SDEs, described here.
3. Exact solution

Seeding the random number generators:

The Random123 library can be used to supply normally distributed samples used for solving the SDEs. The seeds of the random number generators can be supplied to the mechanisms as parameters, however, it is better to think of them as properties of the simulation object.

The user can supply a seed to the simulation object (per cell_group?), the simulator, together with the cell_group can then generate seeds for the needed processes based on the user-supplied seed in combination with (gid, lid) identifiers.

[noraabiakar]

@max9901 Is this a good description of what you need? You mentioned some time ago that what you really want out of this feature is the seed generation. Can you describe what you expect there?

[schmitts]

@noraabiakar Thanks. Looks good. I've also asked @jlubo for comments.

[jlubo]

I just wrote in issue #1643 about a current workaround using the Euler method and the minstd_rand0 random generator of the C++ standard library. @noraabiakar What are the arguments in favor of Random123? Couldn't the minstd_rand0 generator also be used for arbor?

Regarding the solving method, maybe for selected random processes (such as Ornstein-Uhlenbeck) exact solutions could be provided...

[bcumming]

The requirements for random generators for stochastic processes, where there are multiple instances of processes on multiple cells, are that the generator:

1. can be seeded properly
2. can be used to generate independent non-correlated streams
3. produces high quality random streams that pass statistical tests.
4. can be used to generate reproducible results when then same model running on all architectures, CPU and GPU.

minstd_rand0 is not high quality (fails basic statistical tests RNG quality), generates highly correlated random streams, and is more challenging too seed. It is very simple, so a GPU implementation could be hand rolled.

Random123 is a counter-based random number generator, that makes it possible to generate reproducible high quality random number streams in parallel, with efficient CPU and GPU implementations. By virtue of using being counter based, seeding the value according cell gid, process id, etc is easier.

[jlubo]

Thank you for your answer, Ben. For my purposes, minstd_rand0 with careful seeding worked fine, so far. But I see that for general use in arbor Random123 will be much better.

[noraabiakar]

More details about the proposed implementation:

• Allow SDEs of the form d X(t) = f(X,t) dt + g(X,t) d W(t) where W(t) is the Wiener process.

  dX(t) = f(X,t)dt + g(X,t)dW(t)
  
  dX(t)/dt = f(X,t) + g(X,t)dW(t)/dt
  
  // NMODL equiavlent 
  // X’ = f(X,t) + g(X,t)ζ(t)
  // where ζ(t) = dW(t)/dt
  
  // Explicit Euler(-Maruyama)
  X_n+1 = X_n + Δt*f(X_n,t_n) + ΔW_n*g(X_n,t_n)
  
  // ΔW_n are independent and identically distributed
  // normal random variables with μ=0 and σ²=Δt
  

  Using the Euler-Maruyma method, the solution is given by X_n+1 = X_n + Δt*f(X_n,t_n) + ΔW_n*g(X_n,t_n). This can be implemented using a new solver (see here). All of the currently implemented solvers use implicit Euler or an "exact" solution using symbolic differentiation.
  The solver would generate the solution of an SDE in terms of some variable representing ΔW. There would be a different variable ΔWk for each SDE in the system. These variables would then be replaced at the code generation step by values sampled from a normal distributions with μ=0 and σ²=Δt
• On the NMODL side, represent SDEs of the form X’ = f(X,t) + g(X,t)ζ(t), where ζ is a new symbol that represents n-dimensional Gaussian white noise (ζ=dW(t)/dt).
• Use the Random123 library and Box-Muller we to generate normally distributed random numbers. Given 2 uniformally distributed random numbers, Box-muller will generate 2 normally distributed random numbers. Random123 has a GPU implementation which has not been tested in Arbor, which will need to be tested for use in the GPU implementation of SDE mechanisms.
• To seed the random number generators, use a combination (or hash) of a global seed and a value unique to the instance and index of the mechanism, e.g. global_seed ⊕ mechanism_id ⊕ index.
  • global_seed: gid hashed with another seed provided by the user through the simulation object.
  • mechanism_id: unique integer id per mechanism name.
  • index: index per mechanism instance for each mechanism. (not per lid, an lid can contain multiple instances of a gap-junction mechanisms.)
    The seed can also be hashed with the index of the SDE in a system of stochastic ODEs to provide different noise for each SDE.
• The random number generators need to provide different samples at each iteration of the integration loop given the same seed (seed = global_seed ⊕ mechanism_id ⊕ index ⊕ SDE_index). This can be done by incorporating the exact time of the simulation into the seed, or by keeping in the shared_state an external counter per mechanism instance which is incremented after each iteration of the integration loop.

[thorstenhater]

Except #1987 this is done.