We begin with the spatially homogeneous, cancer stem cell model developed by Hillen et al. (2013). By spatial homogeneity, we mean that cell density, cell growth, and the distribution of chemicals are homogeneous throughout the tumor region.

where U(t) is the volume fraction of cancer stem cells (CSCs) with respect to the total domain of interest, which contains both tumor and host cells. Similarly, V(t) is the volume fraction of non-stem tumor cells (TCs) with respect to the total domain of interest. The total volume fraction of the tumor is represented by P(t), that is, P(t) = U(t) + V(t). The parameter δ is the probability that a CSC will give rise to another CSC, when it divides. Thus, 1 − δ is the probability that a CSC will give rise to one CSC and one TC, when it divides. It is assumed that the parent CSC remains (Sell, 2004). The growth rates of the CSCs and TCs are given by m_U and m_V, respectively. The apoptosis rate of the TCs is represented by a_V; we assume that CSCs do not undergo apoptosis since they have unlimited replicative potential. Cell growth and differentiation are tempered by k(P(t)), which is essentially a volume constraint. Hillen et al. (2013) assume that k(P) is monotonically decreasing in P and piecewise differentiable, and they set k(P) > 0 for P ∈ [0,P*) and k(P) = 0 for all P ≥ P*, for some P* > 0. For the purposes of this paper, we adopt the version of k(P) used by Hillen et al. (2013) and assume normalization of P* = 1 limiting P to a maximum volume fraction of one, and k(0) = 1. For simulations we use: (3)kP=max1-P4,0

To match the notation used by Youssefpour et al. (2012) we set δ = 2p − 1. It follows that 1 − δ = 2(1 − p). In this case, p is the probability that a CSC gives rise to two CSCs, rather than two TCs, when it divides. That is, p is the probability that a CSC renews itself, and 1 − p is the probability that a CSC differentiates. While this model of CSC division ignores asymmetric division, it is equivalent to the model in equations (1) and (2), as shown in the Appendix of Hillen et al. (2013). The resulting model is given in equations (4) and (5).

As noted in Hillen et al. (2013), if there are no CSCs, the TC population is governed by the equation V°(t)=mVk(V(t))V(t)-aVV(t) and since k is assumed to be decreasing, the TC population is fated to die out if m_Vk(0) < a_V. Further, if we assume that k is strictly decreasing, then the TC population dies out if (6)mVk0≤aV since either V = 0 and the TCs are already extinct, or V > 0 and so k(V) < k(0) and m_Vk(V) < a_V for all V > 0.

The steady states of the model defined in (1, 2) are discussed in detail in Hillen et al. (2013), where it is assumed that the growth rates m_U and m_V are both one and that the TC apoptosis rate a_V is greater than zero. Here, we give the main results, which also apply to the model as stated in (1, 2) or equivalently in (4, 5). We note that in the untreated tumor, we assume δ ∈ (0, 1), that is p ∈ (0.5, 1), such that (2p − 1) > 0. The steady states of the system are X0=0,0,XV=0,V0,XU=1,0,withkV0=aVmv.

The origin, X₀, has eigenvalues λ₁ = (2p − 1)m_Uk(0) > 0 and λ₂ = m_Vk(0) − a_V. Thus, X₀ is an unstable steady state. The TC only steady state, X_V, occurs where V₀ solves m_Vk(V₀) = a_V and has eigenvalues λ₁ = (2p − 1)m_Uk(V₀) > 0 and λ₂ = m_Vk′(V₀)V₀. Thus, X_V is also unstable. The linearization for the pure CSC steady state, X_U, has negative trace, m_Uk′(1) − a_V, and positive determinant, −a_V(2p − 1)m_Uk′(1), thus both eigenvalues are negative, and X_U is a stable steady state. Hillen et al. (2013) have shown that X_U is globally asymptotically stable in the biologically relevant region where U ∈ [0, 1], V ≥ 0, and U + V ≤ 1.

Hillen et al. (2013) derive the slow manifold of the system defined by (1, 2), where they take m_U = m_V = 1. The slow manifold is a subset of the phase space ((U, V)-plane) which describes the long time dynamics of the system. As shown in Hillen et al. (2013), solutions very quickly converge to the slow manifold, and then they slowly follow this manifold getting closer to the attractor at (1, 0). This same slow manifold applies to the system defined by equations (4) and (5). We simply restore m_U and m_V and set δ = 2p − 1, as described above, giving the slow manifold (see Figure 1A): M:=U,V:(aV-mVk(P))V=mUk(P)U,P=U+V

The main result of Hillen et al. (2013) is the existence of the Tumor Growth Paradox. They show that a tumor with larger death rate a_V grows quicker on the slow manifold. As a consequence, tumors with larger death rate outgrow tumors with lower death rate. The reason is that increased TC death can liberate CSC which were surrounded by TC, and it can allow CSC to replicate and produce more CSCs. As a result, the tumor becomes bigger. See Hillen et al. (2013) for the detailed argumentation using geometric singular perturbation analysis of the system.

Following Youssefpour et al. (2012), we model differentiation therapy through a simple relationship between the average level of the differentiation promoter, which we denote C_F, and the probability of CSC self-renewal, p. Unlike the model of Youssefpour et al. (2012) our model does not include a self-renewal promoter; thus, we use the relationship set forth by Youssefpour et al. (2012) but we omit the self-renewal promoting factor: (7)pt=pmin+pmax-pmin11+ψCFt where p_max is the maximum probability of self-renewal, and p_min is the minimum probability of self-renewal. The value of p_max is attained if no differentiation promoter C_F is present, while p_min is attained for C_F → ∞. Youssefpour et al. (2012) choose p_max = 1 and p_min = 0.2 in their therapy simulations. Unlike Youssefpour and coworkers, we do not model the production of differentiation promoters by tumor cells. Thus, C_F solely represents the level of differentiation promoter prescribed during differentiation therapy. To address this lack of endogenous differentiation promoters, we choose p_max = 0.505, which is equivalent to setting δ = 0.01, as was done by Hillen et al. (2013). Following Youssefpour et al. (2012) we choose p_min = 0.2. The parameter ψ models the sensitivity of the CSCs to the differentiation promoter. The dependence of p(t) on the sensitivity ψ is shown in Figure 1B. Other possible effects of differentiation therapy, such as effects on growth rates, are ignored, as they are by Youssefpour et al. (2012).

To model the average level of differentiation promoter within the spatially homogeneous tumor as a function of time, C_F(t), we assume that the tumor resides in a spherical region of tissue and that the differentiation promoter enters this area through the boundary. The ODE system [equations (4) and (5)] gives the mean tumor behavior in this spherical tissue region. Growth promoter that enters the region from the boundary will diffuse very quickly and attain a steady state distribution over this region. To compute this value of C_F(t) we solve the problem of diffusion over a sphere of radius R and average the solution over the volume of the sphere. A schematic is given in Figure 1C. We use a lower case letter to describe the radial symmetric solution c_F(r, t) of the following boundary value problem ∂cF∂t=ω∂∂r∂cF∂r+2r∂cF∂rcFR,t=CF0t.

Here ω is the effective diffusivity of the differentiation promoter. We set ω = 10⁻⁷ cm²/s throughout our simulations. Before differentiation therapy begins, C_F0(t) is zero. When differentiation therapy begins, the boundary condition on the sphere is set to C_F0(t) = 1 and the promoter diffuses into the sphere. When differentiation therapy ends, the boundary condition is simply set to zero and the promoter diffuses out of the sphere. We then set CFt=3R3∫0RcFr,tr2dr.

To model external beam fractionated radiotherapy, we apply the broadly used linear quadratic (LQ) model. The surviving fraction of cells, S(d), after a single fraction of d grays (Gy) of radiation, is given by (8)Sd=exp-αd-βd2 where α may be interpreted as lethal damage due to a single track of radiation, and β may be interpreted as lethal damage due to the misrepair of DNA damage produced by two separate tracks of radiation (Sachs et al., 2001). As a simple approximation of the radiation resistance of CSCs, we assume that they are better able to repair DNA double strand breaks such that the quadratic interaction term β is zero for CSC. We further assume that there is no interaction between DNA damage produced by separate fractions of radiation, owing to the relatively large time between fractions, typically 1 day, when compared to typical DNA repair times on the order of 1 h (O’Rourke et al., 2009).

Rather than incorporate an appropriate form of the LQ model into the system of ODEs [equations (4) and (5)], we simply apply equation (8) to the CSC volume fraction, U, and to the TC volume fraction, V, at scheduled times during the simulation, using α and β values appropriate for each cell type. For example, if a fraction is scheduled to be delivered at the beginning of the two-hundredth day of the simulation, the simulation is stopped at this time, the LQ model is applied to U and V, using their respective parameter values, and the simulation is continued at 200 days plus the fraction duration, using the surviving fractions given by equation (8) as the new initial conditions. We assume fraction durations of 10 min throughout our simulations.

We use tumor control probability (TCP) to model the probability that the cells remaining after treatment will die out. To reflect the fact that we must eliminate all CSCs for treatment success (Dingli and Michor, 2006), and the fact that TCs are doomed to die out in the absence of CSCs, we calculate TCP based on the number of CSCs remaining after treatment, using the Poisson TCP formula (see Gong et al., 2013 for Poisson TCP and other TCP models). (9)TCP=exp-NU≈exp-U43πR3ρ where N_U stands for the number of CSCs; R is the radius of the spherical region of interest in cm, as described in the section on differentiation therapy; and ρ is the density of cells in the region of interest, which we assume to be 10⁹ cells per cm³, a typical cell density for tumors (for example, see Joiner et al., 2009) The closer TCP is to one, the greater the probability that all CSCs die out and the tumor is controlled.

For all numerical simulations of the tumor model [equations (4) and (5)] we assume the mitosis rates of the CSCs and TCs are equal. That is, m_U = m_V. Further, following Youssefpour et al. (2012), we assume the apoptosis and mitosis rates of the TCs are equal: a_V = m_V. These assumptions imply that the TC populations dies out if k(0) ≤ 1, which is equation (6) for this case. When combined with our earlier assumptions regarding k(P) and with our chosen form for k(P) [equation (3)], we see that in our model TCs are always doomed to die out in the absence of CSCs. Our assumptions regarding the mitosis rates and TC apoptosis rate also simplify the form of the slow manifold to (10)M:=U,V:1-kPV=kPU,P=U+V for all simulations.

Whenever we apply radiation therapy, we assume no difference in the ability of CSCs and TCs to withstand lethal single track damage. Thus, we use the same α value for both cell types. As mentioned previously, we set β = 0 for CSCs to simulate perfect repair of two-track non-lethal damage.

All numerical simulations are carried out in Maple™, using the dsolve ODE solver employing the rfk45 numerical method. For every simulation, the initial conditions are (U₀, V₀) = (0.1, 0.1), and therapy begins on the two-hundredth day. These settings allow the tumor system to hit the slow manifold, M, before treatment begins, in each of our simulations.

To prevent negative volume fractions during numerical simulation, we introduce a simple cutoff function (11)Gx=1,x>λ0,x≤λ where λ is chosen to allow the TCP to approach one before the cutoff is imposed. The system we use for numerical simulation, incorporating the cutoff function is U°(t)=(2p-1)mk(P(t))U(t)G(U(t))V°(t)=21-pmk(P(t))U(t)G(U(t))+mk(P(t))V(t)G(V(t))-mV(t)G(V(t)) where m = m_U = m_V = a_V.

As a measure of treatment success, we calculate the TCP [equation (9)] using the value of U obtained at the end of treatment, which is defined as the latter of: (a) the completion of the final radiation fraction, and (b) the point in time when p(t) reaches 0.5, after differentiation therapy has ended. This second point (b), accounts for the effect of lingering differentiation promoter, after the promoter is no longer being applied.

To simulate ahead and neck tumor, we choose the following parameter values for the LQ model: an α/β ratio of 10 Gy and an α value of 0.35 Gy⁻¹ (Fowler, 2010). We set the mitosis rates of the TCs and CSCs to ln 2/3 day⁻¹, using a cell doubling time of 3 days as per Fowler (2010). We note that this is an estimate of cell doubling time for cells undergoing cytotoxic treatment, which tend to have shorter doubling times than untreated cells (Fowler, 2010). For the radius of the domain of interest, R, we choose 1.5 cm. All simulations of radiation therapy use fraction sizes of 2.53 Gy, delivered once per day, on weekdays only. We used one of the optimized head and neck radiation schedules recommended by Fowler as a constraint on radiotherapy, which is 25 fractions of 2.53 Gy each, for a total of 63.25 Gy delivered over 32 days (Fowler, 2010). This schedule is optimized to satisfy a late tissue constraint of 70 Gy EQD3/2 and an acute mucosal constraint of 51 Gy EQD10/2 while delivering the maximum possible BED to the tumor, given the chosen fraction size and weekday only schedule (Fowler, 2010). We indicate the position of this schedule in Figure 2 using a black plane at 63.25 Gy, which we take as our constraint on radiation therapy.

The only parameter not yet specified is the sensitivity ψ toward the differentiation promoter and we have no experimental data available. In Figure 2, we show four simulations of the tumor control probability for four different sensitivities ψ = 0.5, 2, 5, 50 which covers a wide range of possible values. The x-axis denotes the duration of the differentiation treatment and the y-axis denotes the total radiation dose. The black plane indicates the maximum tolerable radiation dose in this particular treatment. The colored plane is the tumor control probability (TCP). We see in all four Figures that the TCP is 0 near the origin and it rises sharply to values close to one as both treatment modalities are increased. In Figure 2B for example, we see that radiation alone reaches a TCP of about 60% for the maximum dose. In combination with differentiation therapy of 50 days, we observe treatment success already at total dose of 40 Gray. This effect is more pronounced for higher sensitivity parameter ψ. Notice that the curve for 0 DT days is the same in all four figures.

A good quantitative measurement for efficiency of a treatment is the TCP = 50% value. To illustrate how the treatment regimens change for a fixed TCP, we list a few treatment regimens that result in a 50% TCP in Table 2. We see that for large enough sensitivity ψ, the total radiation dose can be drastically reduced if differentiation therapy is applied.

To simulate a brain cancer, we use an average α/β ratio of 12 Gy and α value of 0.3 Gy⁻¹, as estimated by Yuan et al. (2008) for brain cancers (primary tumors as well as brain metastatic cancers). This gives a β value of 0.025 Gy⁻². For the radius of the domain of interest, R, we use 1.9 cm, which is roughly the radius of a sphere of volume 28.8 cm³, the volume of a brain metastatic cancer arising from non-small-cell lung cancer, reported in the same paper (Yuan et al., 2008). For the CSC and TC growth rates, we use ln 2/3.9 day⁻¹, where 3.9 is an estimate of the mean potential doubling time of brain metastatic cancer originating from various primary cancers, as measured by flow cytometry (Struikmans et al., 1997). All simulations involving radiation use a fraction size of 3.8 Gy, delivered once per day on weekdays only. This fraction size is listed by Yuan et al. (2008) as part of a hypofractionated stereotactic radiotherapy regimen involving 15 fractions, and it approaches the radiation tolerance for normal brain tissue. We take the total dose of 57.5 Gy listed by Yuan et al. (2008) as our constraint on radiation therapy, which appears as a black plane in Figure 3.

The results as documented in Figure 3 are very similar to those for the head and neck cancer. One difference is that without any differentiation therapy, the cancer cannot be controlled by radiation alone. At least not within the given parameter values. In Table 3 we list some TCP 50% values for this case.

To simulate the treatment of a small breast tumor, perhaps remaining after the resection of a large tumor, we choose R = 0.25 cm. The CSC and TC growth rates are set to ln 2/8.2 day⁻¹, where 8.2 is the median potential doubling time of human breast tumors measured by Rew et al. (1992) using flow cytometry. Plausible parameter values for the LQ model [equation (8)] are taken from Qi et al. (2011): α/β = 2.88 Gy and α = 0.08 Gy⁻¹. We use a fraction size of 2.26 Gy, delivered once per day on weekdays only. Our radiation constraint, indicated by a black plane in Figure 4, is 65.54 Gy, which corresponds to 29 fractions, the maximum number of fractions that satisfy the late tissue constraint of 70 Gy EQD3/2 and the acute mucosal constraint of 51 Gy EQD10/2 given in Fowler (2010).

Since the breast tumor in this example is late responding (low α/β-ratio), it is very difficult to control the cancer with radiation alone. The maximum tolerable dose of 65.54 Gy is reached much earlier than the TCP shows any growth. Using radiation in combination with differentiation therapy gives some hope that the cancer can be eradicated. Provided, however, that the CSC cells are sensitive enough to the differentiation promoter. In Table 4 we list a few TCP 50% values.