Gregory J. Kimmel, Philip Gerlee, Philipp M. Altrock

Read the publicationInfiltration into nearby tissue is a hallmark of disease progression in cancer. The rate at which cancer cells migrate can then be a useful prognostic tool for understanding disease stages. In addition, finding ways to control this migration rate by altering the tumor micro-environment could be used therapeutically to improve treatment outcomes by minimizing disease spread.

Therefore, understanding the biological system that governs the speed at which cells move into unoccupied space is warranted. From a simplified view, the effect of cancer spreading outward takes the form of a traveling wave. A classic tool for understanding these waves is using the Fisher-Kolmogorov-Petrovsky-Piskunov (F-KPP) equation (wow that is a mouthful!).

Before diving into the rich history of the F-KPP, I would first like to motivate it by appealing to those familiar with population growth. The logistic model $$ \begin{equation} \dot n = r n \left(1 - \frac{n}{K} \right) \end{equation} $$ has been used for over 150 years to describe population dynamics of cells

Equation (1) can be thought of as a well-mixed model. That is, a population whose spatial dynamics can be neglected. This can happen when the time scale on which movement or migration occurs is much more rapid than the time scales involving replication.

A method ubiquitous in science, a complicated process is first broken up into simpler parts and then more complex dynamics are slowly reintroduced. Sometimes this leads to minimal changes. Other times, new phenomenons arise. A natural extension to (1) is to include the spatial effects that were originally neglected. A common way to introduce spatial effects is to include a diffusion term. The result is the famous F-KPP equation: $$ \begin{equation} \dot u = D \nabla^2 u + r u (1 - u), \end{equation} $$ where we have defined $u = n/K$ to simplify the notation. The diffusion term creates random motion that causes the population density $u$ to spread out evenly through the domain over time. Other modifications to include spatial effects (e.g. directed motion, chemotaxis) exist as well

More general forms of (2) can be introduced that are more problem specific. A recent paper written in collaboration with experiments has been helpful in obtaining analytical estimates for the infiltration speed of cancer cells

Around the same time, Kolmogorov, Petrovsky and Piskunov (KPP) also arrived at (2) when discussing the movement of advantageous genes

The F-KPP equation has become ubiquitous to many areas of science, ranging from biology

To see this, we start from equation (2) with the above Ansatz and obtain $$ \begin{equation} 0 = DU'' + \eta U' + r U(1-U). \end{equation} $$ The unstable steady state is $U = 0$ and the stable state is $U = 1$. Linearizing about $U = 0$ results in $0 = DU'' + \eta U' + rU$, from which it can be shown that the minimum achievable speed is $\eta = 2 \sqrt{r D}$. This is speed predicted from linear theory. This turns out to be the speed that waves approach under many conditions

Yet another approach for analyzing the behavior of the wave speed $\eta$ involves a form of dimensional analysis. Appeal to (3), but introduce a rescaling of the derivative by a factor $\alpha$, to be determined. This leads to: $$ \begin{equation*} 0 = D\alpha^2 U'' + \alpha \eta U' + rU(1-U). \end{equation*} $$ Lets choose $\alpha$ so all terms are of the same order. This leads to $D \alpha^2 \sim \alpha \eta \sim r$. The first and second implies that $\eta \sim D \alpha$ and the first and third imply that $\alpha \sim \sqrt{r/D}$. Therefore $\eta \sim \sqrt{r D}$, which is the same as the analytically predicted speed up to a constant.

When the wave profile and speed is properly captured by linear theory, the waves are known as "pulled" waves. Often though, the wave profile and speed is not captured by linear theory, and is dependent on the full nonlinearity. These are known as "pushed" waves and are much more difficult to analyze and so we typically resort to numerical methods.

The method is justified semi-rigorously in our recent work by appealing to differences in behavior at different length scales

Mathematically, the time $\tau$ can be split into three components $$ \begin{equation} \tau = \tau_{\text{wave formation}} + \tau_{\text{wave propagation}}(d) + \tau_{\text{boundary}}(\vec u), \end{equation} $$ where $\vec u$ is the spatially averaged value "close" to the boundary and $d$ is the distance traveled. We also know that the wave propagation term is proportional to the length of the domain. The key insight now is to recognize that we can take two lengths $L_1$ and $L_2$, compute $\tau_1$ and $\tau_2$, the resulting difference will contain only the wave propagation term. This occurs because the behavior of the boundary and the wave formation details are

After rearrangement (and noting that $\tau_{\text{wave propagation}}(d) = d/\eta$), we arrive at a result which bares a striking resemblance to (4) $$ \begin{equation}\label{eq:new method} \eta = \frac{d_2 - d_1}{\tau_2 - \tau_1} = \frac{\Delta d}{\Delta \tau}. \end{equation} $$ However, the similarity ends at form. Equation (6) is a large improvement in the standard method for calculating the wave speed.

It can be shown (see footnote 1) that our method is first-order accurate in time and space, with accuracy improving with finer discretization or faster wave speeds: $$ \begin{equation} \eta_{\text{predicted}} = \eta_{\text{actual}} \left[ 1 + \eta_{\text{actual}}^{-1} \mathcal O(\Delta t,\Delta x) \right]. \end{equation} $$

**Figure 1:** A: The estimated wave speed using (4) (solid,blue) and (6) (dotted, red). As expected, the standard method begins to converge to the true wave speed, but requires a very large domain. In contrast, (6) is accurate using even the smallest two sizes $L = 12,25$. The predicted wave speed $\eta_0 = 2\sqrt{rD} \approx 0.447$ with $r = 0.05$ and $D = 1.0$. B: Our method gives a robust way to distinguish that a transition to nonlinear wave speeds are occurring and that it is not a domain size issue. (adapted from ref. 13).

The size of the domain needed for the standard method to give accurate results depends strongly on the wave shape as well as the speed. Slower waves, will have more time to properly form, but typically involve higher resolution, a shortcoming not observed with the new method. In figure 2, we see that speed approaches the theoretical speed, but only near the end of the simulation where we also see the calculation begin to give false results as it hits the boundary $(L = 50)$. Further, oscillations in the wave speed prediction are common as we are restricted by the discretization of the domain. We took $\Delta t = 0.01$ and $\Delta x = 0.1$. The relative error ranged from 0.7% - 2.3%. The situation is worse for $L = 25$ which ranged from 3.23% - 4.65%. In contrast taking $L = 25$ and $L = 50$, with $r = 1, D = 0.1$ results in a relative error of 1.41%.

**Figure 2**: Wave speed calculation as the wave forms and moves through the domain. The speed asymptotes to the predicted wave speed as it approaches the proper wave form.

From this we could design experiments to test whether altering these biological parameters can slow down cancer cell infiltration and thus slow disease progression.

1. To see this, let $\eta_0$ be the true wave speed in the limit $\Delta t, \Delta x \to 0$. Taylor expanding the terms in (5) with two different lengths with $h, k \ll 1$, we have
$$
\begin{align*}
\tau^{(1)} & = \tau_0^{(1)} + \left(h \frac{\partial \tau^0_{WF}}{\partial h} + k \frac{\partial \tau^0_{WF}}{\partial k} \right) - \frac{d^{(1)}}{\eta_0^2} \left(h \frac{\partial \eta_0}{\partial h} + k \frac{\partial \eta_0}{\partial k} \right) + \left(h \frac{\partial \tau^0_B}{\partial h} + k \frac{\partial \tau^0_B}{\partial k} \right) \\
\tau^{(2)} & = \tau_0^{(2)} + \left(h \frac{\partial \tau^0_{WF}}{\partial h} + k \frac{\partial \tau^0_{WF}}{\partial k} \right) - \frac{d^{(2)}}{\eta_0^2} \left(h \frac{\partial \eta_0}{\partial h} + k \frac{\partial \eta_0}{\partial k} \right) + \left(h \frac{\partial \tau^0_B}{\partial h} + k \frac{\partial \tau^0_B}{\partial k} \right) \\
\Delta \tau & = \Delta \tau_0 - \frac{\Delta d}{\eta_0^2} \left(h \frac{\partial \eta_0}{\partial h} + k \frac{\partial \eta_0}{\partial k} \right) \\
\frac{1}{\eta} & = \frac{1}{\eta_0} \left[ 1 - \frac{1}{\eta_0} \left(h \frac{\partial \eta_0}{\partial h} + k \frac{\partial \eta_0}{\partial k} \right) \right].
\end{align*}
$$
(jump back)

- Daniel A Charlebois and Gabor Balazsi. Modeling cell population dynamics. In silico biology, 13(1-2):21–39, 2019.
- Matthew D Johnston, Carina M Edwards, Walter F Bodmer, Philip K Maini, and S Jonathan Chapman. Mathematical modeling of cell population dynamics in the colonic crypt and in colorectal cancer. Proceedings of the National Academy of Sciences, 104(10):4008–4013, 2007.
- Gregory J Kimmel, Philip Gerlee, Joel S Brown, and Philipp M Altrock. Neighborhood size-effects shape growing population dynamics in evolutionary public goods games. Communications biology, 2(1):1–10, 2019.
- LM Cook. Oscillation in the simple logistic growth model. Nature, 207(4994):316–316, 1965.
- Pierre-Fran ̧cois Verhulst. Notice sur la loi que la population suit dans son accroissement. Corresp. Math. Phys., 10:113–126, 1838.
- James D Murray. Mathematical biology: I. An introduction, volume 17. Springer Science & Business Media, 2007.
- Gregory J Kimmel, Mark Dane, Laura M Heiser, Philipp M Altrock, and Noemi Andor. Integrating mathematical modeling with high-throughput imaging explains how polyploid populations behave in nutrient-sparse environments. Cancer Research, 80(22):5109–5120, 2020.
- Ronald Aylmer Fisher. The wave of advance of advantageous genes. Annals of eugenics, 7(4):355–369, 1937.
- AN Kolmogorov, IG Petrovsky, and NS Piskunov. Investigation of the equation of diffusion combined with increasing of the substance and its application to a biology problem. Bull. Moscow State Univ. Ser. A: Math. Mech, 1(6):1–25, 1937.
- Oskar Hallatschek. The noisy edge of traveling waves. Proceedings of the National Academy of Sciences, 108(5):1783–1787, 2011.
- Gyorgy Bazsa and Irving R Epstein. Traveling waves in the nitric acid-iron (ii) reaction. The Journal of Physical Chemistry, 89(14):3050–3053, 1985.
- Wim Van Saarloos. Front propagation into unstable states. Physics reports, 386(2-6):29–222, 2003.
- Gregory J Kimmel, Philip Gerlee, and Philipp M Altrock. Time scales and wave formation in non-linear spatial public goods games. PLoS computational biology, 15(9):e1007361, 2019.