This enable us dealing with rate models, like Wilson-Cowan models with delay.
I am not quit familiar with packages for converting string equations to proper codes but I know Pydelay package use weave for that.
I can also provide Python and C++ code to solve these type of equations with euler and Runge-Kutta schemes.
Please let me know whether it is possible to add this feature.
Best,
Abolfazl
Hi. This is certainly an interesting topic, and one that we started discussing more than 5 years ago (see this issue: https://github.com/brian-team/brian2/issues/467) 
. I have to admit we did not think about specific algorithms to solve these kind of equations, but rather considered using the standard schemes with access to previous values of variables. But I have to say I don’t know much about DDEs and how to solve them, so pointers/code would definitely appreciated.
I do not so much see the problem in generating the code from the equations. Our code generation machinery takes care of this, and we can probably extend it quite easily to support this use case (if I am not mistaken about what is needed for the numerical solution). In the above-linked issue and similar discussions, we encountered two major roadblocks:
One thing regarding the syntax I am still not sure about: the proposed syntax is super-general in the sense that the delay can be different for different variables, can change over time, etc. This is certainly nice, but I think the most common use case is a single delay per connection. In this case, you could also expect something like
delayed_connection = Synapses(..., '''w : 1
inp_post = rate_pre * w : Hz (summed)''',
delay=5*ms)
to work, but synaptic delays are tied to pathways, i.e. “on_pre” or “on_post” events. But maybe an easier solution than the proposed new syntax for “buffered variables”, would be to somehow integrate things in the Synapses definition. E.g. something like:
delayed_connection = Synapses(..., '''w : 1
inp_post = rate_pre * w : Hz (summed, delay=5*ms)''')
or
delayed_connection = Synapses(..., '''w : 1
inp_post = rate_pre * w : Hz (summed)''')
delayed_connection.inp_post_update.delay = 5*ms
Just to get some idea. What do delays look like in your use case (e.g. do they differ across connections, do they change over time, …?)
And finally, at the risk of losing users to competitors: I think the ANNarchy simulator supports rate-based connections with delays…
Nice.
I usually work on a system of DDE with constant and heterogeneous delays.
Some thing like this for a system of Wilson-Cowan equations:
For a toy example let’s consider the following DDE and step by step we can add more complexity to our solver.
dy0dt = -y0 + y1(t-1) + y1(t-10)
dy1dt = y0 + y1(t-1) - y1
dy2dt = y1 - y1(t-10)
we suppose the delays are integer multiple of time step dt (we can release this restriction later).
For numerical integration we evaluate the state of the system at discrete time grids with equal space (for fixed time step) or unequal ones (adaptive methods).
If delays are integer multiple of time step (the delays need to be larger than the time step)
we can have access to previous states of the system only by using the proper index like y[n-delay_index0][1] at the following equations.
For Euler scheme :
There is no need to touch the standard scheme.
The integrator look like this:
def euler(f, a, b, h):
'''
f: (function) RHS of the DDE
a : [float] start time
b : [float] end time
h : [float] time step
'''
n_iteration = (int)((b-a)/h)
n = len(y)
for i in range(nstart, nstart+n_iteration):
tmp = h * f(y[i], 0.0, i)
y.append(tmp + y[i])
t.append((1+i)*dt-maxdelay)
The equations are :
def dde(x, t, n):
# define dde system of equations
f0p = -y[n][0] * y[n-delay_index0][1] + y[n-delay_index1][1]
f1p = y[n][0] * y[n-delay_index0][1] - y[n][1]
f2p = y[n][1] - y[n-delay_index1][1]
return np.asarray([f0p, f1p, f2p])
where
delays = [1.0, 10.0]
delay_index0 = int(delays[0]/dt)
delay_index1 = int(delays[1]/dt)
maxdelay = max(delays)
we also need to prepare some history for the beginning of the simulation.
at least int(maxdelay/dt) steps. The history values affect the numerical value of the next steps. Usually we know the history or just put a random or fixed value for that.
def set_history(hist, nstart, maxdelay):
N = len(hist)
for i in range(nstart+1):
t.append(-(nstart - i) / float(nstart) * maxdelay)
# x: nstart x N
for i in range(nstart+1):
y.append(hist)
this example use a list of list to record state of the system at each time step which is slow and can be replaced with 2d numpy array.
To make program more memory efficient we probably can store only int(maxdelay/dt) time step.
If we want to have arbitrary delays, then we don’t have the previous points necessarily on the evaluated points on the time grid. So we need to interpolate the state of the system on those points. np.interp is a choice.
The post will be too long to cover other integration schemes. May be It’s better to send you some test and benchmarks. I will be available if there be any ambiguity.
This is the result:
dt = 0.001
y0 = [5.0, 0.1, 1.0]
delays = [1.0, 10.0]
delay_index0 = int(delays[0]/dt)
delay_index1 = int(delays[1]/dt)
maxdelay = max(delays)
nstart = delay_index1
t_sim = 100
y = []
t = []
set_history(y0, nstart, maxdelay)
euler(dde, 0, t_sim, dt)
Reference: look the book appendix
Delay-Coupled Complex Systems and Applications to Lasers
Flunkert , Valentin
I agree with any of suggested syntax.
cheers,
Abolfazl
Many thanks for this detailed reply, that will be very helpful. I will try to look at this more closely next week, I am currently quite busy with other things unfortunately.
Hi,
Dealing with initial discontinuities is what I did not talk about it previously. I think it has been discussed here properly.
