Towards a Field Theoretical Stochastic Model for Description of Tumour Growth

Show more

1. Introduction

In a recent letter [1] , we have shown how the standard field theoretical language based on creation and annihilation operators (building blocks of the second quantization method [2] ) may be used for a straightforward derivation of closed master equations [3] to describe the population dynamics of multivariate stochastic epidemic models. This was mainly motivated by the observation that, as remarked in [4] , for the kinds of model studied in population biology and epidemiology, a field theoretical description is notationally neater and more manageable than standard methods, in often replacing sets of equations with single equations with the same content. Indeed, a single hamiltonian function sums up the dynamics compactly, even when births and deaths allow the population size to change, and it may be easily written down from a verbal description of the transitions presented in these models.

In the present work, we employ the very same methodology established in [1] to develop a field-theory inspired stochastic model for description of tumour growth based on an analogy with an SI epidemic model [5] , where the susceptible individuals (S) would represent the healthy cells and the infected ones (I), the cancer cells. From this model, we were able to obtain a curve to describe the time evolution of the tumour volume, which is then compared to available experimental data for liver cancer. Our main motivation comes from the pioneering ideas about cancer as a phase transition presented in [6] , more specifically, from the following observations: (1) progression of cancer must also involve a population-level shift due to a competition between two co-existing phenotypes: normal and cancerous ; (2) given enough time and resources, cancer cells will usually outcompete healthy cells in the organ or tissue where they coexist, in the competition for space and resources.

The rest of this work is organized as follows. In Section 2, we introduce the basic aspects of our model which allow us to obtain differential equations to describe the time evolution of the mean number of individuals in the interacting populations we are dealing with, i.e., normal (healthy) and cancer cells. The analytical solutions for these differential equations, besides the specific curve which describes the time evolution of the tumour volume and its comparison to experimental data are presented in Section 3. Finally, in Section 4 we present our concluding remarks.

2. Building the Model

We will start by considering interacting populations, whose total sizes are allowed to change, composed of two types of individuals: the normal (healthy) cells and the tumour (cancer) cells. Let us introduce $\mathcal{N}\left(t\right)$ and $\mathcal{C}\left(t\right)$ as random variables which represent, respectively, the number of normal and cancer cells at a given time instant t.

We will then consider a bivariate process ${\left\{\mathcal{N}\left(t\right)\mathrm{;}\mathcal{C}\left(t\right)\right\}}_{t=0}^{\infty}$ with a joint probability function given by

${p}_{\left(n\mathrm{,}m\right)}\left(t\right)=\text{Prob}\left\{\mathcal{N}\left(t\right)=n\mathrm{;}\mathcal{C}\left(t\right)=m\right\}\mathrm{.}$ (1)

Our aim is to compute time-dependent expectation values of the observables $\mathcal{N}\left(t\right)$ and $\mathcal{C}\left(t\right)$ , which may be defined in terms of the configuration probability according to

$\begin{array}{l}\langle \mathcal{N}\left(t\right)\rangle ={\displaystyle \underset{n\mathrm{,}m}{\sum}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}n\text{\hspace{0.05em}}{p}_{\left(n\mathrm{,}m\right)}\left(t\right)\mathrm{,}\\ \langle \mathcal{C}\left(t\right)\rangle ={\displaystyle \underset{n\mathrm{,}m}{\sum}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}m\text{\hspace{0.05em}}{p}_{\left(n\mathrm{,}m\right)}\left(t\right)\mathrm{.}\end{array}$ (2)

Let us represent the probabilistic state of the system by the vector

$|\mu ,\nu \rangle ={\displaystyle \underset{n,m}{\sum}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}{p}_{\left(n\mathrm{,}m\right)}\left(t\right)|n,m\rangle ,$ (3)

with the normalization condition ${\sum}_{n,m}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}{p}_{\left(n,m\right)}\left(t\right)=1$ .

As an example, the vector $\frac{1}{4}\left(|1,1\rangle +|2,1\rangle +|1,2\rangle +|2,2\rangle \right)$ represents the probability distribution where there are 1 or 2 healthy/cancer cells present, each one with probability 1/4, i.e. ${p}_{\left(1,1\right)}={p}_{\left(2,1\right)}={p}_{\left(1,2\right)}={p}_{\left(2,2\right)}=\frac{1}{4}$ .

Since the configurations are given entirely in terms of occupation numbers $\left(n\mathrm{,}m\right)$ , which calls for a representation in terms of second-quantized bosonic operators [7] , we will introduce creation and annihilation operators for the normal cells, respectively, ${h}^{\u2020}$ and $h$ , and for the cancer cells, namely, ${c}^{\u2020}$ and $c$ . These operators must satisfy the following commutation relations

$\left[h,\text{\hspace{0.05em}}{h}^{\u2020}\right]=\left[c,\text{\hspace{0.05em}}{c}^{\u2020}\right]=1,$

$\left[h,\text{\hspace{0.05em}}c\right]=\left[h,\text{\hspace{0.05em}}{c}^{\u2020}\right]=\left[c,\text{\hspace{0.05em}}{h}^{\u2020}\right]=\left[{h}^{\u2020},\text{\hspace{0.05em}}{c}^{\u2020}\right]=0.$ (4)

As usual in the second quantization framework, we say that ${h}^{\u2020}$ and ${c}^{\u2020}$ “create”, respectively, normal and cancer cells when applied over the reference

(vacuum) state $|\mathrm{0,0}\rangle $ . This allows us to build our space from basis vectors of the form $|n,m\rangle ={\left({h}^{\u2020}\right)}^{n}{\left({c}^{\u2020}\right)}^{m}|0,0\rangle $ .

This vacuum state has the following properties: $h|0,0\rangle =c|0,0\rangle =0$ (from which “annihilation” operators) and $\langle 0,0|0,0\rangle =1$ (inner product).

Following the above definitions, we also have

$\begin{array}{l}{h}^{\u2020}\text{\hspace{0.05em}}|n,m\rangle =|n+1,m\rangle ;\text{\hspace{0.05em}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}{c}^{\u2020}\text{\hspace{0.05em}}|n,m\rangle =|n,m+1\rangle ,\\ h|n,m\rangle =n|n-1,m\rangle ;\text{\hspace{0.05em}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}c\text{\hspace{0.05em}}|n,m\rangle =m|n,m-1\rangle .\end{array}$ (5)

At this point it is worth to note that ${h}^{\u2020}h\text{\hspace{0.05em}}|n,m\rangle =n|n,m\rangle $ and ${c}^{\u2020}c\text{\hspace{0.05em}}|n,m\rangle =m|n,m\rangle $ . Thus, the operators $n={h}^{\u2020}h$ and $m={c}^{\u2020}c$ just count the number of cells in a definite state. This is the main reason why they are usually called number operators. The vector state of our system may be then rewritten in terms of creation and annihilation operators as

$|\mu ,\nu \rangle ={\displaystyle \underset{n,m}{\sum}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}{p}_{\left(n,m\right)}\left(t\right){\left({h}^{\u2020}\right)}^{n}{\left({c}^{\u2020}\right)}^{m}|\mathrm{0,0}\rangle .$ (6)

The time evolution of our system will then be generated by a linear operator $\ufffd$ (called hamiltonian) which may be constructed directly from the transition rates present in our model according to Table 1 (cf. [4] , Table 1). Note that, upon summing up the terms presented in Table 1, we may write our hamiltonian as

Table 1. Transition rates presented in our model and corresponding terms in the hamil- tonian $\mathcal{H}$ .

$\mathcal{H}=+\left({b}_{h}+{d}_{h}+\lambda \right)n+\left({b}_{c}+{d}_{c}\right)m-\left({b}_{h}{h}^{\u2020}n+{b}_{c}{c}^{\u2020}m+{d}_{h}h+{d}_{c}c+\lambda {c}^{\u2020}h\right).$ (7)

The notational advantage of this field theoretical description is made clear at this point if we observe that, from the above definitions, the equation which represents the flux of probability between states at rates defined by our model (the so-called master equation or forward Kolmogorov differential equation [5] ) takes the very compact form of an imaginary-time Schrödinger-type linear equation, namely

$\frac{\text{d}}{\text{d}t}|\mu ,\nu \rangle =-\mathcal{H}|\mu ,\nu \rangle .$ (8)

We can get, after some algebra, a more common representation for the master equation by substituting the expressions for the hamiltonian, Equation (7), and the vector state, Equation (3), into Equation (8)

$\begin{array}{c}\frac{\text{d}}{\text{d}t}{p}_{\left(n\mathrm{,}m\right)}=-\left[\left({b}_{h}+{d}_{h}+\lambda \right)n+\left({b}_{c}+{d}_{c}\right)m\right]\text{\hspace{0.05em}}{p}_{\left(n\mathrm{,}m\right)}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+{b}_{h}\left(n-1\right){p}_{\left(n-1,m\right)}+{b}_{c}\left(m-1\right){p}_{\left(n,m-1\right)}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+{d}_{h}\left(n+1\right)\text{\hspace{0.05em}}{p}_{\left(n+1,m\right)}+{d}_{c}\left(m+1\right){p}_{\left(n,m+1\right)}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\lambda \left(n+1\right){p}_{\left(n+1,m-1\right)}.\end{array}$ (9)

In order to compute the time-dependent expectation values of the observables $\mathcal{N}\left(t\right)$ and $\mathcal{C}\left(t\right)$ through the above master equation, we will follow the well- established methodology presented in [5] and introduce the following moment generating function (mgf )

$M\left(\varphi ,\theta ;t\right)=\langle {\text{e}}^{\varphi \mathcal{N}\left(t\right)}{\text{e}}^{\theta \mathcal{C}\left(t\right)}\rangle ={\displaystyle \underset{n,m}{\sum}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}{p}_{\left(n,m\right)}{\text{e}}^{n\varphi +m\theta}.$ (10)

Note that from the above equation we have

$\begin{array}{l}{\left[\frac{\partial M}{\partial \varphi}\right]}_{\varphi ,\theta =0}={\displaystyle \underset{n,m}{\sum}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}n\text{\hspace{0.05em}}{p}_{\left(n\mathrm{,}m\right)}=\langle \mathcal{N}\left(t\right)\rangle ,\\ {\left[\frac{\partial M}{\partial \theta}\right]}_{\varphi ,\theta =0}={\displaystyle \underset{n,m}{\sum}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}m\text{\hspace{0.05em}}{p}_{\left(n\mathrm{,}m\right)}=\langle \mathcal{C}\left(t\right)\rangle \end{array}$ (11)

and, in general

${\left[\frac{{\partial}^{k}M}{\partial {\theta}^{k}}\right]}_{\varphi ,\theta =0}=\langle {\mathcal{C}}^{k}\left(t\right)\rangle ;\text{\hspace{0.05em}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}{\left[\frac{{\partial}^{k}M}{\partial {\varphi}^{k}}\right]}_{\varphi ,\theta =0}=\langle {\mathcal{N}}^{k}\left(t\right)\rangle .$ (12)

After multiplying Equation (9) by $exp\left(n\varphi +m\theta \right)$ and summing on $\left(n\mathrm{,}m\right)$ , we are led, after some algebra, to

$\begin{array}{c}\frac{\partial M}{\partial t}={\displaystyle \underset{n,m=0}{\sum}}\frac{\text{d}\text{\hspace{0.05em}}{p}_{\left(n,m\right)}}{\text{d}t}{\text{e}}^{n\varphi +m\theta}\\ =+\left[{b}_{h}\left({\text{e}}^{\varphi}-1\right)+{d}_{h}\left({\text{e}}^{-\varphi}-1\right)+\lambda \left({\text{e}}^{-\varphi +\theta}-1\right)\right]\frac{\partial M}{\partial \varphi}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\left[{b}_{c}\left({\text{e}}^{\theta}-1\right)+{d}_{c}\left({\text{e}}^{-\theta}-1\right)\right]\frac{\partial M}{\partial \theta}.\end{array}$ (13)

Finally, by differentiating the above equation with respect to $\varphi $ and evaluating the result at $\varphi =\theta =0$ we get the following differential equation for $\langle \mathcal{N}\left(t\right)\rangle $

${\left[\frac{{\partial}^{2}M}{\partial t\text{\hspace{0.05em}}\partial \varphi}\right]}_{\varphi ,\theta =0}=\frac{\text{d}}{\text{d}t}\langle \mathcal{N}\left(t\right)\rangle =\left({b}_{h}-{d}_{h}-\lambda \right)\langle \mathcal{N}\left(t\right)\rangle .$ (14)

On the other hand, if we differentiate Equation (13) with respect to $\theta $ and evaluate at $\varphi =\theta =0$ we get the following differential equation for $\langle \mathcal{C}\left(t\right)\rangle $

${\left[\frac{{\partial}^{2}M}{\partial t\text{\hspace{0.05em}}\partial \theta}\right]}_{\varphi ,\theta =0}=\frac{\text{d}}{\text{d}t}\langle \mathcal{C}\left(t\right)\rangle =\lambda \langle \mathcal{N}\left(t\right)\rangle +\left({b}_{c}-{d}_{c}\right)\langle \mathcal{C}\left(t\right)\rangle .$ (15)

3. Comparison to Liver Cancer Data. Fitting the Parameters in the Analytical Solution for $\langle \mathcal{C}\left(t\right)\rangle $

By defining $\langle \mathcal{N}\left(0\right)\rangle \equiv {N}_{0}$ , $\langle \mathcal{C}\left(0\right)\rangle \equiv {C}_{0}$ , ${\beta}_{h}\equiv {b}_{h}-{d}_{h}$ and ${\beta}_{c}\equiv {b}_{c}-{d}_{c}$ we obtain the following analytical solutions for Equations (14) and (15)

$\langle \mathcal{N}\left(t\right)\rangle ={N}_{0}\text{\hspace{0.05em}}{\text{e}}^{\left({\beta}_{h}-\lambda \right)t}$ (16)

and

$\langle \mathcal{C}\left(t\right)\rangle =\frac{\lambda {N}_{0}}{{\beta}_{h}-\lambda -{\beta}_{c}}\left[{\text{e}}^{\left({\beta}_{h}-\lambda \right)t}-{\text{e}}^{{\beta}_{c}t}\right]+{C}_{0}{\text{e}}^{{\beta}_{c}t}.$ (17)

Finally, if we consider that the number of normal cells is approximately constant $\left({\beta}_{h}\approx 0\right)$ and that the volume of a cancer cell (v) is approximately the same of a normal cell, we may write the following expression for the time evolution of the tumour volume

${V}_{c}\left(t\right)=\frac{\lambda {\left({V}_{h}\right)}_{0}}{-\lambda -{\beta}_{c}}\left[{\text{e}}^{-\lambda t}-{\text{e}}^{{\beta}_{c}t}\right]+{\left({V}_{c}\right)}_{0}{\text{e}}^{{\beta}_{c}t},$ (18)

where ${V}_{c}\left(t\right)\equiv v\langle \mathcal{C}\left(t\right)\rangle $ , ${\left({V}_{h}\right)}_{0}\equiv v{N}_{0}$ is the initial volume of the normal tissue, and ${\left({V}_{c}\right)}_{0}\equiv v{C}_{0}$ is the initial tumour volume.

The above expression is compared to experimental data for liver cancer (average tumour volume) in Figure 1. The corresponding data have been ob- tained through the analysis of computed tomography (CT) scans for a set of 34 patients available in [8] .

4. Concluding Remarks

^{1}A comprehensive review of the most common mathematical models for description of tumour growth may be found in [10] .

We would like to finish this work by making a few comments about Equation (18). As far as we know, this is the first time that an expression describing the time evolution of the volume of tumours^{1} is obtained from basic assumptions about cancer as a phase transition and, most certainly, this is the very first model for description of tumour growth built by using standard field theoretical language commonly found in models describing fundamental interactions of elementary particles. In a future work we are going to present a qualitative

Figure 1. The solid line represents Equation (18) for an initial total volume of the patient liver ${\left({V}_{T}\right)}_{0}={\left({V}_{h}\right)}_{0}+{\left({V}_{c}\right)}_{0}=2153.00\text{\hspace{0.05em}}\text{\hspace{0.17em}}{\text{cm}}^{3}$ ; $\lambda =-7.07\times {10}^{-3}\text{\hspace{0.05em}}{\text{months}}^{-1}$ ; ${\beta}_{c}=5.67\times {10}^{-1}\text{\hspace{0.05em}}{\text{months}}^{-1}$ and ${\left({V}_{c}\right)}_{0}=29.87\text{\hspace{0.05em}}\text{\hspace{0.17em}}{\text{cm}}^{3}$ . This curve was obtained by fitting the experimental data available in [8] (square points in the figure), using the software OriginPro 8 [9] .

analysis of the behaviour of solutions for the system of first order linear differential equations composed by (14) and (15) in the phase plane, which we believe will shed more light on how our model is indeed connected to the idea of cancer as a phase transition.

Last but not least, other possible extensions for the present model should consider the inclusion of other kinds of dependence in the birth/death rates of cells, as temperature and/or concentration gradients of toxic carcinogens, for example.

Acknowledgements

The author would like to thank Prof. Jack A. Tuszynski for the valuable discussions and the very stimulating scientific environment shared at University of Alberta’s Li Ka Shing Applied Virology Institute, where the main results of this work were obtained.

References

[1] Mondaini, L. (2015) Second Quantization Approach to Stochastic Epidemic Models. Biomedical Sciences Today, 2, e8.

http://www.mdtcanada.ca/bmst/v2/e8.pdf

[2] Greiner, W. and Reinhardt, J. (1996) Field Quantization. Springer-Verlag, Berlin.

[3] Gardiner, C.W. (2004) Handbook of Stochastic Methods for Physics, Chemistry and the Natural Sciences. 3rd Edition, Springer, Heidelberg.

[4] Dodd, P.J. and Ferguson, N.M. (2009) A Many-Body Field Theory Approach to Stochastic Models in Population Biology. PLoS ONE, 4, e6855.

https://doi.org/10.1371/journal.pone.0006855

[5] Allen, L.J.S. (2011) An Introduction to Stochastic Processes with Applications to Biology. Chapman and Hall/CRC, Boca Raton.

[6] Davies, P.C.W., Demetrius, L. and Tuszynski, J.A. (2011) Cancer as a Dynamical Phase Transition. Theoretical Biology and Medical Modelling, 8, 30.

https://doi.org/10.1186/1742-4682-8-30

[7] Cardy, J.L. and Täuber, U.C. (1998) Field Theory of Branching and Annihilating Random Walks. Journal of Statistical Physics, 90, 1.

https://doi.org/10.1023/A:1023233431588

[8] Lieffers, J.R., et al. (2009) A Viscerally Driven Cachexia Syndrome in Patients with Advanced Colorectal Cancer: Contributions of Organ and Tumor Mass to Whole-Body energy Demands. The American Journal of Clinical Nutrition, 89, 1173.

https://doi.org/10.3945/ajcn.2008.27273

[9] OriginLab Corp (2007) OriginPro 8. OriginLab Corp., Northampton.

[10] Benzekry, S., et al. (2014) Classical Mathematical Models for Description and Prediction of Experimental Tumor Growth. PLOS Computational Biology, 10, e1003800.

https://doi.org/10.1371/journal.pcbi.1003800