Note

Go to the end to download the full example code

# 1.3. Delay Coupling

In this tutorial, we will go through the different options you have to define delay-coupled dynamical systems in PyRates. As an example, we will use a model of \(N = 2\) interconnected leaky integrator units, the evolution equations of which are given by

with individual rates \(x_i\), a global time scale \(\tau\) and coupling weights \(J_{ij}\). This model does not include any delays in the coupling yet. Below, we will go through the two main options for introducing delay-coupling to this model.

## 1.3.1. References

As preparation, lets import all required packages and load the leaky integrator model into PyRates.

```
# external imports
import numpy as np
import matplotlib.pyplot as plt
# pyrates imports
from pyrates import CircuitTemplate, NodeTemplate, clear
# node definition
li = NodeTemplate.from_yaml("model_templates.base_templates.tanh_node")
```

To get an idea of how the model behaves without delays, lets connect two leaky integrators into a circuit and simulate its dynamics.

```
# define connection weights
J_21 = 5.0
J_12 = -5.0
# define circuit model
net = CircuitTemplate(name="nodelays", nodes={"p1": li, "p2": li},
edges=[("p1/tanh_op/m", "p2/li_op/m_in", None, {"weight": J_21}),
("p2/tanh_op/m", "p1/li_op/m_in", None, {"weight": J_12})])
# define simulation parameters
T = 10.0
dt = 1e-5
dts = 1e-3
times = np.linspace(0, T, int(np.round(T/dt)))
inp = 1.0/(1.0 + np.exp(10.0*np.sin(2.0*np.pi*0.7*times)))
# plot input
plt.plot(times, inp)
plt.title("Input")
plt.show()
# perform simulation
res = net.run(simulation_time=T, step_size=dt, sampling_step_size=dts, vectorize=True, in_place=False,
inputs={"p1/li_op/u": inp}, outputs={"p1": "p1/li_op/r", "p2": "p2/li_op/r"})
# plot the results
plt.plot(res)
plt.legend(res.columns.values)
plt.title("LI Signals: No Delays")
plt.show()
clear(net)
```

We can see that the system dynamics mostly follow the periodic input. Lets see how they are affected by the introduction of delays.

### 1.3.1.1. Option 1: Using Delayed Differential Equations

PyRates allows for the implementation of delayed differential equation (DDE) systems. In the case of our exemplary model, we would like to implement a leaky integrator system with the following evolution equations:

where \(d_{ij}\) are scalar delays specific for the connection from \(j\) to \(i\).
This can be achieved by using the `delay`

keyword to define \(d_{ij}\) for an edge.

```
# define delays
d_21 = 0.2
d_12 = 0.3
# define circuit model with discrete delays
net = CircuitTemplate(name="dde", nodes={"p1": li, "p2": li},
edges=[("p1/tanh_op/m", "p2/li_op/m_in", None, {"weight": J_21, "delay": d_21}),
("p2/tanh_op/m", "p1/li_op/m_in", None, {"weight": J_12, "delay": d_12})]
)
# perform simulation
res = net.run(simulation_time=T, step_size=dt, sampling_step_size=dts, vectorize=True, in_place=False,
inputs={"p1/li_op/u": inp}, outputs={"p1": "p1/li_op/r", "p2": "p2/li_op/r"})
# plot the results
plt.plot(res)
plt.legend(res.columns.values)
plt.title("LI Signals: Scalar Delays")
plt.show()
clear(net)
```

It seems like the addition of discrete delays led to the emergence of a stable periodic solution, the frequency of which differs from the driving frequency. Note that the definition of DDEs does not permit the use of adaptive step-size solvers for most backends at the moment. An exception is the Julia backend, which can be used in combination with the option solver=’julia_dde’ in order to solve a DDE system via any adaptive step-size algorithm available via DifferentialEquations.jl. Alternatively, it is possible to generate the run function for a DDE system via any backend and employ third-party software packages to solve it. An example of that can be found in the model analysis gallery.

### 1.3.1.2. Option 2: Using Distributed Delays

As an alternative, you can define distributed delays in PyRates as convolutions of the source variable of an edge with a gamma kernel. This procedure is explained in detail in [1] and has been applied to a system of coupled QIF neurons in [2]. In our specific example, this would change the model equations as follows

where \(*\) is the convolution operator and \(a_{ij}\) and \(b_{ij}\) are the parameters of the gamma
kernel \(\Gamma_{ij}\). In PyRates, such a convolution can simply be added by specifying an additional keyword
`spread`

for a given edge definition:

```
# define variances
v_21 = 0.1
v_12 = 0.2
# define circuit model with distributed delays
net = CircuitTemplate(name="gamma", nodes={"p1": li, "p2": li},
edges=[("p1/tanh_op/m", "p2/li_op/m_in", None, {"weight": J_21, "delay": d_21, "spread": v_21}),
("p2/tanh_op/m", "p1/li_op/m_in", None, {"weight": J_12, "delay": d_12, "spread": v_12})]
)
```

In that case, `delay`

and `spread`

are interpreted as the mean \(\mu\) and variance \(\sigma^2\)
of the gamma kernel that the source variable should be convoluted with, respectively.
These quantities are related to \(a\) and \(b\) via \(\mu = \frac{a}{b}\) and
\(\sigma^2 = \frac{a}{b^2}\). In addition, PyRates automatically uses the linear chain trick described in [1]
and [2] to translate the convolution operation into a set of coupled ODEs that will be added to the model equations.
Let’s see how these changes to our model affected the dynamics.

```
# perform simulation
res = net.run(simulation_time=T, step_size=dt, sampling_step_size=dts, vectorize=True, in_place=False,
inputs={"p1/li_op/u": inp}, outputs={"p1": "p1/li_op/r", "p2": "p2/li_op/r"}, solver="scipy")
# plot the results
plt.plot(res)
plt.legend(res.columns.values)
plt.title("LI Signals: Distributed Delays")
plt.show()
clear(net)
```

We can see that the period of the periodic solution of the model dynamics were sped up slightly by the addition of distributed delays in comparison to scalar delays.