Hamiltonian systems and closed orbits in replicator dynamics of cancer

June 24, 2016 by Artem Kaznatcheev 2 Comments

Last month, I classified the possible dynamic regimes of our model of acidity and vasculature as linear goods in cancer. In one of those dynamic regimes, there is an internal fixed point and I claimed closed orbits around that point. However, I did not justify or illustrate this claim. In this post, I will sketch how to prove that those orbits are indeed closed, and show some examples. In the process, we’ll see how to transform our replicator dynamics into a Hamiltonian system and use standard tricks from classical mechanics to our advantage. As before, my tricks will draw heavily from Hauert et al. (2002) analysis of the optional public good game. Studying this classic paper closely is useful for us because of an analogy that Robert Vander Velde found between the linear version of our double goods model for the Warburg effect and the optional public good game.

The post will mostly be about the mathematics. However, at the end, I will consider an example of how these sort of cyclic dynamics can matter for treatment. In particular, I will consider what happens if we target aerobic glycolysis with a drug like lonidamine but stop the treatment too early.

Let us recall the reduced and factored dynamics of our system:

$\begin{aligned} \dot{p} & = p(1 - p)(\overbrace{\frac{b_a}{n + 1} - q(b_v - c))}^{\text{gain function for } p} \\ \dot{q} & = q(1 - q)(\underbrace{\frac{b_v}{n + 1}[\sum_{k = 0}^n p^k] - c)}_{\text{gain function for } q} \end{aligned}$

where p is the proportion of glycolytic cells, q is the proportion of VEGF-overproducers among aerobic cells, $b_a$ is the benefit per unit of acidity, $b_v$ is the benefit per unit of vasculature, c is the cost of calling for a unit of vasculature, and n is the interaction group size. Remember that the cost for glycolysis is an opportunity cost in forfeiting the benefits of oxygen from vasculature. In this model, when $\frac{b_a}{n + 1} < b_v - c < cn$ , an internal fixed point exists at (p*,q*) given by the roots of the gain functions for q and p.

CyclesSimplex_marked An example of dynamics in this regime is given at right. The figure shows the proportions of glycolytic cells ( $x_G = p$ ), aerobic VEGF overproducers ( $x_V = (1 - p)q$ ), and aerobic non-producers (defectors, $x_D = (1 - p)(1 - q)$ ). The seven example runs have differing initial tumour composition (p(0), q(0)), but the same micro-environmental parameters: $b_a = 2.5, b_v = 2, c = 1, n = 4$ .

Although the figure illustrates the cycles, it is not proof by itself. However, it does point a way forward. In particular, since the simplex only shows a flow through proportions, one cannot read temporal information — beyond sequence ordering — off the graph. Just from looking at the graph, I can’t tell how long it takes to go around the cycles or if the first half takes the same amount of time as the second half. Most importantly, for the existence of orbits, this information simply doesn’t matter (at least, if we are willing to avoid questions about the temporal extent of orbits). Thus, we can rescale time however we want, and as long as we don’t stop, reverse it, or make it go infinitely fast the resulting orbits would still look the same in this figure.

Let’s do this stretching of time. Since pq(1 – p)(1 – q) is strictly positive for $p,q \in (0,1)$ , dividing the right hand side of our system in equations by pq(1 – p)(1 – q) acts as a dynamic rescaling of time, leaving the orbits unchanged. But it transforms our equations, giving us a new system:

$\begin{aligned} \dot{p} &= \frac{b_a}{(n + 1)q(1 - q)} - \frac{b_v - c}{1 - q} \\ \dot{q} &= \frac{b_v}{(n + 1)p(1 - p)}[\sum_{k = 0}^n p^k] - \frac{c}{p(1 - p)} \end{aligned}$

Note that the right hand side of the top equation depends only on q, and the bottom only on p. This means that we have a Hamiltonian system. In particular, there exists a time-independent Hamiltonian H, such that the above system can be rewritten as $\dot{p} = - \frac{\partial H}{\partial q}$ and $\dot{q} = \frac{\partial H}{\partial p}$ .

We will find this Hamiltonian with a smart guess (or antidifferentiation). To make our life easier, let’s divide $\sum_{k = 0}^n p^k$ by p(1 – p) in our second equation. This yields:

$\dot{q} = \frac{b_v}{n + 1}\Big( \frac{1}{p(1 - p)} + \frac{n}{1 - p} - \sum_{k = 0}^{n - 2}(n - 1 - k)p^k \Big) - \frac{c}{p(1 - p)}$

Then define the following primitives:

$\begin{aligned} R(p) & := \frac{b_v - c(n + 1)}{n + 1}\ln \frac{p}{1 - p} - \frac{b_v}{n + 1} ( n \ln (1 - p) + \sum_{k = 0}^{n - 2} \frac{n - 1 - k}{k + 1}p^{k + 1})\\ S(q) & := \frac{b_a}{n + 1}\ln\frac{q}{1 - q} + (b_v - c) \ln (1 - q) \end{aligned}$

and it is not hard to check that H = R – S works as the Hamiltonian for our system.

As you might remember from classical mechanics and Noether’s theorem, that a Hamiltonian that does not explicitly depend on time conserved energy. In particular if we define an energy E := H(p(0),q(0)) then for all future time t, H(p(t),q(t)) = E. In the case of our system, each energy will pick out an orbit $\{ (p,q) | H(p,q) = E \}$ . To show that the orbits of our replicator dynamic are periodic, we need to show (1) that they remain within $(0,1)^2$ (where the rescaling applies) and that (2) the internal fixed point is unique. Both conditions are met when $b_a < (b_v - c)(n + 1)$ and $b_v < c(n + 1)$ , and can be checked by (1) seeing that H is finite inside $(0,1)^2$ but the limit of H as each of p or q approach 0 or 1 is positive infinity, and by (2) finding that the root of $\frac{\partial H}{\partial p} = 0, \frac{\partial H}{\partial q} = 0$ is unique in $(0,1)^2$ .

Illustration of treatment

TimingGLYSmall Although cycle dynamics are interesting in their own right, they can also be used to illustrate certain dangers in imperfect treatment. If a cell subtype-targeting treatment is strong enough and can be applied long enough to drive that subtype to extinction then that is a viable intervention. But what if it isn’t applied for long enough? I consider this in the figure at right. In both panels, we have the same micro-environmental parameters as in the simplex from earlier: $b_a = 2.5, b_v = 2, c = 1, n = 4$ . And we imagine a highly glycolytic tumour with p(0) = 0.9, q(0) = 0.6. In the top panel, we see the cyclic dynamics that ensue if the tumour is left untreated. For the solid lines, in green we have the proportion of glycolytic cells ( $x_G = p$ ), in red we have the proportion of VEGF over-producers ( $x_V = (1 - p)q$ ) and in blue is the proportion of aerobic non-producers ( $x_D = (1 - p)(1 - q)$ ). For convenience, the dotted red line shows the proportion of VEGF over-producers among aerobic cells ( $q = \frac{x_V}{x_V + x_D}$ ).

In the second panel, we see the effects of treatment. In particular, I introduced a fitness cost of 3 (pretty large number compared to the GLY base fitness) to the glycolytic cells as a caricature of a glycolysis-inhibiting drug like lonidamine. The drug is applied in the grey coloured region of the graph. We can see that by time 3, the drug has greatly reduced the proportion of GLY. Suppose that this reduction is below the detection threshold, but above the extinction threshold (in this case set at $10^{-4}$ ). In this case: not seeing glycolytic cells in his tests, the physician stops administering the drug. Since we are now in a much wider orbit, the stage of competition between the aerobic cells is long enough to drive VOP below the extinction threshold that GLY never reached. With VOP out of the picture, as GLY begin to recover, there is nothing to resist them and they end up taking over the whole population. By time-step 30, the exact opposite of what the physician intended has happened. The glycolytic cells not only rebounded, but they went from being 90% of the population to being the whole tumour.

If, instead, the physician waited until timestep 10 when the tumour was at its lowest proportion of GLY then that same treatment would have been enough to drive them to extinction. Although a near miss there (stopping earlier than 3 time steps), might still be enough to drive VOP extinct in the recovery. To avoid that, one would want to target the tumour when q is at its highest levels. Of course, like most of my models, this is meant as a heuristic and should be used to seed and explore ideas not to guide actual treatment.

As always, if you want to play with this model for yourself then you can find code for it on my GitHub. See the Mathematica notebooks dgLinB_Simplexes.nb to generate the simplexes, and dgLinB_TreatPlayer.nb to target the glycolytic cells.

Hauert, C., De Monte, S., Hofbauer, J., & Sigmund, K. (2002). Replicator dynamics for optional public good games. Journal of Theoretical Biology, 218(2): 187-94.

Filed under Analytic, Models, Preliminary, Technical Tagged with mathematical oncology, public good, replicator dynamics

About Artem Kaznatcheev
From the Department of Computer Science at Oxford University and Department of Translational Hematology & Oncology Research at Cleveland Clinic, I marvel at the world through algorithmic lenses. My mind is drawn to evolutionary dynamics, theoretical computer science, mathematical oncology, computational learning theory, and philosophy of science. Previously I was at the Department of Integrated Mathematical Oncology at Moffitt Cancer Center, and the School of Computer Science and Department of Psychology at McGill University. In a past life, I worried about quantum queries at the Institute for Quantum Computing and Department of Combinatorics & Optimization at University of Waterloo and as a visitor to the Centre for Quantum Technologies at National University of Singapore. Meander with me on Google+ and Twitter.