Mathematical Modelling
in Immunotherapy of Cancer

 

Anna Kraut

  
Biomedical problems have become a popular application for mathematical modelling over the last years. In contrast to economic or climate systems, many biological systems are easy to control and manipulate in the laboratory. Population dynamic models for evolution are particularly suitable in this context and are studied in the group of Anton Bovier at the Institute for Applied Mathematics in Bonn. To do efficient simulations, stochastic-deterministic hybrid algorithms are developed. The interplay of mathematical und medical research leads to interesting questions for both sides.

Inspiration from Biomedicine

The development of tumors is an especially interesting application of mathematical models for evolution. Due to the high mutation rates, evolutionary processes can be observed and influenced within a relatively short time frame. Opposed to, for example, the evolution of mammals, model and reality can thus be directly compared and matched to each other. In cooperation with the group of Michael Holzel at the university hospital Bonn, we are studying the immunotherapy of malignant melanoma, a harmful type of skin tumor. Experiments in mouse models are the foundation of the mathematical model. Here we can observe that there are different types of tumor cells. Usually, they can be identified by a particular protein on their cell surface, a so-called antigen. During the course of the treatment, specific immune cells are injected that identify the tumor cells by this antigen and kill them. However, at the same time, the active immune cells are secreting messenger molecules that cause the tumor cells to change their type and no longer present the antigen. Thus, they cannot be attacked any more and the therapy fails. The goal of our collaborative research is to better understand the processes inside the tumor, to try out different therapeutic approaches, as for example a combination with chemotherapy, and to study survival mechanisms of the tumor.

Deduction of the Model

Our mathematical model is a Markov process, meaning a stochastic process whose future evolution only depends on the current state and not on the past of the system. We are considering a so-called individual-based Markov process, where single players interact with each other. The state of the process is described by a population vector that species the number of individuals of dierent traits. For traits $x \in\mathcal{X} $ and $t \ge 0$, $n_x(t)$ denotes the number of individuals with trait $x$ at time $t$ and $n(t)\in\mathbb{N}^\mathcal{X}$. The dynamics of the process are determined by different events, occurring at random, exponentially distributed times, that change the state of the population. For example, every individual with trait $x$ can produce an offspring of the same trait at rate $b_x$. This changes the state from $n(t)$ to $n(t)+\delta_x$, where $\delta_x$ denotes the corresponding unit vector. Since there exist $n_x(t)$ individuals of trait $x$, there are $b_xn_x(t)$ 'birth' events happening within one unit of time, on average. However, they are not occurring regularly but in varying, random intervals. The dynamics of the process are determined by different events, occurring at random, exponentially distributed times, that change the state of the population. For example, every individual with trait $x$ can produce an offspring of the same trait at rate $b_x$. This changes the state from $n(t)$ to $n(t)+\delta_x$, where $\delta_x$ denotes the corresponding unit vector. Since there exist $n_x(t)$ individuals of trait $x$, there are $b_xn_x(t)$ 'birth' events happening within one unit of time, on average. However, they are not occurring regularly but in varying, random intervals. For a general set $\mathcal{E}$ of of possible events, we denote by $v_e$ the change of the population state and by $r_e(n)$ the overall rate at which it occurs. For the example above ($e=\text{birth of trait }x$), we obtain $v_e=\delta_x$ and $r_e(n)=b_xn_x$, which is composed of the individual rate and the number of individuals. A Markov process is characterized by its so-called generator. In this case, it is of the form

$\mathcal{L}\varphi(n)=\sum_{e\in\mathcal{E}}(\varphi(n+v_e)-\varphi(n))r_e(n)$,

for measurable, bounded $\varphi:\mathbb{N}^\mathcal{X}\to\mathbb{R}$. n recent years, different versions of this process have been studied theoretically as a model for evolution. Primarily, the mechanisms of birth, death, mutation, and competition between individuals have been considered. It is of great interest to approximate these processes by simpler processes to be able to derive their long term behavior. The model that we are using to study immunotherapy of melanoma is of the same type. We distinguish different types of individuals, namely common (Tum$^1$) and 'hidden' (Tum$^2$) tumor cells, as well as immune cells and messenger molecules. In addition to birth/cell division, death (natural or due to competition), and mutation, we have to include a number of different mechanisms. For example, immune cells only reproduce in the presence of their target, the common tumor cells. Thus, their overall birth rate is of the form $b_\text{Imm}n_\text{Imm}n_\text{Tum$^1$}$. Furthermore, several events can be coupled. For instance, $m$ messenger molecules are secreted during the proliferation of immune cells, resulting in a state change from $n$ to $n+\delta_\text{Imm}+m\delta_\text{Msgr}$, and when a common tumor cell changes its type, this corresponds to $v_e=-\delta_\text{Tum$^1$}+\delta_\text{Tum$^2$}$ and $r_e(n)=s_\text{Tum$^1$,Tum$^2$}\ n_\text{Tum$^1$}$. The hypothesis on how individual cells and molecules in the tumor tissue interact is summarized in a diagram.

Interaction diagram of the tumor. Each arrow represents a possible event, the formulae exemplify some of the individual rates at which events occur.

Such a model is applicable whenever different players are interacting at exponential rates. After the experiment has been abstracted and the (presumably) most important players and events have been identified, the precise parameters of the system, for example $b_\text{Tum$^1$}$ und $b_\text{Tum$^2$}$, have to be determined. To do so, we evaluate the experimental data and utilize previous medical knowledge. Next, the stochastic process can be implemented as an algorithm to simulate different scenarios.

Stochastic vs Deterministic

The Markov process described above seems to be more realistic than a deterministic model. For example, cells do not divide exactly every couple of hours, but randomly with a certain average rate. For theoretical considerations, this does not impose a problem. However, simulations of the stochastic model can be very time-consuming. We are using a so-called Gillespie algorithm. It is constructed such that each step simulates a single event, like a birth. First, the waiting time until the next event is determined and afterwards it is decided which of the previously described state changes takes place. The interesting thing is that we do not have to generate an exponential waiting time for every possible event and then choose the smallest one to determine the next event. Alternatively, we can utilize the properties of exponentially distributed random variables. The minimum of a finite set of independent, exponential random variables is again exponentially distributed with the sum of all rates as its parameter. Thus, only a single random time has to be generated. The probability that a particular event realizes this minimum is proportional to its respective rate and therefore can be determined by a second random number. In a large population, many of these events are happening within a short time since many cells divide, die etc. As a consequence, the algorithm has to take many steps to simulate an interesting time span. Our way out is a result from theoretical research on individual-based Markov processes. In the eighties, Ethier and Kurtz have shown that the rescaled stochastic process converges to the solution of a system of ordinary differential equations in the limit of large populations. In the above notation, this system is given by

$\frac{d}{dt}n_x(t)=\sum_{e\in\mathcal{E}}(v_e)_xr_e(n(t))$.

In some kind of law of large numbers, the fluctuations around the mean, given by the differential equations, average out. This deterministic system can be used to approximate the stochastic system even for finite populations in certain cases and, with the help of a classical Runge-Kutta method, it is much faster to simulate. However, the solution is not to only look at the deterministic system from now on. Certain effects are portrayed incorrectly or not at all. If we, for instance, consider a small population, like a tumor at the beginning of its growth, it is very important in which order the different events occur. Even if the reproduction rate of the cells is larger than their death rate, and hence the overall growth rate in the deterministic system would be positive, it can happen that the cells die out before they can divide and therefore the tumor becomes extinct. In the deterministic system, this effect can never be witnessed, which is connected to the fact that an exponentially declining function can never reach zero. This is of particular importance if for example only a few tumor and immune cells are left. In this case, it is critical which cell type dies out first. It decides whether the tumor can be defeated or it regrows since it cannot be fought anymore. In our simulations, we therefore combine both methods. Together with the group of Martin Rumpf at the Institute for Numerical Simulation we developed a hybrid algorithm, which simultaneously approximates frequent events with the deterministic system and simulates rare events as in the stochastic system. Thereby, computing time is saved while at the same time random effects in small populations are preserved. The main question now is, in which situations stochastic effects are important. This is not only the case in small populations or with rare events but in general near 'critical points' of the deterministic system. With this we mean states of the process in which small fluctuations have a large effect on the evolution of the population. This can also happen in the interior of the state space if for example a bifurcation occurs, hence if a continuous change of a parameter causes a sudden qualitative change in the behavior of the solution. Our current research deals with how these critical points can be quantified in the best way and how an algorithm can efficiently detect them during the simulation to adaptively change between the stochastic and the deterministic method. Meanwhile, the computing time should not be largely increased to keep the advantage over the purely stochastic simulation.

Back to the Application

A critical point in the interior of the state space can also be observed in the example of immunotherapy of melanoma. In the figure, we can see how the results of the simulations reflect the experimental data. In both cases, it can be observed how the therapy is only successful for smaller tumors. Above a certain threshold size, the immune system cannot cope anymore and the tumor can outgrow the therapy. This threshold marks a bifurcation of the deterministic system and in the stochastic system, both outcomes can be reached by the same initial values.

Tumor diameter over time for different mice in the experiment (left) and varying initial data in the simulation (right). Initially, the tumor grows untreated, therapy is started on day 18. Red curves correspond to cases where the tumor outgrows the therapy, while gray curves mark successfully treated tumors.

The question remains, what we gain from being able to reproduce data that has already been determined experimentally through simulations. It is not a proof but a strong indicator that the underlying hypothesis on the mechanisms of the biological system was correct. Quite often, the experimental data cannot be matched by the first version of the mathematical model - independently of the choice of parameters. Then, the hypothesis, and thus the model, has to be changed and possibly expanded. This way, important mechanisms can be identified and new research questions can be developed, also for our medical collaborators. Furthermore, we can do simulations beyond the experiments to validate them (naturally occuring mutants in simulations instead of artificially introduced mutation in experiments) and make predictions. The latter are, depending on the underlying experimental data and the resulting accuracy of parameters, not quantitatively reliable. However, if certain effects appear independently of the choice of parameters, qualitative predictions can be made.

Literature

[1] M. Hölzel, A. Bovier, and T. Tüting, Plasticity of tumor and immune cells: a source of heterogeneity and a cause for therapy resistance? Nature Reviews Cancer 13, 365-376  (2013)

[2] M. Baar, L. Coquille, H. Mayer, M. Hölzel, M. Rogava, T. Tüting, and A. Bovier. A stochastic model for immunotherapy of cancer. Scientific Reports 6, 24169 (2016)

First published in German in Mitteilungen der Deutschen Mathematiker-Vereinigung, Band 26, Heft 4, Seiten 175–177, ISSN (Online) 0942-5977, ISSN (Print) 0947-4471, DOI: https://doi.org/10.1515/dmvm-2018-0054

Translation by Carolin Kaffiné and Jessica Theisen.

Anna Kraut works in the stochastics department in the group of Prof. Anton Bovier and is writing her doctoral thesis on population dynamic models and their application in medicine.