Rework delay-coupling implementations

[Richert]

Currently, PyRates allows to define edge delays either within the framework of ordinary differential equations (ODEs) as convolutions with gamma kernels or within a delayed differential equation (DDE) framework via scalar delays. While the former works just fine the way it is currently implemented, the latter comes with some restrictions.

It is currently implemented via the creation of local buffer variables for the delayed state variables. This implementation only allows for the usage of the simple forward Euler algorithm for solving the DDE system, since higher-order methods would require PyRates to add interpolation schemes to the user-supplied equations that have not yet been implemented. Also, the latter would require the definition of 2D interpolation methods for any of the different backends. While that would technically be possible, it would be challenging to ensure that an appropriate interpolation algorithm is used, given the solver algorithm chosen by the user. If this is not ensured, users might think that they solve their DDE system with a 5th order solver algorithm, even though the delay buffer interpolation is only performed at third order, for example.

I thus suggest that we implement an interface to commomly available DDE solvers instead, via the get_run_func method. Such DDE solvers are available via Matlab's dde-biftool, Julia's DifferentialEquations.jl, or jiTCDDE in Python, for example. , and they implement numerical solvers specifically dedicated to solving DDEs. I suggest the following new feature for PyRates that would allow for such an interface:

• add a new backend function past(x,tau) that takes two positional arguments: a state variable x and a delay tau.
• given a choice of backend, translate this function call into the required syntax of a particular DDE solver
• if required, translate the state vector variable y from a 1D vector into a 2D matrix or callable (see requirements of the different DDE solvers listed above)

This implementation would probably work best, if we implement specific backends or backend keyword arguments for each DDE solver, since they all require slightly different implementations of the DDE system. This could be done similarly to the Auto-07p interface of the Fortran backend.

[Richert]

other ideas related to this new feature:

• add the past(x, tau) method automatically to edges that have a state variable as source variable, if an appropriate backend was chosen.
• construct history functions that return the initial states for any time t and add them to the generated function file
• translate past(x,tau_1) calls into something like x_d1 = y_hist[0, tau_1] (will differ from backend to backend) and replace all occurrences in the equations with x_d1

[maedoc]

hi, fellow developer of neural networks simulations at TVB. I was just poking around, noticed

This implementation only allows for the usage of the simple forward Euler algorithm for solving the DDE system

We recommend the 2nd order Heun method because it also avoids interpolation, and it's usually much more efficient for nonlinear models. We don't really bother with higher-order, because there are no convergence results for stochastic delay systems beyond Heun.

[Richert]

Hi @maedoc ,

thank you for the advice! Following that, I did implement Heun's method with the latest PyRates release and made it available with all backends. In addition to that, adaptive step-size solvers for DDEs are available via the Julia backend since PyRates 0.13.0. In general, it is possible to implement delay-coupling via the above mentioned past calls in any backend and generate function equations that can be used with most available DDE solvers. For an example, see the gallery example on readthedocs. This way, higher order methods can be accessed.

We don't really bother with higher-order, because there are no convergence results for stochastic delay systems beyond Heun.

Here, higher order methods can still be of interest, since you can implement mere DDEs without any stochasticity in PyRates. Actually, SDEs are currently barely support by PyRates. The only option users have in this regard is to define noisy time series manually and provide it as extrinsic input to the respective model. By appropriately scaling the noise variance with the integration step-size, this would allow to implement and solve SDEs via Euler's method.