Evolutionary dynamics of cancer in the bone

I don’t know about you, dear reader, but when I was a senior in highschool, I was busy skipping class to play CounterStrike. And I wasn’t even any good at it. Pranav Warman, however, is busy spending his senior year curing cancer. Or at least modeling it. On Friday, David Basanta, Pranav, and I spent much of the evening trying to understand prostate cancer after it has metastasized to the bone. Below, you can see us trying to make sense of some Mathematica calculations.[1]

In this post, I want to sketch some of the ideas that we fooled around with. First is a model of healthy bone. Second is an introduction of the tumour into the system. Third, we will consider a model of a simple chemotherapy as treatment. You might notice some similarities to Warman et al. (2015) and my old discussions of the Basanta et al. (2012) model of tumour-stroma interaction. This is not accidental.

Healthy bone

Cancer breaks your healthy systems. Before you try to understand how it does this, it is good to understand how the system is suppose to work under normal circumstances. Your bones are made mostly of bone (or osseous, for the fancy) tissue. These aren’t cells themselves, but a mineralized connective tissue that is formed by other cells — osteoblasts. When you suffer damage to your bones, not only do you need to form new tissue, but you need general bone remodeling. This combines the activity of two different stromal cells: the tissue creators — osteoblasts — and the tissue destroyers — osteoclasts.

If we let $p_{B}$ be the proportion of osteoblasts in the stroma, and thus $1 - p_{B}$ — the proportion of osteoclasts, then we can write a simple logistic growth equation for the density of bone:

$\dot{B} = B(1 - B)(p_{B} - (1 - p_{B}))$

However, the proportion of osteoblasts and osteoclasts is not static either. They are competing in an evolutionary game against each other. Through a series of intermittent steps, osteoclasts promote osteoblast populations to grow, and osteoblasts inhibit the growth of osteoclast populations. The two strategies are also sensitive to the amount of bone tissue, with bone tissue inhibiting osteoblast growth and encouraging osteoclast growth. From this we can write down the fitness functions for osteoblasts ($w_B$) and osteoclasts ($w_C$):

\begin{aligned} w_B & = (1 - p_B) - B \\ w_C & = - p_B + B \end{aligned}

We can plug these fitnesses into our replicator dynamics to get an ODE for proportion on osteoblasts:

\begin{aligned} \dot{p_B} & = p_B(1 - p_B)(w_B - w_C) \\ & = p_B(1 - p_B)(((1 - p_B) - B) - (B - p_B)) \end{aligned}

This gives us two coupled equations for $p_B$ and B.

The solutions to this system of equations are stable oscillations. If the bone tissue density increases too much, more osteoclasts are recruited to destroy some of the bone, but this recruits osteoblasts to rebuild it which then increases the bone tissue density and the cycle continues. You can see this in the figure at right. In blue is the proportion of osteoblasts in the stroma, and in purple is the density of bone tissue.

The astute evolutionary game theorist will not be surprised by this result. Although the equation for B is a logistic growth, it is mathematically equivalent to replicator dynamics on two strategies. With this in mind, it is easy to see that the system of these two two-strategy replicator dynamics is just a factorization of the Rock-Paper-Scissors game. Just like with RPS, this system can be tuned by adding weights to the interactions between osteoblasts, osteoclasts, and bone tissue. Depending how you set those weights you can then move the center of oscillation, and/or turn the closed orbits into spirals either inwards or outwards. Personally, I would tune the system to have dampened orbits by introducing a small negative self-interaction into the fitness function for osteoblasts.

Introducing the tumour

The next step is to introduce the prostate cancer cells that have metastasized to the bone. These interact with both the osteoblasts and the osteoclasts, but in different ways. The cancer cells (M) encourage the growth of the osteoblast population but do nothing to the osteoclasts, giving us a new fitness function for the osteoblasts:

$w_B = (1 - p_B) - B + aM$,

where a is the strength of the encouragement.

However, as the osteoblasts create bone tissue this does nothing for the tumour, but as the osteoclasts destroy bone, they release products that are useful to the growth of the tumour. Thus, giving us the following logistic growth equation for the tumour:

$\dot{M} = M(1 - M)(r + b(1 - p_B)B)$,

where r is the tumour’s intrinsic growth rate and b is the usefulness of the products that osteoclasts release when destroying bone tissue.

If we add these modifications to our healthy system then we can see the emergence and rapid growth of the tumour and its effects on the bone tissue density and proportion of osteoblasts. You can see an example at right with r = 0.4 and a = b = 0.2. As before, $p_B$ is in blue, B is in purple, and the tumour burden M is in gold. The tumour was introduced into the system at time step 20.

As you can tell from the figure, the result of cancer is to shift the point of oscillation for the bone-blast-clast system to a slightly higher density of bone tissue. In this particular case, the amplitude of oscillations also increased, but by shifting the timing of the first arrival of cancer cells to a different point in the bone-blast-clast cycle, you can also produces decreases in the amplitude. The change in the point around which oscillations happen, however, will depend only on a.

Chemotherapy as time modulation

Of course, the whole point of mathematical oncology is to help treat and prevent cancer. In fact, that is Moffitt’s mission statement: “To contribute to the prevention and cure of cancer”. So no model would be complete without some consideration of intervention. In this case, let us consider a simple non-specific chemotherapy. This sort of chemotherapy targets cells during mitosis. So cancer cells, which are reproducing more rapidly, will feel it more than normal cells which are reproducing at a controlled rate. However, both will be affected. The bone tissue, not being cells, will not be directly affected.

As a simple caricature, let f be the strength of our chemo: when a cell tries to undergo mitosis, it will be killed with probability f. This means that at f = 0.5, there will be no growth in the affected population, since half the time that a cell tries to divide, it succeeds and becomes two cells and half the time it dies and becomes 0 cells. For 0 < f < 0.5, the growth rate is slowed; and for 0.5 < f < 1 the growth rate is reversed. The easiest way to incorporate this into the model is to multiply each cell type’s fitness by a factor of (1 – 2f).

For example, the fitnesses for osteoclasts and osteoblasts become:[2]

\begin{aligned} w_B & = (1 - 2f)((1 - p_B) - B + aM) \\ w_C & = (1 - 2f)(B - p_B) \end{aligned}

It is important to notice that these common factors, can be carried through the equations, giving us the following dynamics for the tumour and the stroma system:

\begin{aligned} \dot{p_B} & = (1 - 2f)p_B(1 - p_B)(((1 - p_B) - B + aM) - (B - p_B)) \\ \dot{M} & = (1 - 2f)M(1 - M)(r + b(1 - p_B)B) \end{aligned}

which are equivalent to our previous equations if we reparametrize time as $\tau = (1 - 2f)t$. In other words, chemotherapy affects cells like time modulation. It slows down time for 0 < f < 0.5, stops it for f = 0.5, and reverses it for 0.5 < f < 1. The intuition is then: since the tumour grew to its size so much faster than the healthy tissue then if we reverse time, it will go back to non-existence faster than the healthy tissue. I guess Doctor Who was an oncologist.

Unfortunately, this interpretation of chemotherapy does not extend to the bone tissue. However, we do want the change in bone tissue to slow down with increased intensity of therapy. The simplest way to implement this is the following:[3]

$\dot{B} = (1 - f)B(1 - B)(p_{B} - (1 - p_{B}))$

Notice how the factor of 2 is missing. Here, the intensity is acting as a proxy for the number or effectiveness of stroma cells. They won’t start doing the opposite of their jobs, but at full strength chemotherapy, the population might be either depleted or the cells stressed to the point of not being able to do their job.

To see an example of these dynamics, consider the tumour from the previous section, with r = 0.4 and a = b = 0.2, appearing at time step 20. Except in this case, Oncologist Who intervened at timestep 33 with a chemotherapy of strength f = 0.6 for 35 time steps. This intervention is coloured in grey.

During the therapy, the tumour burden is drastically decreased. However, since the therapy is stopped before extinction, the tumour quickly rebounds. Unfortunately, the homeostasis of the bone is upset by the chemotherapy. Although I picked an example where the effects are mild. In general, therapy destabilizes the formerly stable fixed point around which we had closed orbits of bone density and osteoblast proportion. B and $p_B$ then spiral away either towards 1 if they were above the fixed point or towards 0 if they were below it. In this case, they were slightly above. When the therapy is stopped, the stable orbits return, but are now of greater amplitude since the population is further away from the fixed point. These oscillations are often further aggravated by the return of the tumour.

Overall, the system is very sensitive to the timing of treatment. If the treatment starts during the peaks of the bone-blast-clast oscillations then the system can quickly spiral towards extremely high or low bone density. However, the system is self-correcting when not under treatment. This means that a staggered regiment of on-off chemo is worth exploring in the model. Of course, it is important to include aspects like chemo-resistant tumour subpopulations. On Friday, we built our model to include this as a replicator dynamic between naive and resistant tumour cells,[4] but I will save that discussion for another day.

Notes and References

1. If you want then you can also puzzle over our Mathematica calculations. The linked script has the equations I describe in this post and provides you with a manipulate environment where you can play with some of the parameters.
2. It is important to note that a global factor of (1 – 2f) is not always a reasonable simplification. Cell fitnesses are affected both by intrinsic reproduction and programmed death or depletion. Chemotherapy creates death during reproduction, but it does not create reproduction during death or depletion. For the tumour cells, this is not unreasonable, since we can assume that they have disabled their programmed death circuits. However, osteoblasts undergo programmed death or depletion during their interaction with bone, so a better fitness function might be:

$w_B = (1 - 2f)((1 - p_B) + aM) - B$.

It is less clear with osteoclasts, if the negative interaction with osteoblasts is due to programmed death or depletion or a decrease in reproduction. The equations used in the post are fine for the latter case, but in the former case, a better fitness function would be:

$w_C = (1 - 2f)B - p_B$.

3. This is the most grievous simplification in the model. A proper analysis would instead model the stroma population explicitly, and couple that to the bone with something like:

\begin{aligned} \dot{B} & = B(1 - B)S(p_{B} - (1 - p_{B})) \\ \dot{S} & = S(1 - s)(w_B + w_C + (1 - 2f)r_S) \end{aligned},

where S is the stroma population and $r_S$ is a net growth rate shared by both osteoblasts and osteoclasts (and thus cancelled out in their relative fitnesses).

4. If you look at the Mathematica code on github then you will see that this consists in introducing an extra replicator dynamic for the proportion of resistent tumour:

$\dot{p_M} = p_M(1 - p_M)(2f(b(1 - p_B) + r) - c)$,

where c is the cost of resistance. The fitness term in the tumour’s logistic growth equation would also have to be modified to account for an average over the two types of tumours, but I will leave the dedicated reader to figure it out on their own or look it up in the linked code.

Basanta, D., Scott, J.G., Fishman, M.N., Ayala, G., Hayward, S.W., & Anderson, A.R. (2012). Investigating prostate cancer tumour-stroma interactions: clinical and biological insights from an evolutionary game. British journal of Cancer, 106 (1), 174-81 PMID: 22134510

Warman, P., Araujo, A., Lynch, C., & Basanta, D. (2015). An evolutionary game theory approach to evolutionary-enlightened application of chemotherapy in bone metastatic prostate cancer. bioRxiv: 030262.