Elsevier

Applied Mathematical Modelling

Modeling continuous levels of resistance to multidrug therapy in cancer

Abstract

Multidrug resistance consists of a series of genetic and epigenetic alternations that involve multifactorial and complex processes, which are a challenge to successful cancer treatments. Accompanied by advances in biotechnology and high-dimensional data analysis techniques that are bringing in new opportunities in modeling biological systems with continuous phenotypic structured models, we investigate multidrug resistance by studying a cancer cell population model that considers a multi-dimensional continuous resistance trait to multiple drugs. We compare our continuous resistance trait model with classical models that assume a discrete resistance state and classify the cases when the continuum and discrete models yield different dynamical patterns in the emerging heterogeneity in response to drugs. We also compute the maximal fitness resistance trait for various continuum models and study the effect of epimutations. Finally, we demonstrate how our approach can be used to study tumor growth with respect to the turnover rate and the proliferating fraction. We show that a continuous resistance level model may result in different dynamics compared with the predictions of other discrete models.

Introduction

The biological mechanisms responsible for the emergence of drug resistance and its propagation often involve a multifactorial and complex process of genetic and epigenetic alternations [1], [2], [3], that arise through a series of genetic and non-genetic changes [4], [5], [6], [7]. Such changes can be due to drug administration (drug induced resistance), or they can emerge independent of therapy due to intrinsic mechanisms [8], [9]. Cancer cells may develop simultaneous resistance to structurally and mechanistically unrelated drugs, leading to multidrug resistance (MDR) [1], [2], [10]. The complex dynamical nature of MDR is one of the most challenging obstacles to successful treatment.

The complexity of the mechanisms underlying drug resistance has led to multiple studies that integrate mathematical modeling with experimental and clinical data. Such models aim at providing quantitative tools for testing therapies that circumvent or at least delay the unfortunate consequences of drug resistance. Examples include the models of Goldie and Coldman [11], [12], [13] that are based on resistance due to point mutations. These works were followed by many studies in which stochastic models (including branching processes and multiple mutations) were used to study MDR and optimal control of drug scheduling [14], [15], [16], [17]. Alternative approaches include continuum deterministic models using ordinary differential equations for modeling kinetic resistance [18] and point mutations [19], and partial differential equations, where spatial heterogeneity and vascularization can be readily incorporated [20], [21], [22]. For additional approaches see [17], [23], [24], [25], [26], [27], [28].

In addition to these modeling approaches, new tools for collecting phenotypical data are generating new opportunities for mathematical modeling of biological systems. Current technology allows cytometry data to be collected up to O(50) dimensions, Methylation profiles in the scale of O(1000), and gene-expression profile in the scale of O(10000) [29], [30], [31], [32], [33], [34]. In particular, recent advances in single cell RNA sequencing technologies enabled a new high-dimensional definition of cell states, that is on the order of 20,000 protein encoding genes that compose the transcriptome [32], [35]. The high-dimensionality of the data makes it practically impossible to consider a meaningful model on the original space in which the data is collected. Thus, various dimensionality reduction techniques, such as, principal component analysis [36], [37], t-distributed stochastic neighbor embedding [30], [38], [39], diffusion maps [40], [41], and machine learning techniques [42], [43], have been employed to reduce the dimensionality and to identify the critical directions. In contrast to classical biology and modeling approaches, where cell types are classified into discrete states and differentiation is considered as a stepwise process of a binary branching decision, the new technologies and data analysis tools enabled considering cell differentiation as a continuous process that can be mapped into a continuum of cellular and molecular phenotypes [31], [33], [40]. In other words, the high-dimensional configuration space is mapped into a continuous trait in a lower-dimensional space. Fig. 1 shows two examples of high-dimensional cell data mapped into a continuous trait in a lower dimensional space using stochastic neighbor embedding (viSNE) [38] and diffusion mapping [41]. This reveals the continuous phenotypic trait space where resistance can be locally characterized. For instance, the left figure shows that relapsed leukemia cells are associated with high expression of CD34, and the ALDH1 in the right figure is related to cancer stem cells in mammary gland and breast cancer [44]. Such tools open the door to mathematical models that assume a continuous trait space [45], [46].

Among continuous phenotypic structured models, recent studies in Lorz et al. [[47], [48], Greene et al. [49], Cho and Levy [50], Lorenzi et al. [51], Cho and Levy [52] consider a continuous trait variable that represents the level of cytotoxic drug resistance. This framework allows to explicitly model the heterogeneous response to drugs and effectively study the selection dynamics under microenvironmental constraints and chemotherapy. The asymptotic distributions on the resistance trait space are obtained in Lorz et al. [47]. The following works in Greene et al. [49], Lorenzi et al. [51] extend it to include mutations and epimutations. The distribution of resistance levels can be then translated to therapeutic recommendations. The effectiveness of a combination of cytotoxic and cytostatic drugs when cytotoxic resistance emerge is studied in Lorz et al. [47]. An optimal combination therapy to eliminate the most resistant clones is proposed in Cho and Levy [52]. Moreover, the authors extend the framework that was restricted to radially symmetric tumors with a fixed boundary [48] to an asymmetric tumor growth model with a moving boundary. However, this framework is limited to a single trait variable corresponding to a cytotoxic drug.

In this paper, we extend the framework of Cho and Levy [52] to a multi-dimensional resistance trait. We compare our approach that allows for a continuous drug response to more traditional approaches that assume a discrete response to drugs. The paper is organized as follows. In Section 2, we introduce a mathematical model for MDR assuming continuous trait variables. We parameterize our model as an extension of a discrete resistance state model in Section 2.1 and compute the maximal fitness trait of resistance in Section 2.2 for different types of continuum models. This allows us to characterize the cases when the solutions of the continuous models are qualitatively different than the corresponding discrete models. Section 2.3 presents simulation results for the different cases of cytotoxic and cytostatic drugs studied in 2.2. The impact of mutations and epimutations is studied in Section 2.4. In Section 3 we simulate tumor growth and resistance dynamics subject to MDR on different types of tumors characterized by turnover rates and the proliferating ratios. Our simulations correspond to the discrete MDR models studied by Komarova and Wodarz [53] and Gardner [54]. We observe that a combination therapy with multiple cytotoxic drugs is also effective in high turnover tumors using relatively high dosages. Increasing the dosage in a low turnover tumor is effective only for certain drug uptake functions. In addition, the drug response function plays a key role in determining the tumor growth dynamics when combination therapy is administered using cell-cycle nonspecific cytotoxic drugs, such as Cyclophosphamide and Doxorubicin. Conclusions and future directions are discussed in Section 4.

Section snippets

Models of multidrug resistance

Let us consider a cancer growth model under multidrug therapy that depends on an M-dimensional phenotype variable θ = ( θ 1 , … , θ M ) ∈ Γ ≐ Π i = 1 M Γ i . The phenotype variable in the ith direction θ i ∈ Γ i = [ 0 , 1 ] characterizes the resistance level to the i-th drug or the i-th drug mechanism, where θ i = 0 and θ i = 1 represents the fully-sensitive cells and fully-resistant cells to drug i, respectively. The value of θi can be obtained by normalizing the expression level of a gene or a gene cluster that is linked to the

Simulating tumor growth under multidrug therapy

In this section, we demonstrate how our continuous phenotype structured modeling framework can be used to study MDR. The impact of the tumor's turnover rate and the proliferating fraction of cancer cells have been studied within a discrete phenotype framework by Komarova and Wodarz [53] and by Gardner [54]. Here, we compare the results obtained with our approach with the conclusions of Komarova and Wodarz [53] and Gardner [54]. In particular, we focus on the relapse after therapy that

Conclusion

In this paper, we propose a mathematical model for multidrug resistance, assuming a continuous resistance phenotype space. The multidrug resistance trait variable represents the level of resistance to various drugs including cell-cycle specific and nonspecific cytotoxic drugs, as well as cytostatic drugs. We classify the proliferation and drug uptake functions and identify the cases where the continuum model results in an intermediate maximal fitness resistance, i.e., the cases in which the

Acknowledgments

The work of DL was supported in part by the National Science Foundation under Grant Number DMS-1713109 and by the Jayne Koskinas Ted Giovanis Foundation.

Cited by (9)

View full text

© 2018 Elsevier Inc. All rights reserved.