Computational Model of Heterogeneity in Melanoma: Designing Therapies and Predicting Outcomes

Cutaneous melanoma is a highly invasive tumor and, despite the development of recent therapies, most patients with advanced metastatic melanoma have a poor clinical outcome. The most frequent mutations in melanoma affect the BRAF oncogene, a protein kinase of the MAPK signaling pathway. Therapies targeting both BRAF and MEK are effective for only 50% of patients and, almost systematically, generate drug resistance. Genetic and non-genetic mechanisms associated with the strong heterogeneity and plasticity of melanoma cells have been suggested to favor drug resistance but are still poorly understood. Recently, we have introduced a novel mathematical formalism allowing the representation of the relation between tumor heterogeneity and drug resistance and proposed several models for the development of resistance of melanoma treated with BRAF/MEK inhibitors. In this paper, we further investigate this relationship by using a new computational model that copes with multiple cell states identified by single cell mRNA sequencing data in melanoma treated with BRAF/MEK inhibitors. We use this model to predict the outcome of different therapeutic strategies. The reference therapy, referred to as “continuous” consists in applying one or several drugs without disruption. In “combination therapy”, several drugs are used sequentially. In “adaptive therapy” drug application is interrupted when the tumor size is below a lower threshold and resumed when the size goes over an upper threshold. We show that, counter-intuitively, the optimal protocol in combination therapy of BRAF/MEK inhibitors with a hypothetical drug targeting cell states that develop later during the tumor response to kinase inhibitors, is to treat first with this hypothetical drug. Also, even though there is little difference in the timing of emergence of the resistance between continuous and adaptive therapies, the spatial distribution of the different melanoma subpopulations is more zonated in the case of adaptive therapy.

Spatial fluxes describe the movement of cells in space. As in previous models of the spatial evolution of cancer Chaplain and Lolas (2005); Hodgkinson et al. (2019), our model considers that cells may move either in a disorganised, or in an directed manner. The disorganised movement corresponds to the space diffusion flux that, according to the Fick law, is proportional to the gradient of cell density and is directed opposite to this gradient. Directed movement corresponds to chemotactic, or haptotactic fluxes that are proportional to chemical gradients, or to the gradient of ECNE, respectively. In both cases cells move uphill gradients, towards higher densities of cells in chemotaxis, or towards higher density of ECNE in haptotaxis. We consider two types of chemotactic gradients: m 1 (x) is a nutrient produced by ECNE and m 2 (x) a chemo-attractant produced by the cancer cells.
The spatial flux term is, therefore, given by: The factor 1 − P c dy has a stabilizing role and was introduced to avoid the Keller-Siegel instability (Keller and Segel, 1970).
Structural fluxes describe the continuous changes of the cell state described as changes of gene expression. Although these changes can be modeled mechanistically using gene and protein interaction networks, in this paper we use a phenomenological model of structural fluxes, based on single cell transcriptomic data. Structural fluxes are also of two types, diffusive and advective, corresponding to random, zero mean, changes of the cell state and deterministic changes of the cell state, respectively. For the diffusive terms, we assume that these vary with the two structural variables y 1 and y 2 as described in the Methods section and Figure 1f in the main text. The advective terms are similarly described in Methods and Figure 1e in the main text and are assumed to have the proliferative and SMC subpopulations, as well as all the points linearly interpolated between these states as stable steady states. The stability of these states was considered higher for larger y 2 , which was modelled as a linear factor (1 − r min )y 2 + r min , 0 < r min < 1, in the advection flux. Thus, we have used the following model for structural fluxes: where (ξ ϕ hyp , ν ϕ hyp ), (ξ ϕmax , ν ϕmax ) are the structural space centroids of the SMC and the proliferative cells subpopulations, respectively.
Cell sources terms are of two types: positive proliferation terms, and negative degradation terms.
Degradation terms are proportional to a combination of the drug concentrations p 1 f 1 (y) + p 2 f 2 (y), where p 1 , p 2 are the drug concentrations and f 1 (y), f 2 (y) are effectiveness functions of the two drugs (see Supplementary Figure 1). The drug effectiveness functions cope with the fact that cells with different states y are eliminated differently by the two drugs.
The proliferation terms take into account competition on resources (logistic growth) and are proportional to the concentration of nutrient m 1 produced by ECNE. The proliferation term also contain a factor depending on the cell state y, that mimics the cell metabolic activity derived in Rambow et al. (2018). For modelling this factor, it is assumed that all cells have a basic proliferation rate, ϕ c ϕ 0 ; that there is a region (centered in (ξ ϕup , ν ϕup ), in between the states defined in the main text, and close to the west part of the structural domain) exhibiting elevated proliferative activities and an elevated rate, ϕ c ϕ up ; and that the cells in the proliferative state exhibit the greatest proliferation rate, ϕ c ϕ max . Meanwhile, the SMC state is assumed to exhibit a significantly lower rate than the remainder of the domain and this is achieved by dividing the entire term by 1, plus a Gaussian function centered in (ξ ϕ hyp , ν ϕ hyp ). We use the following expression for the source terms: Spatial dynamics of other components. As well as the cancer cell population, three types of other species are considered in our model; namely the extra-cellular nutritional environment (ECNE), the chemical species, and the drug species.
As in several previous models of cancer invasion Gatenby and Gawlinski (1996); Chaplain and Lolas (2005); Trucu et al. (2013); Hodgkinson et al. (2019), we assume that the cancer cell population invades the stroma by degrading the ECNE to make space for invading cells. Hence, there is a chemically-mediated degradation, with the rate vector δ v = 0 δ v (only the acidic species m 2 degrade ECNE). The degradation of the ECNE is also achieved through a natural decay term, with rate δ v,0 . Furthermore, the ECNE grows logistically with a rate ϕ v .
The chemical species, themselves, diffuse with a rate vector D m . In particular, the current model considers only two chemical species; the nutritional species m 1 , produced by the ECNE, and the acidic or degradative species m 2 , produced by the cancer cells. The production of these species is given by the rate vector ϕ m . Logistic factors m i (1 − m i ), i ∈ {1, 2} multiply the production terms, in order to keep the chemical species concentrations bounded 0 ≤ m i ≤ 1 for i ∈ {1, 2}. Chemical species are assumed to undergo degradation with rates given by the vector δ m .
We consider two drug species, representing BRAF/MEKi (for simplicity we consider only one species even if two kinase inhibitors are administered) and the hypothetical cancer treatment (HCT), respectively. The drug species exhibit diffusive spatial dynamics, with a rate vector D p . The two drug species, are differentially administered to the patient, on the basis of the required treatment. This differential treatment is mathematically represented by the time-dependent, user-defined, function, Θ p (t). Drug species are assumed to be degraded by natural processes, with rates δ p , and by cellular uptake and degradation, with rates ℓ p .
Full System of Equations. Combining the above equations with the assumptions for the ECNE, chemical, and drug populations, we may then write the following complete set of PDEs as: . The drug acts mainly on cells whose states are located at maxima of the effectiveness function.

Movies
Supplementary