The aim of this assignment is to study one of the following articles:

  • A Validated Mathematical Model of Cell-Mediated Immune Response to Tumor Growth (pdf)
  • A general framework for modeling tumor-immune system competition and immunotherapy: Mathematical analysis and biomedical inferences (pdf)
  • Environmental Brownian noise suppresses explosions in population dynamics (pdf)
  • Oscillatory Regulation of Hes1: Discrete Stochastic Delay Modelling and Simulation (pdf)
  • Nonlinear Dynamics of Immunogenc Tumors: Parameter Estimation and Global Bifurcation Analysis (pdf)
  • Mathematical modelling of stem cell differentiation: the PU.1–GATA-1 interaction (pdf)
  • Cancer cell population growth kinetics at low densities deviate from the exponential growth model and suggest an Allee effect (pdf)

Codes from immunogeneic_tumor_growth on github can be useful.

Exercise 1 - The model

  1. Describe the main question addressed in the article: present the general biological setting, and the key mechanismss that are studied in the article. Define all terms that do not have an obvious meaning, for example tumour dormancy, Allee effect. Use the descriptions and definitions from the article, but also other sources from the scientific literature. Cite all sources with a complete list of refence at the end of the report.

  2. Describe the mathematical model or models. Present the variables, parameters and equationss. Describe the biological mechanism modelled by each term in the equations. Simplifications are often made to reduce the original model to a simpler one. List the hypotheses that were made to achieve the simpler model.

  3. List all parameter values in a table (specific values or ranges)

Exercise 2 - Steady states

We want to find and characterise the steady states of the model. To preserve biological realism of the model, we are looking for non-negative steady states (0 included).

  1. Following the methods used in the article you study, find all stationary, steady states of the system. Present suitable condtions for the existence of non-negative steady states. In particular, we would like to know how many steady states can co-exist for a given parameter set.

Exercise 3 - Phase portrait

We now try to establish the phase portrait of the model. To do this, we need first to determine the steady states, and characterise their types: (saddle, node, focus, ...). and their stability. Second, we need to determine the stable and unstable manifolds. In particular, the stable manifold of a saddle will separate the phase space into two part, so that every trajectory is completely contained in one of those parts.

1 Phase portrait with default parameters. Plot the phase protrait with default model parameters. Follow the steps:

  • Identify (analytically or numerically) all the steady state
  • Linearise the model around each of the steady states. (Determine the Jacobian matrix at each steady state).
  • Compute the eigenvalues and the eigenvectors associated to each steady state.
  • Characterise (asymptotic stability and type) of each steady state.
  • Plot in the phase space a few selected trajectories (try to select trajectories that cover most of the phase spaces are reach the stable steady states).
  • Plot the stable and unstable manifolds of the saddles. (No need to plot the manifolds for the oother type of steady states). To do that, use the information given by the eigenvectors of the Jacoban matrix. Let $p \in \mathbb{R^2}$ a saddle point with $\lambda_1 < 0$, $\lambda_2 > 0$, and $v_1$ $v_2$ the associated eigenvectors. Then, in a neighbourhood of $p$, the stable and unstable manifolds are tangent to the stable and unstable subspaces of the Jacobian, which are spanned by the eigenvectors. To plot the unstable manifold, you can take an initial condition on the corresponding linear subspace: $p \pm \epsilon v_2$, avec $\epsilon > 0$ small. To plot the stable manifold, you can take an initial condition on the corresponding linear subspace: $p \pm \epsilon v_2$, avec $\epsilon > 0$ small, and integrate in negative times.

For the stochastic models, use the determinic ODE versions to plot characterise steady states and the stable and unstable manifolds. Use the stochastic version to compute the trajectories

You can use or adapt the codes immunogeneic_tumor_growth.m on immunogeneic_tumor_growth.

Exercises 4 Bifurcations of steady states

Choose a key parameter of your system as your bifurcation parameter, which we will call $\lambda$. This can be one used in the article, but you can also choose another one.

Plot the bifurcation diagram for the model, with respect to $\lambda$. If your model has several variables, plot one bifurcation diagram for each variable.

Characterise all the bifurcations you see: saddle-node, transcritical, pitchfork, or Hopf. Are there any other types of bifurcations?

Exercise 5 Global bifurcations and two-parameters bifurcations

Depending on your model, global bifurcations may occur, such as heteroclinic bifurcations. Try to characterise such bifurcations by plotting the phase portrait before, at, and after the bifurcation.

It can also be interesting to use second bifurcation parameter to follow saddle-node points and see whether they can form a cusp.

Exercise 6 Go further

If your model is an ODE model, you can re-write it as a birth and death process. To convert the continuous ODE variables $x_1, ..., x_n$ to discrete values, use a reference volume $\Omega$ such that the number each species is $N_i = \Omega x_i$. The volume $\Omega$ is a parameter of the birth and death process you will have to take into account.

  • List all the possible birht death events, and provide the rate at which each event occurs.
  • Implement the process with the Stochastic Simulation Algorithm (SSA).
  • Compute the distribution of the birth/death process by running the SSA many times. In birth/death processes, the initial conditions are also random variables. Try to specify biologically relevant initial conditions for the simulations.

You can use or adapt the codes stoch_tumor_growth.m on immunogeneic_tumor_growth.

As discussed in class, the rule-of-thumb for the variability of species numebrs is one over the square root of the mean species number $1/\sqrt{\langle n \rangle}$. Choose the volume $\Omega$ such that the mean numbers are small (10-100) and large (over 1000).

Another possibility is to expand, or simplify the model. Based on the modellling hypotheses listed in exercise 1 and on the discussion in the original article, find a modification to the model that could provide either a more realistic picture of the biological system or a simpler mathematical analysis. Highlight the main differences between the original and the modifief model.