1. Adams-Bashforth Method

The Adams-Bashforth Method is a two-step technique designed for the integration of a single-variable Ordinary Differential Equation (ODE) of \(N(t)\). At each iteration, the method requires knowledge of two data points, represented as (\(t_n\), \(N_n\)) and (\(t_{n+1}\), \(N_{n+1}\)), to compute (\(t_{n+2}\), \(N_{n+2}\)).

The iterative formula of the Adams-Bashforth Method is expressed as follows:

\[\begin{equation} N_{n+2} = N_{n+1} + \frac{3}{2}hf(t_{n+1}, N_{n+1}) - \frac{1}{2}hf(t_n, N_n) \end{equation}\]

\(h = \Delta t\) is the time step size. Starting from the initial condition \(N_0\), we initiate the process by employing the Euler method for a single step to obtain the next point \(N_1\). Subsequently, the Adams-Bashforth method is applied iteratively to compute \(N_2\), \(N_3\), etc.

(a) Proof the Adams-Bashforth method and show that it is second order (,which means that its accuracy is at \(O(\Delta t^3)\). You will need to use Taylor expansion (for a one-variable function and a two-variable function) and the formula of a total derivative:

\[\frac{df(X,t)}{dt} = \frac{\partial f}{\partial t} + \frac{\partial f}{\partial X}\frac{dX}{dt} \] (b) Provide a pseudocode for the Adams-Bashforth method.

(c) Implement the Adams-Bashforth method. Utilize an existing ODE integrator from numericalR as a template, modifying the code to define a function for the Adams-Bashforth method.

(d) Employ the exponential growth model as an example. Apply the Adams-Bashforth method and compare the ODE simulation with the exact solution obtained through analytical calculation. To assess performance, vary time step sizes and compute:

(1). The Root Mean Square Deviation (RMSD) between the simulation and the exact solution.

(2). The CPU time cost for the method for different time step sizes, for instance, utilizing microbenchmark in R or timeit in Python)

2. Fourth order Runge-Kutta

In part 2B Equation (8), we show a more accurate numerical integrator of ODEs: the fourth order Runge-Kutta method, or RK4. Please proof the method and calculate its accuracy in terms of \(O(\Delta t)\)

3. Hysteresis

Consider a gene circuit with one self-activating gene, described by the following ODE.

\[\frac{dX}{dt} = 10 + 45\frac{X^4}{X^4+200^4} - kX\]

In this problem, we will explore a phenomenon called hysteresis, where the system’s behavior depends on its history.

(1) Simulate, with an ODE integrator of your choice (e.g., RK4), the gene expression dynamics of the circuit by gradually increasing \(k\) from 0.1 to 0.2 for a duration of \(t_{tot} = 1000\), starting from the steady state when \(k = 0.1\). Plot X(t). Hint: to simulate the circuit with uniformly increasing \(k\), you can define the derivative function \(f(X, t) = 10 + 45\frac{X^4}{X^4+200^4} - kX\), where \(k = 0.1 + 0.1*t/t_{tot}\).

(2) Simulate the gene expression dynamics of the same circuit by gradually decreasing \(k\) from 0.2 to 0.1 for a duration of \(t_{tot} = 1000\), starting from the steady state when \(k = 0.2\). Plot X(t). Hint: think what to do for the simulation with uniformly decreasing \(k\).

(3) Plot \(X(k)\) for the above two time trajectories together with the bifurcation diagram obtained from Part 02E. Briefly describe your observations.

4. Another bifurcation

We consider a gene circuit model related to the one discussed in Part 02E, described by the following ODE.

\[\begin{equation} \frac{dX}{dt} = f(X,g) = 10 + g\frac{X^4}{X^4+200^4} - 0.15X \end{equation}\]

Here, \(g\) is the control parameter in the range of (0, 100). Plot the bifurcation diagram of the steady-state \(X\) as the function of the control parameter \(g\).

5. The Spruce Budworm model

Problem 2. The Spruce Budworm model

The spruce budworm, an insect inhabiting spruce-fir forests in the USA and Canada, typically maintains low population levels. However, there are intermittent surges in budworm populations that inflict significant damage on the forests. If we consider the intrinsic growth of budworm populations in isolation, they adhere to the logistic growth model. An additional factor in the dynamics of budworm involves predation by birds. Consequently, we can characterize the dynamics of the budworm using the following equation.

\[ \frac{dN}{dt} = rN(1-\frac{N}{K}) - \frac{\rho N}{N+A}\]

We choose the following parameters: growth rate \(r = 0.2\), carrying capability \(K=1000\), and half saturation population of predation \(A = 250\).

(1) Model implementation. Define a function to compute the derivatives of the Spruce Budworm model. Ensure that the function takes the model parameters \(r\), \(K\), \(A\), and \(\rho\) as arguments. We will utilize the ODE integrator RK4 from numericalR (Python version also available in the folder extra), so the function should also take \(t\) and \(N\) as inputs.

(2) Multistability. We will assess the multistability of the ODE system by sampling different initial conditions \(N(t=0)\) and, for each, we simulate the ODE for a fix time duration of t.total = 100. Explore the multistability of the model for the following three conditions. (a) \(\rho = 40\); (b) \(\rho= 70\); (c) \(\rho = 100\). Based on your simulations, conclude how many steady states are available and their stability.

(3) Bifurcation diagram. Explore the bifurcation diagram of the Spruce Budworm population (\(N\)) by varying the predation parameter \(\rho\) in the range from 0 to 100. Plot the steady-state \(N\) for different \(\rho\) values and annotate the stable/unstable properties of the steady states.

Hint: Aim to obtain a bifurcation diagram similar to the one depicted below. You can use any numerical method discussed in class.

Figure 1
Figure 1

(4) Parameter-varying dynamics. Here, we consider \(\rho\) not at a constant level but instead oscillating over time according to a sinusoidal function.. The following R code generates such \(\rho\) dynamics for an oscillation period of 50, mean signal of 60 and oscillation amplitude of 20. You can directly use the code provided if you use R; otherwise, write your code in Python/MATLAB to achieve the same.

rho_oscillation <- function(t, period, rho_mean, rho_amp){
    return(rho_mean + rho_amp * sin(t/period))
}
period = 50; rho_mean = 60; rho_amp = 20
t_all = seq(0,1000)
rho_all = rho_oscillation(t_all, period = period, rho_mean = rho_mean, rho_amp = rho_amp)
plot(t_all, rho_all, type = "l", col = 2, xlab = "t", ylab = "rho")

Now, devise a new derivative function so that \(\rho\) can be considered as an oscillatory signal. Choose an initial population, e.g., \(X_0 = 100\). Vary the oscillation periods, and for each, simulate the ODE for a long time duration (e.g., t.total = 1000). Explore the differences in the dynamical properties of the system.

Hint: in the new derivative function, first compute \(\rho\) for a given \(t\), and then use the computed \(\rho\) to calculate the derivative.