A Review of Mathematical Models for Tumor Dynamics and Treatment Resistance Evolution of Solid Tumors.

Increasing knowledge of intertumor heterogeneity, intratumor heterogeneity, and cancer evolution has improved the understanding of anticancer treatment resistance. A better characterization of cancer evolution and subsequent use of this knowledge for personalized treatment would increase the chance to overcome cancer treatment resistance. Model-based approaches may help achieve this goal. In this review, we comprehensively summarized mathematical models of tumor dynamics for solid tumors and of drug resistance evolution. Models displayed by ordinary differential equations, algebraic equations, and partial differential equations for characterizing tumor burden dynamics are introduced and discussed. As for tumor resistance evolution, stochastic and deterministic models are introduced and discussed. The results may facilitate a novel model-based analysis on anticancer treatment response and the occurrence of resistance, which incorporates both tumor dynamics and resistance evolution. The opportunities of a model-based approach as discussed in this review can be of great benefit for future optimizing and personalizing anticancer treatment.

Drug resistance is one of the major reasons for patients experiencing treatment failure in the area of oncology.1 Increasing knowledge of intertumor and intratumor heterogeneity that suggests distinct cells exist in different or the same tumors as well as cancer evolution have improved the understanding of anticancer treatment resistance.2 It thereby pushes forward the necessity of precision medicine rather than a one‐size‐fits‐all approach.2 To rationalize the treatment personalization and address treatment failure, the use of modeling and simulation, which can quantitatively characterize and predict the relationships between drug exposure/pharmacokinetics (PK), drug effects/pharmacodynamics (PD), and disease progression, is widely accepted to support drug decision making.3, 4, 5, 6

Mathematical models that characterize the effects of anticancer drug treatment for solid tumors based on tumor size dynamics, which is typically quantified with measurements of tumor diameter and volume, represent one key class of models applied in cancer pharmacology. Various tumor growth modeling strategies have been previously reviewed, including agent‐based models,7 image‐based models,8 multiscale models,9 and PK/PD models.10, 11

Currently, an increasing number of studies concerning the gene sequencing of tumor biopsies in different cancer types have demonstrated the dynamics of cancer evolution.2, 12 Intratumor heterogeneity that results from cancer evolution and an evolving adaption of heterogeneous tumor to treatment are also increasingly acknowledged as key factors related to the development of resistance.2, 12 To better characterize this process and to account for tumor heterogeneity, mathematical models that consider the evolution of tumors have been proposed.13, 14, 15, 16, 17 Potentially, such evolution models in conjunction with tumor growth models could be of benefit to interpret both tumor size change and evolving tumor progression during treatment and thereby ultimately rationalize adaptive treatments for individual patients and overcome treatment resistance.

To identify the challenges and opportunities of characterizing tumor size change and resistance evolution simultaneously with a model‐based approach that can facilitate anticancer treatment optimization and personalized medicine, an overview of the current available model structures is needed. Thus, in the current review, we comprehensively summarized mathematical models for the characterization of tumor growth (inhibition) dynamics in solid tumors and the relevant clonal evolution of drug resistance by a systematic search and study of previous literature. The focus in this review lies particularly on models that are applicable for clinical data.

Literature Search

Studies that characterized tumor growth (inhibition) dynamics and clonal evolution of drug resistance with mathematical models were systematically retrieved and studied from the PubMed database to provide a comprehensive and unbiased review. In total, 274 and 85 publications were obtained, respectively, for studies of tumor dynamics and tumor resistance evolution based on established search terms. Details of the literature search are described in and . Ultimately, 61 and 25 papers, among which 13 and 2 papers were obtained from the publications’ references, which introduced corresponding original models or demonstrated application examples of certain model structures, were included, respectively, for tumor dynamics and resistance evolution modeling. Model structures, cancer types, treatments, and the ways of reporting tumor sizes were extracted from the included papers. The identified model structures were classified by equation types in later sessions and were summarized in Tables  1 and 2. Data input, knowledge requirement, study type, and objectives related to different model structures were summarized in Table  3 to provide a reference of the selection of different model structures. The information of software that was used to perform the corresponding modeling and simulation analysis was also obtained and are summarized in and .

Table 1

Modeling frameworks for characterizing tumor dynamics

Models/assumptions	Equations		Refs.
Ordinary differential equations
Basic functions describing natural tumor growth
Linear growth	dTdt=k_g	Eq. 1	18
	dTdt=k_g-d·T	Eq. 2	21
Exponential growth	dTdt=k_g·T	Eq. 3	20
	dTdt=k_g·T-d·T	Eq. 4	22, 23
Logistic growth	dTdt=k_g·T·1-TT_max	Eq. 5	24, 25
Gompertz growth	dTdt=k_g·T·lnT_maxT	Eq. 6	27, 29
	dlnTdt=a-blnT	Eq. 7	28
Combination of exponential and linear growth	dTdt=λ_0·T1+λ_0λ_1·T^20^120	Eq. 8	31
Model structures integrating tumor heterogeneity
Tumor burden(T) = Proliferative component (P) + Quiescent component (Q)	dPdt=fP-m_1·PdQdt=m_1·P	Eq. 9	25
	dPdt=fP-m_1·P+m_2·QdQdt=m_1·P-m_2·Q	Eq. 10	22, 33
Tumor burden (T) = Sensitive component (S) + Resistant component (R)	dSdt=fSdRdt=fR	Eq. 11	24
	dSdt=fS-m_1·SdRdt=fR+m_1·S	Eq. 12	23, 35
	dSdt=fS-m_1·S+m_2·RdRdt=fR+m_1·S-m_2·R	Eq. 13	34, 36
Model structures integrating tumor biology process
Angiogenesis	dTdt=fT-k·BM_0-BM_tBM_0·T=fT-k·1-BM_tBM_0·T	Eq. 14	31, 39
	dTdt=k_g·T·1-TEdEdt=k_1·T^12	Eq. 15	40
	dTdt=k_g·V·logETdEdt=k_1·T-d·T^23·E	Eq. 16	41, 42
Immune system	dTdt=fT-fI·T·hT+hfI=d_1·I	Eq. 17	44
	dTdt=fT-fI·T·hT+hfI=d_1·I_1+d_2·I_2·I_3I_3+g	Eq. 18	43
	dTdt=fT-d_1·I·T-d_2·N·T	Eq. 19	46
	dTdt=fT-d_1·I·T	Eq. 20	47
Empirical model structures describing therapeutic effect
First‐order treatment effect (“log‐kill” pattern)	dTdt=fT-k_d·T	Eq. 21	18
Exposure‐dependent treatment effect	dTdt=fT-k_d·Exposure·T	Eq. 22	22, 25
Exposure‐dependent treatment effect with resistance (TGI model)	dTdt=fT-k_d·e^-λ·t·Exposure·T	Eq. 23	20, 48, 49
Introducing a damaged cell compartment	dSdt=fS-k_d·Exposure·SdDdt=k_d·Exposure·S-d·DT=S+D	Eq. 24	24, 25
Nonlinear drug exposure–effect relationship	k_g^′=k_g·1-E_max·ExposureIC_50+Exposure	Eq. 25	21
Algebraic equations
Two‐phase model	T=e^-k_d·t+e^k_g·t-1·BASE	Eq. 26	50, 51, 55
	T=e^-k_d·t+e^k_g·t-τ-1·BASE	Eq. 27	50
	T=φ·e^-k_d·t+e^k_g·t-φ·BASE	Eq. 28	55
Model proposed by Wang et al.	T = BASE · e −A·t + B · t	Eq. 29	52, 56
An extension of Eq. 30	T = BASE · e −A·t + B · t + C · t 2	Eq. 30	53
	A=θ_1+θ_2·Dose100mg	Eq. 31	53
Simplified TGI model	T=BASE·e^k_g·t-k_dλ·1-e^-λ·t	Eq. 32	54, 57, 58, 59, 60
Partial differential equations
Proliferation‐invasion model	∂cx,t∂t=Dif·∇^2cx,t+fcx,t	Eq. 33	61, 63, 64, 69, 70
	v=2Dif·ρ	Eq. 34
	∂cx,t∂t=Dif·∇^2cx,t+fcx,t-k_d·cx,t	Eq. 35	67
	Surv=e^-α·Dose+β·Dose^2	Eq. 36	64
	∂cx,t∂t=Dif·∇^2cx,t+fcx,t-1-Surv·fcx,t	Eq. 37
	∂cx,t∂t=Dif·∇^2cx,t+fcx,t-Gx,t	Eq. 38	74

α, β, radio sensitivity parameters; A, exponential shrinkage rate constant as a result of treatment; a, b, constants; B, linear growth rate constant; BASE, baseline of tumor burden; BM0, baseline of biomarkers; BM, biomarker amount at time point t, which could be assumed to remain constant and equal to baseline in the absence of treatment; C, coefficient of quadratic growth term; , tumor cell concentration/density at location x at time t; D, damaged cells; d, death rate constant; d 1, d 2, rate constants; Dif, diffusion coefficient; E, vessel endothelial cells; Emax, maximal fraction of inhibition; f(P), f(S), f(R), f(T), growth function of proliferative cells (P), sensitive cells (S), resistant cells (R), and tumor tissue (T), respectively; G(x,t), surgical term;h, g, constants; I, I 1, I 2, I 3, components in the immune system; IC50, the drug exposure that produces 50% of Emax; k, k 1, rate constants; k d, shrinkage rate constant of tumor as a result of drug treatment; k g, growth rate/growth rate constant; k g ′, tumor growth rate constant under treatment; m 1, m 2, conversion rate constants that can be set as 0; N, normal cells; Surv, the probability of tumor cell survival; T, tumor burden; TGI, tumor growth inhibition; Tmax, carrying capacity; λ, treatment efficacy decay rate constant; λ 0, exponential growth rate; λ 1, linear growth rate; τ, delayed time of tumor regrowth; ϕ, sensitive fraction of the tumor; ρ, growth rate constant; ∇2, a Laplacian operator; , tumor proliferation function.

Table 2

Modeling frameworks for characterizing tumor resistance evolution

Models	Equations		Refs.
Stochastic models
Probability model assuming branching process	Pn+1,m|n,m=b_s·1-u·n·ΔtPn-1,m|n,m=d_s·n·ΔtPn,m+1|n,m=b_r·m·Δt+b_s·u·n·ΔtPn,m-1|n,m=d_r·m·ΔtPn,m|n,m=1-b_s+d_s·n·Δt+b_r+d_r·m·Δt	Eq. 39	76, 81, 83
Stochastic differential equation	dS=k_g·S·1-S+RT_max·dt-u·S·dt-k_d^′·S·dt+σ_1·S·dW_1-q_M·K·S·dN_1k_d^′=k_d·C_DK_D+C_D	Eq. 40	90
Deterministic models
Ordinary differential equation	dSdt=k_g-d-k_d·SdRdt=k_g-d·R+u·S	Eq. 41	91
	dSdt=k_g·1-u-d-k_d·SdRdt=k_g-d·R+k_g·u·S	Eq. 42	88
	dSdt=k_g-k_d·Dose·SdRdt=k_g·R+k_g·u·S	Eq. 43	92
	dSdt=k_g1-u_1-k_d·C_D·S+u_2·RdRdt=k_g2-u_2·R+u_1·R	Eq. 44	93
Game theory	Wi=∑p_j·Payoffij=1-r_i-d_i+1-p_i·X_i	Eq. 45	94
	W¯=∑p_i·Wi	Eq. 46
	dp_idt=p_i·Wi-W¯	Eq. 47
Integral‐differential equation	∂nx,t∂t=rx·1-θ-cx-Gρt·dx·nx,t+θ·∫01ry·My,x·ny,t·dyρt=∫01nx,tdx	Eq. 48	95, 96

n, numbers of sensitive cells; m, numbers of resistant cells; b s, birth rate of sensitive cells; d s, death rate of sensitive cells; u, mutation probability in one cell division; b s, birth rate of resistant cells; d s, death rate of resistant cells; P, probability of cell number changing from current generation to the next; S, sensitive cells, R, resistant cells, k g, , , growth rate constant; d,death rate constant; k d, shrinkage rate constant as a result of drug treatment; C D, drug concentration; K D, drug concentration that produces 50% of maximum treatment effect; dW 1, stochastic cell diffusion in a small time interval (Wiener process); dN 1, stochastic dissemination in a small time interval (Poisson process); σ 1, diffusion rate; q M, dissemination rate; K, angiogenesis; u 1, u 2, mutation rate; , fitness of type i cell; , payoff of type i cells when they meet cell type j; p , p , proportion of cells; r , cost of resistance; d , cost as a result of treatment; X , benefit for resistant cells when interacting with susceptible cells; x, y, resistance levels; , cell density with resistance level x at time t; r (x), r (y), cell division rate; c(x), treatment effect; d (x), cell death rate; , a density dependence term; θ, mutation fraction; , probability that cell y mutates to cell x.

Table 3

Data input and knowledge requirement, study type, and objectives (besides characterizing treatment effect in cancer patients) for applying different model structures

Models	Data input and knowledge requirement	Study types	Study objective
Tumor growth models
Models considering tumor heterogeneity (ODEs)
Proliferative + Quiescent (Eqs. 9 and 10)	Longitudinal TS measurement	Estimation‐based study and simulation study	Optimize treatment
	Mechanism of treatment (cell‐cycle specific or not)	Mixed‐effect model possible
	Applicable for different treatments and cancer types
	Applicable for monotherapy or combination therapy
Sensitive + Resistant (Eqs. 11–13)	Longitudinal TS measurement	Estimation‐based study and simulation study	Identify resistance type (acquired or primary) and mechanism
	Applicable for different treatments and cancer types	Mixed‐effect model possible
	Applicable for monotherapy (or combination therapy)
Model developed by Ideta et al.23 (Eq. 12)	IAS and CAS	Simulation study	Describe and predict PSA change under treatment
	Prostate cancer
Model developed by Hirata et al.36 (Eq. 13)	Longitudinal PSA measurement	Estimation‐based study	Describe and predict PSA change under treatment
	IAS and CAS	Estimate parameters for each subject	Optimize treatment
	Prostate cancer		Individualize treatment
Models incorporating biological factors (ODEs)
Angiogenesis biomarkers (Eq. 14)	Longitudinal TS measurement	Estimation‐based study	Identify clinically relevant outcome predictors and optimal time to measure the biomarkers
	Longitudinal biomarkers measurement or previously reported models for treatment–biomarker interaction	Mixed‐effect model possible
	Mechanism of treatment
	Applied mainly for antiangiogenesis treatment, monotherapy, or combination therapy
	Applicable for different cancer types
Tumor vascular support (Eqs. 15 and 16)	Longitudinal TS measurement	Estimation‐based study and simulation study	Optimize treatment
	Mechanism of treatment
	Applicable for different treatments and cancer types	Mixed‐effect model possible
	Applicable for monotherapy (or combination therapy)
Immune system (Eqs. 17 and 18)	Longitudinal TS/PSA measurement	Estimation‐based study and simulation study	Optimize treatment
	Immunotherapy	Estimate parameters for each subject
	Mechanism of treatment
	Applicable for different cancer types
Immune system (Eqs. 19 and 20)	General cancer	Simulation study	Optimize treatment
	Chemotherapy
Treatment effect model (ODEs)
First‐order treatment effect (Eq. 21)	Drug dependent	Estimation‐based study and simulation study	Characterize treatment effect
	Parameter can be different for different dose group	Mixed‐effect model possible
Exposure‐dependent treatment effect (Eqs. 22–25)	Longitudinal concentration data	Estimation‐based study and simulation study Mixed‐effect model possible	Characterize relationship between tumor‐size change and treatment exposure Optimize treatment
	Or PK model (newly developed or previously published) for simulating drug exposure or dose (as a metric of drug exposure)
	Applicable for different treatments and cancer types
	Applicable for monotherapy or combination therapy
	A compartment for damaged cells (Eq. 24) applied mainly for chemotherapy and/or radiotherapy
Algebraic equation
Two‐phase model (Eqs. 26–28)	Longitudinal TS/PSA data	Estimation‐based study	Investigate the relationship between tumor growth rate and survival
	PK information is not necessary	Estimate parameters for each subject
	Applicable for different treatments and cancer types
	Applicable for monotherapy or combination therapy
Model developed by Wang et al.52 (Eq. 29)	Longitudinal TS data	Estimation‐based study	Elucidate relationship between metrics of tumor size and survival
	PK information is not necessary	Mixed‐effect model possible
	Applicable for different treatments, monotherapy or combination therapy
	Mainly in NSCLC patients
Extension of model developed by Wang et al.52 (Eqs. 30 and 31)	Longitudinal TS data	Estimation‐based study	—
	Mainly in RCC patients treated with pazopanib	Mixed‐effect model possible
	Dose‐depended treatment effect can be incorporated
Simplified TGI model (Eq. 32)	Longitudinal TS data	Estimation‐based study	Elucidate relationship between metrics of tumor size and survival
	PK information is not necessary	Mixed‐effect model possible
	Applicable for different treatments and cancer types
	Applicable for monotherapy or combination therapy
Partial differential equation
Proliferation–invasion model (Eq. 33–38)	Two pair of T1‐Gd and T2 MRI data before treatment	Estimation‐based study and simulation study	Predict patient survival and tumor size after treatment
	Or one pair of T1‐Gd and T2 MRI data before treatment, with available parameters in previous studies	Estimate parameters for each subject	Simulate patient outcome under different treatments
	Glioblastoma		Personalize treatment
	Resection, radiotherapy, or without treatment		Investigate the application of a novel model
	DW‐MRI data, one before and two after treatment, for tumor cell number calculation	Estimate parameters for each subject (the growth rate is the net growth rate considering both tumor growth and death)	Predict tumor burden at the conclusion of treatment
	Breast cancer patients with neoadjuvant therapy
	Consider mass effect
Tumor resistance evolution models
Probability model (Eq. 39)	Parameter values (from previous studies or by estimating clinical or preclinical data)	Simulation study	Elucidate the resistance evolution of cancer
	If available, longitudinal or static ctDNA measurement can be used to estimate parameters or evaluate simulation results	Apply proposed equations in data obtained in clinical study (mainly in lung cancer, colorectal cancer, and leukemia treated with targeted treatment)	Propose equations for estimating and investigating E R, P R, total tumor cells, and treatment success rate
			Estimate the detection time
	General cancer and treatment	No mixed‐effect model applied yet	Predict treatment outcome and optimize treatment
	Single drug or multidrug resistance		Demonstrate if resistance exist at the start of treatment
Stochastic differential equation (Eq. 40)	Parameter values (from previous studies or by estimating clinical or preclinical data)	Simulation study	Connect cellular mechanisms underlying cancer drug resistance to patient survival
	Mechanism of treatment	No mixed‐effect model applied yet
	If available, longitudinal ctDNA measurement can be used to evaluate simulation results
	Mainly in melanoma patients treated with BRAF and MEK inhibitor
ODEs (Eqs. 41–44)	General cancer and treatment	Simulation study	Predict treatment outcome and optimize treatment
	Single drug or multidrug resistance	No mixed‐effect model applied yet	Propose model
Game theory (Eqs. 45–47)	Payoff matrix	Simulation study	To understand experimental results
	Combination treatment	No mixed‐effect model applied
Integral‐differential equation (Eq. 48)	General cancer and treatment	Simulation study	Describing multidrug resistance
	Single drug or multidrug resistance	No mixed‐effect model applied	Demonstrating the evolving resistance under treatment

BRAF, B‐Raf kinase; CAS, continuous androgen suppression; ctDNA, circulating tumor DNA; DW‐MRI, diffusion‐weighted MRI; E R, expected number of resistant cells; IAS, intermittent androgen suppression; MEK, mitogen‐activated protein kinase kinase; MRI, magnetic resonance images; NSCLC, non–small cell lung cancer; ODEs, ordinary differential equations; PK, pharmacokinetics; P R, probability of resistance; PSA, prostate‐specific antigen; RCC, renal cell carcinoma; T1‐Gd, gadolinium‐enhanced T1 weighted; T2, T2 weighted; TS, tumor size.

Tumor Dynamics Modeling

Ordinary differential equation

Basic growth model

A majority of the included studies applied ordinary differential equations (ODEs) to describe tumor burden change. The natural growth of a tumor without treatment is commonly characterized with several basic functions, including linear, exponential, logistic, Gompertz, and combined exponential and linear models (Table  1). The time curves of different models were simulated and are presented in Figure  1. Differential equations were solved with the RxODE package implemented in R software (version 3.4.1; R Foundation for Statistical Computing, Vienna, Austria).

Figure 1

Simulated time curves of tumor burden (T) with tumor natural growth models displayed by Eqs. 1–6 and 8. k g is the tumor growth rate / growth rate constant, d is the tumor death rate constant, Tmax is the carrying capacity, λ 0 is the exponential growth rate, and λ 1 is the linear growth rate. The baseline of tumor burden is 5. Parameter values used for the simulations are as follows: Models 1 and 2 (Eqs. 1 and 2), k  = 2; Model 2 (Eq. 2), d = 0.01; Models 3–6 (Eqs. 3–6), k g = 0.1; Model 4 (Eq. 4), d = 0.01; Models 5 and 6 (Eqs. 5 and 6), Tmax = 120; Model 7 (Eq. 8), λ 0 = 0.1, λ 1 = 2.

The linear tumor growth assumes a constant zero‐order growth rate (Eq. 1; Figure  1).10 It has been applied to describe the natural tumor growth of metastatic renal cell carcinoma18 based on the measurements of sum of longest diameters (SLD) of the target lesions in patients.

The exponential growth assumes the growth rate of a tumor is proportional to tumor burden (first‐order growth; Eq. 3; Figure  1).10, 19 It has been adopted in a widely used tumor growth inhibition (TGI) model developed by Claret et al. to describe nature tumor growth.11, 20

The linear and exponential growth models have also been expanded by introducing a first‐order shrinkage term describing natural tumor death. For example, a model with a linear growth and a first‐order shrinkage (Eq. 2) was applied to describe the natural tumor growth in patients with advanced solid malignancies based on SLD measurements.21 An exponential growth with a first‐order shrinkage (Eq. 4) was also used as part of the model structure to describe the natural growth of pediatric neuroblastoma based on tumor volume measurements.22 The same model structure was also adopted for the description of the change of prostate cancer burden reflected by the level of prostate‐specific antigen (PSA).23

When compared with the unlimited growing pattern of linear and exponential growth models, the logistic and Gompertz growth models provide a biologically realistic change of the growth rate as the tumor burden increases6 (Figure  1). The logistic growth model assumes that the growth is limited by a carrying capacity (Eq. 5)10 whereas the Gompertz model assumes the growth rate of tumor decreases over time (Eqs. 6 and 7).10, 11 Many clinical studies have applied the logistic24, 25, 26 and Gompertz models11, 27 as well as simulation studies.28, 29

Finally, a combination of exponential and linear growth models (Eq. 8) has also been introduced to describe tumor growth in patients, although it was proposed primarily for characterizing xenograft tumor dynamics.30 This combined model structure assumes that an exponential (first‐order) growth switches to a linear (zero‐order) growth after reaching a threshold (Figure  1). It was well used to describe the natural growth of vestibular schwannoma volume in patients with neurofibromatosis type 2.31 Setting the power term as 20 allows the switch between two growth patterns sharply enough.30

Tumor heterogeneity

As a result of the increasing awareness of the relevance of considering tumor heterogeneity, model structures displayed by ODEs that incorporate tumor heterogeneity and mutations have been developed for the characterization of tumor dynamics as was described in a simulation study.32 The general used model structures concerning tumor heterogeneity are shown in Table  1.

Proliferative and quiescent cells

One frequently made assumption when modeling the growth of heterogeneous tumors is to separate total tumor mass into proliferative and quiescent cells.22, 25, 33 The increase of quiescent tumor cells is assumed to result from a first‐order conversion from proliferative tumor cells instead of their own proliferation (Eq. 9). A reversed conversion can also be assumed to be present (Eq. 10). The growth of proliferative cells may follow the patterns as were introduced in the Basic growth model section. Based on these assumptions, the time courses of mean tumor diameter (MTD) in patients with low‐grade glioma25 and that of tumor volume in pediatric neuroblastoma patients were successfully described.22 A similar model structure was also used to predict the effect of different treatment regimens taking tumor cell number as a target.33 Drug treatment effect could work on both kinds of tissues,25 only on the proliferative tissue,22 or on targeted tissues depending on the types of drug.33

Sensitive and resistant cells

Another commonly made assumption is that tumors are composed of drug‐sensitive and drug‐resistant cells.24, 34 These two cell types both proliferate, but drug treatment can only decrease the amount of drug‐sensitive cells. Primary and acquired resistance can both be taken into consideration. For illustrating the acquired resistance, the resistant cells are mostly assumed to mutate from sensitive cells because of the treatment with a first‐order process23, 24, 34, 35 (Eqs. 12–13). By separating tumor mass into sensitive and resistant cells, the dynamics of low‐grade glioma measured with MTD in patients was well described with models assuming that primary resistant cells or both primary and acquired resistant cells are present in the tumor.24 In the study, the natural growth of drug‐sensitive and primary‐resistant cells were described separately without any conversion (Eq. 11). The acquired resistant cells are assumed to emerge exponentially from damaged sensitive cells as a result of treatment. Also, by assuming that resistant cells can also convert back to sensitive cells (Eq. 13), the dynamics of the PSA level in prostate cancer patients was well described, where the rate constants of cell proliferation, apoptosis, and conversion are expressed as functions of intracellular concentration of androgen receptors.34

In addition, the treatment sensitivity of both proliferative and quiescent cells can also be considered when modeling tumor growth, leading to a combination of previous introduced model structures. One example can be seen from a study that assumed proliferative and quiescent cells form a tumor and the proliferative cells could mutate from drug sensitive to drug resistant, which is biologically plausible.33

Androgen‐dependent cells and androgen‐independent cells

Studies regarding prostate cancer often consider prostate tumors consists of androgen‐dependent (AD) and androgen‐independent (AI) cells.23, 36, 37, 38 PSA levels are commonly used to represent tumor burden in this case. Two frequently reported model structures for describing the growth of prostate cancer were proposed by Ideta et al.23 and Hirata et al.36

The former model structure assumes that prostate cancer consists of AD and AI cells, and AD cells can mutate exponentially to AI cells when treatment alters the androgen level. The model structure is shown in Eq. 12. The natural proliferation and apoptosis rate constants of AD and AI cells were expressed as functions of the androgen level.23 The net growth rate of AD decreases when the androgen level decreases because of treatment, whereas that of AI cells increases. When the androgen level is normal, three cases of the net growth rate of AI cells were considered: larger than 0, equal to 0, and smaller than 0. This model was recently extended by accounting for competition between two kinds of cells and the finite carrying capacity environment.35

The latter model structure assumes that besides AD cells, reversible and irreversible AI cells exist. All types of cells are assumed to proliferate and convert to each other exponentially. It is assumed that AD cells convert to both types of AI cells during on‐treatment status and reversible AI cells convert back to AD cells during off‐treatment status. The model structure is expressed with Eq. 13. This model has been applied to adequately describe patient data.37, 38

Integration of biology process

Tumor growth models displayed by ODEs that additionally incorporate biological factors and processes have also been developed,6 such as integration of angiogenesis biomarkers and the dynamics of components in the immune system (Table  1). To apply these methods, apart from tumor burden measurements, knowledge related to the biological processes is also needed.

Angiogenesis

Concentration of vascular endothelial growth factor (VEGF) or soluble VEGF receptor may serve as biomarkers indicating the treatment effect for patients treated with angiogenesis inhibitors.11 Incorporating the dynamics of angiogenesis biomarkers in tumor growth modeling enables better understanding and prediction of tumor progression. A model structure showed as Eq. 14, where the change of biomarkers from baseline affects the tumor decay rate, was applied in two studies.31, 39 One study characterized the time course of SLD in patients with gastrointestinal stromal cancer undergoing sunitinib treatment. The natural growth of the tumor was described with the exponential model, and the model‐predicted relative change of the biomarker's amount was incorporated to affect the shrinkage of the tumor.39 The other study well characterized the dynamics of tumor volume measured in neurofibromatosis patients undergoing bevacizumab and everolimus. The natural tumor growth was described by the combined exponential and linear model (Eq. 8), and the amount of unbound VEGF was considered to affect a first‐order apoptosis of the tumor.31

Another way to account for angiogenesis effect on tumor growth is by assuming the carrying capacity of the tumor is determined by the effective tumor vascular support that is in turn affected by the tumor volume (Eqs. 15 and 16).40, 41 Logistic and Gompertz model structures were applied under this assumption. A model structure displayed by Eq. 15 was applied to well characterize the tumor growth in renal cell carcinoma (RCC) patients based on SLD measurements.40 The carrying capacity in this study was assumed to expand because of proangiogenic factors. Another similar model structure is shown by Eq. 16. Although as far as we know there is no clinical study that utilized this model framework, it has been used to perform simulations to optimize the delivery of therapeutic agents for enhancing targeted therapies for liver cancer41 and to investigate the optimization of antiangiogenic treatment.42

Immune system

Apart from angiogenesis, the effect of the immune system has also been incorporated in the tumor growth model when patients were undergoing immunotherapy.43, 44 The proposed model structure is presented in Eqs. 17 and 18, where the rate of first‐order decline of tumor burden was assumed to depend on the amount of immune component and decrease while tumor burden was increasing. This model structure was adopted to characterize the growth of prostate cancer by accounting for the dynamics of the immune system. Tumor cells were assumed to proliferate exponentially, and the amount of cytotoxic T lymphocytes affected the cell decline rate (Eq. 17).44 The applicability of this model was validated by the results of a clinical trial where PSA measurements were obtained from prostate cancer patients treated with a vaccine. Considering the effect of more than one immune component, another study developed a model structure to simulate the growth of bladder cancer undergoing immunotherapy.43 The growth of tumor cells was described with a logistic model, and the cell decline rate was set to be linearly or nonlinearly related to the amount of immune components (Eq. 18).

Another concept model structure described tumor burden dynamics by a logistic growth, a first‐order damage resulting from immune cells, and a first‐order competition with normal cells (Eq. 19).45 This model structure was recently adopted to obtain an optimal dosing regimen for cancer patients based on simulation.46 A model structure that omits the competition with normal cells (Eq. 20) was also proposed to investigate treatment optimization.47

Treatment effect

Empirical method

Tumor shrinkage resulting from drug treatment is typically quantified with an empirical drug‐induced shrinkage term as has previously been summarized.10 Commonly used equations identified from included papers are presented in Table  1. The time curves of these equations were simulated with R and are shown in Figure  2, assuming an exponential growth with the growth rate constant k g = 0.1.

Figure 2

Simulated time curves of total tumor burden (T) with tumor dynamic models incorporating treatment effect with Eqs. 21–25 and assuming an exponential growth (growth rate constant k g = 0.1). k d is the tumor shrinkage rate constant due to drug treatment, λ is the treatment efficacy decay rate constant, S is the drug sensitive cells, D represents the damaged cells, d is the death rate constant, Emax is the maximal fraction of inhibition, and IC50 is the drug exposure that produces 50% of Emax. The baseline of total tumor burden is 30. Parameter values used for the simulations are as follows: Model 1 (Eq. 21), k d = 0.4; Models 2–4 (Eqs. 22–24), k d = 0.04; Model 3 (Eq. 23), λ = 0.1; Model 4 (Eq. 24), d = 0.1; Model 5 (Eq. 25), Emax = 2, IC50 = 5. Drug exposure was simulated with Hill's equation: , where Epmax represents the maximum exposure at steady state and Ept50 represents the time when the exposure reaches half maximum value.

A log‐kill pattern is commonly used for modeling treatment effect, which assumes that the shrinkage rate of the tumor as a result of drug treatment is proportional to tumor burden.6 The simplest way to adopt this pattern is using Eq. 21, where k d is the drug‐induced tumor shrinkage rate constant. Such an equation has been used to well described the treatment effect of everolimus on metastatic RCC patients.18 The estimates of k d in that study were different between two dose groups.

The rate of drug‐induced shrinkage can also be considered to depend on drug exposure, i.e., drug concentration and area under the concentration‐time curve or drug dose. A linear drug exposure–effect relationship can be quantified using Eq. 22.22, 25 Meanwhile, drug resistance can also be taken into consideration by introducing a e ‐ term on the basis of Eq. 22 to quantify the decline of drug effect over time (Eq. 23; Figure  2). This model structure has been applied to characterize the effect of pazopanib on RCC patients.40 Setting , an exposure‐driven TGI model was developed based on SLD measurements from colorectal cancer patients receiving capecitabine and fluorouracil.20 It has then been widely applied to various cancer types and drugs as was reviewed previously.11 Two more recent studies also adopted this model structure to characterize the tumor SLD change in metastatic breast cancer patients treated with eribulin48 and in metastatic ovarian cancer patients receiving carboplatin or gemcitabine plus carboplatin,49 respectively.

In addition, a damaged cell compartment (D) has also been introduced in studies to account for the damage on cell DNA as a result of the treatment, as is displayed by Eq. 24, which can result in a delay on drug onset (Figure  2). This model structure was used in two studies that characterized the MTD change in low‐grade glioma patients treated with chemotherapy or radiotherapy.24, 25 In these two studies, the damaged cell compartment was used to characterize the treatment effects on drug‐sensitive cells24 and quiescent cells25 respectively. Part of the damaged cells eventually died, and the rest were assumed to become drug‐resistant cells24 and proliferative cells25 respectively.

Apart from the linear drug exposure–effect relationship, a nonlinear drug exposure–effect relationship can also be considered to characterize treatment effect particularly for targeted anticancer treatment.21 An Emax model is commonly used in this circumstance. An example equation is showed as Eq. 25, which was derived from a model where the studied medicine was assumed to inhibit the zero‐order growth rate of advanced solid malignancies following the nonlinear drug exposure–effect relationship.21

Considering biomarkers

When biomarkers that represent the drug‐targeting system are incorporated in the tumor dynamic models, treatment effect can be added on the dynamics of biomarkers according to corresponding mechanisms.

In the study where neurofibromatosis patients were treated with bevacizumab and everolimus, the decrease of the unbound VEGF amount because of the binding with bevacizumab was considered in the model.31 Meanwhile, the inhibition of the zero‐order production rate of total VEGF because of everolimus was described with a nonlinear exposure–effect relationship: , where IC50 is the drug exposure that produces 50% of the maximal inhibition effect. As a result of the quantity decrease of biomarkers, the shrinkage rate of tumor burden increased (Eq. 14). The delayed activation of tumor proliferation result from the continuous use of everolimus was also integrated in their model structure.31 In the study where gastrointestinal stromal cancer patients were treated with sunitinib, the effect of sunitinib was described by a nonlinear inhibition on the zero‐order production rate or first‐order decline rate of biomarkers using , where Imax is the maximal fraction of inhibition.39 The negative item in Eq. 23 was also included to quantify the treatment effect and resistance.39

In addition, the effect of angiogenesis inhibition treatment can also be incorporated by introducing a first‐order drug exposure dependent decline term (Eq. 22) on the dynamics of tumor vascular support,40, 41 when the vascular support was assumed to determine the carrying capacity of tumor (Eqs. 15 and 16).

Studies where patients were treated with immunotherapy have also considered drug interaction with the immune system. The presence of immunotherapeutic agents is frequently assumed to affect the dynamics of components in the immune system, and the amount of those components can affect the decrease rate of tumor burden (Eqs. 17 and 18).43, 44 For example, the model structure proposed to describe PSA change in prostate cancer patients treated with a vaccine assumed that the presence of the vaccine upregulated the zero‐order production rate of mature dendritic cells and therefore increased the number of cytotoxic T lymphocytes, which increased the decay of tumor tissue.44

Algebraic equation

Besides using ODEs, model structures displayed by algebraic equations have also been developed to characterize the dynamics of tumor directly as is summarized in Table  1 . 50, 51, 52, 53, 54 The simulated time curves of tumor dynamics given by these models are shown in Figure  3. Although these equations could be treated as analytical solutions of ODEs, they provided different shapes of time curves when compared with what was introduced previously.

Figure 3

Simulated time curves of tumor burden (T) with tumor dynamic models displayed by algebraic equations that describe both tumor natural growth and treatment effect (Eqs. 26–30 and 32). k g is the tumor growth rate constant, k d is the tumor shrinkage rate constant due to drug treatment, τ is the delayed time of tumor regrowth, ϕ is the sensitive fraction of the tumor, A is the exponential shrinkage rate constant due to treatment, B is the linear growth rate constant, C is the coefficient of quadratic growth term, BASE is the baseline of tumor burden, and λ is the treatment efficacy decay rate constant. Parameter values used for the simulations are as follows: Models 1–3 (Eqs. 26–28), k g = 0.1, k d = 0.4, BASE = 30; Model 2 (Eq. 27), τ = 10; Model 3 (Eq. 28), ϕ = 0.6; Models 4 and 5 (Eqs. 29 and 30), A = 0.4, B = 2, C = 0.05, BASE = 30; Model 6 (Eq. 32), k g = 0.1, k d = 0.4, λ = 0.1.

A novel two‐phase model that combines exponential tumor regrowth and regression was developed to interpret serial PSA measurements from AI prostate cancer patients50 and metastatic castration–resistant prostate carcinoma patients undergoing combination therapy.51 The corresponding model equation is shown in Eq. 26, where k g is the tumor regrowth rate constant and k d is the drug‐dependent tumor regression rate constant. The same model structure was also utilized to assess the therapeutic efficacy of bevacizumab in patients with RCC using the sum of perpendicular diameter measurements.55 On the bases of this model structure, an extra parameter τ has been introduced to account for the delayed tumor regrowth as presented in Figure  3 (Eq. 27).50 In addition, a parameter ϕ has also been introduced on the basis of Eq. 26 to differentiate the sensitive and resistant part of the tumor (Eq. 28),55 which results in a less degree of tumor shrinkage at the early phase (Figure  3). This model structure was found to be applicable when sufficient data points were available, and the estimation of growth rate constant was similar to what was obtained by the original equation (Eq. 26).

Another model structure was proposed by Wang et al. to describe the time courses of tumor SLD data of non–small cell lung cancer (NSCLC) patients from four clinical trials treated with eight treatments/placebos,11, 52 as shown by Eq. 29. A and B represent the rate constants of exponential shrinkage as a result of treatment and linear growth, respectively. The treatment effect was also characterized as a drug‐dependent manner. This model structure has been successfully applied afterward11 and was recently applied to analyze SLD measurements collected from NSCLC patients from three clinical studies to identify the obstacles to wider use of quantitative measures.56

A quadratic growth term with a coefficient C was later introduced to this model structure as is shown in Eq. 30.53 This model structure was demonstrated to have the best performance on characterizing the SLD measurements in RCC patients receiving pazopanib or placebo, and predictive patient‐specific covariates were also identified.53 Treatment effect, which is reflected by parameter A, was described in a dose‐depended manner for one group of the patients in this case (Eq. 31).

In addition, a simplified version of the previously introduced TGI model, which was displayed by an algebraic equation, was also developed (Eq. 32).54 This model structure also assumes an exponential tumor growth with growth rate (k g) while the treatment effect is described in a drug‐dependent manner with parameters account for tumor growth inhibition (k d) and drug resistance (λ). By applying this model structure, the tumor size change in metastatic colorectal cancer patient treated with bevacizumab and chemotherapy was described satisfactorily.54 This model structure has been well applied to describe tumor size change in metastatic RCC patients treated with cytokine, mammalian target of rapamycin inhibitor, and VEGF receptor inhibitors;57 in NSCLC patients undergoing treatment of carboplatin/paclitaxel combining motesanib or not;58 in NSCLC patients treated with bevacizumab and erlotinib;59 and in gastric cancer patients treated with bevacizumab and chemotherapy.60

Partial differential equation

Natural growth

Partial differential equations (PDEs), which take the change of a dependent variable in time and space into consideration, have also been adopted in the modeling of solid tumor dynamics in clinical research. One common application is known as a proliferation–invasion model or a reaction–diffusion model, which hypothesize that it is the net proliferation and invasion that contribute to the growth of cancer.61 This model formation has been typically used in studies where imaging observations of tumor, especially brain tumors, were available to describe and predict tumor expansion.8 The equation of this model structure is shown as Eq. 33 in Table  1, where the dynamics of tumor cell concentration/density at location x at time t (c(x,t)) is described.8, 61 The tumor proliferation in this model can be expressed by exponential, logistic, or Gompertz functions.8, 61 Moreover, this model mathematically regards the expansion of imaging detectable tumor edge as a “traveling wave,” and the velocity of tumor expansion is a constant that is determined by the diffusion coefficient (Dif) and growth rate constant ρ (Eq. 34).61 This linear radius/diameter expansion was confirmed in a group of grade II gliomas patients with magnetic resonance image (MRI) measurements before any oncological treatment.62

Studies applying the proliferation–invasion model to characterize tumor dynamics typically have interest in estimating the rate constants of net proliferation and invasion. An application of this model structure can be found in a study where the tumor volumes obtained from the MRI imaging were available for 70 patients with previously untreated glioblastoma.61 The tumor proliferation was described by a logistic function (Eq. 5) with a growth rate constant ρ. The ratio Dif/ρ was estimated for each patient based on MRI observations. Subsequently, setting ρ as a reported mean value and estimating Dif, the velocity of tumor radial expansion was estimated, and the survival time of patients underwent tumor resection were satisfactorily predicted by the estimated time of reaching a target radius. The same model structure was also applied on serial available MRI data from 32 glioblastoma patients before treatment.63 The net proliferation and invasion rates they quantified were significantly associated with the survival of patients. Another study characterized tumor natural growth for nine patients with glioblastoma with the same model.64 This study demonstrated that the parameter estimated based on pretreatment MRIs had high prediction accuracy for responses after treatment for these patients. Using the same model structure, the correlation between proliferation rate and hypoxic volumes based on imaging data from newly diagnosed glioblastomas patients was demonstrated.65 This model structure was also recently used to investigate the personalization of radiotherapy strategy for brain cancer patients.66

Setting , a similar model structure was also used to simulate the growth of glioblastoma based on previous reported parameters estimated from patients and estimated the survival times of patients under different parameter settings.67

Likewise, the proliferation–invasion model with logistic growth function was also successfully applied in breast cancer patients to characterize and predict their tumor burden.68 The model developed based on MRI data that were available from the early treatment phase was demonstrated to be able to predict patient response at the end of treatment.69, 70 In these studies, an apparent diffusion coefficient was estimated based on diffusion‐weighted MRI data and was then transformed to an estimate of tumor cell number, which was the dependent variable in the model. Moreover, the inhibitory effect of tumor diffusivity resulting from the stress and the deformation of surrounding tissue forced by the tumor cells were also considered in these studies,69, 70 which is called “mass effect".8 More examples of the application of the proliferation–invasion model can be found in a previous review.8

Apart from taking the diffusion coefficient as a constant, the difference between diffusion rates in gray and white matter can also be considered, such as setting Dif as two different constants for the cells in gray and white matter, respectively.71 The proportions of white and gray matter (i.e., ) have also been taken into account when computing the diffusion coefficient with the following equation: . The prediction of the model was validated with clinical imaging data from one glioma patient case.72

Recently, a threshold and a necrosis rate were also introduced into the proliferation–invasion model structure, which assumes an exponential decay will occur once the tumor cell amount exceeds the threshold.73

When using the proliferation–invasion model, the treatment effect can also be expressed by subtracting an extra term (corresponding equations are shown in Table  1). The effect of chemotherapy can be expressed with Eq. 35, where k d is the drug effect rate constant.67 For radiotherapy, a linear‐quadratic equation has been used to estimate the probability of tumor cell survival (Surv) after the administration of radiation with dose Dose (Eq. 36). The effect of radiotherapy can thus be incorporated as presented by Eq. 37.64 In addition, it is also possible to incorporate the effect of resection in the proliferation–invasion model to describe tumor growth after surgery. The resection can be simulated by setting the cell concentration in the resected region as zero at the time point of surgery.61 Subtracting a surgical term (Eq. 38) was also found to be applied to simulate the resection of tumor.74

Tumor Resistance Evolution Modeling

Tumor clonal evolution

Theoretically, three models of tumor evolution have been reported. One is a selective sweep model, which is also known as “linear” model.14, 75 It holds that during cancer initiation, mutations with fitness advantage are raised and then selectively take over the whole population sequentially.14, 75 However, because intratumor heterogeneity was identified and evidence of branching growth was found from multibiopsy and genomewide studies, a branching evolution theory where multiple subclones are considered to present and compete was developed.14, 75 Another “big bang” model of tumor evolution was observed in colorectal tumors, which suggests that advantage mutations arise and cumulate during the early phase of cancer development and the tumor then grows as a neutral single clonal.14, 75

Mathematical models that characterize tumor initiation and progression as an evolving process, including stochastic models and deterministic models, were sufficiently introduced in previous reviews.13, 14 A well‐mixed cell population is typically assumed.13 Modeling strategies that focus on describing the evolution of cancer resistance have also been discussed.15, 17 In the following sections, we will mainly give an introduction about different mathematical modeling strategies that were used to characterize cancer resistance with the tumor evolution principle.

Stochastic model

Probability model assuming the branching process

The branching process, which is also called the birth–death process, is a commonly adopted stochastic process that is used to characterize the evolving dynamics of cancer resistance.13, 15, 17 The Markov property is adopted in this model. Normally, at least two cell types, i.e., sensitive cells and resistant cells, are considered. It assumes that a tumor grows exponentially and that each sensitive cell has a certain birth rate, death rate, and a mutation probability in one cell division, and each resistant cell also has a certain birth rate and death rate. The probability of cell number change from current generation to the next could therefore be expressed with these parameters, as is shown in Eq. 39 (Table  2 ). n and m represent the numbers of sensitive cells and resistant cells, respectively. Substantially, stochastic simulation could be performed and the probability of resistance (the probability of at least one resistant cell is present; P R) and the expected number of resistant cells (E R) could be calculated with probability‐generating function.

Resistance evolution before treatment

By applying the branching process, the resistance evolution before treatment can be investigated. One study estimated the P R and E R of a cell population reached a certain size through the branching process starting with one sensitive cell.76 The fitness of the resistant cells that is relative to sensitive cells was also taken into consideration.76 The derived equations were later adopted to estimate the resistance probability of colorectal cancer prior to endothelial growth factor receptor (EGFR) antibody treatment, where the parameters were estimated based on longitudinal KRAS mutation amount measurements.77 The results indicated that the resistant mutation was highly likely to be present prior to the initiation of treatment. The same process has also been applied to investigate the evolution of drug resistance in chronic lymphocytic leukemia before treatment,78 where the growth and death rates of cancer cells were set based on patient results. In this case, besides estimating P R and E R at the time of treatment start, a time needed for the resistant population to reach a detectable level after treatment was also estimated based on which disease progression was analyzed and compared with real patient data.

Another study proposed functions for estimating the expected and median cell numbers for each resistant subclone in a metastatic lesion containing a certain number of cells with the branching process starting with a single sensitive cell.79 The predictions of relative cell numbers of resistant subclones assuming resistant cells were neutral were demonstrated to be in agreement with what was estimated based on the mutation concentrations in circulation tumor DNA (ctDNA) obtained from colorectal cancer patients treated with an EGFR blockade.79

Resistance evolution during treatment

The branching process has also been applied to simulate the evolving resistance during treatment. Regarding treatment initiation as the starting point, the dynamics of resistance evolution has been investigated with branching stochastic processes. Starting with a group of drug‐sensitive cells, Foo and Michor80 proposed functions of P R and E R during treatment depending on the length of treatment on and break time for continuous and pulsed dosing strategies. Treatment effect was incorporated by setting different birth and death rates for sensitive and resistant cells, if considering partial resistance, at on‐treatment and off‐treatment periods, respectively. They also estimated P R, E R, and variance of resistance cell number during treatment as functions of time considering with or without preexisting resistant cells.81 Treatment effect in this study was incorporated by making the birth and death rates of both sensitive and (partial) resistant cells affected by drug concentration. The treatment schedule could therefore be optimized by minimizing resistance risk or limiting the size of resistant clones. Corresponding equations were later adopted to simulate the time curve of E R and P R, and thereby to identify a relatively best treatment strategy for EGFR‐mutant NSCLC patients receiving erlotinib.82 In that study, the birth and death rates of different types of cells were obtained from in vitro experiments, and the birth rates were affected by drug concentration.82 Three cases of mutation rate change because of drug dose were also considered in the study.

Cancer progression under combination therapies has also be investigated with evolution models to predict the outcome of multiple treatment strategies in EGFR‐mutant lung cancer patients treated with two drugs.83 Tumor evolution after treatment initiation was modeled as a branching process with at least three types of cells considered: one type of sensitive cell and two types of preexisting resistant cells that are resistant to only one of the two drugs, respectively. The expected numbers of each type of cells were thereby estimated and the sum of which was the total expected cancer cell number (treatment outcome). The treatment effect was described by decreasing the birth rates of cells depending on drug concentration, and drug interaction was also taken into consideration.83

Besides separating tumor cells as being sensitive and resistant to treatment, one study also separated cells (subclones) according to resistant status and the number of accumulated drivers.84 In the stochastic branching process of tumor progression, subclones were assumed to have probabilities of raising a driver mutation and a resistant mutation during division. The accumulation of driver mutations resulted in an increase in the fitness of cells, whereas resistance was related to a fitness cost, and the fitness of nonresistant cells decreased because of treatment. By modeling the probability change of each cell type, the expected tumor size and the average frequency of resistant cells were estimated as functions of time. Subsequently, tumor detection time was calculated and used to compare the effect of prevention and postdiagnostic interventions.84

Tumor eradication

Considering that resistant mutations may die out as a result of stochastic drift during branching evolution, tumor eradication (treatment success) probability has also been investigated. One study modeled tumor progression as the following three phases: expansion with decreasing division rate until steady state, maintaining steady state, and treatment phase, starting with a single sensitive cell.85 Treatment was assumed to decrease the division rate and increase the death rate of sensitive cells. A formula of the probability of resistant cells arising but becoming extinct by the end of the treatment in each phase was then proposed, and the overall probability of treatment success was estimated as the product of the three probabilities.85

Multidrug resistance

The evolution of multidrug resistance has also been elucidated by a stochastic model where drug‐sensitive and drug‐resistant cells can divide, die (naturally and as a result of treatment), and mutate with certain probabilities.86, 87 In this model, cells accumulate one mutation that leads to resistance to one drug each time, and all mutations must be accumulated to make a cell resistant to all drugs. The treatment success probability (probability of extinction) as well as the probabilities of resistance when resistant cells generated exclusively before and during treatment were estimated, respectively. Based on the derived equations, the tumor size at which a certain percentage of patients were treated successfully were investigated under various numbers of drugs, mutation rates, and the turnover rates of cancer cells.86, 87 This model structure and the derived equation of treatment success probability were later utilized to optimize cyclic treatment scheduling.88 Moreover, taking the contribution of quiescent tumor cells into consideration by incorporating the branching process of both cycling cells and quiescent cells, the effect of quiescent cells on the treatment outcome, such as the resistance probability, of chronic myelogenous leukemia patients has also been investigated.89

Stochastic differential equation

In addition of the probability models, another stochastic modeling strategy that has been applied to characterize the development of resistance during treatment is by using stochastic differential equations. An example can be found in a study on melanoma cancer patients.90 Three types of cancer cells, including sensitive, resistant, and metastasis cells, and angiogenetic cells were considered. The dynamics of the number of drug‐sensitive cells is described by Eq. 40 (Table  2). In this differential equation, cell growth, mutation, and death were described deterministically, whereas cell diffusion and dissemination were considered as stochastic processes. Logistic growth function was used to describe the growth of cells, and the mutation from sensitive to resistant cells is described with a first‐order process. The death of sensitive cells was caused by drug treatment, and the nonlinear drug exposure–effect relationships was adopted (Eq. 40). Wiener process and Poisson process were incorporated to account for stochastic cell diffusion and dissemination, respectively. The effect of angiogenesis was also included. A drug‐induced resistance factor, which depends on drug concentrations, was integrated to increase the growth and dissemination rates. The model predictions of the progression‐free survival and number of metastasis cells were demonstrated to be, respectively, comparable with the observed progression‐free survival and ctDNA level obtained from melanoma patients treated with B‐Raf kinase and mitogen‐activated protein kinase kinase inhibitors.90

Deterministic model

ODEs

Other than stochastic models, deterministic differential equations have also been used to study the evolution toward drug resistance, especially for a population with a large size that often behaves nearly deterministically.13 The dynamics of sensitive cells and resistant cells can be modeled with ODEs similar to what were introduced in the “Tumor Heterogeneity” section, but the transition from resistant to sensitive cells is often neglected. The model structures that have been identified are shown in Table  2.

One model of resistance evolution displayed by ODEs is shown as Eq. 41, where drug resistance is considered to raise due‐to‐point mutations.91 When considering multiple drug resistance, multiresistant cells were assumed to only be mutated from single‐resistant cells. Starting with a certain number of sensitive cells, the resistance amount by the time of treatment initiation and during treatment was estimated under different conditions. The authors demonstrated that the simpler ODE model provided comparable results to previous models that were derived from more complicated stochastic models.91 Another example can be seen in Eq. 42. This model was used to investigate the preferable treatment by controlling the total amount of fully resistant mutants, which can be acquired from sensitive cells and single‐resistant cells.88 In addition, a model with treatment effect being proportional to drug dose has also been used to model evolving tumor resistance (Eq. 43).92 Multiresistant cells were also considered and were assumed to mutate only from single‐resistant cells. Based on this model structure, the survival of patients undergoing different treatment strategies, such as the strategy of minimizing the total cell population or minimizing the multiresistant population, was investigated.92 Another model structure of resistance evolution that includes the transition from resistant to sensitive cells (Eq. 44) has also been adopted to investigate the optimization of treatment.93

Game theory

Evolutionary game theory has also been used to investigate the evolution of cancer resistance, especially under combination therapy.94 It assumes the fitness of one type of cell, which can be understood as the growth rate, changes when the cells interact with different types of other cells. This can be expressed with a payoff matrix, and the final fitness of one type of cell is their expected payoff of this “game".13 An example was found from a study where a well‐mixed population and a deterministic dynamic of the evolving process were considered.94 The evolutionary game theory was adopted to investigate and understand the evolving resistance for small cell lung cancer patients under a combination of chemotherapy and tumor suppressor p53 vaccine treatment.94 Three cell populations, including cells that are sensitive to both treatments and cells that are resistant to one of the treatments but sensitive to the other, were considered to constitute the total tumor population. As presented in Table  2, the fitness of type i cell can be expressed as a sum of the product of the payoff of type i cell interacting with type j cell and the proportion of type j cell (Eq. 45), where a cost of resistance and a cost as a result of treatment was considered.94 In addition, to account for the influence of cell interaction on cell sensitivity and fitness, an extra benefit for resistant cells when interacting with susceptible cells under treatment was also introduced (Eq. 45).94 The average fitness was expressed with Eq. 46, where p is the proportion of each type of cells. The dynamics of each cell type under sequencing treatment was described using a replicator equation (Eq. 47), and the time curve of the proportion and fitness of each cell type are two main outcomes of the simulations in this study.

Integral‐differential equation

An integral‐differential equation, where the states of cancer resistance are described in a continuous way ranging from complete sensitivity to complete resistance, has also been used to characterize the evolution of cancer resistance.95, 96 A model structure shown as Eq. 48 has been used to describe the dynamics of cancer cell density with resistance level x at time t,95, 96 where cell division, cell death, treatment effect, and cell mutation were all incorporated (Table  2). Simulations were performed in these studies to illustrate the evolution of resistant level during treatment, but it has not yet been applied in clinical studies.

Model Selection

Applying different model structures to characterize tumor dynamics and tumor resistance evolution may achieve different objectives and require different data input and knowledge (Table  3). The target cancer type and treatment option may also influence the selection of model structure (Table  3).

As for the tumor dynamics models displayed by ODEs and algebraic equations, most models are applicable to describe tumor size change in patients with various kinds of solid tumors and under different kinds of treatment (monotherapy or combination therapy). However, the models specifically developed for prostate cancer are mainly suitable to describe PSA level change, and the models incorporating angiogenesis biomarkers or immune components are normally considered when patients are treated with antiangiogenesis treatment or immunotherapy, respectively.

Longitudinal tumor size data, such as the SLD of target lesions, MTD, or tumor volumes, or PSA measurements are required to estimate model parameters. A mixed‐effect modeling approach has been applied to most model structures that are displayed by ODEs and algebraic equations to account for interindividual variability, whereas the parameters of other structures, such as the two‐phase model, were normally estimated for each subject separately. In the former case, each subject in a group is normally required to contribute at least one measurement before treatment and one thereafter. More data points are preferred to enable the better estimate of all parameters. However, the latter method may require each subject to contribute enough data points to enable parameter estimates. In addition, if a study aims at developing a model incorporating biomarkers, longitudinal biomarker observations or previously reported models for treatment‐biomarkers interaction are required. If no specific biological process is considered, the selection of model structures can also depend on the model fit to the data as long as the model is physiologically or biologically plausible.

Among the functions of the natural tumor growth (Eqs. 1–8), which are always part of the tumor dynamics models, the exponential growth model has been the most frequently selected in clinical studies. The logistic growth model was normally satisfactorily applied when the maximum tumor capacity was fixed. The selection of the basic functions could also depend on the model fit to the data. More than one available pretreatment tumor size measurement would be helpful to find the best fit natural growth model and would enable a more accurate estimate of the tumor natural growth rate.

The treatment effect can be characterized in a drug‐dependent manner or exposure–dependent manner. If a study does not focus on investigating the exposure–effect relationship, using a model with drug‐dependent tumor shrinkage will be enough and drug‐exposure information is not required. For studies aiming at characterizing the relationship between drug exposure/dose and tumor response and/or optimizing treatment regimens for patients based on simulations, the exposure‐dependent (or dose‐dependent) treatment effect structure should be applied. To estimate drug exposure, longitudinal concentration data for PK model development or a previously reported PK model are  needed. In addition, the previous knowledge of the treatment mechanism may also be required to appropriately characterize the treatment effect, especially when applying models considering biological factors.

The proliferation–invasion model that is displayed by PDE has mainly been applied to investigate glioblastoma or breast cancer based on available MRI measurements. The required parameters can be estimated for each patient separately based on two sets of pretreatment MRI data or one before treatment and one thereafter. Simulations can then be performed to predict patient outcome with the model or with the velocity function of tumor radius expansion (Eq. 34). The mixed‐effect modeling approach has not been found to be applied in these studies yet.

The model structures of tumor resistance evolution have been mainly applied to perform simulations to understand evolving resistance and optimize the treatment. The equations derived from the branching process can be applied to answer clinical questions. Available longitudinal or static ctDNA measurements can be utilized to determine the parameter values and to evaluate the simulation results. Although no mixed‐effect modeling approach has been applied in these studies yet, the model structures displayed by ODEs, which can provide comparable results to stochastic models, are considered to be potentially able to account for interindividual variability.

Discussion

Overcoming treatment resistance with a better understanding of cancer evolution and personalizing treatment brings opportunities to treat cancer as a chronic disease and has been increasingly studied in the oncology field. Model‐based approaches incorporating tumor growth and resistance evolution may help achieve this goal. By applying mathematical models, prior knowledge derived from clinical trials and routine patients care can be utilized to quantitatively understand drug PK profiles, the drug–response relationship, and evolving resistance in cancer patients. These profiles can be predicted accordingly for future patients, which could be beneficial for identifying optimized therapeutic regimens. Furthermore, by accounting for interindividual variability with a mixed‐effect modeling approach, treatment individualization can also be designed and guided rationally.97

In the current review, feasible model structures that have been used to describe and predict tumor dynamics and resistance evolution during treatment for patients with solid tumors are discussed. Models concerning tumor evolution in leukemia were included because they provide reference value for solid tumors. Apart from what has been introduced, more extensive models have also been found in the literature search, such as agent‐based models and the cellular automata approach. The agent‐based models often include components from two or more spatial or temporal scales, ranging from molecular to tissue,7 and the cellular automata approach adopts a discrete dynamical system of time and space.9 Although tumor growth can be simulated in silico realistically with these approaches, because they require infeasible information input (e.g., cell location, nutrition distribution, and/or oxygen amount) from clinical patients, they were excluded from the current review. Studies applying the proliferation–invasion model, which are expressed with PDE, were not excluded, although tumor cell location is also one of the variables. It is because two main parameters in this model structure, the diffusion coefficient Dif and growth rate constant ρ, can be estimated directly based on MRI results obtained from patients, and the velocity of tumor radius expansion can then be estimated and utilized for prediction.

Models displayed by ODEs, algebraic equations, and PDE are commonly reported for the modeling of tumor size change and, in the case of prostate cancer, PSA amount change. Five main basic natural tumor growth model structures were frequently reported. The diversity in model selection can be explained by the difficulties of assessing real long‐term natural tumor growth pattern in patients.11 Although setting the maximum boundaries of tumor growth is more biologically plausible, the models without such limits, especially the exponential growth models, have also been used extensively. The concept of linear growth is also reflected in the studies that applied the proliferation–invasion model, as the expansion of tumor radius has a constant velocity under such a model, and this concept has been used to predict tumor radius.63, 64, 65

For characterizing treatment effect, empirical methods are relatively simple to apply for describing the effect of various kinds of drugs and are therefore more generally applicable. The shrinkage rate of tumor burden caused by treatment can be described to be proportional to drug exposure/dose or by utilizing drug‐dependent parameters, although the latter method does not allow differentiation among different dosing regimens. In addition, when the dynamics of biomarkers are available and are incorporated in the tumor dynamics models, the treatment effect on the production of biomarkers can be integrated according to drug mechanism.31, 39 Furthermore, the regrowth of a tumor during treatment can be considered in several ways. Studies applying algebraic equations generally characterize the decline and regrowth of a tumor by a single equation. For studies that used ODEs, tumor regrowth was mainly characterized by separating the tumor in two parts consisting of drug‐sensitive cells and drug‐resistant cells or by adding the e ‐ term.

The resistance evolution of cancer has been mainly characterized by stochastic models within which the branching process is reported most frequently. However, in studies applying the branching process, the focus was mainly on the expected outcome of tumor evolution, such as the P R and E R. Therefore, relatively simpler deterministic models are considered to be good alternative choices. It has already been demonstrated that ODE models can provide comparable results to those that are derived from stochastic models.91 Given that the goal is to characterize evolving tumor resistance based on clinical data, applying deterministic models might be more suitable given clinical available data generally represents the apparent response of each patient.

Among the studies included in this review, the detailed data of resistance evolution have not yet been incorporated in tumor size–based modeling of anticancer treatment response. However, genetic biomarkers that represent tumor heterogeneity and resistance evolution become increasingly available as a result of novel technologies. For example, in a clinical setting, a feasible genetic biomarker that is also correlated with tumor burden has been identified as ctDNA.98 Three of the included studies have already utilized the available ctDNA data to support the estimation of parameters in the tumor evolution model or to evaluate the model simulation results.77, 79, 90 It has also been demonstrated that the mutation in ctDNA, which represents treatment resistance, is detectable before disease progression,99 suggesting the predictive value of ctDNA to the development of drug resistance. By applying longitudinal monitoring of ctDNA, an adaptive treatment for individual patients may be achieved by selecting drugs that target emerging actionable mutations.98 Therefore, it is feasible to obtain the information of evolving cancer resistance and, to increase the chance to overcome treatment resistance, it would be helpful if such information could be incorporated in future model‐based studies.

Based on what was learned from previous reported studies, as is introduced in this review, model structures displayed by ODEs are considered to be feasible for the characterization of both tumor size change and resistance evolution in cancer patients. A mathematical model can be developed based on the input data of tumor size, mutation load of ctDNA, and treatment information over time. The emergence and dynamics of mutations in ctDNA can provide insight of the occurrence, growth, decay, and mutation for different tumor subclones. External data sets, if available, can be used to further evaluate the developed model structure. Subsequently, the effect of sequential treatment regimens with different dose levels or starting times of therapies can be explored with simulation and thereby to facilitate the identification of an optimal regimen. Moreover, because the parameter values can be estimated for each individual and the variability of which can be partially explained by patient characteristics, the treatment personalization can also be rationally guided based on the modeling and simulation results. These will be the ultimate output of the model‐based study.

However, challenges remain beyond what is already stated. First, in terms of data collection, previous knowledge of the mutations that represent resistant subclones is required. Second, if sequencing data of the subclones (ctDNA) over time are available, efforts need to be made to handle the vast amount of genetic data in a quantitative manner in relation to tumor size dynamics. Third, the optimal method on how to predict a newly acquired mutation that has not yet occurred in the data needs to be further explored. Finally, because in‐depth knowledge is required from multiple aspects of tumor and clone dynamics as well as complex modeling and simulation, a multidisciplinary collaboration is essential to enable the achievement of the ultimate goal of optimizing and personalizing anticancer treatment.

In conclusion, based on a systematic search of studies from the literature, mathematical models that have been used to describe and predict tumor size change, drug effect, and resistance evolution based on clinically available data were introduced in this review. The results may facilitate the model‐based anticancer treatment response analysis that accounts for both tumor growth inhibition and resistance evolution, although important challenges still need to be overcome. An ultimate model structure handling all of these aspects would be of great benefit for optimizing and personalizing anticancer treatment.

Funding

No funding was received for this work.

Conflict of Interest

The authors declared no competing interests for this work.

Author Contributions

A.Y., D.J.A.R.M., J.G.C.v.H., J.J.S., and H.‐J.G. wrote the manuscript. A.Y. performed the research. D.J.A.R.M. and H.‐J.G. designed the research.

Figure S1. Diagram of literature scanning for (a) tumor dynamics and (b) tumor resistance evolution.

Click here for additional data file.

Table S1. Software that was applied in studies concerning tumor dynamics (TD) and tumor evolution (TE) to perform parameter estimation and data simulation.

Click here for additional data file.

Supplementary Material S1. Literature searching method and software summarization.

Click here for additional data file.