Tumour Control Probability in Cancer Stem Cells Hypothesis

The tumour control probability (TCP) is a formalism derived to compare various treatment regimens of radiation therapy, defined as the probability that given a prescribed dose of radiation, a tumour has been eradicated or controlled. In the traditional view of cancer, all cells share the ability to divide without limit and thus have the potential to generate a malignant tumour. However, an emerging notion is that only a sub-population of cells, the so-called cancer stem cells (CSCs), are responsible for the initiation and maintenance of the tumour. A key implication of the CSC hypothesis is that these cells must be eradicated to achieve cures, thus we define TCPS as the probability of eradicating CSCs for a given dose of radiation. A cell surface protein expression profile, such as CD44high/CD24low for breast cancer or CD133 for glioma, is often used as a biomarker to monitor CSCs enrichment. However, it is increasingly recognized that not all cells bearing this expression profile are necessarily CSCs, and in particular early generations of progenitor cells may share the same phenotype. Thus, due to the lack of a perfect biomarker for CSCs, we also define a novel measurable TCPCD+, that is the probability of eliminating or controlling biomarker positive cells. Based on these definitions, we use stochastic methods and numerical simulations parameterized for the case of gliomas, to compare the theoretical TCPS and the measurable TCPCD+. We also use the measurable TCP to compare the effect of various radiation protocols.


Mathematical Model
With the goal of developing a fully stochastic model for tumor growth and treatment by radiation, we first define a hierarchy of the heterogenous cell populations within a tumor that are central for further analysis. Following Refs. [1,2], we consider three populations, stem S, progenitor P , and mature M cells. Fundamentally, stem cells differentiate into progenitor cells, which differentiate into mature cells, or S → P → M . However, we note that while stem cells have the capacity for unlimited division, progenitor cells divide only a limited number of times, and mature cells do not divide. Thus, we assume that a progenitor cell divides exactly N times before finally differentiating into a mature cell M . That is, we have the modified hierarchy S → P 1 → . . . → P N → M [1,2]. Additionally, we note that any type of cell may be killed due to radiation, and we assume that it occurs with rate Γ i , for i = s, p, m, representing stem, progenitor, and mature cells, respectively. The model includes the following division pathways (note that r 1 + r 2 + r 3 = 1, and i = 1, . . . , N ): Here ρ s and ρ p i refer to the proliferation rates for stem and progenitor cells, respectively, and r 1 , r 2 , r 3 refer to the probabilities of each type of stem cell division (symmetric self-renewal, asymmetric self-renewal, or symmetric commitment). For the purposes of this model, as well as the subsequent analysis, we assume that the cell deaths are independent of one another.
In order to compute the TCP, we first define the joint probability function for the system. That is, we define the probability that the system contains a given number of each type of cell at the time t [3]. We make the assumption that at the initial time t 0 , the number of each type of cell is known and these values are denoted n 0 S , n 0 P 1 , . . . , n 0 P N , n 0 M , for stem cells, each of the generations of progenitor cells, and mature cells, respectively. We denote this probability function by p n S ,n P 1 ,n P 2 ,...,n P N ,n M (t). The corresponding master equation [3] is then given by: We have omitted indices of p(t) that remain unchanged for brevity. The initial condition is given by (δ i,j is the Kronecker delta function): p n S ,n P 1 ,n P 2 ,...,n P N ,n M (t 0 ) = δ n S ,n 0 S δ n P 1 ,n 0 From this, we can obtain master equations for the probability functions for the number of stem cells only, as well as for the number of marker-positive cells (S, P 1 , P 2 , . . ., P k ). Denoting the probability function for the number of stem cells, as u n S (t), we observe that, by definition: u n S (t) = n P 1 ,...,n P N ,n M ≥0 p n S ,n P 1 ,n P 2 ,...,n P N ,n M (t) Thus we obtain with the initial condition u n S (t 0 ) = δ n S ,n 0 S . Following Van Kampen [3], and Zaider and Minerbo [5], we solve for u n S (t) analytically by introducing the generating function U (s, t) = ∞ i=0 u i (t)s i . We then obtain the following partial differential equation(PDE) for U (s, t): This PDE can be solved via the method of characteristics [4], along with the initial condition U (s, 0) = s n 0 S . Assuming that all parameters are constant, we obtain U (s, t) = (s − 1)(ρ s r 3 + Γ s )e (ρs(r 1 −r 3 )−Γs)t − ρ s r 1 s + ρ s r 3 + Γ s (s − 1)(ρ s r 1 )e (ρs(r 1 −r 3 )−Γs)t − ρ s r 1 s + ρ s r 3 + Γ s n 0 S .
From this, we observe that the probability that there are no stem cells remaining at time t is given by We thus define this quantity to be the theoretical TCP, denoted TCP S (t).
Using the same technique, we can obtain the probability function for the number of marker positive cells, denoted v n S ,n P 1 ,...,n P k (t). We have: v n S ,n P 1 ,...,n P k (t) = n P k+1 ,...,n P N ,n M ≥0 p n S ,n P 1 ,n P 2 ,...,n P N ,n M (t) .
This gives (where again we have omitted the indices which remain unchanged, for brevity): The initial condition is v n S ,n P 1 ,...,n P k (t 0 ) = δ n S ,n 0 S δ n P 1 ,n 0 P 1 . . . δ n P k ,n 0 P k . To solve for v(t), similar to the previous case, we introduce the generating function, We can obtain the following PDE for V (x S , x 1 , . . . , x k , t):

Results
The model defining the effects of radiation on the cell populations is essentially the same as that defined in [7]. Other model parameters were taken from [8], except for β s , β p , which were calculated from data collected in [7], and ω, which was estimated specifically for glioblastoma treatment. The parameters used are as follows: ρ s = 0.6931 (1/day), ρ p = ρ p i = 0.6931 (1/day), r 1 = 0.15, r 2 = 0.7, r 3 = 0.15, α s = 0.2 Gy −1 , β s = 0.02 Gy −2 , α p = 0.2 Gy −1 , β p = 0.02 Gy −2 , n 0 S = 100, n 0 P 1 = 100, n 0 P 2 = 100, n 0 P 3 = 100, ω = 15 min, and N = 7. The hazard function, a function describing radiation-induced cell death, is based on the linear-quadratic (LQ) model for cell survival for each dose fraction [5]. Thus, the hazard function for an individual dose fraction of radiation, with dose d and irradiation period ω, takes value: during the period of irradiation (0 ≤ t ≤ ω), and for all t > ω or t < 0, is identically 0 (with j = s, p). Using this, we define the overall function describing the rate of cell death for the entire treatment, consisting of all dose fractions as: In the above expressions, α s , α p , β s , β p are parameters describing the radiosensitivities of stem (s) and progenitor (p) cells in the LQ model, t i is the time of the i th fraction, d i is the dose administered on the i th fraction, and the sum is over all fractions administered.