Integration of stochastic differential equation

[Michel Besserve]

[Michel Besserve]

Dear all,
I am new to PyDSTool and try to simulate a differential equation which contains a noise term. I am looking for the simplest way to do this with pyDStool. While there is a tutorial on simulating a noisy channel, I could not follow what was stochastic in the code. In principle, I would like to represent a Wiener process in the differential equation, but this cannot be written as a simple function of the parameters and state variables: not only a sample should be drawn from a normal distribution at each step, but the noise amplitude depends on the integration step. Is it possible to integrate such a process at all?
As I do not need to be very rigorous about the way the noise is integrated, I could also just simulate a stochastic process independently on a grid with another library, and use this signal as an input to the system enforcing interpolation between grid points. In case this second option is the easiest, could you point me to some related example?
Thanks in advance,
Michel

[Rob Clewley]

Hi Michel, yes the second approach is the only one available. You'll
have to do some kind of analytic process around the results of
multiple runs with different instances of the stochastic process. Any
of the examples that use the 'inputs' parameter to the (usually named)
DSargs struct are basically the same idea. For instance,
https://github.com/robclewley/pydstool/blob/master/examples/test_hybrid_extinputs.py

Let us know if you have further questions.

On Tue, Apr 11, 2017 at 3:59 PM, Michel Besserve mbesserve@users.sf.net wrote:

Dear all,
I am new to PyDSTool and try to simulate a differential equation which
contains a noise term. I am looking for the simplest way to do this with
pyDStool. While there is a tutorial on simulating a noisy channel, I could
not follow what was stochastic in the code. In principle, I would like to
represent a Wiener process in the differential equation, but this cannot be
written as a simple function of the parameters and state variables: not only
a sample should be drawn from a normal fistribution at each step, but the
noise amplitude depends on the integration step. Is this possible?
As I do not need to be very rigorous about the way the noise is integrated,
I could also just simulate a stochastic process independently on a grid with
another library, and use this signal as an input to the system enforcing
interpolation between grid points. In case this second option is the
easiest, could you point me to some related example?
Thanks in advance,
Michel

---

Integration of stochastic differential equation

---

Sent from sourceforge.net because you indicated interest in
https://sourceforge.net/p/pydstool/discussion/472291/

To unsubscribe from further messages, please visit
https://sourceforge.net/auth/subscriptions/

[Michel Besserve]

Hi Rob,
Thanks for the quick answer. I followed your advice and this works using the Vode integrator, but it crashes when in try to use Dopri ('No Trajectory Created' message...). Vode is really slow (I think also due to the added external interpolated input), so I would like to try whether Dopri could do better. Otherwise would you have some tips to optimize the introduction of the external input in the model in order to make the integration faster? I provide the full code below.
Thanks,
Michel

# -*- coding: utf-8 -*-
"""
Spyder Editor

This is a temporary script file.
"""

import PyDSTool as dst
import numpy as np
from matplotlib import pyplot as plt

# define interpolated noise to emulate stochasticity
tdomain = [0,12000]
timeData = dst.linspace(tdomain[0], tdomain[1], 100000)
noise_std = .3
sindata = noise_std*np.random.randn(len(timeData)) 
xData = {'in': sindata}

my_input = dst.InterpolateTable({'tdata': timeData,
                              'ics': xData,
                              'name': 'interp'
                              }).compute('interp')

# we must give a name
DSargs = dst.args(name='IsolatedThalamusModel')
# parameters
DSargs.pars = {'Qmax': 400e-3,# max firing rate
               'theta': -58.5,# firing threshold
               'sigmaT': 6.0, # neural gain of TC neurons
               'sigmaR': 6.0, # neural gain of RT neurons
               'tau': 20.0, # membrane time constant
               'EL': -70.0, # leak current reversal potential
               'Eampa': 0.0, # reversal potential of excitatroy synapse
               'Egaba': -70.0, # reversal potential of inhibitory synapse
               'Ek': -100, # potassium
               'Eh': -40, # for the rectifier current
               'Eca': 120.0,
               'gamma_e': 70e-3, 
               'gamma_r': 100e-3,
               'Ntr': 5.0,
               'Nrt': 3.0,
               'Nrr': 25.0,
               'gLK': 0.018, # leak conductance, setting for isolated thalamus
               'gTt': 3.0, # T channel conductance
               'gTr': 2.3, #
               'ginc': 2.0,
               'gh': .06, # conductance of the anomalous rectifier, setting for isolated thalamus
               'Cm': 1.0, # coupling capacitance
               'k1': 2.5e7,
               'k2': 4e-4,
               'k3': 1e-1,
               'k4': 1e-3,
               'alpha_Ca': 51.8e-6,
               'tau_Ca': 10.0,
               'Ca0': 2.4e-4
               }
DSargs.inputs = {'in': my_input.variables['in']}
# auxiliary helper function(s) -- function name: ([func signature], definition)
DSargs.fnspecs  = {'Qt': (['v'], 'Qmax / (1 + exp( -(v-theta)/sigmaT ))'),
                  'Qr': (['v'], 'Qmax / (1 + exp( -(v-theta)/sigmaR ))'),
                  'm_inf_t': (['v'], '1 / (1 + exp( -(v+59)/6.2 ))'),
                  'm_inf_r': (['v'], '1 / (1 + exp( -(v+52)/7.4 ))'),
                  'hinfT_t': (['v'], '1 / (1 + exp( (v+81)/4 ))'),
                  'hinfT_r': (['v'], '1 / (1 + exp( (v+80)/5 ))'),
                  'tau_ht': (['v'],'(30.8+(211.4+exp((v+115.2)/5))/(1+exp((v+86)/3.2)))/exp(1.2*log(3))'),
                  'tau_hr': (['v'],'(85+1/(exp((v+48)/4)+exp(-(v+407)/50)))/exp(1.2*log(3))'),
                  'mhinf': (['v'],'1/(1+exp((v+75)/5.5))'),
                  'tau_hm': (['v'],'(20+1000/(exp((v+71.5)/14.2)+exp(-(v+89)/11.6)))'),
                  'Ph':(['C'],'k1*C*C*C*C/(k1*C*C*C*C+k2)')
                  }

# rhs of the differential equation, including dummy variable w
DSargs.varspecs = {'Vt': '1/tau*(-(Vt-EL) -S_et*(Vt-Eampa) -S_rt*(Vt-Egaba))\
                           -1/Cm*(gLK*(Vt-Ek) +gTt*m_inf_t(Vt)*m_inf_t(Vt)*hTt*(Vt-Eca) +gh*(mh1 + ginc*mh2)*(Vt-Eh))', # membrane equation
                   'Vr': '1/tau*(-(Vr-EL) -S_er*(Vr-Eampa) -S_rr*(Vr-Egaba))\
                           -1/Cm*(gLK*(Vr-Ek) +gTr*m_inf_r(Vr)*m_inf_r(Vr)*hTr*(Vr-Eca) ) ', # membrane equation
                   'S_et': 'dS_et' ,# synaptic filtering
                   'dS_et': 'gamma_e*gamma_e * (in - S_et)-2*gamma_e*dS_et', # ampa synapse of TC pop
                   'S_rt': 'dS_rt' ,# synaptic filtering
                   'dS_rt': 'gamma_r*gamma_r * (Ntr*Qr(Vr) - S_rt)-2*gamma_r*dS_rt', # gaba synapse of TC pop
                   'S_rr': 'dS_rr' ,# synaptic filtering
                   'dS_rr': 'gamma_r*gamma_r * (Nrr*Qr(Vr) - S_rr)-2*gamma_r*dS_rr', # gaba synapse of RT pop
                   'S_er': 'dS_er' ,# synaptic filtering
                   'dS_er': 'gamma_e*gamma_e * ( Nrt*Qt(Vt) - S_er)-2*gamma_e*dS_er', # gaba synapse of RT pop
                   'hTt': '(hinfT_t(Vt)-hTt)/tau_ht(Vt)',
                   'hTr': '(hinfT_r(Vr)-hTr)/tau_hr(Vr)',
                   'mh1': '(mhinf(Vt)*(1-mh2)-mh1)/tau_hm(Vt) -k3*Ph(Ca)*mh1 +k4*mh2',
                   'mh2': 'k3*Ph(Ca)*mh1 - k4*mh2',
                   'Ca': '-alpha_Ca*gTt*m_inf_t(Vt)*m_inf_t(Vt)*hTt*(Vt-Eca)  - (Ca-Ca0)/tau_Ca'
                   }
# initial conditions
DSargs.ics      = {'Vt': -60,
                   'Vr': -60,
                   'S_et': 10e-3,
                   'S_rt': 10e-3,
                   'S_er': 10e-3,
                   'S_rr': 10e-3,
                   'dS_et': .0,
                   'dS_rt': .0,
                   'dS_er': .0,
                   'dS_rr': .0,
                   'hTt': .05,
                   'hTr': .05,
                   'mh1': .05,
                   'mh2': .02,
                   'Ca': 3e-3}

DSargs.tdomain = [0,12000]                         # set the range of integration.
DSargs.pdomain = {'gh': [.01, .07], 'gLK': [.005, .03] }
#ode  = dst.Generator.Dopri_ODEsystem(DSargs)     # an instance of the 'Generator' class.
ode  = dst.Generator.Vode_ODEsystem(DSargs)
traj = ode.compute('essai2')              # integrate ODE
pts  = traj.sample(dt=.1)                      # Data for plotting

# PyPlot commands
plt.plot(pts['t'], pts['Vt'])
plt.plot(pts['t'], pts['Vr'])
plt.xlabel('time')                              # Axes labels
plt.ylabel('voltage')                           # ...
plt.title(ode.name)                             # Figure title from model name
plt.show()

plt.clf()
plt.plot(pts['t'], 100*pts['Ca'])
plt.plot(pts['t'], pts['mh1'])
plt.plot(pts['t'], pts['mh2'])
plt.plot(pts['t'], pts['hTt'])
plt.plot(pts['t'], pts['hTr'])
plt.xlabel('time')                              # Axes labels
plt.ylabel('voltage')                           # ...
plt.title(ode.name)                             # Figure title from model name
plt.show()

[Rob Clewley]

The diagnostics returned with that error suggest that it's struggling
to compute steps, and possibly the step size is shrinking to zero. I
suspect it might be because you have to be more careful when adding a
noisy input to what's meant to be a stationary, deterministic RHS
definition! The deterministic ODE solver sees this "input" as a
parameter that's changing, which can cause a lack of convergence if
the amplitude is large enough or the frequency of change is high
enough. Can you try running this without the input and check that it
works, first. Then, I'd suggest using a smoother noise than "white
noise".

Does that help?

On Wed, Apr 12, 2017 at 12:32 AM, Michel Besserve
mbesserve@users.sf.net wrote:

Hi Rob,
Thanks for the quick answer. I followed your advice and this works using the
Vode integrator, but it crashes when in try to use Dopri ('No Trajectory
Created' message...). Vode is really slow (I think also due to the added
external interpolated input), so I would like to try whether Dopri could do
better. Otherwise would you have some tips to optimize the introduction of
the external input in the model in order to make the integration faster? I
provide the full code below.
Thanks,
Michel

-- coding: utf-8 --

"""
Spyder Editor

This is a temporary script file.
"""

import PyDSTool as dst
import numpy as np
from matplotlib import pyplot as plt

define interpolated noise to emulate stochasticity

tdomain = [0,12000]
timeData = dst.linspace(tdomain[0], tdomain[1], 100000)
noise_std = .3
sindata = noise_std*np.random.randn(len(timeData))
xData = {'in': sindata}

my_input = dst.InterpolateTable({'tdata': timeData,
'ics': xData,
'name': 'interp'
}).compute('interp')

we must give a name

DSargs = dst.args(name='IsolatedThalamusModel')

parameters

DSargs.pars = {'Qmax': 400e-3,# max firing rate
'theta': -58.5,# firing threshold
'sigmaT': 6.0, # neural gain of TC neurons
'sigmaR': 6.0, # neural gain of RT neurons
'tau': 20.0, # membrane time constant
'EL': -70.0, # leak current reversal potential
'Eampa': 0.0, # reversal potential of excitatroy synapse
'Egaba': -70.0, # reversal potential of inhibitory synapse
'Ek': -100, # potassium
'Eh': -40, # for the rectifier current
'Eca': 120.0,
'gamma_e': 70e-3,
'gamma_r': 100e-3,
'Ntr': 5.0,
'Nrt': 3.0,
'Nrr': 25.0,
'gLK': 0.018, # leak conductance, setting for isolated
thalamus
'gTt': 3.0, # T channel conductance
'gTr': 2.3, #
'ginc': 2.0,
'gh': .06, # conductance of the anomalous rectifier, setting
for isolated thalamus
'Cm': 1.0, # coupling capacitance
'k1': 2.5e7,
'k2': 4e-4,
'k3': 1e-1,
'k4': 1e-3,
'alpha_Ca': 51.8e-6,
'tau_Ca': 10.0,
'Ca0': 2.4e-4
}
DSargs.inputs = {'in': my_input.variables['in']}

auxiliary helper function(s) -- function name: ([func signature],

definition)
DSargs.fnspecs = {'Qt': (['v'], 'Qmax / (1 + exp( -(v-theta)/sigmaT ))'),
'Qr': (['v'], 'Qmax / (1 + exp( -(v-theta)/sigmaR ))'),
'm_inf_t': (['v'], '1 / (1 + exp( -(v+59)/6.2 ))'),
'm_inf_r': (['v'], '1 / (1 + exp( -(v+52)/7.4 ))'),
'hinfT_t': (['v'], '1 / (1 + exp( (v+81)/4 ))'),
'hinfT_r': (['v'], '1 / (1 + exp( (v+80)/5 ))'),
'tau_ht':
(['v'],'(30.8+(211.4+exp((v+115.2)/5))/(1+exp((v+86)/3.2)))/exp(1.2log(3))'),
'tau_hr':
(['v'],'(85+1/(exp((v+48)/4)+exp(-(v+407)/50)))/exp(1.2log(3))'),
'mhinf': (['v'],'1/(1+exp((v+75)/5.5))'),
'tau_hm':
(['v'],'(20+1000/(exp((v+71.5)/14.2)+exp(-(v+89)/11.6)))'),
'Ph':(['C'],'k1CCCC/(k1CCCC+k2)')
}

rhs of the differential equation, including dummy variable w

DSargs.varspecs = {'Vt': '1/tau(-(Vt-EL) -S_et(Vt-Eampa)
-S_rt(Vt-Egaba))\
-1/Cm(gLK(Vt-Ek)
+gTtm_inf_t(Vt)m_inf_t(Vt)hTt(Vt-Eca) +gh(mh1 + gincmh2)(Vt-Eh))', #
membrane equation
'Vr': '1/tau(-(Vr-EL) -S_er(Vr-Eampa)
-S_rr(Vr-Egaba))\
-1/Cm(gLK(Vr-Ek)
+gTrm_inf_r(Vr)m_inf_r(Vr)hTr(Vr-Eca) ) ', # membrane equation
'S_et': 'dS_et' ,# synaptic filtering
'dS_et': 'gamma_egamma_e * (in - S_et)-2gamma_edS_et',

ampa synapse of TC pop

               'S_rt': 'dS_rt' ,# synaptic filtering
               'dS_rt': 'gamma_r*gamma_r * (Ntr*Qr(Vr) -

S_rt)-2gamma_rdS_rt', # gaba synapse of TC pop
'S_rr': 'dS_rr' ,# synaptic filtering
'dS_rr': 'gamma_rgamma_r * (NrrQr(Vr) -
S_rr)-2gamma_rdS_rr', # gaba synapse of RT pop
'S_er': 'dS_er' ,# synaptic filtering
'dS_er': 'gamma_egamma_e * ( NrtQt(Vt) -
S_er)-2gamma_edS_er', # gaba synapse of RT pop
'hTt': '(hinfT_t(Vt)-hTt)/tau_ht(Vt)',
'hTr': '(hinfT_r(Vr)-hTr)/tau_hr(Vr)',
'mh1': '(mhinf(Vt)(1-mh2)-mh1)/tau_hm(Vt) -k3Ph(Ca)mh1
+k4mh2',
'mh2': 'k3Ph(Ca)mh1 - k4mh2',
'Ca': '-alpha_CagTtm_inf_t(Vt)m_inf_t(Vt)hTt(Vt-Eca)
- (Ca-Ca0)/tau_Ca'
}

initial conditions

DSargs.ics = {'Vt': -60,
'Vr': -60,
'S_et': 10e-3,
'S_rt': 10e-3,
'S_er': 10e-3,
'S_rr': 10e-3,
'dS_et': .0,
'dS_rt': .0,
'dS_er': .0,
'dS_rr': .0,
'hTt': .05,
'hTr': .05,
'mh1': .05,
'mh2': .02,
'Ca': 3e-3}

DSargs.tdomain = [0,12000] # set the range of
integration.
DSargs.pdomain = {'gh': [.01, .07], 'gLK': [.005, .03] }

ode = dst.Generator.Dopri_ODEsystem(DSargs) # an instance of the

'Generator' class.
ode = dst.Generator.Vode_ODEsystem(DSargs)
traj = ode.compute('essai2') # integrate ODE
pts = traj.sample(dt=.1) # Data for plotting

PyPlot commands

plt.plot(pts['t'], pts['Vt'])
plt.plot(pts['t'], pts['Vr'])
plt.xlabel('time') # Axes labels
plt.ylabel('voltage') # ...
plt.title(ode.name) # Figure title from model
name
plt.show()

plt.clf()
plt.plot(pts['t'], 100*pts['Ca'])
plt.plot(pts['t'], pts['mh1'])
plt.plot(pts['t'], pts['mh2'])
plt.plot(pts['t'], pts['hTt'])
plt.plot(pts['t'], pts['hTr'])
plt.xlabel('time') # Axes labels
plt.ylabel('voltage') # ...
plt.title(ode.name) # Figure title from model
name
plt.show()

---

Integration of stochastic differential equation

---

Sent from sourceforge.net because you indicated interest in
https://sourceforge.net/p/pydstool/discussion/472291/

To unsubscribe from further messages, please visit
https://sourceforge.net/auth/subscriptions/

[Michel Besserve]

Hi Rob,
Yes it works with a smoother noise! I ll try to optimize the way to set the noise, Thanks for the insight, I did not realize that the solver considers it as a (smoothly varying) parameter.
Best,
Michel