R script: dde.R

$$ \newcommand{\esp}[1]{\mathbb{E}\left(#1\right)} \newcommand{\var}[1]{\mbox{Var}\left(#1\right)} \newcommand{\deriv}[1]{\dot{#1}(t)} \newcommand{\prob}[1]{ \mathbb{P}\!(#1)} \newcommand{\eqdef}{\mathop{=}\limits^{\mathrm{def}}} \newcommand{\by}{\boldsymbol{y}} \newcommand{\bc}{\boldsymbol{c}} \newcommand{\bpsi}{\boldsymbol{\psi}} \def\pmacro{\texttt{p}} \def\like{{\cal L}} \def\llike{{\cal LL}} \def\logit{{\rm logit}} \def\probit{{\rm probit}} \def\one{{\rm 1\!I}} \def\iid{\mathop{\sim}_{\rm i.i.d.}} \def\simh0{\mathop{\sim}_{H_0}} \def\df{\texttt{df}} \def\res{e} \def\xomega{x} \newcommand{\argmin}[1]{{\rm arg}\min_{#1}} \newcommand{\argmax}[1]{{\rm arg}\max_{#1}} \newcommand{\Rset}{\mbox{$\mathbb{R}$}} \def\param{\theta} \def\setparam{\Theta} \def\xnew{x_{\rm new}} \def\fnew{f_{\rm new}} \def\ynew{y_{\rm new}} \def\nnew{n_{\rm new}} \def\enew{e_{\rm new}} \def\Xnew{X_{\rm new}} \def\hfnew{\widehat{\fnew}} \def\degree{m} \def\nbeta{d} \newcommand{\limite}[1]{\mathop{\longrightarrow}\limits_{#1}} \def\ka{k{\scriptstyle a}} \def\ska{k{\scriptscriptstyle a}} \def\kel{k{\scriptstyle e}} \def\skel{k{\scriptscriptstyle e}} \def\cl{C{\small l}} \def\Tlag{T\hspace{-0.1em}{\scriptstyle lag}} \def\sTlag{T\hspace{-0.07em}{\scriptscriptstyle lag}} \def\Tk{T\hspace{-0.1em}{\scriptstyle k0}} \def\sTk{T\hspace{-0.07em}{\scriptscriptstyle k0}} \def\thalf{t{\scriptstyle 1/2}} \newcommand{\Dphi}[1]{\partial_\pphi #1} \def\asigma{a} \def\pphi{\psi} \newcommand{\stheta}{{\theta^\star}} \newcommand{\htheta}{{\widehat{\theta}}}$$

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("
input = {a, b, tau1, tau2}

  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)