How to implement stochastic differential equation solvers?

Particularly useful chunks of hoc and/or NMODL code. May be pedestrian or stunningly brilliant, may make you gasp or bring tears to your eyes, but always makes you think "I wish I had written that; I'm sure going to steal it."
Post Reply
nicola
Posts: 3
Joined: Tue May 17, 2016 12:18 pm

How to implement stochastic differential equation solvers?

Post by nicola »

I am using some data and multi-compartment models implemented in NEURON by the Blue Brain Project (BBP, http://www.cell.com/cell/fulltext/S0092-8674(15)01191-5) in order to test some parameter estimation methods. Since I am not very familiar with the HOC programming language, in my coding I mostly use the python interface (NEURON+Python).

For technical reasons regarding the methods I am using, I would like to turn the ODE models I have into their stochastic counterparts. Such a stochastic differential equation (SDE) model would essentially result from adding some Brownian-noise perturbation in the membrane potential and activation variables. However, given the numerical solvers available in NEURON (implicit Euler, Crank-Nicolson, and Adams-Bashforth, fixed-step or adaptive, if I correctly got the whole picture), I can't find any intuitive "hack" to implement any SDE solver. For instance, I could easily implement a Euler-Maruyama scheme if a simple explicit Euler method for ODEs was available.

Is there any way I could modify the standard NEURON's solvers in order to adapt a given numerical method to my needs? Thanks.
hines
Site Admin
Posts: 1682
Joined: Wed May 18, 2005 3:32 pm

Re: How to implement stochastic differential equation solver

Post by hines »

I'm afraid the only fixed step integrator for the current balance equations is the backward euler method (or crank-nicolson when secondorder=2) and, when playing a random normal variable into a current, I
am don't know whether that is substantially numerically equivalent to the euler-maruyama method. I played for a while with what I thought was a reasonable implementation

Code: Select all

from neuron import h, gui

s = h.Section()
s.insert("pas")
s.diam = 10
s.L = 100/(h.PI * s.diam) # area 100 um2   (1nA same as 1mA/cm2)
s.g_pas = .001
s.e_pas = 10  
h.v_init = -65

stim = h.IClamp(s(.5))
stim.delay = 0
stim.dur = 1e9
stim.amp = 0 # will play random variable into this

r = h.Random()
r.Random123(0,0,0)
sigma = 0.001
r.play(stim._ref_amp)           

def set_r():
 var = sigma/h.sqrt(h.dt))        
 print ("set_r dt=%g var=%g" % (h.dt, var))
 r.normal(0, var)

fih = h.FInitializeHandler(1, set_r)


#h.load_file("temp.ses") # pop up RunControl and voltage graph

But I was disappointed that the statistics for 40 points/ms was very different from the statistics for 400 and 4000 points per ms.
I seemed to get more similar voltage excursions using
var = sigma/h.dt
nicola
Posts: 3
Joined: Tue May 17, 2016 12:18 pm

Re: How to implement stochastic differential equation solver

Post by nicola »

Dear Michael Hines,
thank you very much for your reply. I think that analyzing the code you provided I found some kind of solution to this problem.

First, I realized that it does exist a implicit version of the standard Euler-Maruyama explicit scheme for SDEs (see for instance section 12.2 of Kloeden PE, Platen E (1999) Numerical Solution of Stochastic Differential Equations, Springer), so that the code you proposed should indeed be an implementation of a linear SDE with stochastic Brownian-motion-like input current.

However, from my understanding of SDEs, I reckon that the variance of the Brownian-motion term should actually be proportional to the numerical time step (i.e. var=sigma*h.dt). I am not sure what you refer to when mentioning the statistics for 40-400-4000 steps/ms, but in case you are talking about comparing the mean and standard deviation of the different SDE solutions, such adjustment should do the trick.
In fact, running the slightly modified code provided below, I get much more consistent results in that respect.

Code: Select all

from neuron import h, gui

s = h.Section()
s.insert("pas")
s.diam = 10
s.L = 100/(h.PI * s.diam) # area 100 um2   (1nA same as 1mA/cm2)
s.g_pas = .001
s.e_pas = 10  
h.v_init = -65

stim = h.IClamp(s(.5))
stim.delay = 0
stim.dur = 1e9
stim.amp = 0 # will play random variable into this

h.tstop = 100
h.steps_per_ms=40
h.dt = 1/h.steps_per_ms

r = h.Random()
r.Random123(0,0,0)
sigma = 0.1
r.play(stim._ref_amp)     #At the beginning of every call to :func:`fadvance` and :func:`finitialize` var is set to a new value equivalent to var=r.repick()

def set_r():
    var = sigma*h.dt #var = sigma/h.dt  #var = sigma/h.sqrt(h.dt)
    print ("set_r dt=%g var=%g" % (h.dt, var))
    r.normal(0, var)
fih = h.FInitializeHandler(1, set_r)

v_vec = h.Vector()   # Membrane potential vector at soma
v_vec.record( s(0.5)._ref_v )

print "tstop =", h.tstop, "dt =", h.dt
h.run()
print "run complete up to t =", h.t

print "solution mean =", v_vec.mean(), "and std ", v_vec.stdev()

Now, I would also like to be able to add a stochastic term in the equations for the activation variables in a Hodgkin-Huxeley-type neuron model. Unfortunately, I guess there is not such a thing as an IClamp stimulator acting on hidden variables (e.g. m,h, and n in HH model) in NEURON, correct?
hines
Site Admin
Posts: 1682
Joined: Wed May 18, 2005 3:32 pm

Re: How to implement stochastic differential equation solver

Post by hines »

All I can say is that with the fixed step method, every step can be thought of as a new initial value problem. So you are allowed to change any/every state
of the mode arbitrarily on return from or before entry to fadvance()
nicola
Posts: 3
Joined: Tue May 17, 2016 12:18 pm

Re: How to implement stochastic differential equation solver

Post by nicola »

Right, this is what I was thinking of doing actually.

Thank you very much for your very helpful support.
Post Reply