In this section, we will illustrate the application of DDEs to some classic models what we have previously discussed. We will investigate how time delays affect the dynamical behaviors of these systems. In the following we will use the Heun method for DDE integration.

# Heun's method for DDEs with a constant time delay
DDE_Heun <- function(derivs, t0, X0, t.total, dt, tau, ...){
  # derivs: the function of the derivatives 
  # t0: initial time
  # X0: initial condition
  # t.total: total simulation time
  # dt: time step size 
  # tau: a constant time delay 
  # (this illustrate will not work for cases using multiple time delays)
  t_all = seq(t0-tau, t.total, by=dt)   # save data for all previous time points too
  n_all = length(t_all)
  X_all = numeric(n_all)
  
  n_delay = as.integer(tau/dt) # number of delay time steps
  X_all[1:(n_delay+1)] = X0   # a constant initial condition between t in [-tau, 0]
  
  for (i in (n_delay+1):(n_all-1)) {
    k1 = dt * derivs(t_all[i],X_all[i],X_all[i-n_delay],...)
    X_all[i+1] = X_all[i] + k1   # X_all_i+1 is first assigned here, so that the function works for tau = 0
    k2 = dt * derivs(t_all[i+1],X_all[i+1], X_all[i+1-n_delay],...)
    X_all[i+1] = X_all[i] + (k1+k2)/2
  }
  return(cbind(t_all, X_all))   # the output is a matrix of t & N(t) for all time steps
}

Delay logistic growth model

We now consider a delay logistic growth model following

\[\frac{dN}{dt} = r N(t) [1 - \frac{N(t- \tau)}{B}] \tag{4}\] \(B\) is the carrying capacity.

# derivative function for delay logistic growth model
# t: time
# N: population size of the current time N(t)
# N_old: population size of the previous time N(t-tau). 
#        We will use the DDE integrator to specify the constant time delay.
# r: parameter, growth rate
# B: parameter, carrying capacity
dN_logistic_delays <- function(t, N, N_old, r, B) {
   return(r*N*(1-N_old/B))
}

We select \(B = 100\), \(\tau = 1\). When \(r = 0.3\), the system exhibits a normal logistic growth.

r1 = 0.3
results_4 = DDE_Heun(derivs = dN_logistic_delays, t0 = 0, X0 = 1, t.total = 40, 
                       dt = 0.01, tau = 1, r = r1, B = 100)

plot(results_4[,1], results_4[,2], xlab="t", ylab="N", xlim=c(-1,40), ylim=c(0,120),
     type = "l", col=2)
legend("topleft", inset=0.02, legend = c("r = 0.3"), col=2,lty=1, cex=0.8)

When \(r = 2\), the system quickly converges to a stable oscillation.

r2 = 1.7 
results_5 = DDE_Heun(derivs = dN_logistic_delays, t0 = 0, X0 = 1, t.total = 80, 
                       dt = 0.01, tau = 1, r = r2, B = 100)

plot(results_5[,1], results_5[,2], xlab="t", ylab="N", xlim=c(-1,80), ylim=c(0,240),
     type = "l", col=2)
legend("topright", inset=0.02, legend = c("r = 1.7"), col=2,lty=1, cex=0.8)

When \(r = \pi/2\) (a critical point), the system slowly decays to a stable oscillation.

r3 = pi/2 
results_6 = DDE_Heun(derivs = dN_logistic_delays, t0 = 0, X0 = 1, t.total = 200, 
                       dt = 0.01, tau = 1, r = r3, B = 100)

plot(results_6[,1], results_6[,2], xlab="t", ylab="N", xlim=c(-1,200), ylim=c(0,220),
     type = "l", col=2)
legend("topright", inset=0.02, legend = c("r = pi/2"), col=2,lty=1, cex=0.8)

When \(r = 1.5\), the system has a damping oscillation to reach the carrying capacity.

r4 = 1.5
results_7 = DDE_Heun(derivs = dN_logistic_delays, t0 = 0, X0 = 1, t.total = 200, 
                       dt = 0.01, tau = 1, r = r4, B = 100)

plot(results_7[,1], results_7[,2], xlab="t", ylab="N", xlim=c(-1,200), ylim=c(0,220),
     type = "l", col=2)
legend("topright", inset=0.02, legend = c("r = 1.5"), col=2,lty=1, cex=0.8)

One-gene circuit with self-inhibition

In the third example, we consider a one-gene circuit with self-inhibition. We use a fixed time delay \(\tau\) to represent the time cost to induce transcription. The rate equation of the gene circuit is

\[ \frac{dX(t)}{dt} = g_0 + g_1\frac{1}{1+ (X(t-\tau)/X_{th})^n} - kX(t) \tag {5}\]

\(g_0 + g_1\) represents the maximum transcriptional rate (why?), \(g_0\) is the leakage transcriptional rate, \(X_{th}\) is the Hill threshold level of \(X\), \(n\) is the Hill coefficient, and \(k\) is the degradation rate.

# Definition of the inhibitory Hill function
hill_inh <- function(X,X_th,n) {
  a = (X/X_th)**n
  return(1/(1+a))
}

# DDE derivative function for a one-gene circuit with self-inhibition
# t: current time point
# X: current X level (for t)
# X_old: previous X level (for t - tau)
# g0: parameter, leakage transcriptional rate
# g1: parameter, g0 + g1 is the maximum transcriptional rate
# X_th: parameter, Hill threshold 
# n: parameter, Hill coefficient 
# k: parameter, degradation rate
f_1g_self_inh_delays <- function(t, X, X_old, g0, g1, X_th, n, k) {
   return(g0 + g1 * hill_inh(X_old, X_th, n) - k * X)
}

We set \(g_0 = 10\), \(g_1 = 60\), \(X_{th} = 200\), \(n = 4\), and \(k = 0.1\). For different time delays, the circuit exhibits different types of time dynamics.

plot(NULL, xlab="t", ylab="X",
      xlim=c(-20,200), ylim=c(150,450))

tau_all = seq(5,20,5)
label_all = paste0("tau = ", tau_all)
ntau = length(tau_all)

ind = 1
for(tau in tau_all){
  results = DDE_Heun(derivs = f_1g_self_inh_delays, t0 = 0, X0 = 200, t.total = 200, 
                       dt = 0.01, tau = tau, g0 = 10, g1 = 60, X_th = 200,
                       n = 4, k = 0.1)
  lines(results[,1], results[,2], type = "l", col = ind)
  ind = ind + 1
}
legend("topleft", inset=0.02, legend = label_all,
       col=1:ntau, lty=1, cex=0.8)

When \(\tau = 5\), the system converges to the steady state at around 250. When \(\tau\) increases, the time dynamics change to damping oscillation, and eventually stable oscillation. Therefore, with different delays, a self-inhibitory gene may have either a stable steady state or an oscillatory state.

Lotka-Volterra model with delays

In the last example, we consider the Lotka-Volterra model with time delays.

\[\begin{equation} \begin{cases} \frac{dN}{dt} = N(t) [a - bP(t-\tau_1)] \\ \frac{dP}{dt} = P(t) [cN(t-\tau_2) -d] \tag {6}\end{cases} \end{equation}\]

# derivatives of the Lotka-Volterra model with two time delays
# Xs_old inputs the previous system levels
derivs_LV_delays <- function(t, Xs, Xs_old, a, b, c, d) {
  N = Xs[1]
  P = Xs[2]
  N_old = Xs_old[1]
  P_old = Xs_old[2]
  dNdt = N * (a- b*P_old)
  dPdt = P * (c*N_old - d)
  return(c(dNdt, dPdt))
}

For simplification, we consider \(\tau_1 = \tau_2 = \tau\). Here is the modified DDE_Heun function for a multi-variable system.

# Heun's method for DDEs with a constant time delay (for a generic multi-variable system)
DDE_Heun_generic <- function(derivs, t0, X0, t.total, dt, tau, ...){
  # derivs: the function of the derivatives 
  # t0: initial time
  # X0: initial condition (vector)
  # t.total: total simulation time
  # dt: time step size 
  # tau: a constant time delay 
  # (this illustrate will not work for cases using multiple time delays)
  t_all = seq(t0-tau, t.total, by=dt)   # save data for all previous time points too
  nt_all = length(t_all)
  nx = length(X0)
  X_all = matrix(0, nrow = nt_all, ncol = nx)
  
  n_delay = as.integer(tau/dt) # number of delay time steps
  for (i in 1:(n_delay + 1)){
    X_all[i,] = X0   # a constant initial condition between t in [-tau, 0]
  }
  
  for (i in (n_delay+1):(nt_all-1)) {
    k1 = dt * derivs(t_all[i],X_all[i,],X_all[i-n_delay,],...)
    X_all[i+1,] = X_all[i,] + k1   # X_all_i+1 is first assigned here, so that the function works for tau = 0
    k2 = dt * derivs(t_all[i+1],X_all[i+1,], X_all[i+1-n_delay,],...)
    X_all[i+1,] = X_all[i,] + (k1+k2)/2
  }
  return(cbind(t_all, X_all))   # the output is a matrix of t & X(t) for all time steps
}

We set \(a = d = 1\), \(b = 0.03\), \(c = 0.02\), and the initial condition \({X(t = 0)} = (30, 10)\). Without time delays, the system exhibits an oscillation.

a = 1; b = 0.03; c = 0.02; d = 1
X0 = c(30, 10)
# without time delay tau = 0
results_LV_nodelay = DDE_Heun_generic(derivs = derivs_LV_delays, t0 = 0, X0 = X0, t.total = 50, 
                       dt = 0.01, tau = 0, a, b, c, d)

plot(results_LV_nodelay[,1], results_LV_nodelay[,2], xlab="t", ylab="Population", xlim=c(-1,50), ylim=c(0,160),
     type = "l", col=2)
lines(results_LV_nodelay[,1], results_LV_nodelay[,3], col = 3)
legend("topright", inset=0.02, legend = c("N", "P"), col=c(2,3),lty=1, cex=0.8)

Even with very small time delays (\(\tau = 0.01\), just a delay by a single time step), the system starts to deviate from the stable oscillation. The deviation becomes even larger for larger \(\tau\).

# with time delay tau = 0.01
results_LV_delay = DDE_Heun_generic(derivs = derivs_LV_delays, t0 = 0, X0 = X0, t.total = 50, 
                       dt = 0.01, tau = 0.01, a, b, c, d)

plot(results_LV_nodelay[,2], results_LV_nodelay[,3], xlab="N", ylab="P", xlim=c(0,180), ylim=c(0,150),
     type = "l", col=2)
lines(results_LV_delay[,2],results_LV_delay[,3], col =4)
legend("topright", inset=0.02, legend = c("No delay", "with delay"), col=c(2,4),lty=1, cex=0.8)

Note that this approach can be further modified for cases where \(\tau_1\) and \(\tau_2\) are different (think how to make it work).