Article Contents
Article Contents

# SDE-driven modeling of phenotypically heterogeneous tumors: The influence of cancer cell stemness

• * Corresponding author: Christina Surulescu

This work is dedicated to Peter Kloeden on the occasion of his 70th birthday

• We deduce cell population models describing the evolution of a tumor (possibly interacting with its environment of healthy cells) with the aid of differential equations. Thereby, different subpopulations of cancer cells allow accounting for the tumor heterogeneity. In our settings these include cancer stem cells known to be less sensitive to treatment and differentiated cancer cells having a higher sensitivity towards chemo- and radiotherapy. Our approach relies on stochastic differential equations in order to account for randomness in the system, arising e.g., due to the therapy-induced decreasing number of clonogens, which renders a pure deterministic model arguable. The equations are deduced relying on transition probabilities characterizing innovations of the two cancer cell subpopulations, and similarly extended to also account for the evolution of normal tissue. Several therapy approaches are introduced and compared by way of tumor control probability (TCP) and uncomplicated tumor control probability (UTCP). A PDE approach allows to assess the evolution of tumor and normal tissue with respect to time and to cell population densities which can vary continuously in a given set of states. Analytical approximations of solutions to the obtained PDE system are provided as well.

Mathematics Subject Classification: Primary: 92C50, 60H10, 35Q92; Secondary: 34F05, 35Q84.

 Citation:

• Figure 1.  Interaction of differentiated cancer cells $c$ and cancer stem cells $s$

Figure 2.  Several trajectories of the solution components of system (1) and averages (black lines) of $10^3$ trajectories solved with the Euler-Maruyama algorithm and using parameter values as in Table 2. Initial condition: $(c,s)^t = (0.9,0.1)^t$

Figure 3.  Different trajectories of the SDE system (1) and their averages (thick black lines) for $10^3$ trajectories when applying chemotherapy. Parameter values as in Table 2, dose-dependent probability $a_3$ of symmetric splitting of a cancer stem cell into differentiated cells as given in (3). Upper row: $\tilde \alpha = 0.5$, lower row: $\tilde \alpha = 10$

Figure 4.  Different trajectories of the solutions components of system (1) with the application of radiation therapy and their averages (thick black lines) over $10^3$ trajectories. Parameter values as in Tables 2 and 4, follow up to $70$ days

Figure 5.  Trajectories of the SDE model (1) with the simultaneous application of radio- and chemotherapy (differentiation promoter) and their averages (thick black lines) over $10^3$ simulated trajectories. Parameter values as in Tables 2 and 4, with $\tilde{\alpha} = 10$, follow up to $70$ days

Figure 6.  Different trajectories of the SDE model (1) with the alternating application of chemotherapy and radiation therapy and their mean values (thick black lines) for $10^3$ simulated trajectories. Parameter values as in Tables 2 and 4 with $\tilde{\alpha} = 10$, follow up to $70$ days

Figure 7.  Trajectories of (1) and their averages (thick black lines) over $10^3$ simulated trajectories with a time-discrete radiotherapy solved with the Euler-Maruyama scheme, parameter values as in Table 2. Subfigures 7(a) and 7(b) show the results of the simulation for the two subpopulations of cancer cells and 7(c) shows the processes $L_t^{(s)}$, $L_t^{(c)}$

Figure 8.  TCP for the different therapy strategies based on (1), with $10^3$ simulations

Figure 9.  Trajectories of (12) for normal tissue density and of CSCs and DCs as solutions to (1). Their averages over $10^3$ trajectories are shown as thick black lines. Parameter values are given in Table 5. The dotted grey line represents the threshold value $\gamma$ for the normal cell density

Figure 10.  Trajectories of (12) for normal tissue density and of CSCs and DCs as solutions to (1) for the differentiation therapy approach in Section 3.1, with $\widetilde \alpha = 10$. Their averages over $10^3$ trajectories are shown as thick black lines. Parameter values are given in Table 5. The dotted grey line represents the threshold value $\gamma$ for the normal cell density

Figure 11.  Trajectories of the SDE (12) and of system (1) and their averages (black lines) over $10^3$ trajectories. Parameters as in Tables 2, 4, 5 and 6. Plots 11(a)-11(c): effect of radiation therapy only. Plots 11(d)-11(f): concurrent chemo and radiotherapy. Plots 11(g)-11(i): alternating chemo and radiotherapy

Figure 12.  Trajectories of the SDE (12) and of processes $c(t)$, $s(t)$ obtained by applying time discrete therapy: radiotherapy alone (plots 12(a)-12(c)); in a concurrent (plots 12(d)- 12(f)) or an alternating way (plots 12(g)-12(i)). Averages over $10^3$ trajectories are shown in thick black lines. Parameters as in Tables 2, 4, 5 and 6

Figure 13.  NTCP for the different therapy strategies based on SDEs, $10^3$ simulations

Figure 14.  UTCP for the different treatment strategies based on $10^3$ simulations of (1) and (12)

Figure 15.  Solution components of (27) and (29) at various time moments, no therapy applied. Left column: expected total tumor burden, middle: expected density of cancer stem cells, right: expected deviations of normal cells from critical threshold for which $\gamma = 0.7$. Plots 15(a)-15(c): initial state, 15(d)-15(f): 10 days, 15(g)-15(i): 30 days, 15(j)-15(l): 70 days

Figure 16.  Solution of (29) at different times, with no therapy applied. Shown are expected deviations of normal cells from the critical threshold for which $\gamma = 0.85$. Plot 16(a): 10 days, 16(b): 30 days, 16(c): 70 days

Table 1.  Transition probabilities for ${\bf X} = (c,s)^t$, according to the scheme in Figure 1

 Changes in ${\bf X}$ Probability $\Delta X^{(1)} = (1,0)^t$ $p_1=a_2(c) b_ss \Delta t$ $\Delta X^{(2)} = (-1,0)^t$ $p_2=d_c c \Delta t$ $\Delta X^{(3)} = (0,1)^t$ $p_3= a_1(c)b_s s \Delta t$ $\Delta X^{(4)} = (0,-1)^t$ $p_4= d_s s\Delta t$ $\Delta X^{(5)} = (2,-1)^t$ $p_5= a_3(c)b_s s \Delta t$ $\Delta X^{(6)} = (-1,1)^t$ $p_6=d_cca_1(c)b_ss\Delta t$ $\Delta X^{(7)} = (2,0)^t$ $p_7=0$ $\Delta X^{(8)} = (-1,-1)^t$ $p_8=d_cc d_ss\Delta t$ $\Delta X^{(9)} = (1,-1)^t$ $p_9=b_cc d_ss\Delta t=0$ $\Delta X^{(10)} = (1,1)^t$ $p_{10}=b_cc a_1(c)b_ss\Delta t=0$ $\Delta X^{(11)} = (3,-1)^t$ $p_{11}=0$ $\Delta X^{(12)} = (0,0)^t$ $p_{12}=1-\sum_{i=1}^{11} p_i$

Table 2.  Summary of the model parameters for cancer cells based on [61]

 Parameter $b_s$ $d_c$ $d_s$ $\widetilde{a_1}$ $\widetilde{a_2}$ $\widetilde{a_3}$ Value $\log 2/55$ $\log 2/60$ $\log 2/200$ $0.5$ $0.05$ $0.1$

Table 3.  Persistence times for the whole tumor $c+s$, averaged over 1000 trajectories, for different compositions $(c,s)$ of the initial tumor

 Initial values $(c_0,s_0)^t$ Persistence time (in days) $(0.01, 0.01)$ 16.313 $(0.1,0.1)$ 56.052 $(0.99,0.01)$ 69.463 $(0.8, 0.2)$ 69.790 $(0.95,0.05)$ 69.791 $(0.7, 0.3)$ 69.799 $(0.4,0.6)$ 69.804 $(0.3,0.7)$ 69.807 $(0.6, 0.4)$ 69.871 $(0.9,0.1)$ 69.892 $(0.5,0.5)$ 69.909 $(1, 1)$ 70

Table 4.  Model parameters for radiation therapy based on [53,66], where no difference was made between various phenotypes of tumor cells. The parameter $\beta _s$ was chosen here slightly larger than $\beta _c$ to account for the reduced sensitivity of CSCs towards therapy

 Parameter $\alpha\ (\text{in Gy}^{-1})$ $\beta_c \ (\text{in Gy}^{-2})$ $\beta_s \ (\text{in Gy}^{-2})$ $d$ (in Gy) $\nu$ (in days) Value $0.0906$ $0.006$ $0.00605$ $2$ $25$

Table 5.  Model parameters for normal tissue (without therapy)

 Parameter Description Value $\mu_n$ constant maximum birth rate (in $\text{day}^{-1}$) $0.01$ $\delta_n$ interaction factor $0.002$ $\gamma$ threshold of functionality $0.7$ $n_0$ initial density of normal tissue $1$ $M$ actual carrying capacity (normalized) $1$

Table 6.  Model parameters for radiotherapy affecting normal tissue, based on [21]

 Parameter Description Value $\alpha _n$ sensitivity parameter (in Gy$^{-1}$) $0.0025$ $\beta _n$ sensitivity parameter (in Gy$^{-2}$) $0.003$

Figures(16)

Tables(6)