R script: dde.R


# 1 Introduction

A system of delay differential equations (DDEs) can be implemented in a block EQUATION of the section [LONGITUDINAL] of a script Mlxtran.

Mlxtran provides the command delay(x,T) where x is a one-dimensional component and T is the explicit delay. Therefore, DDEs with a nonconstant past of the form \begin{align} \deriv{x} &= f(x(t),x(t-T_1), x(t-T_2), \ldots), \quad \text{for } \ \ t \geq 0\\ x(t) &= x_0(t) \quad \text{for } \ \ \min(T_k) \leq t \leq 0 \end{align} can be solved.

# 2 Example

We consider in this example a model where two functions $$f_1$$ and $$f_2$$ depending on a vector of parameters $$\phi=(a, b, \tau_1, \tau_2)$$ are defined by a system of DDEs. Here, $$\tau_1$$ and $$\tau_2$$ are delays.

\begin{align} \deriv{f_1} &= a\, f_2(t-\tau_2) - b\, f_1(t-\tau_1) \\ \deriv{f_2} &= b\, f_1(t-\tau_1) - a\, f_2(t-\tau_2) \end{align}

with the following initial condition: \begin{align} f_1(t) &= \max(10 + t, 0) \quad {\rm for } \ t\leq 0 \\ f_2(t) &= 3 \ \ \quad {\rm for } \ t\leq 0 \end{align}

This model can easily be implemented with Mlxtran:

dde.model <- inlineModel("
[LONGITUDINAL]
input = {a, b, tau1, tau2}

EQUATION:
t0 = 0
f1_0 = 10+t
f2_0 = 3
ddt_f1 = a*delay(f2,tau2) - b*delay(f1,tau1)
ddt_f2 = b*delay(f1,tau1) - a*delay(f2,tau2)
")

We will use Simulx for evaluating $$f_1$$ and $$f_2$$ every hour between time -10h and time 100h, with $$a=0.08$$, $$b=0.05$$, $$\tau_1=5$$ and $$\tau_2=10$$:

p <- c(a = 0.08, b  = 0.05, tau1 = 5, tau2 = 10)

out <- list(name=c('f1', 'f2'), time=seq(-10,100,by=1))

res <- simulx(model     = dde.model,
parameter = p,
output    = out)

plot1=ggplot(data=res$f1, aes(x=time, y=f1)) + geom_line() plot2=ggplot(data=res$f2, aes(x=time, y=f2)) + geom_line()
gridExtra::grid.arrange(plot1, plot2, ncol=2)