The development of drug resistance in metastatic tumours under

chemotherapy: an evolutionary perspective

Federica Padovano1,* , Luis Almeida1 , and Chiara Villa2,*
1 Sorbonne Université, CNRS, Université de Paris, Laboratoire Jacques-Louis Lions UMR 7598, 75005 Paris, France
2 Sorbonne Université, CNRS, Université de Paris, Inria, Laboratoire Jacques-Louis Lions UMR 7598, 75005 Paris, France
* Correspondance to: federica.padovano@sorbonne-universite.fr (FP), chiara.villa.1@sorbonne-universite.fr (CV)

May 31, 2024

arXiv:2405.20203v1 [q-bio.CB] 30 May 2024

Abstract
We present a mathematical model of the evolutionary dynamics of a metastatic tumour under
chemotherapy, comprising non-local partial differential equations for the phenotype-structured cell pop-
ulations in the primary tumour and its metastasis. These equations are coupled with a physiologically-
based pharmacokinetic model of drug delivery, implementing a realistic delivery schedule. The model is
carefully calibrated from the literature, focusing on BRAF-mutated melanoma treated with Dabrafenib
as a case study. By means of long-time asymptotic analysis, global sensitivity analysis and numerical
simulations, we explore the impact of cell migration from the primary to the metastatic site, physiological
aspects of the tumour sites and drug dose on the development of drug resistance and treatment efficacy.
Our findings provide a possible explanation for empirical evidence indicating that chemotherapy may
foster metastatic spread and that metastatic sites may be less impacted by chemotherapy.

1 Introduction
1.1 Biological context and motivation
In 2022 around 9.7 million of people died for cancer worldwide [103], accounting for approximately 15% of
total deaths, making it one of the primary health problems globally. In particular, metastases are responsible
for approximately 90% of cancer-related mortality [45, 54]. These form following a multi-step process com-
prising the migration of cancer cells from the original site to regional or distant organs and lymph nodes, e.g.
by accessing the lymphatic or blood vessels. In this latter case the metastatic steps include local invasion of
regions surrounding the primary tumour, the cancer-induced formation of new blood vessels to access more
nutrients and foster cancer growth, process known as angiogenesis, intravasation, survival in the circulatory
system and extravasation through vascular walls of distant sites [54]. The process of primary tumour cells
colonization of other tissues is known as metastatic spread, or primary seeding, meanwhile when cells return
to the original site or an existing metastasis they contribute to self-seeding [74, 79]. Moreover, if the popu-
lation in the metastatic site grows and acquires itself the ability to metastasise, it can spread to other sites,
phenomenon known as secondary seeding [79].
Chemotherapy, i.e. the use of cytotoxic drugs to kill cancer cells, is to this day considered the most
effective, and thus most widely used, modality of cancer treatment. The development of drug resistance
confers a selective advantage upon the cancer cell population, as cancer cells exhibit reduced sensitivity to
cytotoxic compounds, and represents a significant challenge in cancer therapies as it often contributes to
disease relapse [27, 42, 44, 97, 101]. Despite the efforts put into the development of treatment strategies
evading or reverting the development of drug resistance [11], this process becomes significantly more com-
plex in metastatic cancers, given the exposure of tumours in distinct organs to disparate concentrations
of therapeutic agents and environmental factors [104]. In fact, it has been reported that even if standard
chemotherapy effectively controls disease progression at the primary tumour site, it often fails to influence
metastatic populations [34, 80]. Moreover, many studies suggest that there is a strong correlation between
treatment resistance and metastatic ability, i.e. cancer cells with lower sensitivity to cytotoxic agents usu-
ally also manifest enhanced disseminating properties [34, 52, 54, 63, 93]. Therefore the development of drug

1
resistance in the primary tumour during chemotherapy, even in reversible cases, may favour the persistence
of its metastases. Nevertheless, it is still unclear whether metastases are intrinsically more resistant than
the primary tumour or if their reduced sensitivity arises because they originate from particularly aggressive
cell subpopulations or due to subsequent evolution after dissemination [54]. It may therefore be beneficial to
test existing hypotheses on the development of drug resistance in metastatic tumours under chemotherapy
using proof-of-concept theoretical frameworks.

1.2 Mathematical modelling background

Mathematical modelling is a valuable tool for gaining insights into the mechanisms driving cancer evolution
and the emergence of intratumour and intertumour heterogeneity, for simulating complex and long-term dy-
namics that may take years to observe in real-life settings, and for designing optimised therapeutic strategies.
All three objectives share a common aim of paving the way towards precision medicine, i.e. the innovative
idea of targeting therapies on patients in order to increase their efficacy.
Many mathematical models have been previously employed to describe the different steps of the metastatic
cascade, see the review article [90] and references therein. These may focus, for instance, on the develop-
ment of metastatic phenotypes [66, 67], their intravasation and extravasation [9, 39, 86], or extensive cell
dissemination and potential evolution of cells from the primary to the metastatic site [12, 30, 47, 72, 73, 89].
Moreover, some of them also investigate the effect of therapeutic treatments, such as chemotherapy, radio-
therapy and immunotherapy, alone or combined, on metastatic cancers [41, 77, 88, 94]. Among these, Sun et
al. [94] accounted for the existence of a sensitive and a drug-resistant subpopulation in the primary site. Nev-
ertheless, drug-resistance levels may be better captured on a continuum [56], and phenotypic heterogeneity
in the metastatic site may also play an important role in the failure of treatment.
Many studies in literature investigate the adaptive processes that underlie the development of intra-
tumour phenotypic heterogeneity and the emergence of resistance to chemotherapeutic agents, see for in-
stance [24, 106] and references therein. Among these, many works consider models comprising nonlocal
partial differential equations (PDEs) modelling the adaptive dynamics of cancer cell populations structured
by a continuous phenotypic trait linked with the cell sensitivity to cytotoxic agents [15, 17, 18, 19, 20, 21,
23, 28, 55, 60, 99]. These models enabled the study of the progressive development of increasing levels
of drug resistance during treatment and the role of intratumour heterogeneity in cancer persistence after
therapy. This was also possible thanks to their higher amenability to analytical investigations compared to
their stochastic counterparts and the mathematical theory supporting such analysis [5, 31, 36, 38, 81]. Some
of these studies further applied optimal control methods in order to explore the most effective treatment
strategies [3, 78, 83]. While this mathematical framework had not been applied to metastatic tumours prior
to this manuscript, Mirrahimi [68] considered the adaptive dynamics of phenotype-structured populations
in communicating patches, further developed with applications to speciation in ecology or host-pathogen
interactions [2, 7, 46, 58, 69], which is easily comparable with the scenario of connected primary tumour and
metastatic sites.
Pharmacokinetic (PK) modelling effectively captures the time course of the drug concentration according
to various administration modes and allows to quantify the relationship between the dose and the in situ
drug concentration [1, 51, 109]. In these models, the body is divided into compartments, also referred to as
building blocks, each of which is associated to a variable, or quantity, describing the drug concentration in the
block. The compartments are then linked to each other through the drug exchange. This kind of models are
usually coupled with a pharmacodynamic (PD) one, which captures the effect of the drug on the considered
disease, and employed to find suitable drug dosing schedules to achieve optimal drug concentrations in the
target tissues [40, 71, 107]. Moreover, they can also be expanded to include physiological aspects of the
considered tissues, such as as organ blood flow and size, in order to explore their influence on the drug
exposure [40, 48].

1.3 Synopsis and paper structure

We propose a mathematical model for the evolutionary dynamics of metastatic tumours under chemotherapy,
where physiological differences of the different tumour sites are integrated within a PK model of drug delivery.
In particular, the model comprises a system of coupled nonlocal PDEs for the phenotypic distribution of
the cancer cells in the primary tumour and the metastatic site, structured by their level of resistance to
the chemotherapeutic agent, and ordinary differential equations (ODEs) for the pharmacokinetics. We focus

2
on a biological framework including a primary tumour that already faced angiogenesis and a newly-formed
yet growing metastasis, where cancer cells characterised by higher levels of drug resistance are assumed to
be more aggressive and thus able to migrate to distant sites at higher rates. We restrict our attention to
a highly perfused primary tumour and a metastatic site where cells are proliferating but not yet able to
disseminate, as they chose in [94], exploring different levels of tumour vascularisation for the metastasis.
We specifically consider BRAF-mutated melanoma, which is a form of skin cancer that develops in
melanocytes, i.e. the cells responsible for melanin production, from the uncontrolled proliferation of cells
induced by a mutation of the BRAF gene. Despite its lower incidence, BRAF-mutated melanoma is the
most aggressive and lethal among the skin cancers [13, 95], particularly due to its high metastatic rate [10].
These tumours can spread locally, regionally and distantly, with the most common metastatic sites being
skin and subcutaneous tissue, followed by lungs, liver, bones, and brain [95]. Metastatic melanomas are
often treated with the chemotherapeutic agent Dabrafenib, a kinase inhibitor of mutated BRAF. Despite the
rapid response, with a median time around 6 weeks, and short-term increase in patient survival, resistance
to Dabrafenib persists with a median progression-free survival of approximately 6–8 months [8]. For these
reasons, a metastatic BRAF-mutated melanoma under Dabrafenib treatment constitutes the ideal case study
to adopt for our model, which we carefully calibrate from the literature employing PK parameter values
estimated from in vivo and ex vivo data from patients.
The paper is organised as follows. In Section 2 we introduce the model assumptions and equations. In
Section 3 we carry out a formal asymptotic analysis of evolutionary dynamics. In Section 4 we perform
global sensitivity analysis and conduct further numerical investigations, to check the analytical results and
explore the role of evolutionary and physiological parameters on the outcome of treatment and the timescale
of development of drug resistance. Section 5 concludes the paper and provides a brief overview of the model
limitations and possible research perspectives.

2 Description of the model

We present a mathematical model of evolutionary dynamics of a metastatic tumour under chemotherapy,
comprising non-local phenotype-structured PDEs for the primary tumour and its metastasis. In order
to effectively capture the time course of the drug concentration in each tumour site depending on the
administration mode and dose, as well as physiological aspects of each tumour site, we complement the
evolutionary dynamics model with a physiologically-based PK model for the drug delivery.

2.1 Evolutionary dynamics model of metastatic cancer

We model the evolution of two tumour cell populations, i.e. the primary tumour and its metastasis, exposed
to a chemotherapeutic agent. In order to consider a metastasised tumour, we assume that the primary tumour
is vascularised and the sites are sufficiently well connected so that the cancer cells in the primary tumour can
intravasate and, subsequently, extravasate in the secondary site. The opposite process, known as secondary
self-seeding [89] is also allowed. In particular, the metastatisation process is modeled by allowing tumour cells
to transition from one site to another [70, 89]. Building upon the ideas presented in [3, 68, 99], we introduce
the variable y ∈ [0, 1], which represents the cell phenotypic state linked to its level of chemoresistance. We
assume that the phenotypic variant y = 1 endows cells with the highest level of cytotoxic-drug resistance
and the greatest migratory abilities. This is motivated by evidence suggesting phenotypes with higher drug
resistance are more aggressive, and they are often associated with a higher invasive potential [93]. On the
contrary, the state y = 0 corresponds to cells with the lowest level of cytotoxic-drug resistance and the lowest
migratory abilities.
From now on we will make use of the index i ∈ {1, 2}, to represent each tumour site. In particular, i = 1
corresponds to the primary tumour, while i = 2 is the metastasis. The phenotypic distribution of tumour
cells at time t ∈ [0, ∞) and site i is described by the function ni (t, y). Moreover, at each time t, we define
the density of tumour cells in site i as Z 1
Ii (t) := ni (t, y)dy , (1)
0
with the corresponding local mean phenotypic state and related variance defined, respectively, as
Z 1 Z 1
1 1
µi (t) := yni (t, y)dy and σi2 (t) := y 2 ni (t, y)dy − µ2i (t) . (2)
Ii (t) 0 Ii (t) 0

3
The phenotypic distribution of tumour cells in each site, ni (t, y), is governed by the following non-local PDE
with given initial and Neumann boundary conditions:

 2
 ∂t ni − βi ∂yy ni = R(y, Ii , Ci )ni + νj,i (y)nj − νi,j (y)ni , i ̸= j, in (0, ∞) × [0, 1]



 Z 1

I (t) :=
i ni (t, y)dy,
0 i = 1, 2 . (3)



 ni (0, y) = ni,0 (y),



∂ n (t, 0) = ∂ n (t, 1) = 0 ,
y i y i

The diffusion term in (3)1 models the effects of spontaneous epimutations, which occur at rate βi ∈ R>0
[17, 60]. The non-local reaction term takes into account the effects of cell proliferation, natural death,
death due to competition for resources, and the cytotoxic action of the drug. The functional Ri (y, Ii , Ci ) ≡
Ri (y, Ii (t), Ci (t)), models the fitness of tumour cells in site i in the phenotypic state y and under the local
environmental conditions at time t, characterised by the cell density Ii ≡ Ii (t), and the chemotherapeutic
agent concentration Ci ≡ Ci (t). Building upon the modelling strategies presented in [3, 99], the fitness
function Ri is defined as follow

Ri (y, Ii , Ci ) := pi (y) − di Ii − ki (y, Ci ). (4)

In definition (4) the term di Ii translates into mathematical terms the idea that a higher total cell number
corresponds to less available resources and space in the system, hence a higher rate of death due to intrapop-
ulation competition. The parameter di ∈ R>0 is related to the local carrying capacity of the tumour. The
function pi (y), also referred to as the background fitness in the absence of treatment, stands for the net
proliferation rate of cancer cells in the phenotypic state y and site i, and based on the ideas proposed in [99]
we define it as
pi (y) := δi (1 − y 2 ) + φi (1 − (1 − y)2 ) , (5)
where δi ∈ R>0 represents the maximum background fitness for highly proliferating and drug sensitive cells
in site i, while φi ∈ R>0 represents the maximum background fitness for weakly proliferating and fully
chemoresistant cells in tumour i, and we further assume δi ≫ φi . Under definition (5), pi (y) reaches its
minimum at y = 1, i.e. the phenotypic trait characterised by the highest level of cytotoxic-drug resistance
and the slowest proliferation in the absence of a drug, since we assume that slowly proliferating cells are less
susceptible to chemotherapy and thus more likely to develop resistance [22, 25]. Meanwhile ki (y, Ci ) is the
rate of death induced by the cytotoxic drug, and based on the ideas proposed in [99], we define

ηi · Ci
ki (y, Ci ) := (1 − y)2 , (6)
αi + Ci
where ηi ∈ R>0 is the maximal death rate of highly drug sensitive phenotypic variants due to the cytotoxic
action of the chemotherapeutic agent, and αi ∈ R>0 is the Michaelis-Menten constant of the chemothera-
peutic agent. Under definition (6), ki (y, Ci ) is a decreasing function of y, i.e. the rate of drug-induced death
decreases as the level of chemoresistance of the cells increases, and it is null for y = 1, consistently with
the assumption that such a phenotype is completely resistant to the chemotherapeutic agent. With these
definitions, after a little algebra, the fitness function in (4) can be rewritten as

2
Ri (y, Ii , Ci ) := ai (Ci ) − bi (Ci ) y − hi (Ci ) − di Ii , (7)

where
·Ci 2

ηi ·Ci
ηi · Ci φi + αηii+C ηi · Ci φi + αi +Ci
ai (Ci ) = δi − + i
·Ci
, bi (Ci ) = δi +φi + and hi (Ci ) = ·Ci
. (8)
αi + Ci δi + φi + αηi+C αi + Ci δi + φi + αηii+C
i i i

Here ai ≡ ai (Ci ) is the maximum fitness, bi ≡ bi (Ci ) is the nonlinear selection gradient, and hi ≡ hi (Ci )
is the phenotypic trait associated to the maximum fitness corresponding to the chemotherapeutic agent
concentration Ci (t). We observe that a higher drug concentration Ci results in a lower maximum fitness ai
because of the greater cytotoxic activity due to the compound, a greater fittest phenotypic trait hi and a

4
stronger selective pressure on tumour cells, i.e. a greater bi .

The terms νi,j (y) and νj,i (y) in equation (3) are non-negative functions representing the rate of transition
of cells in the phenotypic state y from site i to site j ̸= i, and vice versa, which depend on the ability of
the phenotype to intravasate, survive in the circulation and extravasate, given the connectivity of sites i and
j. Motivated by [89], we assume that all cells with phenotypic state y > 0 can access the bloodstream and
reach other sites, and we assume that when cancer cells migrate from one site to the other they maintain
their original phenotypic trait [54]. Specifically, we define the migration rates as

νj,i (y) = ν̂j,i y 2 , i ̸= j, (9)

where ν̂j,i ∈ R>0 is the maximum transition rate from site j to site i, with i ̸= j. Under definition (9),
migration is almost absent as y → 0, while the migration rate approaches its maximum velocity ν̂j,i as y → 1,
i.e. cells with higher levels of cytotoxic-drug resistance migrate at more significant rates, consistently with
our assumptions [93].

2.2 PK model for anticancer drug

We employ a PK model in order to better predict the drug concentration in two physiologically different
tumour sites, given distinct administration methods and doses. To achieve this, we take into account a five-
compartments structure, as shown in Figure 1, comprising the administration site, the central and peripheral
compartments, and the primary tumour and metastasis blocks. The administration site allows us to take into
account the scenarios of extravascular drug injection, i.e. when the compound does not enter directly into
the circulating system, e.g. per os. The central compartment represents plasma and well-perfused tissues,
i.e. tissues that receive a rich blood supply relative to their organ weight [26], such as kidneys, liver, heart,
and brain. Meanwhile, the peripheral compartment represents poorly-perfused tissues, such as muscle and
skin. The primary tumour and metastasis blocks are included in the model structure for two main purposes:
to capture the delay between the drug injection and its therapeutic effect, and to introduce in the model
some of the physiological differences that can affect the compound distribution in different tumour sites,
specifically the tumour size and its vascularization.
We denote by Ca (t) ≥ 0 the quantity of administered anticancer drug in the administration site, and let
Cc (t), Cp (t), C1 (t), C2 (t) ≥ 0 denote the drug concentration at time t in the central compartment, peripheral
block, primary tumour and metastatic site, respectively. Building upon the ideas of [109] we adopt a first-
order kinetics strategy to describe the drug exchange between compartments. Then the mass balance for
the first compartment in the model reads as
(
dCa
dt = rF − ka Ca , (10)
Ca (0) = 0 ,

where F ∈ (0, 1] is the bioavailability, i.e. the fraction of the administered dose that reaches the systemic
circulation, ka ∈ R>0 is the absorption rate, and r ≡ r(t) ≥ 0 is the quantity of administered drug over time.
The evolution in time of the drug concentration in the peripheral compartment is given by
(
dCp
dt = kin,p Cc − kout,p Cp , (11)
Cp (0) = 0 ,

where kin,p ∈ R>0 and kout,p ∈ R>0 are the first-order constant rates of drug concentration respectively
entering and exiting the peripheral compartment. The two physiological aspects we want to consider, i.e.
tumour size and vascularization, are captured by quantities such as tumour volume and in situ blood flow.
We thus introduce the parameters V1 , V2 ∈ R>0 representing the volumes of the primary tumour and the
metastasis, respectively, and Q1 , Q2 ∈ R>0 representing the blood flows through the primary tumour and
metastasis, respectively. We assume the primary tumour site to be more vascularised and thus take Q1 /V1 >
Q2 /V2 . In order to set the dependency of the first-order distribution constant rates of the two tumour sites
on their physiological features, we take inspiration from the physiologically-based PK models and adopt the
perfusion rate-limited strategy [51]. This latter assumes that the drug distributes freely and instantly across
membranes, hence blood perfusion becomes the only limiting process to the drug distribution in the two

5
Figure 1: Schematic of the model. The five boxes represent the building blocks of the PK model, while
the two circles depict the local environment in primary tumour (blue) and metastatic (orange) sites. As
illustrated by the overlapping of the boxes and circles, the local environmental conditions in each tumour
site will be affected by the drug concentration in the respective PK compartment. The red pills and the grey
syringe represent oral administration and intravenous injection, respectively. The straight arrows represent
the drug exchange between compartments, while the dashed ones indicate migration of cancer cells. The
variables used in the model for the drug concentration in each compartment (Ci ) and the phenotypic distri-
bution in each tumour site (ni ) are indicated in the respective box/circle. The parameters and/or functions
used in the model to represent the rate of drug or cell flow are indicated next to the respective arrow.

tumour sites. Therefore, the ODEs describing the evolution in time of the compound concentration in the
primary tumour and its metastasis, are defined as follows
( (
dC1 Q1 R Q1 R dC2 Q2 R Q2 R
dt = V1 Cc − V1 K1 C1 − Ψ1 , and dt = V2 Cc − V2 K2 C2 − Ψ2 , (12)
C1 (0) = 0 , C2 (0) = 0 ,

where the constant R ∈ (0, 1] is the blood-to-plasma drug partition coefficient, K1 , K2 ∈ (0, 1] are the
tumour-to-plasma partition coefficients for primary tumour and metastasis, respectively, and will affect how
the chemical distributes throughout the tissues, w.r.t. the plasma concentration, and are an important part
of any pharmacokinetic study. The terms Ψ1 ≡ Ψ1 (C1 , n1 ) and Ψ2 ≡ Ψ2 (C2 , n2 ) model consumption of the
the cytotoxic drug by the cancer populations and, following with the ideas proposed in [99], are defined as
Z 1
ηi · Ci (t)
Ψi (Ci , ni ) := ψi · (1 − y)2 ni (t, y)dy i = 1, 2 , (13)
αi + Ci (t) 0
where αi and ηi were introduced in (6), while ψi ∈ R>0 is a conversion factor. Finally, the mass balance
equation describing the evolution in time of the drug concentration in the central compartment is given by
  
 dCc = ka C + Q1 C + Q2 C + k C − Cl
+ Q1 +Q2
+ k Cc ,
dt Vc a Vb K1 1 Vb K2 2 out,p p Vc Vb in,p
(14)

Cc (0) = 0

where Vc ∈ R>0 is the volume of the central compartment and Vb ∈ R>0 is the blood volume. This latter
parameter is introduced in the ODE because the tumour sites uptake drug from the circulating system and

6
not from the entire central compartment. Moreover, Cl ∈ R>0 represents the clearance, i.e. the volume of
plasma cleared of a drug over a specified time period. Under the scenario of drug infusion, the factor ka Ca
in equation (14)1 will be replaced by the infusion rate.

3 Analysis of evolutionary dynamics

In this section, we study the long-time behaviour of the system (3), i.e. its solution for t → ∞, with Ri in
the form of (7). To achieve this, we focus on a scenario where the concentration of chemotherapeutic agent
is constant in time, i.e. the function Ci (t) is given and satisfies the following assumption

Ci (t) ≡ ci ≥ 0 , i ∈ {1, 2}.

Consequently, the time dependency of the fitness function is no longer mediated by the drug concentration,
and we here make use of the simplified notation Ri (y, Ii ) ≡ Ri (y, Ii , ci ), along with ai ≡ ai (ci ), bi ≡ bi (ci )
and hi ≡ hi (ci ) for the factors appearing in (7) and defined in (8).

3.1 Assumptions
Fitness functions Ri . We assume there exist positive constants Im , IM such that the following hold:

∂Ri
(y, Ii ) < 0 ∀ y ∈ [0, 1] , Ii ∈ [Im , IM ] , i = 1, 2 ; (15)
∂Ii

∂ 2 Ri
(y, Ii ) < 0 ∀ y ∈ [0, 1] , Ii ∈ [Im , IM ] , i = 1, 2 . (16)
∂y 2
We have introduced the natural assumption that growth is saturated by an overall higher population density,
due to competition for space and resources, which translates mathematically into (15), i.e. each fitness
function is a monotonically decreasing function of the local cell density. In addition, we have assumed that
in each given local environment there is only one fittest phenotypic trait, which translates mathematically
into (16), i.e. each fitness function is strictly concave in y and thus presents only one maximum in [0, 1].
The fittest function definition we give in (7) satisfies assumptions (15) and (16).

Migration rates νi,j (y). We consider generic phenotype-dependent migration rates νi,j (y) ≥ 0, and as-
sume there exist positive constants νm , νM ∈ R≥0 so that νi,j (y) satisfy the following properties:

0 ≤ νm ≤ νi,j (y) ≤ νM < ∞ ∀y ∈ [0, 1], i, j = 1, 2 i ̸= j, (17)

′
νi,j (y) ≥ 0 ∀y ∈ [0, 1], i, j = 1, 2 i ̸= j. (18)
Property (17) expresses the natural assumption that migration rates from one site to the other are non-
negative and bounded for all phenotypic variants. Moreover, we assume that cells with higher levels of
cytotoxic-drug resistance migrate to different sites at higher rates [93], which translates mathematically
into (18), i.e. the migration rates are increasing functions of the phenotypic state y. Definition (9) for the
migration rates satisfies assumptions (17) and (18).

Effective fitness. Notice that we may rewrite (3)1 in terms of the effective fitness function Ri (y, Ii )−νij (y)
for i, j = 1, 2 (i ̸= j). Building on the ideas proposed in [68], we assume this is such that
 
arg max Ri (y, Ii ) − νij (y) is a singleton for i, j = 1, 2 i ̸= j , (19)
y∈[0,1]

   
arg max R1 (y, I1 ) − ν12 (y) ∩ arg max R2 (y, I2 ) − ν21 (y) = ∅ , (20)
y∈[0,1] y∈[0,1]

that is the effective fitness function of each population allows only one maximum and the traits corresponding
to these maxima are distinct.

7
Additional technical assumptions. Following the ideas proposed in [68], and given νm and νM intro-
duced in (17), the analysis will also rely on the additional technical assumption that there exists a constant
δ > 0 such that
  
νm  
δ ≤ min Ri y, Im , Ri y, Im ∀ y ∈ [0, 1] , i = 1, 2 , (21)
νM

    
νm
max Ri y, IM , Ri y, IM ≤ −δ ∀ y ∈ [0, 1] , i = 1, 2 . (22)
νM

3.2 Hyperbolic time scaling and steady state problem

From literature we know that the rate of spontaneous phenotypic variations βi occur on a slower time scale
compared to cell division and death [32, 33]. We therefore introduce a small parameter ε > 0 and assume
both βi := ε2 (i = 1, 2). Following previous studies on the long-time behaviour of non-local PDEs and
integro-differential equations modelling the dynamics of continuously structured populations [6, 16, 29, 31,
50, 62, 70, 82], we use the hyperbolic time scaling t 7→ εt in the conservation

equation (3), and obtain the
following system of non-local PDEs for the phenotypic distributions ni εt , y = niε (t, y) (i = 1, 2):




ε∂t niε − ε2 ∂yy2
niε = Ri y, Iiε niε + νji (y) njε − νij (y) niε


 Z 1
I iε (t) := niε (t, y) dy i = 1, 2 . (23)



 0

∂ n (t, 0) = ∂ n (t, 1) = 0
y iε y iε

As we are interested in the equilibria of (3) in the case of rare phenotypic variations, we follow the strategy
adopted in [68] and investigate the equilibria of (23) in the asymptotic regime ε → 0. Assume that as t → ∞
we have that niε (t, y) → n∞ ∞
iε (y) and Iiε (t) → Iiε (i = 1, 2). Then the equilibria of (23) satisfy the following
system of non-local ODEs

∞
∞ ′′


 Ri y, Iiε niε + ε2 (n∞ ∞ ∞
iε ) + νji (y) njε − νij (y) niε = 0


 Z 1
∞
 Iiε = n∞iε (y) dy
i = 1, 2 . (24)

 0


n∞ (0) = n∞ (1) = 0
iε iε

From now on we will make use of the notation niε (y) and Iiε (i = 1, 2) to refer to n∞ ∞
iε (y) and Iiε , respectively.

3.3 Results of formal analysis

We here summarise the results of the formal analysis, detailed in appendix A, extending the results of
Mirrahimi [68] to the case of phenotype-dependent migration rates between the sites, for our model of
evolutionary dynamics of connected metastatic tumours.

Bounds on I. It can be shown from (24) that for all ε ≤ ε0 , with ε0 small enough, under assumptions
(15), (17), (21) and (22), we have
Im ≤ Iiε ≤ IM i = 1, 2 , (25)

for all ε, where we recall 0 < Im < IM . Proof of this follows analogous steps of that in [68, Lemma 2.1], and
can be found in appendix A.1. As a result, in the asymptotic regime ε → 0 we have Iiε → Ii with

Im ≤ Ii ≤ IM i = 1, 2 . (26)

8
Asymptotic regime ε → 0. Building on the strategies adopted in [68], we introduce the Hopf-Cole
transformation
niε (y) = euiε (y)/ε i = 1, 2 , (27)
2
with uiε (y) semi-convex (i.e. ≥ −E, for some constant E > 0). Then, we expect that I1ε → I1 ,
∂yy uiε
I2ε → I2 , u1ε → u1 and u2ε → u2 in the asymptotic regime ε → 0, where I1 , I2 , u1 and u2 are the leading-
order terms of the asymptotic expansions for I1ε , I2ε , u1ε and u2ε , respectively. Assuming that u1 = u2 = u
and that the semi-convexity of uiε (i = 1, 2) is preserved in the limit such that u = u(y) is also semi-convex,
we have that u is a viscosity solution to the following constrained Hamilton-Jacobi equation
  2
− du

= H(y, I1 , I2 ) ,
dy (28)

 max u(y) = 0 ,
y∈[0,1]

with H(y, I1 , I2 ) being the largest eigenvalue of the matrix



R1 y, I1 − ν12 (y) ν21
(y)
A= .
ν12 (y) R2 y, I2 − ν21 (y)

Given the assumptions introduced in Section 3.1, from (25) and (26) we deduce that n1ε and n2ε converge
weakly to measures n1 and n2 . Moreover, from (28), we have that
K
X
∗
niε (y) −−−⇀ ρik δ(y − yk ) i = 1, 2 , (29)
ε→0
k=1

i.e. the measures n1 and n2 to which n1ε and n2ε converge weakly in the asymptotic regime ε → 0 concentrate
as K Dirac masses centered at yk (k = 1, ..., K), where the weights ρik ≥ 0 must be such that
K
X
Ii = ρik i = 1, 2 , (30)
k=1
 
and where the finite number of points yk ∈ Ω ∩ Γ, with Ω := y ∈ [0, 1] : u(y) = 0 and Γ := y ∈ [0, 1] :
H(y, I1 , I2 ) = 0 . Specifically, this yields


 R1 (yk , I1 ) − ν12 (yk ) ρ1k + ν21 (yk )ρ2k = 0
k = 1, .., K . (31)
 R (y , I ) − ν (y )
ρ + ν (y )ρ = 0
2 k 2 21 k 2k 12 k 1k

The proof of the formal analysis in the asymptotic regime ε → 0 can be found in appendix A.2.

Analytical results for the metastatic spread case: ν1,2 ̸≡ 0 and ν2,1 ≡ 0. The long-time solution of
our system with the fitness function Ri (y, Ii ) defined as in (7) and the migration rate νi,j (y) defined as in
(9), and obtained by solving (31), are the following:

n1 (y) = I1 δ(y − y1 ) and n2 (y) = ρ21 δ(y − y1 ) + ρ22 δ(y − y2 ) , (32)

with
b1
y1 = h1 , y2 = h2 , (33)
b1 + ν̂12
and
   
1 ν̂12 b1 h21 a2 ν̂12 y12
I1 = a1 − , I2 = , ρ21 = min I 1 , I2 , ρ22 = I2 − ρ21 . (34)
d1 b1 + ν̂12 d2 b2 (y1 − h2 )2

The proof can be found in the appendix A.3, while the results for the localised tumours (ν1,2 ≡ 0 and
ν2,1 ≡ 0) and secondary seeding (ν1,2 ̸≡ 0 and ν2,1 ̸≡ 0) scenarios can be found in the supplementary sections
S.1.1 and S.1.2, respectively.

9
3.4 Biological insights
The biological interpretation of the results of the formal analysis summarised in Section 3.3 for the metastatic
spread case yield interesting biological insights. Figure 2 displays the possible model outcomes, based on
the constraints imposed on the migration rates, obtained simulating a non-dimensional version of the model.
We particularly focus on the metastatic spread case here, and further illustrate the dependency of the
solution on certain parameters in Figure 3, but provide more details on the other biological scenarios in the
supplementary material.

Different biological scenarios. From the analysis, and as can be seen in Figure 2, we note that:
(i) In the case of a localised tumour, the primary tumour evolves into a monomorphic population, with
the selected trait being the fittest locally (y1 = h1 ) and the total cell density reaching the carrying
capacity of the site (I1 = ad11 ), and no metastasis forms (cf. Figure 2A).
(ii) In the case of metastatic spread, from (32) we have that the population in the primary tumour will
evolve into a monomorphic population, while the population in the metastatic site will evolve into
a dimorphic population (or a monomorphic one under certain conditions, explained later), as also
illustrated in Figure 2B.
(iii) In the case of metastatic spread with secondary seeding, dimorphism may be observed in both tumours
(cf. Figure 2C) or, in the case of highly connected sites, both tumours may evolve into a monomorphic
population adapted to both environments (cf. Figure 2D), although this scenario is unlikely to be
observed in vivo.

Figure 2: Model outcome dependency on the migration rates. Possible model outcomes under
different proof-of-concept non-dimensional parameter sets. We simulate system (3), under definition (7)
for the fitness functions Ri and definition (9) for the migration rates νij (y), under the initial conditions
n0,1 (y) > 0 and n0,2 (y) = 0 ∀y ∈ [0, 1]. Each plot displays the equilibrium solution in a different biological
scenario: (A) Localised tumour, under the parameter set β1 = β2 = 10−7 , a1 = 8, a2 = 0.1, b1 = 1, b2 = 0.8,
h1 = 0.2, h2 = 0.6, d1 = d2 = 0.2, and ν̂1,2 = ν̂2,1 = 0; (B) Metastatic spread, under the parameter set of
A, except ν̂1,2 (y) = 0.007; (C) Secondary seeding, under the parameter set β1 = β2 = 10−7 , a1 = 6, a2 = 5,
b1 = 1, b2 = 0.6, h1 = 0.2, h2 = 0.6, d1 = d2 = 0.2, ν̂1,2 = 0.1 and ν̂2,1 = 0.05. (D) Secondary seeding in
highly connected sites, under the parameter set of C except ν̂1,2 = 0.6 and ν̂2,1 = 0.3.

Metastatic spread: the primary site. From (32)-(34) we have that the composition and size of the
population in the primary tumour not only depend on the local environment but also on the ability of cells
to intravasate and eventually metastasise. In particular, we observe the following:
(iv) The selected trait in the primary tumour y1 is smaller than the locally fittest one h1 , unlike in the
localised tumour case. Hence, if cells have different metastatic abilities depending on their phenotypic

10
 state, at equilibrium cells less prone to metastasis are selected in the primary tumour, probably as the
more metastatic phenotypes have left the site.
(v) In particular, y1 is a decreasing function of ν̂12 , indicating that the magnitude of this phenotypic shift
is proportional to the migration rate.
(vi) The total population at equilibrium in the primary tumour is lower than the carrying capacity of the
site, i.e. I1 < ad11 , unlike in the case of a localised tumour. Recall however that this comes with
the trade off of having established a dimorphic population in the metastatic site, which increases the
chances of surviving environmental changes, overall resulting in greater evolutionary advantage.

Metastatic spread: the metastatic site. From (32)-(34), we observe that the total cell density of the
metastatic tumour reaches carrying capacity, i.e. I2 = ad22 , but the contribution of cells coming in from the
primary tumour will affect the composition of the local population. In particular, we observe the following:
(vii) We have that ρ21 is an increasing function of I1 and ν̂12 , suggesting that a larger population in the
primary site or a higher migration rate will result in a larger subpopulation in the secondary site
presenting traits selected in the primary, as also illustrated in Figure 3 (third column).
(viii) Moreover, if ν̂12 is significantly high, this might lead to ρ21 = I2 and ρ22 = 0 (cf. bottom plot in the
third column of Figure 3). This means that if the migration rate of cells between sites is particularly high
then the subpopulation in the secondary site presenting traits selected in the primary may outnumber
the subpopulation with traits adapted to the local environment and drive it to extinction. This scenario
can also occur in the case the population of the primary tumour I1 is significantly big w.r.t. to the one
of the metastasis I2 .
(ix) We have that ρ21 is a decreasing function of the nonlinear selection gradient b2 , suggesting that a
stronger selective pressure from the environment in the secondary site will result in a smaller local
subpopulation presenting traits selected in the primary, as also in Figure 3 (first column).
(x) Moreover, in the limit b2 → ∞ we have ρ21 → 0 (cf. bottom plot in the first column of Figure 3).
This means that in the case of extremely strong selective pressure, assuming cells from the primary
site already manage to create a local niche in the metastatic site, this will eventually only evolve into
a monomorphic population, as the trait selected in the primary site will not be fit enough to survive.
(xi) Finally, ρ21 is a decreasing function of the distance between y1 and h2 , indicating that the closer the
selected trait in the primary site is to the fittest trait of the metastatic site, the higher the subpopulation
density in the metastatic site of cells presenting the trait selected in the primary site is. This may
depend on the shift in the phenotypic trait selected in the primary site compared to the fittest one in
the local environment h1 , but also on the similarity of environmental condition sin the two sites, i.e.
the distance |h1 − h2 |, as also illustrated in Figure 3 (second column).

4 Numerical results
We complement the analytical results presented in Section 3 with numerical solutions of the model equations,
focusing on the metastatic spread case and under dynamic drug concentrations as predicted by the PK model,
i.e. we solve system (3), under definitions (4) and (7)-(9), coupled with equations (10)-(14).
In Section 4.1 we report the set-up and parameter values employed for the numerical simulations, along
with the numerical method employed to simulate the model. In Section 4.2 we perform global sensitivity
analysis of the model equations, while in Section 4.3 we investigate how the parameters impact the model
outcome, that is the equilibrium solution of the system of equations, and the time it takes for the system to
reach steadiness, in view of the relative impact this may have during the course of treatment.

4.1 Set-up of numerical simulations and numerical methods.

Set-up and model parametrization. As introduced in Section 1, we consider a metastatic BRAF-
mutated melanoma, assuming no self-seeding takes place (i.e. ν̂2,1 = 0), treated with Dabrafenib orally
administered twice a day in the form of a 150 mg tablet [63]. To carry out numerical simulations of the
system, we calibrate the model with parameters drawn from the literature. For the PK parameters we rely on
the reports on Tafinlar, i.e. brand name for Dabrafenib, realised by the European Medicine Agency (EMA)
[35] and the US Food and Drug Administration (FDA) [98], as reported in Table 1. For the parameters
appearing in the cancer evolutionary dynamic model, we mainly rely on the values reported in [61, 99].

11
Figure 3: Model outcome dependency on input factors in the metastatic spread scenario. Illus-
trative example showing how the selection gradient b2 , the fittest trait h1 and the maximum migration rate
ν̂1,2 may affect the equilibrium distributions of the cancer cell populations of each site, in the metastatic
spread case. We simulate system (3), under definition (7) for the fitness functions Ri and definition (9) for
the migration rates νij (y), under the initial conditions n0,1 (y) > 0 and n0,2 (y) = 0 ∀y ∈ [0, 1]. Each column
displays the equilibrium solutions of three simulations obtained by progressively varying the parameters b2
(first column), h1 (second column) and ν̂1,2 (third column) form their baseline values bB B
2 = 0.8, h1 = 0.2
−7
and ν̂1,2 = 0.007. The remaining parameters are set to β1 = β2 = 10 , a1 = 8, a2 = 0.1, b1 = 1, b2 = bB
B
2 ,
h1 = hB B
1 , h2 = 0.6, d1 = d2 = 0.2, ν̂1,2 = ν̂1,2 and ν̂2,1 = 0.

However, we set the conversion factors ψi (i = 1, 2) in equations (12) to zero, since we consider PK parameters
obtained from a model that was calibrated with ex vivo data and thus the value of the clearance Cl accounts
for the whole process of elimination, including drug consumption by tumour cells [4]. For the migration
rates we consider the ratio ν̂i,j = (Intravasation rate) × (Survival in the circulation) × (Extravasation rate),
where the intravasation rate takes value in the range [10−11 , 10−2 ] cells/day, the probability of survival in
the circulation is in the interval [5 · 10−4 , 2.5 · 10−2 ] and the extravasation probabilities are in the range
[0.1986, 0.5461], as reported in [39]. Building on these works, we identified upper and lower bounds that each
parameter should take, selecting its reference value in this range, as detailed in Table 2.

Numerical methods. Numerical solutions are constructed using a uniform discretisation of the interval
[0, 1], consisting of 101 grid points, as the computational domain of the independent variable y. We consider
t ∈ [0, T ], with T > 0 being the final time of simulations, chosen sufficiently large to reach steady state or
to mimic the duration of therapy. We discretize the interval [0, T ] with a uniform step, sufficiently small to
ensure numerical stability under each parameter set used. We construct numerical solutions employing an
explicit scheme, based on a first order forward difference approximation in time and a second order central
difference approximation in space, applied also to the zero-flux boundary conditions. Following [3, 99], we
set the initial phenotypic distribution in each site ni,0 (y) to be a normal distribution, truncated in [0, 1] and
2
with mean µi,0 and variance σi,0 (i = 1, 2). For both tumour sites i we consider a mean value of µi,0 = 0,
consistently with the assumption that most cancer cells are sensitive to the drug prior to treatment, and
2
a variance of σi,0 = 1/500, and we set as initial population size Ii,0 = dδii . All numerical computations are
performed in Matlab.

12
 PK parameters
Parameter Value Unit Ref
−1
ka 1.8 h [4]
F 0.95 - [35, 98, 84]
Cl 17 h−1 [4, 98]
R 0.54 h−1 [35, 98]
Vc 37.525 l [4, 84]
Vb 5 l [91]
kin,p 0.0974 h−1 [4]
kout,p 0.196 h−1 [4]

Table 1: Pharmacokinetics parameter values.

Tumour site specific parameters

Parameter LB UB RV i = 1 RV i = 2 Unit Ref
−1
Qi - - 0.3 0.01 lh [51]
Vi - - 0.5 0.05 l [40]
Ki - - 0.8 0.5 - [35, 98]
βi 10−13 10−8 10−9 10−9 s−1 [17]
δi 10−5 10−3 10−4 10−4 s−1 [102]
φi 10−7 10−5 10−5 10−5 s−1 [43]
ηi 10−5 10−3 1.8 · 10−4 1.8 · 10−4 s−1 [102]
αi 10−7 10−5 2 · 10−6 2 · 10−6 g cm−3 [75]
di 10−14 10−12 2 · 10−13 2 · 10−13 cm3 s−1 cells−1 [57]
νi,j 10−12 10−8 1.3 · 10−10 0 s−1 [39]

Table 2: Parameter lower bounds (LB), upper bounds (UB) and reference values (RV), for the primary
tumour site (i = 1) and its metastasis (i = 2).

Steady-state criterion. We check the time Tss at which the phenotypic distribution ni (t, y) reaches
steadiness. In order to identify this numerically, we consider a sufficiently small tolerance, tol > 0, below
which the relative difference in the numerical solution at two consecutive time steps is considered to no longer
be significant. Let nki,j denote the numerical approximation of the density of cancer cells in site i (i = 1, 2),
endowed with phenotypic trait yj (j = 0, . . . , J), at time step tk (k = 0, . . . , N ). Then we define the relative
difference dif f i,k between the phenotypic distribution at site i at two consecutive time steps tk−1 and tk as:

J k−1
1 X nki,j − ni,j
dif f i,k = for 1≤k≤N, (35)
J + 1 j=0 nki,j

and denote Tss as the first time step tk after which the inequality dif fi,k < tol is always satisfied for i = 1, 2.

4.2 Global sensitivity analysis

We here make use of Global sensibility analysis (GSA) to verify that the insights gained from the analytical
results of Section 3 still hold outside the asymptotic regime of rare spontaneous phenotypic changes. GSA is
the study of how uncertainty in the output of a model can be apportioned to different sources of uncertainty
in the model input and, specifically, it offers a wide overview of how the parameter interactions impact
the model output. The most suitable GSA techniques for systems modelling highly nonlinear dynamics,
generally the case in biological applications and certainly the case in this work, are the elementary effect
(EE) and Sobol methods [53, 85].

EE method. The EE method is of screening type, i.e. it aims at identifying the parameters that have
negligible impact on the output variability, and is intended for use when dealing with a large number of

13
input parameters, since it has a lower computational cost w.r.t. other methods. Inspired by the radial
design strategy proposed in [14], we consider a vector Z of p independent input parameters Zi (i = 1, ..., p),
varying in their input space. Given two sample points of Z, e.g. a = (a1 , . . . , ap ) and b = (b1 , . . . , bp ),
respectively known as the baseline and the auxiliary points, and given the model output Y = f (Z), we
define the elementary effect EEi for a given (a, b) associated to the i-th input factor as

f ((a1 , . . . , ai , . . . , ap )) − f ((a1 , . . . , bi , . . . , ap ))
EEi = , (36)
ai − bi

and the corresponding measures as

r r  Xr 1/2
1X ∗ 1X 1
EE i = EEik , EE i = |EEik | and SDi = k
(EEi − EE i )2
, (37)
r r r
k=1 k=1 k=1

where r is the number of evaluated EEs. In particular:

∗
• EE i (sensitivity measure) represents the overall influence of the i-th parameter on the output. EE i has
the same meaning but is adopted for avoiding the cancellation effect when dealing with non-monotonic
models.
• SDi (interactions measure) gives information of the non-linearity and/or interaction effects of the i-th
input. If it is small it means that the EEi are similar all along the sample space, suggesting a linear
relationship between the i-th input and the output. On the contrary, if it is large it means that the
EEi are strongly affected by the sample point at which they are computed, thus the parameter effects
are considered to be non-linear and/or due to interaction with other factors.

Sobol method. The Sobol method is a variance-based technique, and is both of screening and ranking
type, i.e. it aims at ordering the inputs based on their impact on the output variability. This approach
is computationally demanding, since it requires a large number of model simulations. For this reason, it is
usually adopted for a small set of input parameters. Given the input vector Z and the model output Y , the
Sobol indices are given by

V ar(E[Y |Zi ]) X V ar(E[Y |Zi , Zu ])

Si = , and ST i = Si + , (38)
V ar(Y ) V ar(Y )
u⊆{1,...,p}\{i}

which are both non-negative. In particular:

• Si provides the first-order contribution of the i-th input to the output variance and measures the main
effect of Zi .
• ST i measures the total effect of the parameter Zi on the output. ST i = 0 implies that the parameter
Zi is non-influential, while a significant difference between Si and ST i suggests that the factor Zi is
involved in important interactions.
We compute the primary and total indices adopting the Monte Carlo method, and in particular we make
use of the estimators for Si and ST i suggested in [87].

Model inputs and outputs. We consider Z = (β1 , β2 , φ1 , φ2 , δ1 , δ2 , α1 , α2 , η1 , η2 , d1 , d2 , ν̂1,2 ) and the

input vector of the GSA, while the other parameters are fixed to the reference values of Tables 1 and 2.
We consider each parameter a random variable, independent from the others, distributed with a log-uniform
distribution in the ranges defined in Table 2. We adopt this distribution in order to better explore the
different orders of magnitude of the input parameters. Furthermore, we take into account a drug schedule
of 150 mg oral tablets twice a day, and set a final time T = 91 days. We consider two model outputs:

• YI = I1 (T ) + I2 (T ), i.e. the total tumour mass at time T ;

µ1 (T )+µ2 (T )
• Yµ = 2 , i.e. the mean phenotypic trait average between the two tumour sites at time T .

We have used Sobol’ quasi-random sequences in order to generate our sets of quasi-random points [14].

14
 (a.1) (b.1)

(a.2) (b.2)

(a.3) . (b.3)

Figure 4: Results of the Elementary effects method for GSA. EE measures associated with the model
output YI (a.1-3) and Yµ (b.1-3), i.e. the total tumour mass and average mean phenotypic state of the sites.
The estimates are computed with r = 500 evaluations of the elementary effect for each parameter. The
elementary effects EEi , defined in (36), calculated for each parameter are displayed in plots (a.1) and (b.1).
∗ ∗
The associated measures EE i , |EE i | and SDi , defined in (37), are shown in plots (a.2) and (b.2). Plots
∗
(a.3) and (b.3) present the ratio between SD and EE i .

15
4.2.1 Results of the EE method for GSA
Figure 4 displays the results of EE method for GSA, from which we remark the following:
∗
1) Figures 4a.3 and 4b.3 display the ratios between SD and EE for the two model outputs. They
both display non-linearity and/or interactions between the parameters, consistently with the nonlinear
nature of the model. In particular, for both outputs the migration rate ν̂1,2 seems to have a highly
non-linear effect and/or to strongly interact with the other factors, consistently with its key role in
connecting they dynamics in the two sites.
∗
2) For YI the indices EE and EE are in agreement (Figures 4a.2 and 4a.1), indicating a lack of cancel-
lation effects and thus a monotonic nature of the input-output response. This is not the case for Yµ
(Figures 4b.2 and 4b.1) for which many parameters show both negative and positive EEi , suggesting
a non monotonic behaviour.

3) Consistently with the analytical results of Section 3.3, in Figure 4a.2 we can observe that parameters
δi , di and ηi (i = 1, 2) have the largest impact on the total cell number YI . Figure 4a.1 shows that
the EEi evaluations for di and ηi are negative, while those for δi are positive, in agreement with
equations (8) and (34).

4) In Figure 4b.2 we can observe that parameters δi , βi , and ηi (i = 1, 2) exhibit the largest impact on
the mean phenotypic trait average Yµ , which is consistent with the analytical results of Section 3.3. In
particular, in Figure 4b.1 we observe that most of the EEi evaluations for ηi are positive, in agreement
with equations (8) and (33) in which yi is proportional to the maximum cytotoxic death rate.

5) It is interesting to notice that βi (i = 1, 2) emerges as a significant input parameter for Yµ , with its
EEi evaluations being all positive (cf. Figure 4b.1). This finding, which could not emerge from the
analytical results due to the asymptotic regime considered in Section 3, reasons with the fact that
higher values of βi should correlate with a larger variance of the phenotypic distributions – e.g. as seen
in [99] – which may shift the mean of the distribution in the bounded domain [0, 1].

6) We note that the migration rate ν̂1,2 has negligible impact on both outputs, likely due to the small
value it may take in the admissible range for the GSA. Considering a wider range of values, i.e. ν̂1,2 ∈
[10−12 , 10−5 ], we observed an increase in the importance of this parameter, as expected, particularly
for Yµ (cf. Figure S.4). Nonetheless, a highly non-linear effect is still observed in both parameter
ranges.

7) The Michaelis-Menten coefficients αi (i = 1, 2) exhibit a really low impact on both outputs, likely due
to the large dose of drug injected. In fact, considering a lower drug dose leads to an increase in the
importance of αi , along with a decrease in the importance of the maximum death rate ηi (cf. Figure
S.5). This suggests that at low drug doses, increasing the cytotoxic efficiency of the compound will not
yield significant therapeutic improvements.

4.2.2 Results of Sobol method for GSA

We conduct further GSA following the Sobol strategy on the most influential parameters highlighted by the
results obtained with the EE method. We thus consider a shorter input vector Z = (β1 , β2 , δ1 , δ2 , η1 , η2 , d1 , d2 )
and fix the rest of the parameters to the values of Table 2. The obtained results, illustrated in Figure 5, are
generally consistent with the findings of the GSA following the EE method. In addition, severable notable
observations can be made:

8) In contrast to the results obtained from the EE method, βi (i = 1, 2) appear to have relatively low
importance for Yµ . Moreover their ST i are significantly bigger than the Si , implying that they mainly
act in interaction with the other parameters. This suggests that higher rates of phenotypic changes,
and thus intra-population heterogeneity, imply greater complexity and non-linearity.

9) Parameters associated with the primary tumour exhibit higher values of Si . This may be due to the
fact that we consider a metastatic spread scenario, where the primary tumour affects the metastasis,
but not the other way around.

16
 (a) Model output YI . (b) Model output Yµ .

Figure 5: Results of the Sobol method for GSA. Sobol measures Si and ST i , defined in (38), associated
with the model output YI (a) and Yµ (b), i.e. the total tumour mass and average mean phenotypic state of
the sites. The estimates are computed from 5000 sample points in the input space.

10) For YI , the significantly higher values of ST i compared to those of Si provide further evidence of
the importance of interactions among these parameters. Meanwhile, for Yµ the Si values occupy a
considerable proportion of ST i , showing that the main effects of ηi and δi are more important than
their interactions.

4.3 Numerical simulations

Given the findings of the GSA, we further explore how the most impactful parameters affect the model
outcome. Specifically, we investigate how different tumour environments may influence possible treatment
outcomes, and how different factors or biological processes may influence the evolutionary timeline of the
tumours under treatment.

4.3.1 Highly communicating sites with different tissue-to-plasma partition coefficients

In this section we take into account two different scenarios, referred to as the baseline and non-baseline
scenarios, which take different values of the maximum migration rate ν̂1,2 and the partition coefficients of
the two sites Ki (i = 1, 2), as summarised in Table 3. Specifically, the non-baseline one is an extreme scenario

Baseline scenario Non-baseline scenario

Parameter i=1 i=2 i=1 i=2 Unit Ref
Ki 0.8 0.5 1 10−4 - [35, 98]
ν̂i,j 1.3 · 10−10 0 1.3 · 10−5 0 s−1 [39]

Table 3: Cancer physiological and evolutionary dynamics parameter values for the baseline and non-baseline
scenarios. Ki represents the partition coefficient of site i, and ν̂i,j the maximum migration rate from site i
to site j, where i = 1 is the primary tumour and i = 2 the metastatic site.

in which we significantly increase the migration rate from primary tumour to the metastasis to a value out of
the range reported in Table 2. We do this for illustrative purposes, as from the results of the GSA we expect
to be able to more easily observe the potential impact of the communication between sites by considering
parameter values outside this range. In addition, we vary the partition coefficients of the tissues in the two
sites, in order to consider how significantly different tumour environments in different sites may affect the
treatment outcome in metastatic tumours. In fact, the parameter Ki indicates the degree of tissue drug
accumulation, attributed to phenomena such as protein binding, lysosomal trapping, and lipid dissolution
[51], and a smaller value of K2 compared to that of K1 implies a lower drug concentration will accumulate
in the metastatic site compared to the primary. In order to mainly focus on the impact of these biological

17
factors on the treatment outcome, the other parameters are fixed to their reference values of Tables 1 and
2, and are thus the same in both scenarios.

Different in situ drug concentrations. Figure 6 displays the drug concentration in four compartments,
excluding the administration site for visualisation purposes. First of all, we remark that the peak of the
plasma compound concentration (Cc ) is reached at the time tmax = 0.93 h, and with a value of Cc,max =
0.0022 g/l, of the same order of magnitude of the values found in literature [35, 98]. Moreover, the drug
concentration in the metastatic site (C2 ) is significantly lower in the non-baseline scenario compared to the
baseline one, reaching a maximum value of the order of 10−7 g/l while in the baseline scenario it reaches higher
peaks comparable to those reached by the drug concentration in the primary site (C1 ). This is consistent
with the biological meaning of Ki , and highlights how different tissue-to-plasma partition coefficients may
result in different drug concentrations in the two tumour sites, creating substantial discrepancies between
the two tumour environments.

(a) Baseline scenario. (b) Non-baseline scenario.

Figure 6: Drug concentrations under baseline and non-baseline scenarios. Pharmacokinetics results
of the numerical simulations under the baseline (a) and non-baseline (b) scenarios with the drug schedule
of 150 mg orally administered twice a day, and a final time T = 4 days. For each scenario we plot the drug
concentration in the central (blue), peripheral (orange), primary tumour (yellow) and metastatic tumour
(purple) compartments. More details on the simulation set-up and numerical methods can be found in
Section 4.1.

Different evolutionary outcomes of the tumours. Such a difference in drug concentration in the
two sites results in different evolutionary dynamics of the two tumour populations, unlike in the baseline
scenario where no significant difference is observed between the phenotypic distributions in the two sites, as
illustrated in Figure 7 comparing the long-term evolution of the tumour populations in the two scenarios.
In particular, in the non-baseline scenario we observe that the lower drug concentration accumulating in the
metastatic site results in a population mostly characterised by lower levels of drug resistance compared to
that in the primary site, exposed to a higher compound concentration. Nonetheless, the high migration rate
results in a small subpopulation in the metastatic site presenting levels of resistance comparable to those
developed in the primary tumour. In the baseline scenario, cell migration has no meaningful impact on the
phenotypic distribution given the similarity of the tumour environments, consistently with the analytical
results illustrated in Figure 3 (second column).

Consistency with analytical results. The obtained results are consistent with the analytical study of
Section 3. In fact, from definitions (8), by significantly decreasing C2 we obtain a lower fittest phenotypic
trait h2 and selective gradient b2 , and a higher maximum fitness a2 . Following (33), this results in a lower
selected trait y2 and a higher metastatic population size I2 , which is exactly what we detect in Figure 7.
Given the considerably different selected traits y1 and y2 , and the increased migration rate, the population
density of the metastatic site shows polymorphism, in accord with equations (32) and (34). Moreover, we
detect a slight decrease in both the primary tumour population size and selected trait, in agreement with
equation (34) considering that the order of magnitude of the migration rate is lower than that of b1 .

18
Sensitivity of cancer evolutionary dynamics to temporal oscillations in drug concentrations.
Indeed temporal oscillations in the in situ compound concentration, due to the drug administration schedule,
imply temporal oscillations in the moments of the phenotypic distribution during treatment, as can be better
observed by plotting the solutions of Figure 7 over just two days (cf. Figure S.6). As already clear from
Figure 7, the great variability of the I2 curve in the non-baseline scenario indicated that the metastatic
population size is considerably susceptible to the drug concentration variation. In fact, it can be detected
from the zoomed-in plots that under both scenarios Ii is the most sensitive to the drug concentration
variations, while both µi and σi2 are substantially less susceptible. This may be explained by noticing that
while the pharmacokinetics are in the scale of hours, the cancer evolutionary dynamics, i.e. the phenotypic
adaptation of the cancer cell populations to the local environments, are in the scale of days. Nonetheless, we
remark that the fluctuations of these quantities are negligible in the baseline scenario and the primary site
in the non-baseline scenario, where their amplitude is small compared to the average value of the respective
quantity. This is not as much the case in the metastatic site of the non-baseline scenario, characterised by
a lower average drug concentration and thus higher sensitivity to small variations of this quantity. We thus
conclude that a higher administered drug dose may result in a lower evolutionary dynamics sensibility to
the compound variability.

4.3.2 Oral administration vs intravenous injection

In the previous section we remarked that under the drug administration regimen of a 150 mg oral dose twice a
day, in the baseline scenario, the impact of the drug concentration oscillations on the evolutionary dynamics
and overall therapeutic outcome is negligible. We thus wish to conduct further numerical investigations
considering the equilibrium phenotypic distribution to which the system converges at the end of treatment,
as well as the time it takes to reach this. However, considering oral drug administration the system does
not reach numerical steadiness for either scenario, as can be seen from the lower-right panels of Figure 7
displaying the step-difference dif fi,k introduced in equation (35), which becomes approximately periodic
with a 12 h period (cf. zoom-in plot in Figure S.6), i.e. the time interval between two consecutive drug
administrations. This would not be the case if we considered a constant intravenous drug injection.
Therefore, we aim at finding the intravenously-injected constant drug dose leading to cancer evolutionary
dynamics analogous to those obtained with an orally-administered drug dose of 150 mg twice a day, which we
will use in the next section for convenience, acknowledging that Dabrafenib is administered per os to patients.
Given the final distribution obtained with the dose administered twice a day n∗i (T ) and the final distribution
obtained with a constant intravenously-injected dose ndose
i (T ), we define the relative error between the two
outcomes as
J
1 X n∗i,j (T ) − ndose
i,j (T )
erri,dose = . (39)
J j=1 n∗i,j (T )

Minimizing the function f (dose) = 0.5 (err1,dose + err2,dose ) in the baseline scenario, we obtain that the
intravenously-injected drug dose that best reproduces the evolutionary outcomes of the oral administration
regimen is of 2.6915 µg/s. We verified that this optimal dose also provides a good approximation of the
dynamics obtained in the non-baseline scenario, as displayed in Figure 8. Moreover, we report the times Tss,i
at which the phenotypic distributions in the two sites reach steadiness, which are larger than the median
time of response and below the median Dabrafenib treatment durations reported in [64].

Tss,1 Tss,2
Baseline 66.98 days 67.02 days
Non-baseline 65.24 days 59.16 days

Table 4: Steady state time of the phenotypic distribution of the primary tumour (Tss,1 ) and the metastasis
(Tss,2 ) in the baseline and non-baseline scenarios under an intravenously-injected constant drug dose of
2.6915 µg/s. The steady state time is the first time at which dif fi,k < tol, with dif fi,k defined in (35) and
a tolerance of tol = 10−6 .

19
 (a) Baseline scenario

(b) Non-baseline scenario

Figure 7: Cancer evolutionary dynamics under baseline and non-baseline scenarios. Cancer
evolutionary dynamics results of the numerical simulations under the baseline (a) and non-baseline (b)
scenarios with the drug schedule of 150 mg orally administered twice a day, and a final time T = 105 days.
For each scenario we plot in the first row the phenotypic distribution over time in the primary tumour (first
panel) and the metastatic tumour (second panel), and compare them at the final time T (third panel), and
we plot in the second row the population size (first panel), mean phenotypic trait (second panel) and step
difference dif fi,k (third panel), as respectively defined in equations (1), (2) and (35). More details on the
simulation set-up and numerical methods can be found in Section 4.1.
20
Figure 8: Oral administration vs intravenous injection under non-baseline scenario. Cancer evolu-
tionary dynamics results under the non-baseline scenario and a drug schedule of 150 mg orally administered
twice a day (VD) or constant intravenous injection of 2.6915 µg/s (CD), and final time T = 105 days. From
left to right, the panels display the cancer population sizes, the mean phenotypic traits and the corresponding
variances over time for the two tumour cell populations, respectively defined in equations (1), and (2), and
the phenotypic distributions at time T . More details on the simulation set-up and numerical methods can
be found in Section 4.1.

4.3.3 Steady-state study of treatment outcomes

Having identified the dose of constant intravenous drug injection that yields analogous therapeutic outcomes
to the oral administration regimen used in the clinic, we now investigate how the therapeutic outcome may be
affected by some key factors. Specifically, we consider the level of vascularisation of the metastatic tumour,
i.e. a physiological factor which will regulate the distribution of the compound in the metastatic site, and
the most important parameters emerged from the GSA focussing on Yµ as the model output , i.e. δi , ηi
(i = 1, 2) and βi . This choices is motivated by the fact that the mean phenotypic trait µi in the two tumours
gives us an indication of the development of drug resistance during treatment. Moreover, as highlighted
in the previous section, µi is the least affected by the drug concentration variability, and we thus expect
conclusions drawn in this section to hold in both drug delivery regiments with higher certainty. We also
remark that from the analytical results of Section 3 we expect that the total population in the sites (Ii )
depends on the selected traits (y1 and y2 ), but not viceversa.
To test the influence of these factors on the treatment outcome, for each parameter we vary its value in
the range defined in Table 2, and compare the model outputs, i.e. the cancer population sizes Ii and mean
phenotypic traits µi in the two sites, after 30 weeks of continuous intravenous constant drug administration,
i.e. the average length of Dabrafenib treatment [64]. For each simulation we also report on when the drug
concentration and the phenotypic distributions in the two sites reach steady state, as it gives insights in
how the different input factors affect the speed of the cancer cell adaptation to the drug regimen. This also
suggests reasonable timeframes in which to assess the advancement of the therapeutic intervention.

Fitness-related parameters. Figure 9a and Figure 9b display the model sensitivity to the maximum
background fitness for highly proliferating yet sensitive cells, δi (i = 1, 2), and the maximal death rate of
sensitive cells, ηi (i = 1, 2), respectively. We notice that incrementing δi speeds up the adaptive dynamics,
increases the cancer cell population sizes, while decreases the mean phenotypic traits. Meanwhile, consis-
tently with its biological meaning, we detect that incrementing ηi decreases the cancer cell population sizes,
while increases the mean phenotypic traits, but has no consistent impact on the adaptation speed. All these
findings are in line with previous analytical results in the literature [3, 99].

Epimutation rate. Since the effects of the epimutation rates βi (i = 1, 2) could not emerge from the
asymptotic analysis of Section 3, and the results of the GSA yield contrasting information on the importance
of these parameters, we consider the results of the numerical steady-state investigation, reported in Figure 9c.
We immediately notice that as βi take greater values the phenotypic distributions ni reach steadiness faster,
implying a greater speed of adaptation. On the other hand, in the parameter range considered βi do not
appear to affect Ii and µi , as suggested by the results of the Sobol method for GSA, although increasing the
epimutation rates results in a higher variance of the phenotypic distribution (cf. Figure S.7), resulting in a

21
 (a) Maximal background fitness

(b) Maximum cytotoxic death rate

(c) Epimutation rate

Figure 9: Treatment outcome dependency on δi , ηi and βi . Treatment outcomes obtained from

numerical simulations of the model under the baseline scenario, for an intravenously-injected drug dose of
2.6915 µg/s and a final time of T = 210 days. The graphs display the steady-state times (first panel), mean
phenotypic traits (second panel) and tumour sizes (third panel) at steady state, varying: a) the maximal
background fitness δi , and ϕi to ensure the assumption δi ≫ ϕi still holds; b) the maximum cytotoxic death
rate ηi ; c) the epimutation rate βi . The steady state time is the first time at which dif fi,k < tol, with dif fi,k
defined in (35) and a tolerance of tol = 10−6 . More details on the simulation set-up and numerical methods
can be found in Section 4.1.

22
more heterogeneous population. These findings are consistent with analytical results in the literature [3, 99],
highlighting that higher epimutation rates increase phenotypic diversity, in turn speeding up natural selection.

Vascularisation of the metastatic site. Lastly, we study the effect of the level of vascularization of
the metastatic site on the evolutionary timeline during treatment by varying the blood flow-to-volume ratio,
i.e. Q2 /V2 . From Figure 10 it is evident that a lower blood flow-to-volume ratio delays the time at which
the in situ compound concentration reaches its equilibrium value. This can be seen in Figure 10a for C2
compared to the drug concentration in the other model compartments, and it is consistent with the biological
interpretation of these factors: given that there is less blood entering the metastatic site, and consequently
less compound, it takes much more time for the tissue to fill up. We remark that for Q2 /V2 ≥ 2 · 10−2
1/h, i.e. for particularly well vascularised metastatic tumours, the steady-state time of the local phenotypic
distribution remains constant, as can be observed in Figure 10b. Nonetheless, for lower values of the ratio
Q2 /V2 we observe that the steady-state time for n2 significantly increases, indicating that for particularly
low levels of tumour vascularisation the slower drug perfusion delays the cancer evolutionary dynamics.

(a) Drug concentrations for Q2 /V2 = 2 · 10−3 1/h. (b) Steady-state times for varying Q2 /V2 .

Figure 10: Evolutionary timescale dependency on the vascularisation of the metastatic site.
Time evolution of the drug concentration and steady-state times obtained from numerical simulations of the
model under the baseline scenario, for an intravenously-injected drug dose of 2.6915 µg/s, and final time
T = 210 days, varying the blood flow-to-volume ratio in the metastatic site Q2 /V2 . In (a) we set Q2 /V2 =
2 · 10−3 1/h, leaving the other parameters at their reference value, and display the drug concentration in the
central (blue), peripheral (orange), primary tumour (yellow) and metastatic tumour (purple) compartments,
up to day 105. In (b) we display the steady-state times for n1 (blue), n2 (red) and C2 (yellow) for different
values of the ratio Q2 /V2 , specifically Q2 /V2 ∈ {2 · 10−5 , 2 · 10−4 , 2 · 10−3 , 2 · 10−2 , 2 · 10−1 , 1} 1/h. The steady
state time is the first time at which dif fi,k < tol, with dif fi,k defined in (35) and a tolerance of tol = 10−6 .
More details on the simulation set-up and numerical methods can be found in Section 4.1.

5 Conclusions
5.1 Summary of major results
In this paper, we conducted a mathematical study of the evolutionary dynamics of metastatic tumours under
chemotherapy. Specifically, our approach involves formal asymptotic analysis and numerical simulations of
a system of non-local PDEs describing the phenotypic evolution of the tumour cells coupled with a system
of ODEs modelling the cytotoxic drug pharmacokinetics.
Focussing on a biological scenario allowing for metastatic spread but no secondary or self-seeding, the
analytical results indicate that, while the primary tumour will evolve into a monomorphic population with
trait selected according to local environmental conditions, the metastasis may display polymorphism if local
environmental conditions differ from those in the primary site, with a subpopulation of cells presenting traits
selected in the primary tumour. Moreover, the size of this subpopulation increases as the connectivity of the

23
sites and the migratory ability of cells selected in the primary tumour increases, or as the selective pressure
in the metastatic site decreases.
We performed global sensitivity analysis and long-term numerical simulations of the full system, to verify
the analytical results outside the asymptotic regime of rare phenotypic changes due to epimutations, and
to better understand the consequences of cancer adaptive dynamics on the outcome of treatment. To do
so we adopted a case study of metastatic melanoma under Dabrafenib monotherapy, through which we
could gain deeper insights into the impact of tissue-specific physiological parameters and drug kinetics on
the evolution of the metastatic cancer. We found that the tumour stage and its location may influence
drug delivery [108], as the tumour angiogenic switch and the organ characteristics may affect the overall
level of tumour vascularisation, ultimately regulating the total amount of in situ cytotoxic compound. As a
consequence, discrepancies in drug concentrations between the two tumour sites may arise during treatment,
resulting in differing local environments to which cells will adapt. In line with the results of the asymptotic
analysis, this fosters phenotypic diversity among tumours at different locations and the emergence of distinct
subclones [96, 105].
In particular, assuming that the primary tumour has significantly increased its vascular supply prior
to undergoing metastatic spread, the higher local concentration of the cytotoxic drug will result in the
selection of more resistant phenotypic traits in the primary tumour. Conversely, assuming the metastasis
is yet to undergo substantial angiogenesis for further dispersal, the expected lower vascular supply of the
secondary tumour should result in a lower drug concentration at this site. This correlates with a lower
selective pressure in the metastatic site which, together with the selection of cells with higher migratory
abilities in the primary tumour, favours the emergence of a subpopulation of highly drug-resistant cells
in the metastasis, in agreement with the results of the asymptotic analysis. Overall, this confers further
evolutionary advantage to the tumour population, as it will have a higher chance of survival to sudden
environmental changes, a big trade off feature of biological dispersal despite leaving the primary site may
be a risk for tumour cells. This is coherent with empirical studies suggesting that chemotherapy fosters
metastatisation and may provide an explanation for the evidence indicating that the metastatic site is often
less impacted by chemotherapy [34, 65, 76, 80].
Finally, investigating the impact of drug kinetics on cancer evolution under the twice-a-day oral adminis-
tration protocol used in the clinic, we observed that the pharmacokinetics and the cancer adaptive dynamics
evolve on different timescales, i.e. hours vs. days. This discrepancy implies that the long-term evolutionary
outcome of treatment may not be significantly influenced by short-term PK dynamics and fluctuations of the
in situ drug concentration, but rather by average cytotoxic levels throughout the duration of chemotherapy,
especially when the drug dose is particularly high.

5.2 Limitations and perspectives

Cell migration and metastatic spread. We chose the migration rate ν12 (y) to significantly increase only
as the phenotypic state of a cell gets closer to the highest levels of drug resistance, but it would be relevant to
consider alternative definitions, in particular non-smooth ones, as they might better reflect the switch in the
intravasation ability of cancer cells as they become more aggressive. The use of different definitions for ν12 (y)
in the PDE for n1 and that for n2 would also be relevant in view of the fact that successful intravasation of
cancer cells does not necessarily correlate with survival in the circulation and successful extravasation to the
metastatic site. In fact, circulating tumour cells in the bloodstream may be attacked by the immune system
and face the challenge of adhesion to the vascular walls to exit the circulation, and may potentially endow
themselves with more resistant traits [37]. The addition of a PDE describing the evolutionary dynamics
of circulating tumour cells may therefore provide a more realistic extension of our modelling framework,
which currently relies on the limiting assumption that cells exiting the primary tumour will immediately
appear in the metastatic site maintaining their original phenotypic state. In addition, the extension of this
model to the one including motion between several patches, e.g. as in [68], could be used to investigate more
complex dynamics arising in cancers presenting multiple metastatic sites. This could allow for the theoretical
exploration of the seed and soil hypothesis [37] of metastatic spread, accounting for different extravasation
probabilities to different sites, e.g. as in [39]. Moreover, this mathematical framework may help explain why
secondary seeding is much more efficient at spreading cells than dispersal from the primary site [89, 100],
thanks to promising results of previous studies on host-pathogen interactions [2, 46] analysing how pathogen
persistence is mediated by the presence of a middle host.

24
Cell proliferation and tumour burden. Building upon the strategies presented in [3, 99], in our model
we took into account a one-dimensional phenotypic trait whereby higher levels of drug resistance correlate
with lower proliferation rates of cancer cells [22, 25]. Nevertheless, subclones simultaneously exhibiting high
proliferation rates and high levels of drug resistance may exist [92]. It would therefore be biologically signif-
icant to consider a two-dimensional phenotypic state, as previously done in [17], to capture the phenotypic
evolution of tumours without any a priori assumption of a link between the level of drug tolerance and
proliferation potential of the cells. This may affect the relation, emergent from the model, between the de-
velopment of resistance and the evolution of the tumour burden during therapy. On this note, it is important
to remark that we included the tumour volume Vi as a physiological parameter in the PK model, but kept
this quantity constant throughout the simulations, for the sake of simplicity. Indeed the tumour volume is
dependent on the number of cancer cells composing it, and it would therefore be significant to allow Vi to
vary in time as a function of the cell number Ii in future work. This would allow to account for the impact
on the drug delivery of changes in the tumour burden expected to occur during chemotherapy.

Therapeutic strategy and personalised treatment. Given the novelty of this work, we here restricted
our attention to a case study of BRAF-mutated melanoma under Dabrafenib monotherapy. In the clinic,
however, this is usually combined with Trametinib [59], a kinase inhibitor of MEK protein adopted to
prevent the development of resistance to the BRAF inhibitor [63]. Therefore, future work investigating
the development of drug resistance in metastatic melanoma under treatment should extend this modelling
framework to include a PK model for Trametinib, as they did in [94], and its drug resistance-inhibiting effect
on cancer cells, e.g. modelled via a drift in phenotypic space, being aware that compounds interact with each
other when administered simultaneously. Moreover, while in this theoretical study we adopted parameter
values collected from the literature, in order to employ this model for treatment optimisation purposes, it
should rely on patient-specific data to improve reliability. To achieve this it is necessary to collect ex vivo
data on plasma and tissue compound concentrations, and use them for model calibration using a two-stage
approach, as done in [48]. First, one should estimate the systemic PK parameters, i.e. the ones related
to central and peripheral compartments, based on the plasma concentration-time profile. Secondly, tissue
concentration–time profiles should be used to estimate the remaining tissue-specific parameters. We note,
however, that experimental methods to measure the partition coefficient parameters, independently of the
mathematical model, do exist [49].

Acknowledgements
This project has received funding from the European Union’s Horizon 2020 research and innovation pro-
gramme under the Marie Sklodowska-Curie grant agreement No 945298-ParisRegionFP. C.V. is a Fellow of
the Paris Region Fellowship Programme, supported by the Paris Region.

References
[1] T. A. Ahmed, Pharmacokinetics of drugs following IV bolus, IV infusion, and oral administration,
vol. 10, IntechOpen Rijeka, 2015.

[2] M. Alfaro, F. Hamel, F. Patout, and L. Roques, Adaptation in a heterogeneous environment

II: To be three or not to be, Journal of Mathematical Biology, 87 (2023), p. 68.

[3] L. Almeida, P. Bagnerini, G. Fabrini, B. D. Hughes, and T. Lorenzi, Evolution of cancer cell
populations under cytotoxic therapy and treatment optimisation: insight from a phenotype-structured
model, ESAIM: Mathematical Modelling and Numerical Analysis, 53 (2019), pp. 1157–1190.

[4] D. Balakirouchenane, S. Guégan, C. Csajka, A. Jouinot, V. Heidelberger, A. Puszkiel,

O. Zehou, N. Khoudour, P. Courlet, N. Kramkimel, et al., Population pharmacoki-
netics/pharmacodynamics of dabrafenib plus trametinib in patients with BRAF-mutated metastatic
melanoma, Cancers, 12 (2020), p. 931.

[5] G. Barles, L. C. Evans, and P. E. Souganidis, Wavefront propagation for reaction-diffusion

systems of PDE, Duke Mathematical Journal, 61 (1990), pp. 835 – 858.

25
 [6] G. Barles, S. Mirrahimi, and B. Perthame, Concentration in Lotka-Volterra parabolic or integral
equations: a general convergence result, Methods and Applications of Analysis, 16 (2009), pp. 321–340.

[7] V. Boussange and L. Pellissier, Eco-evolutionary model on spatial graphs reveals how habitat
structure affects phenotypic differentiation, Communications Biology, 5 (2022), p. 668.

[8] S. Bowyer, R. Lee, A. Fusi, and P. Lorigan, Dabrafenib and its use in the treatment of metastatic
melanoma, Melanoma Management, 2 (2015), pp. 199–208.

[9] G. W. Brodland and J. H. Veldhuis, The mechanics of metastasis: insights from a computational
model, (2012).

[10] J. Budczies, M. von Winterfeld, F. Klauschen, M. Bockmayr, J. K. Lennerz, C. Denkert,

T. Wolf, A. Warth, M. Dietel, I. Anagnostopoulos, et al., The landscape of metastatic
progression patterns across major human cancers, Oncotarget, 6 (2015), p. 570.

[11] K. Bukowski, M. Kciuk, and R. Kontek, Mechanisms of multidrug resistance in cancer chemother-
apy, International journal of molecular sciences, 21 (2020), p. 3233.

[12] I. M. Bulai, M. C. De Bonis, C. Laurita, and V. Sagaria, Modeling metastatic tumor evolution,
numerical resolution and growth prediction, Mathematics and Computers in Simulation, 203 (2023),
pp. 721–740.

[13] L. C. Cabaco, A. Tomas, M. Pojo, and D. C. Barral, The dark side of melanin secretion in
cutaneous melanoma aggressiveness, Frontiers in Oncology, 12 (2022), p. 887366.

[14] F. Campolongo, A. Saltelli, and J. Cariboni, From screening to quantitative sensitivity analysis.
a unified approach, Computer physics communications, 182 (2011), pp. 978–988.

[15] R. H. Chisholm, T. Lorenzi, and J. Clairambault, Cell population heterogeneity and evolu-
tion towards drug resistance in cancer: biological and mathematical assessment, theoretical treatment
optimisation, Biochimica et Biophysica Acta (BBA)-General Subjects, 1860 (2016), pp. 2627–2645.

[16] R. H. Chisholm, T. Lorenzi, and A. Lorz, Effects of an advection term in nonlocal lotka–volterra
equations, Communications in mathematical sciences, 14 (2016), pp. 1181–1188.

[17] R. H. Chisholm, T. Lorenzi, A. Lorz, A. K. Larsen, L. N. d. Almeida, A. Escargueil, and

J. Clairambault, Emergence of drug tolerance in cancer cell populations: an evolutionary outcome
of selection, nongenetic instability, and stress-induced adaptation, Cancer research, 75 (2015), pp. 930–
939.

[18] H. Cho and D. Levy, Modeling the dynamics of heterogeneity of solid tumors in response to
chemotherapy, Bulletin of Mathematical Biology, 79 (2017), pp. 2986–3012.

[19] , Modeling continuous levels of resistance to multidrug therapy in cancer, Applied Mathematical
Modelling, 64 (2018), pp. 733–751.

[20] , Modeling the chemotherapy-induced selection of drug-resistant traits during tumor growth, Jour-
nal of theoretical biology, 436 (2018), pp. 120–134.

[21] , The impact of competition between cancer cells and healthy cells on optimal drug delivery, Math-
ematical Modelling of Natural Phenomena, 15 (2020), p. 42.

[22] E. Chu and A. Sartorelli, Cancer chemotherapy, Basic & Clinical Pharmacology, (2004).

[23] J. Clairambault, An evolutionary perspective on cancer, with applications to anticancer drug re-
sistance modelling and perspectives in therapeutic control, Journal of Mathematical Study, (2019),
pp. 1–21.

[24] J. Clairambault and C. Pouchol, A survey of adaptive cell population dynamics models of emer-
gence of drug resistance in cancer, and open questions about evolution and cancer, Biomath, 8 (2019),
p. 23.

26
[25] P. G. Corrie, Cytotoxic chemotherapy: clinical aspects, Medicine, 39 (2011), pp. 717–722.

[26] C. Cuenod and D. Balvay, Perfusion and vascular permeability: basic concepts and measurement
in DCE-CT and DCE-MRI, Diagnostic and interventional imaging, 94 (2013), pp. 1187–1204.

[27] I. Dagogo-Jack and A. T. Shaw, Tumour heterogeneity and resistance to cancer therapies, Nature
reviews Clinical oncology, 15 (2018), pp. 81–94.

[28] M. Delitala and T. Lorenzi, A mathematical model for progression and heterogeneity in colorectal
cancer dynamics, Theoretical population biology, 79 (2011), pp. 130–138.

[29] L. Desvillettes, P. E. Jabin, S. Mischler, and G. Raoul, On selection dynamics for continuous
structured populations, Communications in Mathematical Sciences, 6 (2008), pp. 729–747.

[30] D. Diego, G. F. Calvo, and V. M. Pérez-Garcı́a, Modeling the connection between primary and
metastatic tumors, Journal of mathematical biology, 67 (2013), pp. 657–692.

[31] O. Diekmann, P.-E. Jabin, S. Mischler, and B. Perthame, The dynamics of adaptation: an illu-
minating example and a Hamilton–Jacobi approach, Theoretical population biology, 67 (2005), pp. 257–
271.

[32] W. Doerfler and P. Böhm, DNA methylation: development, genetic disease and cancer, vol. 310,
Springer Science & Business Media, 2006.

[33] P. Duesberg, R. Stindl, and R. Hehlmann, Explaining the high mutation rates of cancer cells
to drug and multidrug resistance by chromosome reassortments that are catalyzed by aneuploidy, Pro-
ceedings of the National Academy of Sciences, 97 (2000), pp. 14295–14300.

[34] C. D’Alterio, S. Scala, G. Sozzi, L. Roz, and G. Bertolini, Paradoxical effects of chemotherapy
on tumor relapse and metastasis promotion, in Seminars in cancer biology, vol. 60, Elsevier, 2020,
pp. 351–361.

[35] European Medicine Agency, EMA/CHMP/242419/2013/corr 1 - CHMP assessment

report: Tafinlar. https://www.ema.europa.eu/en/documents/product-information/
tafinlar-epar-product-information_en.pdf, 2013.

[36] L. C. Evans, Partial differential equations, vol. 19, American Mathematical Society, 2022.

[37] I. J. Fidler, The pathogenesis of cancer metastasis: the ‘seed and soil’ hypothesis revisited, Nature
reviews cancer, 3 (2003), pp. 453–458.

[38] W. H. Fleming and P. E. Souganidis, PDE-viscosity solution approach to some problems of large
deviations, Annali della Scuola Normale Superiore di Pisa-Classe di Scienze, 13 (1986), pp. 171–192.

[39] L. C. Franssen, T. Lorenzi, A. E. Burgess, and M. A. Chaplain, A mathematical framework

for modelling the metastatic spread of cancer, Bulletin of mathematical biology, 81 (2019), pp. 1965–
2010.

[40] J. M. Gallo, P. Vicini, A. Orlansky, S. Li, F. Zhou, J. Ma, S. Pulfer, M. A. Bookman, and
P. Guo, Pharmacokinetic model-predicted anticancer drug concentrations in human tumors, Clinical
cancer research, 10 (2004), pp. 8048–8058.

[41] A. Ghaffari, B. Bahmaie, and M. Nazari, A mixed radiotherapy and chemotherapy model for
treatment of cancer with metastasis, Mathematical methods in the applied sciences, 39 (2016), pp. 4603–
4617.

[42] R. J. Gillies, D. Verduzco, and R. A. Gatenby, Evolutionary dynamics of carcinogenesis and

why targeted therapy does not work, Nature Reviews Cancer, 12 (2012), pp. 487–493.

[43] J. D. Gordan, J. A. Bertout, C.-J. Hu, J. A. Diehl, and M. C. Simon, HIF-2α promotes
hypoxic cell proliferation by enhancing c-myc transcriptional activity, Cancer cell, 11 (2007), pp. 335–
347.

27
[44] M. M. Gottesman, Mechanisms of cancer drug resistance, Annual review of medicine, 53 (2002),
pp. 615–627.
[45] G. P. Gupta and J. Massagué, Cancer metastasis: building a framework, Cell, 127 (2006), pp. 679–
695.
[46] F. Hamel, F. Lavigne, and L. Roques, Adaptation in a heterogeneous environment I: persistence
versus extinction, Journal of Mathematical Biology, 83 (2021), pp. 1–42.
[47] N. Hartung, S. Mollard, D. Barbolosi, A. Benabdallah, G. Chapuisat, G. Henry, S. Gi-
acometti, A. Iliadis, J. Ciccolini, C. Faivre, et al., Mathematical modeling of tumor growth
and metastatic spreading: validation in tumor-bearing mice, Cancer research, 74 (2014), pp. 6397–6407.
[48] A. Himstedt, C. Braun, S. G. Wicha, and J. M. Borghardt, Towards a quantitative mechanistic
understanding of localized pulmonary tissue retention—A combined in vivo/in silico approach based on
four model drugs, Pharmaceutics, 12 (2020), p. 408.
[49] K. Holt, M. Ye, S. Nagar, and K. Korzekwa, Prediction of tissue-plasma partition coefficients
using microsomal partitioning: Incorporation into physiologically based pharmacokinetic models and
steady-state volume of distribution predictions, Drug metabolism and disposition, 47 (2019), pp. 1050–
1060.
[50] P.-E. Jabin and R. S. Schram, Selection-Mutation dynamics with spatial dependence, Journal de
Mathématiques Pures et Appliquées, 176 (2016).
[51] H. Jones and K. Rowland-Yeo, Basic concepts in physiologically based pharmacokinetic modeling
in drug discovery and development, CPT: pharmacometrics & systems pharmacology, 2 (2013), pp. 1–
12.
[52] G. S. Karagiannis, J. S. Condeelis, and M. H. Oktay, Chemotherapy-induced metastasis: mech-
anisms and translational opportunities, Clinical & experimental metastasis, 35 (2018), pp. 269–284.
[53] A. Kiparissides, S. Kucherenko, A. Mantalaris, and E. Pistikopoulos, Global sensitivity
analysis challenges in biological systems modeling, Industrial & Engineering Chemistry Research, 48
(2009), pp. 7168–7180.
[54] A. W. Lambert, D. R. Pattabiraman, and R. A. Weinberg, Emerging biological principles of
metastasis, Cell, 168 (2017), pp. 670–691.
[55] O. Lavi, J. M. Greene, D. Levy, and M. M. Gottesman, The role of cell density and intratumoral
heterogeneity in multidrug resistance, Cancer research, 73 (2013), pp. 7168–7175.
[56] , Simplifying the complexity of resistance heterogeneity in metastasis, Trends in molecular
medicine, 20 (2014), pp. 129–136.
[57] C. K. Li, The glucose distribution in 9L rat brain multicell tumor spheroids and its effect on cell
necrosis, cancer, 50 (1982), pp. 2066–2073.
[58] S. Lion, M. Boots, and A. Sasaki, Multimorph eco-evolutionary dynamics in structured populations,
The American Naturalist, 200 (2022), pp. 345–372.
[59] G. Long, K. Flaherty, D. Stroyakovskiy, H. Gogas, E. Levchenko, F. De Braud,
J. Larkin, C. Garbe, T. Jouary, A. Hauschild, et al., Dabrafenib plus trametinib versus
dabrafenib monotherapy in patients with metastatic BRAF V600E/K-mutant melanoma: long-term
survival and safety analysis of a phase 3 study, Annals of Oncology, 28 (2017), pp. 1631–1639.
[60] T. Lorenzi, R. H. Chisholm, and J. Clairambault, Tracking the evolution of cancer cell pop-
ulations through the mathematical lens of phenotype-structured equations, Biology direct, 11 (2016),
pp. 1–17.
[61] T. Lorenzi, C. Venkataraman, A. Lorz, and M. A. Chaplain, The role of spatial variations of
abiotic factors in mediating intratumour phenotypic heterogeneity, Journal of theoretical biology, 451
(2018), pp. 101–110.

28
[62] A. Lorz, S. Mirrahimi, and B. Perthame, Dirac mass dynamics in multidimensional nonlocal
parabolic equations, Communications in Partial Differential Equations, 36 (2011), pp. 1071–1098.

[63] S. A. Luebker and S. A. Koepsell, Diverse mechanisms of BRAF inhibitor resistance in melanoma
identified in clinical and preclinical studies, Frontiers in oncology, 9 (2019), p. 268.

[64] S. Martin-Algarra, R. Hinshelwood, S. Mesnage, J. Cebon, P. F. Ferrucci, M. Agli-

etta, B. Neyns, V. Chiarion-Sileni, C. R. Lindsay, M. Del Vecchio, et al., Effectiveness of
dabrafenib in the treatment of patients with BRAF V600–mutated metastatic melanoma in a Named
Patient Program, Melanoma research, 29 (2019), pp. 527–532.

[65] J. Massagué and A. C. Obenauf, Metastatic colonization by circulating tumour cells, Nature, 529
(2016), pp. 298–306.

[66] F. Michor and Y. Iwasa, Dynamics of metastasis suppressor gene inactivation, Journal of theoretical
biology, 241 (2006), pp. 676–689.

[67] F. Michor, M. A. Nowak, and Y. Iwasa, Stochastic dynamics of metastasis formation, Journal of
Theoretical Biology, 240 (2006), pp. 521–530.

[68] S. Mirrahimi, Adaptation and migration of a population between patches, Discrete and Continuous
Dynamical Systems-B, 18 (2012), pp. 753–768.

[69] S. Mirrahimi and S. Gandon, Evolution of specialization in heterogeneous environments: equilib-

rium between selection, mutation and migration, Genetics, 214 (2020), pp. 479–491.

[70] S. Mirrahimi and B. Perthame, Asymptotic analysis of a selection model with space, Journal de
Mathématiques Pures et Appliquées, 104 (2015), pp. 1108–1118.

[71] D. M. Moss and M. Siccardi, Optimizing nanomedicine pharmacokinetics using physiologically

based pharmacokinetics modelling, British journal of pharmacology, 171 (2014), pp. 3963–3979.

[72] P. K. Newton, J. Mason, K. Bethel, L. Bazhenova, J. Nieva, L. Norton, and P. Kuhn,

Spreaders and sponges define metastasis in lung cancer: a markov chain monte carlo mathematical
model, Cancer research, 73 (2013), pp. 2760–2769.

[73] P. K. Newton, J. Mason, N. Venkatappa, M. S. Jochelson, B. Hurt, J. Nieva, E. Comen,

L. Norton, and P. Kuhn, Spatiotemporal progression of metastatic breast cancer: a markov chain
model highlighting the role of early metastatic sites, NPJ breast cancer, 1 (2015), pp. 1–9.

[74] D. X. Nguyen and J. Massagué, Genetic determinants of cancer metastasis, Nature Reviews
Genetics, 8 (2007), pp. 341–352.

[75] E. Norris, J. R. King, and H. M. Byrne, Modelling the response of spatially structured tumours
to chemotherapy: Drug kinetics, Mathematical and Computer Modelling, 43 (2006), pp. 820–837.

[76] A. C. Obenauf, Y. Zou, A. L. Ji, S. Vanharanta, W. Shu, H. Shi, X. Kong, M. C. Bosen-

berg, T. Wiesner, N. Rosen, et al., Therapy-induced tumour secretomes promote resistance and
tumour progression, Nature, 520 (2015), pp. 368–372.

[77] B. Ocaña-Tienda and V. M. Pérez-Garcı́a, Mathematical modeling of brain metastases growth

and response to therapies: A review, Mathematical Biosciences, (2024), p. 109207.

[78] A. Olivier and C. Pouchol, Combination of direct methods and homotopy in numerical optimal
control: application to the optimization of chemotherapy in cancer, Journal of Optimization Theory
and Applications, 181 (2019), pp. 479–503.

[79] K. Pantel and M. Speicher, The biology of circulating tumor cells, Oncogene, 35 (2016), pp. 1216–
1224.

29
[80] A. L. Parker, M. Benguigui, J. Fornetti, E. Goddard, S. Lucotti, J. Insua-Rodrı́guez,
A. P. Wiegmans, and E. C. L. C. of the Metastasis Research Society, Current challenges in
metastasis research and future innovation for clinical translation, Clinical & Experimental Metastasis,
39 (2022), pp. 263–277.

[81] B. Perthame, Transport equations in biology, Springer Science & Business Media, 2006.

[82] B. Perthame and G. Barles, Dirac concentrations in lotka-volterra parabolic pdes, Indiana Univer-
sity Mathematics Journal, (2008), pp. 3275–3301.

[83] C. Pouchol, J. Clairambault, A. Lorz, and E. Trélat, Asymptotic analysis and optimal control
of an integro-differential system modelling healthy and cancer cells exposed to chemotherapy, Journal
de Mathématiques Pures et Appliquées, 116 (2018), pp. 268–308.

[84] A. Puszkiel, G. Noé, A. Bellesoeur, N. Kramkimel, M.-N. Paludetto, A. Thomas-

Schoemann, M. Vidal, F. Goldwasser, E. Chatelut, and B. Blanchet, Clinical pharma-
cokinetics and pharmacodynamics of dabrafenib, Clinical Pharmacokinetics, 58 (2019), pp. 451–467.

[85] G. Qian and A. Mahdi, Sensitivity analysis methods in the biomedical sciences, Mathematical bio-
sciences, 323 (2020), p. 108306.

[86] I. Ramis-Conde, M. A. Chaplain, A. R. Anderson, and D. Drasdo, Multi-scale modelling of

cancer cell intravasation: the role of cadherins in metastasis, Physical biology, 6 (2009), p. 016008.

[87] A. Saltelli, Making best use of model evaluations to compute sensitivity indices, Computer physics
communications, 145 (2002), pp. 280–297.

[88] P. Schlicke, C. Kuttler, and C. Schumann, How mathematical modeling could contribute to the
quantification of metastatic tumor burden under therapy: insights in immunotherapeutic treatment of
non-small cell lung cancer, Theoretical Biology and Medical Modelling, 18 (2021), p. 11.

[89] J. G. Scott, D. Basanta, A. R. Anderson, and P. Gerlee, A mathematical model of tumour

self-seeding reveals secondary metastatic deposits as drivers of primary tumour growth, Journal of The
Royal Society Interface, 10 (2013), p. 20130011.

[90] J. G. Scott, P. Gerlee, D. Basanta, A. G. Fletcher, P. K. Maini, and A. R. Anderson,

Mathematical modeling of the metastatic process, Experimental Metastasis: Modeling and Analysis,
(2013), pp. 189–208.

[91] R. Sharma and S. Sharma, Physiology, Blood Volume., StatPearls, StatPearls Publishing, 2022.

[92] S. V. Sharma, D. Y. Lee, B. Li, M. P. Quinlan, F. Takahashi, S. Maheswaran, U. McDer-

mott, N. Azizian, L. Zou, M. A. Fischbach, et al., A chromatin-mediated reversible drug-tolerant
state in cancer cell subpopulations, Cell, 141 (2010), pp. 69–80.

[93] T. Shibue and R. A. Weinberg, EMT, CSCs, and drug resistance: the mechanistic link and clinical
implications, Nature reviews Clinical oncology, 14 (2017), pp. 611–629.

[94] X. Sun, J. Bao, and Y. Shao, Mathematical modeling of therapy-induced cancer drug resistance:
connecting cancer mechanisms to population survival rates, Scientific reports, 6 (2016), p. 22498.

[95] S. Sundararajan, A. M. Thida, S. Yadlapati, and S. Koya, Metastatic Melanoma, StatPearls,

2022.

[96] T. Tsuruo and I. J. Fidler, Differences in drug sensitivity among tumor cells from parental tumors,
selected variants, and spontaneous metastases, Cancer Research, 41 (1981), pp. 3058–3064.

[97] N. C. Turner and J. S. Reis-Filho, Genetic heterogeneity and cancer drug resistance, The lancet
oncology, 13 (2012), pp. e178–e185.

[98] US Food and Drug Administration, Dabrafenib (Tafinlar) Pharmacology Review. https://www.
accessdata.fda.gov/drugsatfda_docs/label/2022/202806s022lbl.pdf, 2013.

30
 [99] C. Villa, M. A. Chaplain, and T. Lorenzi, Evolutionary dynamics in vascularised tumours un-
der chemotherapy: Mathematical modelling, asymptotic analysis and numerical simulations, Vietnam
Journal of Mathematics, 49 (2021), pp. 143–167.
[100] N. Vitos and P. Gerlee, Model-based inference of metastatic seeding rates in de novo metastatic
breast cancer reveals the impact of secondary seeding and molecular subtype, Scientific Reports, 12
(2022), p. 9455.
[101] X. Wang, H. Zhang, and X. Chen, Drug resistance and combating drug resistance in cancer, Cancer
drug resistance, 2 (2019), p. 141.
[102] J. P. Ward and J. R. King, Mathematical modelling of avascular-tumour growth, Mathematical
Medicine and Biology: A Journal of the IMA, 14 (1997), pp. 39–69.
[103] World Health Organization, Global Cancer Observatory. https://gco.iarc.fr/, 2022.
[104] A. Wu, K. Loutherback, G. Lambert, L. Estévez-Salmerón, T. D. Tlsty, R. H. Austin,
and J. C. Sturm, Cell motility and drug gradients in the emergence of resistance to chemotherapy,
Proceedings of the National Academy of Sciences, 110 (2013), pp. 16103–16108.
[105] M. Yancovitz, A. Litterman, J. Yoon, E. Ng, R. L. Shapiro, R. S. Berman, A. C. Pavlick,
F. Darvishian, P. Christos, M. Mazumdar, et al., Intra-and inter-tumor heterogeneity of
BRAFV600E mutations in primary and metastatic melanoma, PloS one, 7 (2012), p. e29336.
[106] A. Yin, D. J. A. Moes, J. G. van Hasselt, J. J. Swen, and H.-J. Guchelaar, A review of
mathematical models for tumor dynamics and treatment resistance evolution of solid tumors, CPT:
pharmacometrics & systems pharmacology, 8 (2019), pp. 720–737.
[107] Q. Zhou, P. Guo, G. D. Kruh, P. Vicini, X. Wang, and J. M. Gallo, Predicting human
tumor drug concentrations from a preclinical pharmacokinetic model of temozolomide brain disposition,
Clinical Cancer Research, 13 (2007), pp. 4271–4279.

[108] A. Ziemys, K. Yokoi, M. Kai, Y. Liu, M. Kojic, V. Simic, M. Milosevic, A. Holder, and
M. Ferrari, Progression-dependent transport heterogeneity of breast cancer liver metastases as a
factor in therapeutic resistance, Journal of controlled release, 291 (2018), pp. 99–105.
[109] H. Zou, P. Banerjee, S. S. Y. Leung, and X. Yan, Application of pharmacokinetic-
pharmacodynamic modeling in drug delivery: development and challenges, Frontiers in pharmacology,
11 (2020), p. 543082.

31
A Proofs of analytical results
A.1 Proof of the bounds on Iiε
We prove the result by contradiction. Suppose that 0 < IM < I1ε . We integrate (24)1 and adopt the
boundary conditions (24)3 to obtain
Z 1 Z 1 Z 1
R1 (y, I1ε )n1ε (y)dy + ν2,1 (y)n1ε (y)dy − ν1,2 (y)n1ε (y)dy = 0 . (A.1)
0 0 0

For the first term of equation (A.1) we employ assumptions (15) and (22), along with the inequality IM < I1ε ,
to deduce


R1 (y, I1ε ) < R1 (y, IM ) ≤ max R1 (y, IM ) ≤ −δ .
y∈[0,1]

Multiplying both sides of this inequality by n1ε (y) and integrating, we obtain
Z 1 Z 1
R1 (y, I1ε )n1ε (y)dy < −δn1ε (y)dy = −δI1ε .
0 0

For what concerns the second and third terms of equation (A.1), from assumption (17) we can deduce
Z 1 Z 1
ν2,1 (y)n2ε (y)dy ≤ νM I2ε and − ν1,2 (y)n1ε (y)dy ≤ −νm I1ε .
0 0

Combining these results we obtain

Z 1 Z 1 Z 1
0= R1 (y, I1ε )n1ε (y)dy + ν2,1 (y)n1ε (y)dy − ν1,2 (y)n1ε (y)dy < −δI1ε + νM I2ε − νm I1ε ,
0 0 0

and thus
(νm + δ) νm νm
I2ε > I1ε > I1ε > IM . (A.2)
νM νM νM
Now we integrate equation (24)1 , both for i = 1 and i = 2, apply boundary conditions (24)3 , and then sum
them to obtain
Z 1 Z 1
R1 (y, I1ε )n1ε (y)dy + R2 (y, I2ε )n2ε (y)dy = 0. (A.3)
0 0

From assumption (15), inequality (A.2) and the technical assumption (22), we have that
   
νm

R1 (y, I1ε ) < max R1 (y, IM ) ≤ −δ and R2 (y, I2ε ) < max R2 y, IM ≤ −δ .
y∈[0,1] y∈[0,1] νM

Together with equation (A.3) and, again, the previous result (A.2), this gives
Z 1 Z 1
νm

0= R1 (y, I1ε )n1ε (y)dy + R2 (y, I2ε )n2ε (y)dy < −δ(I1ε + I2ε ) < −δ IM + IM < 0 ,
0 0 νM

a contradiction. We thus conclude that I1ε ≤ IM . The proof for the case with 0 < IM < I2ε , and for the
lower bounds of Iiε (i = 1, 2) follow analogous arguments, relying on assumption (21).

A.2 Proof of the results for the asymptotic regime ε → 0

We formally extend the analysis carried out in [68] for the case of constant migration rates to the case of
phenotype-dependent migration rates.

32
Hopf-Cole transformation. Following the work of [68], we introduce the Hopf-Cole transformation
niε (y) = euiε (y)/ε [5, 36, 38], with uiε (y) semi-convex (i.e. ∂yy2
uiε ≥ −E, for some constant E > 0), for
i = 1, 2. Then we can evaluate
d2 niε  2  2 
−1 d uiε −2 duiε
= ε + ε niε ,
dy 2 dy 2 dy
and equation (24)1 becomes

 d2 uiε  du 2 
iε
Ri y, Iiε niε + ε2 ε−1 2
+ ε−2 niε + νji (y) njε − νij (y) niε = 0, y ∈ (0, 1).
dy dy
We can rewrite the system (24)1 for niε (i = 1, 2) as the following steady-state system for uiε (i = 1, 2) in
matrix-vector form:
Aε Nε = Lε Nε (A.4)
where


R1 y, I1ε − ν12 (y) ν21
(y)
Aε = , (A.5)
ν12 (y) R2 y, I2ε − ν21 (y)
   u /ε 
n1ε e 1ε
Nε = = u2ε /ε , (A.6)
n2ε e
 d2 u  2  2
1ε du1ε d2 u2ε du2ε 
Lε = diag −ε − , −ε − . (A.7)
dy 2 dy dy 2 dy
Then given the eigenvalue Hε = Hε (y, I1ε , I2ε ) of Aε with corresponding eigenvector Nε , from (A.4) we have
Lε Nε = Hε Nε , (A.8)
which corresponds to the following equations for uiε
d2 uiε  duiε 2
−ε − = Hε (y, I1ε , I2ε ) i = 1, 2 . (A.9)
dy 2 dy

Hamilton-Jacobi equation. Firstly, given the assumptions introduced in Section 3.1, from (25) and (26)
we deduce that n1ε and n2ε converge weakly to measures n1 and n2 , i.e.
∗
niε (y) −−−⇀ ni (y) i = 1, 2 . (A.10)
ε→0

On the basis of the analysis carried out in [68] for the case of constant migration rates, we expect that in the
asymptotic regime ε → 0 both sequences u1ε and u2ε converge uniformly in [0, 1] to a continuous function
u ∈ C([0, 1]) and (I1ε , I2ε ) converges to (I1 , I2 ). We further assume that the semi-convexity of uiε (i = 1, 2)
is preserved in the limit such that u = u(y) is also semi-convex. If that is the case, then in the limit ε → 0
equation (A.4) yields
AN = LN (A.11)
where


R1 y, I1 − ν12 (y) ν21
(y)
A= , (A.12)
ν12 (y) R2 y, I2 − ν21 (y)
 
n1
N = , (A.13)
n2
  du 2  2 
du
L = diag − ,− . (A.14)
dy dy
Then as ε → 0 equation (A.9) and the bounds (26) lead to the following constrained Hamilton-Jacobi
equation for u(y):
  2
− du

= H(y, I1 , I2 )
dy , (A.15)

 max u(y) = 0
y∈[0,1]

33
where the constraint is required to ensure bounds on Ii (i = 1, 2) in (26) are satisfied. Note that the right-
hand-side of (A.15)1 is zero at points yk ∈ arg max u(y) and negative for all y ∈ [0, 1] that are not stationary
points of u(y). This implies

H(yk , I1 , I2 ) = 0 ∀ yk ∈ arg max u(y),
(A.16)
 max H(y, I1 , I2 ) = 0,
y∈[0,1]

where H = H(y, I1 , I2 ) is the largest eigenvalue of A with corresponding eigenvector N . Hence H is given
by p
H(y, I1 , I2 ) = F + F 2 − 4G , (A.17)
with
F = R1 (y, I1 ) − ν12 (y) + R2 (y, I2 ) − ν21 (y) , (A.18)



G = R1 (y, I1 ) − ν12 (y) R2 (y, I2 ) − ν21 (y) − ν12 (y)ν21 (y) . (A.19)
Note from (A.17) that for (A.16)1 to be satisfied we must have that F is negative at each yk ∈ arg max u(y).
In such case, (A.16) implies

G(yk , I1 , I2 ) = 0 ∀ yk ∈ arg max u(y)
(A.20)
 min G(y, I1 , I2 ) = 0 .
y∈[0,1]

Then from (A.15) and (A.16) we have

supp ni ⊂ Γ i = 1, 2 , (A.21)
where 
Γ := arg max H(y) ≡ y ∈ [0, 1] : H(y, I1 , I2 ) = 0 , (A.22)
y∈[0,1]

which – under assumption (19) – we expect to be discrete and finite, with cardinality K ∈ N.

Remark 1: As we have shown that (A.16), together with definitions (A.17) and (A.18), implies (A.20), we
could alternatively define Γ in (A.21) as

Γ := arg min G(y) ≡ y ∈ [0, 1] : G(y, I1 , I2 ) = 0 .
y∈[0,1]

Remark 2: The support of ni is only a proper subset of Γ as defined in (A.22), since Γ includes both the
points at which u attains its maxima and minima (from (A.15)1 ), while ni will be zero at the points in
arg min u (from (A.15)2 ). Hence (A.21) may be more precisely written as


supp ni ⊆ Ω ⊂ Γ i = 1, 2 ,
where Γ is defined as in (A.22) (or as in Remark 1) and Ω is defined as

Ω := arg max u(y) ≡ y ∈ [0, 1] : u(y) = 0 .
y∈[0,1]

Concentration as Dirac masses. From (25), we deduce that, along subsequences and for i = 1, 2, niε
converges weakly in the asymptotic regime ε → 0 to a measure ni . In view of the results so far obtained on
the support of ni , we have that the measures n1 and n2 to which n1ε and n2ε converge concentrate as Dirac
masses, i.e.
K
X
∗
niε (y) −−−⇀ ρik δ(y − yk ) i = 1, 2 , (A.23)
ε→0
k=1

where the points yk ∈ [0, 1] can be found by solving (A.16) – or alternatively (A.20) – while the weights
ρik ≥ 0 must be such that
XK
Ii = ρik i = 1, 2 , (A.24)
k=1

34
 
and can be found by integrating the system (A.11) over a ball Br (yk ) := y ∈ [0, 1] : |y − yk | ≤ r centred
at yk , with radius r small enough that no other point in the support of ni (i = 1, 2) is contained in Br (yk ).
This yields 

 R1 (yk , I1 ) − ν12 (yk ) ρ1k + ν21 (yk )ρ2k = 0
k = 1, .., K (A.25)
 R (y , I ) − ν (y )
ρ + ν (y )ρ = 0.
2 k 2 21 k 2k 12 k 1k

A.3 Proof the analytical results for the metastatic spread case
If ν12 (y) ≥ 0, with ν12 (y) > 0 for some y ∈ [0, 1], and ν21 (y) ≡ 0 for all y ∈ [0, 1], then solving (A.16) leads
to the following systems:
 
R1 (y1 , I1 ) − ν12 (y1 ) = 0 R2 (y2 , I2 ) = 0
and (A.26)
∂ R (y , I ) − ∂ ν (y ) = 0 ∂ R (y , I ) = 0 .
y 1 1 1 y 12 1 y 2 2 2

Notice that (A.26) together with assumption (20) ensure that F defined in (A.18) is negative at y1 and y2 ,
as required by (A.16) and (A.17). From the second equation of each system in (A.26) we obtain y1 and y2 ,
while from the first one of each system we get I1 and I2 . Moreover, equations (A.24) and (A.25) result in
 
ν12 (y1 )
ρ11 = I1 , ρ12 = 0 ρ21 = min − ρ11 , I2 , ρ22 = I2 − ρ21 . (A.27)
R2 (y1 , I2 )

Solving (A.26) for Ri (y, Ii ) (i = 1, 2) given by (7) and for ν12 (y) given by (9), we have that (A.10) and (A.23)
corresponds to the following asymptotic solution:

n1 (y) = I1 δ(y − y1 ) and n2 (y) = ρ21 δ(y − y1 ) + ρ22 δ(y − y2 ) , (A.28)

with
b1
y1 = h1 , y2 = h2 , (A.29)
b1 + ν̂12
and
   
1 ν̂12 b1 h21 a2 ν̂12 y12
I1 = a1 − , I2 = , ρ21 = min I1 , I2 , ρ22 = I2 − ρ21 . (A.30)
d1 b1 + ν̂12 d2 b2 (y1 − h2 )2

35
 The development of drug resistance in metastatic tumours under
chemotherapy: an evolutionary perspective
Supplementary Material
Federica Padovano∗ Luis Almeida Chiara Villa†
May 2024
arXiv:2405.20203v1 [q-bio.CB] 30 May 2024

S.1 Supplementary analytical results

The following analytical results extend those in Section 3.3 of the Main Manuscript to the case of localised
tumour and secondary seeding. Therefore we refer to equations appearing in Section 3.3 and Appendix A.2
of the Main Manuscript.

S.1.1 Analytical results for the localised tumors case: ν1,2 ≡ 0 and ν2,1 ≡ 0
If both ν12 (y) ≡ 0 and ν21 (y) ≡ 0 for all y ∈ [0, 1], then solving (A.16) leads to the following systems:

Ri (yi , Ii ) = 0
i = 1, 2 . (S.1)
∂ R (y , I ) = 0
y i i i

Notice that (S.1) together with assumption (20) ensure that F defined in (A.18) is negative at y1 and y2 , as
required by (A.16) and (A.17). From (S.1)2 we obtain y1 and y2 , and from (S.1)1 we get I1 and I2 , while
(A.25) leads to having ρ21 = ρ12 = 0 which together with (A.24) imply ρ11 = I1 and ρ22 = I2 . Doing this
for Ri (y, Ii ) (i = 1, 2) given by (7) leads to the following solution:
ai
ni (y) = Ii δ(y − yi ) with yi = hi , Ii = for i = 1, 2 . (S.2)
di

These results, in agreement with previous results in the literature focusing on a single localised tumour,
provide a formal description of the fact that populations living in separate, non-communicating sites evolve
independently of each other. Each population presents monomorphism at equilibrium, and the trait selected
corresponds to the fittest trait hi (i = 1, 2) determined by the local environment. That is, of course, if cells
were present in both sites initially: if the secondary site is representing a metastatic site, in this scenario no
cells ever reach it and we expect n2 ≡ 0 at equilibrium.

S.1.2 Analytical results for the secondary seeding case: ν1,2 ̸≡ 0 and ν2,1 ̸≡ 0
Consider now the case in which ν12 ̸≡ 0 and ν21 ̸≡ 0. In particular we focus on the case in which the selected
traits yk are such that νij (yk ) > 0 (i, j = 1, 2 i ̸= j). In this case, for each k = 1, .., K we have from (A.25)1
that ρ2k > 0 implies ρ1k > 0 and R1 (yk , I1 ) − ν12 (yk ) < 0. Similarly, from (A.25)2 we have that ρ1k > 0
implies ρ2k > 0 and R2 (yk , I2 ) − ν21 (yk ) < 0. Therefore we either have both ρ1k = ρ2k = 0 (trivial case) or
ρ1k > 0, ρ2k > 0 and
Ri (yk , Ii ) − νij (yk ) < 0 , i, j = 1, 2 i ̸= j . (S.3)
∗ Corresponding author address: federica.padovano@sorbonne-universite.fr
† Corresponding author address: chiara.villa1@sorbonne-universite.fr

1
Under definition (A.18), (S.3) ensures that F is negative at every yk (k = 1, .., K) as required by (A.16) and
(A.17). In addition, (S.3) and (A.25) imply
!
!
ν21 (yk ) ν21 (yk ) − R2 yk , I2
ρ1k =
ρ2k = ρ2k for k = 1, ..., K . (S.4)
ν12 (yk ) − R1 yk , I1 ν12 (yk )

For a more precise characterisation of the phenotypic distributions at equilibrium, we must solve (A.20),
2
which is equivalent to simultaneously solving G = 0, ∂y G = 0 and ∂yy G ≥ 0. Due to analytical complexity,
we here do not solve this explicitly. Instead we investigate the number K of discrete points that we might find
solving (A.20), under the specific fitness functions and transition rates, to shed light on the conditions under
which we may expect the population to go extinct, present monomorphism or polymorphism at equilibrium.
Let Ri (y, Ii ) (i = 1, 2) be given by (7) and νij (y) (i, j = 1, 2 i ̸= j) be given by (9). Note first of all
2
that, under these specific definitions, G(y, I1 , I2 ) in (A.19) is a polynomial of order 4 in y. Therefore ∂yy G
2
is quadratic in y, suggesting that G can have at most 2 inflection points, determined by solving ∂yy G = 0 in
y. Figure S.1 illustrates why in order for the population to be dimorphic we must have 2 inflections points:
this will give us a necessary condition for dimorphism and, in the opposite case (i.e. if we have one or no
inflection points), a sufficient condition for monomorphism – always assuming that the population does not
2
go extinct. Said condition is determined by applying the well-known binomial formula to ∂yy G = 0 and
investigating its discriminant D(I1 , I2 ). We here omit simple calculations to derive this and focus on the
obtained necessary condition for dimorphism, which is given by

D(I1 , I2 ) > 0 . (S.5)

Figure S.1: Four different cases for the function G. Possible plots of G(y, I1 , I2 ) as a function of y
displaying how G having two inflection points is necessary – cf. (a) and (b) – but not sufficient – vid (c) and
(d) – for the population to present dimorphism at equilibrium.

Under definitions (7) and (9), condition (S.5) can be rewritten to give the following condition:
 
a1 − d1 I1 a2 − d2 I2
(B1 h1 − B2 h2 )2 + 8B1 BV1 V2 h1 h2 − 2B1 B2 V1 V2 (h21 + h22 ) + 2(1 − V1 V2 ) + > 0,
b1 + ν̂12 b2 + ν̂21

bi ν̂ij
with Bi := , Vi := i, j = 1, 2 i ̸= j .
bi + ν̂ij bi + ν̂ij

While this condition is quite complex, we may still consider the limit as b1 , b2 → ∞, under which we have
Bi → 1, Vi → 0 (i = 1, 2) and therefore the condition reduces to (h1 − h2 )2 > 0, which is satisfied as long
as h1 ̸= h2 . Therefore under extremely strong selective pressure from the two environments we can – and
we might expect to – have dimorphism in both sites, which is in line with the results obtained in the case of
constant transition rates. Figure S.2 displays the possible model outcomes changing some of the parameters,
and obtained simulating a nondimensional version of the model.

2
S.2 Supplementary figures

Figure S.2: Outcome dependency on input factors in the secondary seeding scenario. Illustrative
example showing how the selection gradients b1 and b2 , the fittest trait h1 and the maximum migration
rates ν̂1,2 and ν̂2,1 may affect the equilibrium distributions of the cancer cell populations of each site,
in the secondary seeding case. We simulate system (3), under definition (7) for the fitness functions Ri
and definition (9) for the migration rates νij (y), under the initial conditions n0,1 (y) > 0 and n0,2 (y) = 0
∀y ∈ [0, 1]. Each column displays the equilibrium solutions of three simulations obtained by progressively
varying the parameters b1 and b2 (first column), h1 (second column) and ν̂1,2 and ν̂1,,12 (third column) form
their baseline values bB B B B B
2 = 1, b2 = 0.6, h1 = 0.2, ν̂1,2 = 0.1 and ν̂2,1 = 0.05. The remaining parameters are
set to β1 = β2 = 10−7 , a1 = 6, a2 = 5, b1 = bB B B B
1 , b2 = b2 , h1 = h1 , h2 = 0.6, d1 = d2 = 0.2, ν̂1,2 = ν̂1,2 and
B
ν̂2,1 = ν̂2,1 .

3
 (a)

(b)

Figure S.3: Results of the Elementary effects method for GSA: EE measures comparison. Com-
∗
parison of the EE measures EE i (x-axis) and SDi (y-axis), defined in (37), associated with the model output
YI (a) and Yµ (b). The estimates are computed with r = 500 evaluations of the elementary effect for each
parameter.

4
 (a.1) (b.1)

(a.2) (b.2)

(a.3) (b.3)

Figure S.4: Results of the Elementary effects method for GSA under a wider range of values for
the maximum migration rates. EE measures associated with the model output YI (a.1-3) and Yµ (b.1-3),
i.e. the total tumour mass and average mean phenotypic state of the sites. The estimates are computed
with r = 500 evaluations of the elementary effect for each parameter, and migration rates taking values in
the range ν1,2 ∈ [10−13 , 10−5 ]. The elementary effects EEi , defined in (36), calculated for each parameter
∗ ∗
are displayed in plots (a.1) and (b.1). The associated measures EE i , |EE i | and SDi , defined in (37), are
∗
shown in plots (a.2) and (b.2). Plots (a.3) and (b.3) present the ratio between SD and EE i .

5
 (a.1) (b.1)

(a.2) (b.2)

(a.3) (b.3)

Figure S.5: Results of the Elementary effects method for GSA under a lower drug dose. EE
measures associated with the model output YI (a.1-3) and Yµ (b.1-3), i.e. the total tumour mass and average
mean phenotypic state of the sites. The estimates are computed with r = 500 evaluations of the elementary
effect for each parameter, and an intravenously-injected drug dose of 3.75·10−4 µg/s. The elementary effects
EEi , defined in (36), calculated for each parameter are displayed in plots (a.1) and (b.1). The associated
∗ ∗
measures EE i , |EE i | and SDi , defined in (37), are shown in plots (a.1) and (b.2). Plots (a.3) and (b.3)
∗
present the ratio between SD and EE i .

6
 (a) Baseline scenario (b) Non baseline scenario

Figure S.6: Cancer evolutionary dynamics under baseline and non-baseline scenarios: zoom-in
analysis of the last two days of the simulations. Cancer evolutionary dynamics results of the numerical
simulations under the baseline (a) and non-baseline (b) scenarios with the drug schedule of 150 mg orally
administered twice a day, and time period t ∈ [103, 105] days, i.e. the last 26301 time steps of the simulation
results shown in Figure 7 of the Main Manuscript. For each scenario we plot the population size (top left
panel), mean phenotypic trait (top right panel), the related variance of the phenotypic distribution (bottom
left panel) and the step difference dif fi,k (bottom right panel), respectively defined in equations (1), (2) and
(35). More details on the simulation set-up and numerical methods can be found in Section 4.1.

Figure S.7: Variance dependency on βi . Treatment outcomes obtained from numerical simulations of
the model under the baseline scenario, for an intravenously-injected drug dose of 2.6915 µg/s and a final
time of T = 210 days. The graph displays the variance σi2 , defined in equation (2), of the phenotypic
distribution ni at steady state, varying the epimutation rate βi . The steady state time is the first time at
which dif fi,k < tol, with dif fi,k defined in (35) and a tolerance of tol = 10−6 . More details on the simulation
set-up and numerical methods can be found in Section 4.1.