# American Institute of Mathematical Sciences

August  2019, 24(8): 4629-4663. doi: 10.3934/dcdsb.2019157

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

 1 Basque Center for Applied Mathematics, Alameda Mazarredo 14, 48009 Bilbao, Spain 2 Technische Universität Darmstadt, Schlossgartenstr. 7, 64289 Darmstadt, Germany 3 Technische Universität Kaiserslautern, Paul-Ehrlich-Str. 31, 67663 Kaiserslautern, Germany 4 German Breast Group Forschungs GmbH, 63263 Neu-Isenburg, Germany

* Corresponding author: Christina Surulescu

Received  October 2018 Revised  May 2019 Published  August 2019 Early access  June 2019

Fund Project: 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.

Citation: Julia M. Kroos, Christian Stinner, Christina Surulescu, Nico Surulescu. SDE-driven modeling of phenotypically heterogeneous tumors: The influence of cancer cell stemness. Discrete & Continuous Dynamical Systems - B, 2019, 24 (8) : 4629-4663. doi: 10.3934/dcdsb.2019157
##### References:

show all references

##### References:
Interaction of differentiated cancer cells $c$ and cancer stem cells $s$
. Initial condition: $(c,s)^t = (0.9,0.1)^t$">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$
, 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 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$
and 4, follow up to $70$ days">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
and 4, with $\tilde{\alpha} = 10$, 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
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
. 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 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)}$
TCP for the different therapy strategies based on (1), with $10^3$ simulations
. The dotted grey line represents the threshold value $\gamma$ for the normal cell density">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
. 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
, 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 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
, 4, 5 and 6">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
NTCP for the different therapy strategies based on SDEs, $10^3$ simulations
UTCP for the different treatment strategies based on $10^3$ simulations of (1) and (12)
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
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
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$
 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$
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$
 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$
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
 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
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$
 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$
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$
 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$
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$
 Parameter Description Value $\alpha _n$ sensitivity parameter (in Gy$^{-1}$) $0.0025$ $\beta _n$ sensitivity parameter (in Gy$^{-2}$) $0.003$
 [1] Baba Issa Camara, Houda Mokrani, Evans K. Afenya. Mathematical modeling of glioma therapy using oncolytic viruses. Mathematical Biosciences & Engineering, 2013, 10 (3) : 565-578. doi: 10.3934/mbe.2013.10.565 [2] Christian Engwer, Markus Knappitsch, Christina Surulescu. A multiscale model for glioma spread including cell-tissue interactions and proliferation. Mathematical Biosciences & Engineering, 2016, 13 (2) : 443-460. doi: 10.3934/mbe.2015011 [3] Dominik Wodarz. Computational modeling approaches to studying the dynamics of oncolytic viruses. Mathematical Biosciences & Engineering, 2013, 10 (3) : 939-957. doi: 10.3934/mbe.2013.10.939 [4] Harsh Vardhan Jain, Avner Friedman. Modeling prostate cancer response to continuous versus intermittent androgen ablation therapy. Discrete & Continuous Dynamical Systems - B, 2013, 18 (4) : 945-967. doi: 10.3934/dcdsb.2013.18.945 [5] Eduardo Casas, Boris Vexler, Enrique Zuazua. Sparse initial data identification for parabolic PDE and its finite element approximations. Mathematical Control & Related Fields, 2015, 5 (3) : 377-399. doi: 10.3934/mcrf.2015.5.377 [6] Chun Liu, Huan Sun. On energetic variational approaches in modeling the nematic liquid crystal flows. Discrete & Continuous Dynamical Systems, 2009, 23 (1&2) : 455-475. doi: 10.3934/dcds.2009.23.455 [7] Pedro M. Jordan, Barbara Kaltenbacher. Introduction to the special volume Mathematics of nonlinear acoustics: New approaches in analysis and modeling''. Evolution Equations & Control Theory, 2016, 5 (3) : i-ii. doi: 10.3934/eect.201603i [8] Ziheng Chen, Siqing Gan, Xiaojie Wang. Mean-square approximations of Lévy noise driven SDEs with super-linearly growing diffusion and jump coefficients. Discrete & Continuous Dynamical Systems - B, 2019, 24 (8) : 4513-4545. doi: 10.3934/dcdsb.2019154 [9] Yu Fu, Weidong Zhao, Tao Zhou. Efficient spectral sparse grid approximations for solving multi-dimensional forward backward SDEs. Discrete & Continuous Dynamical Systems - B, 2017, 22 (9) : 3439-3458. doi: 10.3934/dcdsb.2017174 [10] Lisette dePillis, Trevor Caldwell, Elizabeth Sarapata, Heather Williams. Mathematical modeling of regulatory T cell effects on renal cell carcinoma treatment. Discrete & Continuous Dynamical Systems - B, 2013, 18 (4) : 915-943. doi: 10.3934/dcdsb.2013.18.915 [11] Michael Leguèbe. Cell scale modeling of electropermeabilization by periodic pulses. Mathematical Biosciences & Engineering, 2015, 12 (3) : 537-554. doi: 10.3934/mbe.2015.12.537 [12] A. Chauviere, T. Hillen, L. Preziosi. Modeling cell movement in anisotropic and heterogeneous network tissues. Networks & Heterogeneous Media, 2007, 2 (2) : 333-357. doi: 10.3934/nhm.2007.2.333 [13] Shinji Nakaoka, Hisashi Inaba. Demographic modeling of transient amplifying cell population growth. Mathematical Biosciences & Engineering, 2014, 11 (2) : 363-384. doi: 10.3934/mbe.2014.11.363 [14] A. Chauviere, L. Preziosi, T. Hillen. Modeling the motion of a cell population in the extracellular matrix. Conference Publications, 2007, 2007 (Special) : 250-259. doi: 10.3934/proc.2007.2007.250 [15] Reihaneh Mostolizadeh, Zahra Afsharnezhad, Anna Marciniak-Czochra. Mathematical model of Chimeric Anti-gene Receptor (CAR) T cell therapy with presence of cytokine. Numerical Algebra, Control & Optimization, 2018, 8 (1) : 63-80. doi: 10.3934/naco.2018004 [16] Konstantinos Chrysafinos, Efthimios N. Karatzas. Symmetric error estimates for discontinuous Galerkin approximations for an optimal control problem associated to semilinear parabolic PDE's. Discrete & Continuous Dynamical Systems - B, 2012, 17 (5) : 1473-1506. doi: 10.3934/dcdsb.2012.17.1473 [17] Christina Surulescu, Nicolae Surulescu. Modeling and simulation of some cell dispersion problems by a nonparametric method. Mathematical Biosciences & Engineering, 2011, 8 (2) : 263-277. doi: 10.3934/mbe.2011.8.263 [18] Katarzyna A. Rejniak. A Single-Cell Approach in Modeling the Dynamics of Tumor Microregions. Mathematical Biosciences & Engineering, 2005, 2 (3) : 643-655. doi: 10.3934/mbe.2005.2.643 [19] Chien Hsun Tseng. Analytical modeling of laminated composite plates using Kirchhoff circuit and wave digital filters. Journal of Industrial & Management Optimization, 2020, 16 (5) : 2213-2252. doi: 10.3934/jimo.2019051 [20] J. A. López Molina, M. J. Rivera, E. Berjano. Electrical-thermal analytical modeling of monopolar RF thermal ablation of biological tissues: determining the circumstances under which tissue temperature reaches a steady state. Mathematical Biosciences & Engineering, 2016, 13 (2) : 281-301. doi: 10.3934/mbe.2015003

2020 Impact Factor: 1.327