Anni S. H
alkola
A
I 674
A
N
N
A
LES U
N
IV
ERSITATIS TU
RK
U
EN
SIS
ISBN 978-951-29-9012-2 (PRINT)
ISBN 978-951-29-9013-9 (PDF)
ISSN 0082-7002 (PRINT)
ISSN 2343-3175 (ONLINE)
Pa
in
os
al
am
a 
O
y, 
Tu
rk
u,
 F
in
la
nd
 2
02
2
TURUN YLIOPISTON JULKAISUJA – ANNALES UNIVERSITATIS TURKUENSIS
SARJA – SER. AI OSA – TOM. 674 | ASTRONOMICA – CHEMICA – PHYSICA – MATHEMATICA | TURKU 2022
MATHEMATICAL 
MODELLING AND SURVIVAL 
PREDICTION IN CANCER
Anni S. Halkola
MATHEMATICAL MODELLING
AND SURVIVAL PREDICTION IN
CANCER
Anni S. Halkola
TURUN YLIOPISTON JULKAISUJA – ANNALES UNIVERSITATIS TURKUENSIS
SARJA – SER. AI OSA – TOM. 674 | ASTRONOMICA – CHEMICA – PHYSICA – MATHEMATICA | TURKU 2022
University of Turku
Faculty of Science
Department of Mathematics and Statistics
Applied mathematics
Doctoral Programme in Exact Sciences (EXACTUS)
Supervised by
Associate prof. Kalle Parvinen, PhD
Department of Mathematics and Statis-
tics
University of Turku
Turku, Finland
Advancing Systems Analysis Program
International Institute for Applied Sys-
tems Analysis (IIASA)
Laxenburg, Austria
Prof. Tero Aittokallio, PhD
Institute for Molecular Medicine Finland
(FIMM), University of Helsinki
Helsinki, Finland
Institute for Cancer Research
Oslo University Hospital
Oslo, Norway
Oslo Centre for Biostatistics and Epi-
demiology (OCBE), University of Oslo
Oslo, Norway
Reviewed by
Associate prof. Harri La¨hdesma¨ki, D.Sc.
(Tech)
Department of Computer Science
Aalto University
Espoo, Finland
Research Associate Alvaro Ko¨hn-
Luque, PhD
Institute of Basic Medical Sciences
University of Oslo
Oslo, Norway
Opponent
Prof. Joel S. Brown, PhD
Moffitt cancer center
Tampa, FL, USA
University of Illinois
Chicago, USA
The originality of this publication has been checked in accordance with the University
of Turku quality assurance system using the Turnitin OriginalityCheck service.
ISBN 978-951-29-9012-2 (PRINT)
ISBN 978-951-29-9013-9 (PDF)
ISSN 0082-7002 (PRINT)
ISSN 2343-3175 (ONLINE)
Painosalama Oy, Turku, Finland, 2022
iii
To anyone affected by cancer
(so basically everyone)
UNIVERSITY OF TURKU
Faculty of Science
Department of Mathematics and Statistics
Applied mathematics
HALKOLA, ANNI S.: Mathematical modelling and survival prediction in can-
cer
Doctoral dissertation, 144 pp.
Doctoral Programme in Exact Sciences (EXACTUS)
November 2022
ABSTRACT
Cancer is one of the leading causes of death, thus opening a vast need for exten-
sive research and insights. The survival prospects, along with treatment benefits and
costs (economical or health-related), can be predicted with tools from mathematical
modelling and regression analysis. Promising results have been gained on suggest-
ing mono- and combination therapies, that could potentially improve the treatment
strategies. Furthermore, multiple biological features have been recognized as impor-
tant predictors of treatment outcome. However, since cancer remains a challenging
and unpredictable enemy, the need for more effective and personalized predictions
and treatment suggestions remains.
In this dissertation, various modelling approaches were used to predict cancer
behavior, treatment outcomes and patient survival. An ordinary differential equa-
tion model was developed to investigate the changes in the cancer cell density as
different treatment regimen were applied. In addition, we included the immune sys-
tem along with immunotherapy, since the immune response is an important part of
cancer development and has a potential to eradicate tumors. It was noted, that an
adaptive treatment resulted in lower cancer burden and less time in treatment. In
addition, combination treatments (immunotherapy with either chemo- or targeted
therapy) generally resulted in smaller cancer burden than monotherapies, however,
the potential additional side effects of two therapies have to be considered.
A metapopulation model was developed for the cancer development, in which
the focus was on emergence of angiogenesis and cancer cell emigration. We inves-
tigated, in which conditions cancer cells would become angiogenic with or with-
out treatments (anti-angiogenic, cytotoxic or combination). In general, angiogenesis
contribution was desired quality for cancer cells, if no anti-angiogenic treatment was
administrated. With anti-angiogenic treatment, angiogenesis diminished, however
the risk of resistance against anti-angiogenic treatment also increased.
Two new regression methods were developed with focus on survival prediction.
A greedy budget-constrained Cox regression (Greedy Cox) utilizes 𝐿2-penalty and
considers the cost of selected parameters. It was also compared to LASSO selec-
tion (𝐿1). Optimal Subset CArdinality Regression (OSCAR) method was developed
with 𝐿0-pseudonorm penalty to provide sparse models. The costs of measuring the
selected model features were also considered in comparison to prediction accuracy.
The methods were validated on clinical prostate cancer data and it was noted that a
comparable level of prediction accuracy was already reached with a few parameters,
v
Anni S. Halkola
resulting in relatively low costs. All of the investigated methods also selected rea-
sonable, cancer-related parameters such as prostate specific antigen (PSA).
Taken together, this dissertation provides a comprehensive research of novel tools
for modelling and predicting cancer behavior and patient survival. Important hall-
marks of cancer development, such as immune response and angiogenic switch have
been included along with corresponding treatments that have potential to change the
traditional treatment regimens.
KEYWORDS: cancer, killer T-cells, mathematical modelling, personalized medicine,
immunotherapy, anti-angiogenic treatment, cytotoxic treatment, targeted treatment,
combination therapy, side-effects, tumor microenvironment, metapopulation mod-
elling, cost optimization, survival prediction, Cox regression, prostate cancer, feature
selection, 𝐿0-pseudonorm
vi
TURUN YLIOPISTO
Matemaattis-luonnontieteellinen tiedekunta
Matematiikan ja tilastotieteen laitos
Sovellettu matematiikka
HALKOLA, ANNI S.: Mathematical modelling and survival prediction in can-
cer
Va¨ito¨skirja, 144 s.
Eksaktien tieteiden tohtoriohjelma (EXACTUS)
Marraskuu 2022
TIIVISTELMA¨
Syo¨pa¨ on yksi yleisimmista¨ kuolinsyista¨ aiheuttaen tarpeen laajalle ja kattavalle tut-
kimukselle. Selviytymismahdollisuuksia seka¨ hoitojen hyo¨tyja¨ ja haittoja (ekonomi-
sia tai terveydellisia¨) voidaan ennustaa hyo¨dynta¨en matemaattista mallintamista ja
regressioanalyysia¨. Lupaavia teoreettisia ehdotuksia onkin jo saatu esimerkiksi mono-
ja kombinaatiohoidoista, jotka voisivat potentiaalisesti parantaa hoitostrategioita. Li-
sa¨ksi analyysien avulla on tunnistettu useita biomarkkereita, jotka ennustavat hoito-
vastetta. Ta¨sta¨ huolimatta syo¨pa¨ on edelleen hankala ja ennalta-arvaamaton vastus-
taja, jota vastaan tarvitaan yha¨ tehokkaampia ja yksilo¨llisempia¨ vaste-ennusteita ja
hoitoja.
Ta¨ssa¨ va¨ito¨skirjassa on ka¨ytetty erilaisia la¨hestymistapoja syo¨va¨n ka¨ytta¨ytymisen,
hoitovasteiden ja potilaan selviytymisen ennustamiseen. Differentiaaliyhta¨lo¨ihin pe-
rustuvalla mallilla tutkittiin, kuinka syo¨va¨n ma¨a¨ra¨ muuttuu erilaisten hoitojen ja
hoitoaikataulujen seurauksena. Lisa¨ksi malliin sisa¨llytettiin mekanismit immuuni-
vasteelle ja immunohoidolle, koska immuunipuolustuksella on ta¨rkea¨ rooli syo¨va¨n
kehittymisessa¨ ja tuhoamisessa. Havaittiin, etta¨ adaptiivinen hoito tuotti pienemma¨t
syo¨pa¨ma¨a¨ra¨t ja hoitoa tarvittiin va¨hemma¨n. Lisa¨ksi kombinaatiohoidot (immuno-
hoito joko kemoterapian tai targetoidun hoidon kanssa) aiheuttivat yleensa¨ pienem-
ma¨t syo¨pa¨ma¨a¨ra¨t kuin hoidot erikseen annettuina. Kombinaatiohoitojen yhteydessa¨
on kuitenkin huomioitava, etta¨ kahdesta pa¨a¨llekka¨isesta¨ hoidosta voi aiheutua ylima¨a¨-
ra¨isia¨ sivuvaikutuksia.
Toisessa la¨hestymistavassa syo¨pa¨a¨ mallinnettiin metapopulaatiomallin avulla.
Ta¨ssa¨ keskityttiin angiogeneesin kehittymiseen ja solujen liikkumiseen paikasta toi-
seen (emigraatio). Mallin avulla tutkittiin, missa¨ olosuhteissa syo¨pa¨solut muuttuvat
angiogeenisiksi ja miten hoidot vaikuttavat (antiangiogeeninen, sytotoksinen tai kom-
binaatiohoito). Yleisesti syo¨pa¨solujen kannalta oli hyo¨dyllista¨ edista¨a¨ angiogeneesia¨,
jos antiangiogeenista¨ hoitoa ei ollut. Hoidon kanssa angiogeneesi ha¨visi, mutta
toisaalta myo¨s riski hoitoresistenttiydelle kasvoi.
Lisa¨ksi kehitettiin uusia regressiometodeita erityistesti potilaan selviytymistoden-
na¨ko¨isyyden ennustamiseen. Greedy Cox (a greedy budget-constrained Cox re-
gression) metodi tuottaa regressiomalleja, joissa muuttujien ma¨a¨ra¨a¨ rajoitetaan 𝐿2-
normin avulla. Greedy Cox huomioi, kuinka paljon piirteiden mukaan ottaminen
maksaisi ja metodille voidaan antaa ka¨ytetta¨vissa¨ oleva budjetti. Greedy Coxia myo¨s
verrattiin LASSO-pohjaiseen mallinvalintaan (𝐿1-normi rajoittajana). Uusi metodi
OSCAR (Optimal Subset CArdinality Regression) kehitettiin hyo¨dynta¨en 𝐿0-pseu-
vii
Anni S. Halkola
donormia muuttujien ma¨a¨ra¨n rajoittamisessa. Lisa¨ksi valittujen muuttujien aiheut-
tama hinta laskettiin ja verrattiin, miten se suhteutui mallin ennustustarkkuuteen.
Kaikki esitetyt metodit validoitiin kliinisella¨ eturauhassyo¨pa¨aineistolla. Metodit tuot-
tivat hyva¨n ennustetarkkuuden jo pienella¨ ma¨a¨ra¨lla¨ muuttujia, jolloin myo¨s hinta
pysyi matalampana. Metodit myo¨s valitsivat ja¨rkevia¨, (eturauhas)syo¨pa¨a¨n yhdis-
tetta¨via¨ parametreja kuten prostataspesifinen antigeeni (PSA).
Ta¨ma¨ va¨ito¨skirja tarjoaa siis kattavan valikoiman uusia va¨lineita¨ syo¨va¨n ka¨yto¨k-
sen ja potilaan selviytymistodenna¨ko¨isyyden ennustamiseen. Tutkimuksessa huomi-
oitiin ta¨rkeita¨ syo¨va¨n osa-alueita, kuten immuunivaste ja angiogeneesin kehittyminen
seka¨ na¨ihin liittyva¨t hoidot, joilla on mahdollisuus muuttaa perinteisia¨ hoitoka¨yta¨n-
to¨ja¨.
ASIASANAT: syo¨pa¨, tappaja T-solut, matemaattinen mallintaminen, yksilo¨llinen
la¨a¨kehoito, immunohoito, anti-angiogeeninen hoito, sytotoksinen hoito, targetoitu
hoito, kombinaatiohoito, sivuvaikutukset, syo¨va¨n mikroympa¨risto¨, metapopulaatio-
malli, hintaoptimointi, selvia¨misen ennustaminen, Cox-regressio, eturauhassyo¨pa¨,
piirteiden valinta, 𝐿0-pseudonormi
viii
Acknowledgements
First, I warmly thank my supervisors Associate professor Kalle Parvinen and Profes-
sor Tero Aittokallio for their valuable advice, constructive comments and excellent
support during the preparation of this dissertation. I received guidance whenever I
needed it. Especially the weekly meetings in the beginning were crucial for keeping
me on track and on schedule. I thank Kalle for all the courses about adaptive dy-
namics, mathematical modelling and metapopulation models that I took during my
studies. I thank Kalle for supervising my bachelor’s and master’s theses. I am also
grateful to Kalle for general discussions we have had during the past years. I thank
Tero for coming up with interesting research ideas and introducing me to many pos-
sible collaborators. I know most of my co-authors because of Tero. In addition, I
have Tero to thank for the funding that I got for over two years and for travel ex-
penses.
I am grateful to Professor Marko Ma¨kela¨, who was the research director during the
process of this dissertation. I have also attended multiple Marko’s courses and es-
pecially I am thankful for the Modelling Project course (Mallinnusprojekti), that
connected the real world to the theoretical math studies. I thank the Department of
Mathematics and Statistics for accepting me as a student back in 2011 and again as
a PhD student in 2018, as well as providing my education during these years. I am
also thankful for all the work equipment and work facilities. I thank the head of
department Professor Iiro Honkala and other staff at the Department of Mathematics
and Statistics as well as at the Faculty of Science. I am thankful for the financial sup-
port of the Doctoral Programme in Mathematics and Computer Sciences (MATTI)
of University of Turku, the Academy of Finland, Cancer Society of Finland and the
Sigrid Juse´lius Foundation.
I thank all of my co-authors Kalle Parvinen, Tero Aittokallio, Satu Mustjoki, Henna
Kasanen, Teemu D. Laajala, Kaisa Joki, Mika Murtoja¨rvi, Marko Ma¨kela¨, Tuomas
Mirtti, Tapio Pahikkala and Antti Airola. I am also grateful to all the anonymous re-
viewers that commented on the articles. I thank Harri La¨hdesma¨ki and Alvaro Ko¨hn-
Luque for reviewing this dissertation. In addition, I thank Professor Joel Brown for
agreeing to act as my opponent.
ix
Anni S. Halkola
I thank Arho Virkki for being the link between us and TYKS. He provided the re-
mote connection to TYKS server and answered all my questions when I did not know
how to use that connection. I am also grateful to Mika Murtoja¨rvi for conducting the
original data selection and answering numerous questions about the data and the
selection criteria. I thank Teemu D. Laajala for advice on R-coding, PhD studies
and research in general. I have no idea how he had time to help me on top of his
own work but he did it anyway. I thank Satu Mustjoki, Henna Kasanen and Jukka
Westermarck for their knowledge and advice when I needed to know more about
cancer and its behavior.
I thank Leo Lahti, Aaro Salosensaari, Ville Laitinen and Noora Kartiosuo for let-
ting me join them in their machine learning seminar. I appreciate the hard work
that Eva Kisdi, Stefan Geritz and Mats Gyllenberg, among others, did for numerous
biomathematics events, such as seminars, biomathematics days and summer school
that I were able to attend.
I thank my friends and colleagues Kaisa Joki and Teemu D. Laajala for discussions
during and after (remote) meetings. General, non-work-related discussions were
greatly appreciated, especially during the pandemic. I am also grateful to my friend
and former colleague Pauliina Ma¨kinen for our discussions, as well as all the mu-
tual complaining about the PhD research. Finally, I want to thank my other friends
and family. I am grateful to HybridiSpeksi and all the baes for providing me with
a hobby and mental breaks that have kept me somewhat sane (or not?) since 2014.
I thank Eevastiina for our over 11-year-old friendship that was born back in 2011
when we both started our math studies. I thank Elina for our friendship that is old
enough to buy hard liquor. I thank my relatives, especially my late grandparents who
were proud of my education. Last but not least, I thank my parents and brothers for
imagination and everything else. Kiitos!
November 2022
Anni S. Halkola
x
Table of Contents
Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix
Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi
Abbreviations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii
List of Original Publications . . . . . . . . . . . . . . . . . . . . . . xv
1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
2 Cancer cell population model with immune system . . . . 3
2.1 Literature review . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.1.1 Immune system and immunotherapy in modelling . 5
2.1.2 Treatment modelling . . . . . . . . . . . . . . . . . . 6
2.2 Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.2.1 Resources . . . . . . . . . . . . . . . . . . . . . . . . 6
2.2.2 Cancer cells . . . . . . . . . . . . . . . . . . . . . . . 7
2.2.3 Immune system . . . . . . . . . . . . . . . . . . . . . 8
2.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.3.1 No treatment . . . . . . . . . . . . . . . . . . . . . . . 11
2.3.2 Mono-immunotherapy . . . . . . . . . . . . . . . . . 12
2.3.3 Combination therapy . . . . . . . . . . . . . . . . . . 16
3 Cancer metapopulation model with angiogenesis . . . . . . 19
3.1 Literature review . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.1.1 Angiogenesis modelling . . . . . . . . . . . . . . . . 20
3.1.2 Treatment modelling . . . . . . . . . . . . . . . . . . 21
3.2 Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.2.1 Within-patch dynamics . . . . . . . . . . . . . . . . . 22
3.2.2 Metapopulation dynamics . . . . . . . . . . . . . . . 24
3.2.3 Evolutionary dynamics . . . . . . . . . . . . . . . . . 25
3.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.3.1 Single strategies . . . . . . . . . . . . . . . . . . . . 28
xi
Anni S. Halkola
3.3.2 Joint evolution of angiogenesis and emigration strate-
gies . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.3.3 Joint evolution with treatment resistance . . . . . . . 31
4 Cancer survival prediction based on clinical data . . . . . . 34
4.1 Literature review . . . . . . . . . . . . . . . . . . . . . . . . . 34
4.2 Materials and methods . . . . . . . . . . . . . . . . . . . . . 36
4.2.1 Data . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
4.2.2 LASSO selection and Greedy Cox methods . . . . . 37
4.2.3 OSCAR method . . . . . . . . . . . . . . . . . . . . . 38
4.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
4.3.1 LASSO selection and Greedy Cox . . . . . . . . . . 40
4.3.2 OSCAR . . . . . . . . . . . . . . . . . . . . . . . . . . 42
5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
List of References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
Original Publications . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
xii
Abbreviations
AGEGRP Age group
ALB Albumin
ALP Alkaline phosphatase
APM-𝐿0 Augmented penalized minimization-𝐿0
AST Aspartate aminotransferase
BMI Body mass index
CA Calcium
CARs Chimeric antigen receptors
CHF Congestive heart failure
CRAN Central R Archive Network
CREAT Creatinine
CTLA4 Cytotoxic T-lymphocyte antigen-4
B-PVKT Perusverenkuva ja trombosyytit (Basic blood count and throm-
bocytes, here includes HB, HEMAT, RBC, WBC and PLT)
DBDC The double bundle method
DC function Difference of two Convex functions
FGFs Fibroblast growth factors
HB Hemoglobin
HEMAT Hematocrit
LASSO Least absolute shrinkage and selection operator
LDH Lactate dehydrogenase
LMBM The limited memory bundle method
LYMperLEU Ratio of lymphocytes and leukocytes
mCRPC Metastatic castration-resistant prostate cancer
MG Magnesium
NA Sodium
NEU Neutrophils
NEUperLEU Ratio of neutrophils and leukocytes
NP hard Nondeterministic polynomial hard
ODE Ordinary differential equation
OSCAR Optimal Subset CArdinality Regression
PDGFs Platelet derived growth factors
PD-1 Programmed cell death protein 1
xiii
Anni S. Halkola
PD-L1 Programmed death ligand
PLT Platelets
PSA Prostate-specific antigen
TYKS Turun yliopistollinen keskussairaala (Turku University Hospital)
VEGF Vascular endothelial growth factor
WBC White blood cells
xiv
List of Original Publications
This dissertation is based on the following original publications, which are referred
to in the text by their Roman numerals:
I Anni S. Halkola, Kalle Parvinen, Henna Kasanen, Satu Mustjoki, Tero Ait-
tokallio. Modelling of killer T-cell and cancer cell subpopulation dynamics
under immuno- and chemotherapies. Journal of Theoretical Biology, 2020;
488: 110136.
II Anni S. Halkola, Tero Aittokallio, Kalle Parvinen. Tumor microenviron-
ment as a metapopulation model: The effects of angiogenesis, emigra-
tion and treatment modalities. Journal of Theoretical Biology, 2022; 545:
111147.
III Mika Murtoja¨rvi, Anni S. Halkola, Antti Airola, Teemu D. Laajala, Tuo-
mas Mirtti, Tero Aittokallio, Tapio Pahikkala. Cost-effective survival pre-
diction for patients with advanced prostate cancer using clinical trial and
real-world hospital registry datasets. International Journal of Medical In-
formatics, 2020; 133: 104014.
IV Anni S. Halkola, Kaisa Joki, Tuomas Mirtti, Marko M. Ma¨kela¨, Tero Ait-
tokallio, Teemu D. Laajala. OSCAR: Optimal subset cardinality regres-
sion using the L0-pseudonorm with applications to prognostic modelling of
prostate cancer. Submitted for publication.
The original publications have been reproduced with the permission of the copyright
holders.
xv
1 Introduction
Cancer is one of the main causes of death with more than 19 million new cases per
year [1]. Mathematical and prognostic survival modelling offer insights into individ-
ualized patient care and treatment outcomes.
Tumors have their individual, genetically, functionally or epigenetically distinct
cancer cell populations that develop throughout the disease progression and in re-
sponse to the treatments. Understanding the dynamic evolutionary processes that
undergo before and during treatments, enables the prediction of treatment outcomes
and eases the optimization and personalization of cancer treatment. Mathematical
modelling has already provided some insights into treatment sensitivity and resis-
tance along with suggestions of mono- or combination therapies [2; 3; 4; 5; 6; 7; 8].
Modelling aids the understanding of underlying dynamics and investigation of
the qualitative possibilities of different treatment regimens. In addition, data-based
prediction is crucial in the cancer care to predict the risk of death or cancer recur-
rence, and to identify patients that would benefit from a specific treatment. The
treatment outcome is predicted by utilizing laboratory tests. However, depending
on the cancer type, the process may involve numerous tests and the healthcare costs
often rise. The prognostic prediction by survival models (here Cox’s regression [9]),
especially penalized regression, allows the selection of most prominent predictors in
the assessment of treatment outcome with respect to survival (e.g. overall survival or
progress free survival). Thus, the model sparsity enables pruning unnecessary pre-
dictors along with their costs.
This dissertation includes three chapters that focus on the different approaches that
were used in the cancer modelling and outcome prediction. The chapters summarize
the results and related research of this work. In Chapter 2 the focus is on Publication
I, in which cancer was modelled with ordinary differential equations. The aim was
to model the dynamics between immune system and cancer, as well as to predict
possible qualitative outcomes of different treatment modalities. Thus, the model in-
cludes the immune system (active killer T-cells) and an immunotherapy mechanism
that boosts the immune response in case cancer has started to prevent the T-cells’
efficacy. The Chapter 2 is further divided into sections that offer a review on the sub-
ject (Section 2.1), an overview to the presented model (Section 2.2) and the results
(Section 2.3).
1
Anni S. Halkola
Chapter 3 focuses on the metapopulation modelling approach presented in Pub-
lication II. The aim was to gain insights into underlying evolutionary dynamics that
lead to cancer development. The cancer cells had three strategies (contribution to
angiogenesis, emigration and treatment resistance) that evolved based on the current
living conditions that were determined by the cells, resource availability and treat-
ments. Angiogenesis mechanism was included along with anti-angiogenic treatment
and treatment resistance. Chapter 3 is divided into literature review (Section 3.1),
the model and results overviews (Sections 3.2 and 3.3 respectively).
Chapter 4 summarizes both Publication III and IV as both were based on the Cox
regression [9] and predicting the treatment outcome on prostate cancer data cohorts.
In Publication III, the aim was to offer a method to efficiently predict treatment out-
comes (e.g. overall survival or progress-free survival), while considering also the
costs that the model variables would have in practice. Two algorithms were tested:
LASSO-based selection and a greedy budget-constrained Cox regression (Greedy
Cox). In the latter, the budget was considered within the algorithm. In Publica-
tion IV, the aim was to produce prognostic models while offering different levels
of model sparsity. The proposed OSCAR method implemented penalization with
𝐿0-pseudonorm, which is less studied in penalized regression due to its discontin-
uous and nonconvex nature. Chapter 4 is divided into three sections, starting with
literature review (Section 4.1). The second section (Section 4.2) is an overview on
materials and methods, which is further divided into introduction of data and meth-
ods in Publications III and IV. The third section (Section 4.3) gives an overview on
the results of Publications III and IV.
In Chapter 5 the results from all publications are discussed with ideas for future
research. Finally, Chapter 6 concludes the dissertation with a short summary.
2
2 Cancer cell population model with
immune system
In Publication I, cancer has been modelled with an ordinary differential equation
(ODE) model. The aim was to model the dynamic competition between the immune
system and cancer, and to gain insights into qualitative behavior of the cells under
and after different treatment choices. Thus, in addition to cancer cells, the immune
system was included in the form of active killer T-cells. The inclusion of T-cells
enabled the inclusion of immunotherapy. We also modelled the effects of targeted
therapy and chemotherapy and investigated the treatments as monotherapies or as
a combination. See Figure 1 for the components and their interactions. The cancer
cell type was determined by the cell qualities like resource consumption rate, aggres-
siveness against T-cells and the effectiveness of treatments. We investigated a few
different cancer cell types and predicted how different treatment regimens affected
the cell densities.
- Inactivation
+ Activation
- Elimination
- Death - Death
+ T-cell division
+ Inflow - Outflow
- Consumption+ Cell division
Figure 1. The interactions of (cancer) cells (C), active killer T-cells (T) and resources (R) (e.g.
glucose). The ellipses point the affected entity with the sign indicating the change in the density.
Lightning bolt arrows indicate the treatments: chemotherapy affects the cells divisions, targeted
therapy affects the cancer cell death and immunotherapy affects the inactivation of T-cells. (Adopted
with permission from Publication I: Fig. 1.)
2.1 Literature review
As cancer remains one of the most common causes of death, it is also a solid interest
for mathematical modelling. In addition to dynamic modelling of the cancer cells
3
Anni S. Halkola
or density, the models may also suggest theoretical therapy regimens and predict
possible treatment outcomes. Similarly to our model in Publication I, many mod-
els assume that cancer has survived the early stages and has reached a considerable,
possibly detectable, size. This assumption justifies that the tumor can be modelled
as a continuous variable, enabling the use of differential equations in the modelling
of the cancer dynamics. For example, Fassoni et al. [10] proposed an ODE model
for chronic myeloid leukemia. They consider multiple different sub-models with
different functional interactions between cells. Letellier et al. [11] included four
ODEs in their model (host, immune, tumor and endothelial cells). Also others, such
as [12; 13; 7] utilized differential equations in their cancer modelling. Pinho et al.
[12] had five differential equations for normal, cancer and endothelial cells along
with chemotherapy and anti-angiogenic agents. Also Yonucu et al. [13] focused on
the effects of anti-angiogenic and chemotherapies with ODEs for tumor, vasculature,
interstitial fluid pressure and therapy agents. Zhang et al. [7] considered abiraterone
therapy in prostate cancer model with three cancer cell types that differ in their an-
drogen relation. They modelled the system as Lotka-Volterra equations including an
interaction matrix of the different cell types. Louzoun et al. [5] modelled pancreatic
cancer with four ODEs including cancer cells, host cells, T-cells and the ratio of two
macrophage phenotypes.
Many of the above mentioned models included some form of normal or host
cells along with the cancer cells. In our ODE model, the cell type is general, though
mainly thought to be cancerous in the presented analyses. However, the cell type
may be modified by selecting cell type specific parameter values indicating a normal
cell behavior. For example, the proliferation rate may be decreased or the normal
apoptosis rate increased. Additionally, the interaction with the active killer T-cells
may be selected to be small or non-existent.
In addition to differential equations, other approaches have been considered in re-
search [6; 8]. Opposite to the assumption that cancer has reached a considerable
size, other models have considered only a single cell or a few cells. These early
stages may be modelled by branching processes. Avanzini and Antal [14] proposed
a branching process model to investigate the formation and evolution of metastatic
lesions based on the size of the primary tumor. Also Bozic and Nowak [3] focused
on metastatic cancer and model a lesion as a branching process with possibility of
treatment resistance. Bozic et al. [2] investigated the need of targeted combination
therapy with a model based on multitype branching process. Kozłowska et al. [15]
utilized a stochastic multitype branching process in which cancer cells can acquire
multiple resistance mechanisms.
Various other approaches have been used as well, for example, cellular automata
by Lai et al. [16] and Alarco´n et al. [17]. Komarova, Burger and Wodarz [18] used
a stochastic birth-death process model to predict the relapse and emergence of treat-
4
Cancer cell population model with immune system
ment resistance against a targeted therapy (ibrutinib). Fischer, Va´zquez-Garcı´a and
Mustonen [4] used a Wright-Fisher model with stochastic optimal control to investi-
gate adaptive therapies in case of treatment resistance. A metapopulation modelling
approach has been used in Publication II (see Chapter 3).
2.1.1 Immune system and immunotherapy in modelling
Many cancer models do not include the immune system, which has a significant role
in the development and life cycle of cancer. In our model in Publication I, we have
included the immune system in the form of active killer T-cells, that are an impor-
tant part of the immune response. Fassoni et al. [10] included the immune system
cells along with quiescent and proliferating cancer cells. Similarly to our model,
they considered one type of immune cells (e.g., natural killer cells). However, in
their model, cancer cells only affect the immune cell density positively by recruiting,
whereas in our model the cancer cells may also decrease the T-cell density. This neg-
ative effect on the immune system justifies the use of immunotherapy and enables its
inclusion in our model. Fassoni et al. considered targeted therapy (tyrosine kinase in-
hibitors) that affects only proliferating cancer cells. Letellier et al. [11] modelled the
immune cells, however, similarly to Fassoni et al. they did not include immunother-
apy. Instead, their focus was on the combination of anti-angiogenic treatment and
chemotherapy.
In our model in Publication I the active killer T-cells are produced based on the
antigen delivery from dying cancer cells [19; 20; 21]. The active killer T-cells are
in contact with the cancer cells causing the death of the latter. However, we in-
cluded a mechanism in which the cancer cells may cause inactivation of the T-cells
leading to ineffective T-cells. This model mechanism is based on the binding be-
tween programmed death-ligand (PD-L1) on the cancer cell surface and its receptor
programmed cell death protein 1 (PD-1) on the T-cell surface [19]. The binding pre-
vents the T-cell function, however, this may be prevented by a form of immunother-
apy, immune-checkpoint inhibitors (e.g., nivolumab and pembrolizumab). Such a
immune-checkpoint inhibitor is, for example, anti-PD-L1, which binds to PD-L1
preventing the connection with PD-1, or anti-PD-1 which binds to PD-1 leading to
similar results [22; 23]. We included this mechanism in our model in Publication I.
Since the immune response is a complex system, the PD-L1/PD-1 interaction is not
the only mechanism that may be targeted by immunotherapy. The immune response
requires the antigen delivery, antigen presentation to the immune system, priming
and activation of T-cells, trafficking and infiltration of T-cells into tumors and recog-
nition and killing of cancer cells [19].
Chemotherapy, radiation and targeted therapies might enhance the antigen deliv-
ery by increasing the cancer cell death. Antigen presentation can be promoted by im-
5
Anni S. Halkola
munotherapy vaccines. Cancer cells may affect the activation of T-cells, which may
be inhibited by CTLA4 (a negative regulator of T-cells). This interaction is targeted
by anti-CTLA4 antibodies (e.g. ipilimumab) [24]. The trafficking and infiltration of
T-cells into tumor microenvironment might also be compromised. Treatments like
anti-VEGF (bevacizumab) and B-Raf inhibitors (vemurafenib) may enhance the tu-
mor infiltration [25; 26]. In CAR-T immunotherapy, CARs (chimeric antigen recep-
tors) are produced from patient’s T-cells and modified to adhere and attack cancer
cells [27]. Also combinations of multiple immunotherapies have been considered
[28; 29; 30; 31].
2.1.2 Treatment modelling
In addition to immunotherapy, in Publication I, we considered also chemotherapies
(e.g. dacarbazine or temozolomide) and targeted therapies (e.g. BRAF or C-KIT
inhibitors). Since the cytostatic treatment damages the cell cycle benchmarks [32],
rapidly dividing (cancer) cells have less time to recover, thus making them more
vulnerable for the treatment. Targeted treatment recognizes specific proteins, thus
allowing targeted elimination of only cells with the target protein. These are more
traditional treatment modalities and as such are more widely considered in the cancer
models. For example, Fassoni et al. [10] modelled targeted therapies, and Pinho et
al. [12] and Yonucu et al. [13] modelled chemotherapies. In our model in Publica-
tion I, we used these therapies also in a combination with immunotherapy.
Other mono- or combination therapies can also be considered in mathematical
modelling. Depending on the scenario and cancer type, one could consider, for
example, hormone therapy like abiraterone therapy in prostate cancer [7] or radi-
ation therapy. Also anti-angiogenic treatment has been used in cancer treatment
[33; 34; 35] and we have modelled such treatment mechanism in Publication II (see
Chapter 3).
2.2 Model
The minimum model consists of three ordinary differential equations: resources (e.g.
glucose), active killer T-cells and cancer cells. More cell populations could be in-
cluded by adding a differential equation for each with cell type specific parameter
values.
2.2.1 Resources
Since the cells require resources for their proliferation, we included the amount of
resources as a part of the competitive situation. Resources flow in to and out of the
tumor microenvironment following chemostat dynamics in the absence of resource
6
Cancer cell population model with immune system
consuming cells [36]. However, in addition to normal flow through the region, the
cells consume the resources following the law of mass action [37]. Thus the differ-
ential equation for resources is
𝑑𝑅
𝑑𝑡
= 𝜆(?^?−𝑅)⏟  ⏞  
Flow
−
⎛⎜⎝1− 𝑤
∑︀
𝑗
𝐶𝑗
1 +
∑︀
𝑗
𝐶𝑗
⎞⎟⎠∑︁
𝑗
𝛼(𝑠𝐶𝑗 )𝑅𝐶𝑗
⏟  ⏞  
Consumption by cells
, (1)
where 𝜆 is the relative flow speed, ?^? is the baseline resource inflow and𝐶𝑗 is the den-
sity of type 𝑗 cells. Cells are assumed to restrict each other’s resource consumption,
for example, by increasing the distance between veins [38], and 𝑤 is the maximum
proportion at which cells limit each other’s resource consumption.
Figure 2. The form of the resource consump-
tion rate 𝛼(𝑠𝐶𝑗 ) with the lower and upper lim-
its for the parameter values presented in the
Publication I.
Here 𝑠𝐶𝑗 is the proliferation strat-
egy of the cell type 𝐶𝑗 and 𝛼(𝑠𝐶𝑗 ) is
the corresponding resource consumption
rate. The rate 𝛼(𝑠𝐶𝑗 ) is an increas-
ing function with lower and upper lim-
its (see Figure 2). The lower limit
is assumed since cells require at least
some resources to maintain themselves,
and the upper limit is assumed since the
cells cannot consume resources infinitely
fast.
2.2.2 Cancer cells
Cancer cell dynamics are determined by the
proliferation proportional to the resource
consumption and death by normal apoptosis, drugs or active killer T-cells. The cells
attempt to divide based on their resource consumption:
𝛾𝛼(𝑠𝐶𝑖)
⎛⎜⎝1− 𝑤
∑︀
𝑗
𝐶𝑗
1 +
∑︀
𝑗
𝐶𝑗
⎞⎟⎠𝑅𝐶𝑖, (2)
where 𝛾 is the conversion coefficient. However, the division may not be successful
in the presence of cytostatic drugs that interfere with the mitosis. Thus, the division
fails and the cells face death with the probability 𝑝𝑖. Accordingly, the division pro-
ceeds with the probability (1− 𝑝𝑖). The cytostatic treatments affect rapidly dividing
7
Anni S. Halkola
cells more than cells that hardly divide, therefore, the probability 𝑝𝑖 depends on the
proliferation strategy 𝑠𝑖:
𝑝𝑖 =
𝑠𝑖
1 + 𝑠𝑖
𝐻𝑝, (3)
where 𝐻𝑝 is the dose-response effect (Hill equation [39]) for the cytostatic drug. It
is assumed for simplicity that the drug concentration is a positive constant during the
treatment periods and zero afterwards.
In addition to cytostatic drugs, the division may be prevented by the microenvi-
ronment carrying capacity, assumed to depend on, e.g., limitations of space. Thus,
aggregating the above factors, a successful division happens with the likelihood
(1− 𝑝𝑖)
⎛⎜⎝1−
∑︀
𝑗
𝐶𝑗
𝐾
⎞⎟⎠ , (4)
where 𝐾 is the carrying capacity.
Cancer cells of the type 𝑖 die by natural apoptosis with the death rate 𝜇𝑖. This
rate is increased by the rate 𝛽𝑖 in case a targeted drug is administered. Now 𝛽𝑖 = 0, if
the cell type 𝑖 is resistant to the targeted treatment. Similarly to cytostatic treatment,
the targeted drug concentration is considered to be a constant during the treatment
period, and thus the rate 𝛽𝑖 is also considered to be a constant. In addition to drugs,
active killer T-cells may cause cancer death as the proportion 𝜃 of the T-cells enters
the cancer microenvironment and proceeds to attack the cancer cells with the rate 𝜉𝑖.
However, if a cancer cell type lacks the specific antigen, the T-cells ignore the cells
and 𝜉𝑖 = 0. Thus, we obtain the following differential equation for the cancer cells:
𝑑𝐶𝑖
𝑑𝑡
= 𝛾𝛼(𝑠𝐶𝑖)
⎛⎜⎝1− 𝑤
∑︀
𝑗
𝐶𝑗
1 +
∑︀
𝑗
𝐶𝑗
⎞⎟⎠𝑅𝐶𝑖
⏟  ⏞  
Proliferation attempt
⎛⎜⎝ (1− 𝑝𝑖)
⎛⎜⎝1−
∑︀
𝑗
𝐶𝑗
𝐾
⎞⎟⎠
⏟  ⏞  
Proliferation succeeds
−𝑝𝑖⏟ ⏞ 
Failure and death
⎞⎟⎠
−𝜇𝑖𝐶𝑖 − 𝛽𝑖𝐶𝑖 − 𝜉𝑖𝐶𝑖𝜃𝑇⏟  ⏞  
Death
.
(5)
2.2.3 Immune system
The killer T-cells are activated in relation to the antigens delivered by the dying
cancer cells. The amount of antigen delivery rate 𝑑𝐶 depends on the cancer cell
death as the dying cancer cells deliver antigens to the blood stream [19]. The rate
𝑑𝐶 considers the normal cancer cell apoptosis, death by T-cells and death due to the
8
Cancer cell population model with immune system
targeted or cytostatic drugs. The killer T-cell activation rate ℎ(𝑑𝐶) is assumed to
be sigmoid function (Fig. 3a) since the activation may not be initiated properly with
small amount of antigen delivery due to the insufficient signal [40]. Furthermore, the
T-cell activation rate ℎ(𝑑𝐶) has the maximum value of 𝑚.
a) T-cell activation b) The immunotherapy
Figure 3. a) The form of the killer T-cell activation ℎ(𝑑𝐶). b) The form of the overall effect of the
immunotherapy 𝜌𝑖.
The activated T-cells attempt division with the rate 𝛼𝑇 . However, cytostatic drugs
may cause division failure and subsequent cell death with the probability 𝑝𝑇 . Thus
the dividing T-cells die with the rate 𝛼𝑇 𝑝𝑇 and succeed to divide with the rate 𝛼𝑇 (1−
𝑝𝑇 ). The activated T-cells may also die due to natural apoptosis (𝜇𝑇 ). The proportion
𝜃 of the active killer T-cells enters the microenvironment. Those T-cells then proceed
to kill the cancer cells that have the specific antigen. The cancer cells may also make
the T-cells ineffective with the cell type specific rate 𝜙𝑖, while the immunotherapy
(𝜌𝑖) counteracts the action. The effect function of immunotherapy 𝜌𝑖 is sigmoidal,
following the general form of dose-response curves [39] (Fig. 3b). For the effect of
immunotherapy, 𝜌𝑖 = 1 is desired as it would lead to full prevention of PD-1/PD-L1
interaction between type 𝑖 cancer cells and T-cells ((1 − 𝜌𝑖)𝜙𝑖𝐶𝑖𝜃𝑇 = 0). As with
the other treatments, the concentration of the immunotherapy drug is set to a positive
constant or zero depending on whether the patient is on a treatment period or on a
treatment break (drug holiday). Finally, the immune system has self-regulation [41]
with the rate 𝛿. Thus, we obtain the following differential equation for the active
killer T-cells:
𝑑𝑇
𝑑𝑡
= ℎ(𝑑𝐶)⏟  ⏞  
Activation
+𝛼𝑇 (1− 𝑝𝑇 )𝑇⏟  ⏞  
Proliferation
−𝜇𝑇𝑇 − 𝛼𝑇 𝑝𝑇𝑇⏟  ⏞  
Death
−𝛿
2
𝑇 2⏟  ⏞  
Self-regulation
−
∑︁
𝑖
(1− 𝜌𝑖)𝜙𝑖𝐶𝑖𝜃𝑇⏟  ⏞  
PD-1/PD-L1 interaction
.
(6)
9
Anni S. Halkola
2.3 Results
We investigated several case studies, i.e. virtual patients characterized by patient
specific model parameters. We focused on the immunotherapy effect either as a
monotherapy or in a combination with targeted and chemotherapies. The immune-
checkpoint inhibitors generally used in practice include, for example nivolumab ad-
ministered fortnightly or pembrolizumab which is given every 3rd week. However,
these schedules may not be optimal and thus the timing and duration of treatment
remains under investigation.
Since the main focus is on the immune system, the virtual patients differed in
their T-cell functions. We altered the maximum activation rate of killer T-cells as
well as the birth rate of the activated T-cells along with the interaction rates between
T-cells and cancer cells. See Table 1 for the altered parameters and their definitions.
Table 1. The main parameters defining a virtual patient. Two cancer cell types were included, thus
index 𝑖 = 1, 2.
Symbol Meaning
𝛼𝑇 Division attempt rate of active killer T-cells.
𝜉𝑖 Rate at which T-cells kill cancer cell subpopulation 𝑖.
𝜙𝑖
Rate at which cancer cell subpopulation 𝑖
makes killer T-cells ineffective.
𝑚 Maximum rate of killer T-cell activation.
In addition to parameters, different treatment options were considered for different
virtual patients. See Table 2 for the treatment parameters and in which case studies
they were considered. We investigated also a few virtual patients without any treat-
ment to monitor the dynamics in an untreated situation.
Table 2. The treatments altered for each case study. Here index 𝑖 = 1, 2.
Symbol Meaning Case study
𝜌𝑖 Effect of immunotherapy on cell type 𝑖. 1, 2, 3 and 4
𝛽𝑖 Rate at which cell type 𝑖 dies due to targeted treatment. 3
𝑐𝑝 Concentration of the cytostatic drug. 4
We calculated numerically (Runge-Kutta method, RK45) the dynamics of ODEs for
each virtual patient. We compared the numerical solutions qualitatively to observed
clinical outcomes in cancer [42; 43; 44; 45; 46; 47] to justify the majority of selected
parameter values. The virtual patients and their parameter values were selected to
10
Cancer cell population model with immune system
represent different qualitative results, such as progressive, chronic or dormant dis-
ease, that have been observed in real patients. Furthermore, the cancer cell doubling
times were calculated for each case study and they were found reasonable compared
to experimental observations in melanoma [48]. The value of maximum proportion at
which cancer cells limit each other’s resource consumption (𝑤) was selected so that,
multiplied by the saturation of cancer cells, it fell into the limit of 4–56% relative
volume of necrotic tissue in melanoma [49]. For the selection of T-cell parameters,
we assumed that if the active killer T-cells have existed in the system, a small density
remains even if the cancer cell density becomes zero (memory in the immune sys-
tem).
Generally estimated values have to be used, since in clinical practice it is still im-
practical or even impossible to estimate all parameters individually. However, based
on the additional sensitivity analysis, the most important individual parameters relate
to the cancer cell and T-cell interactions (𝜉𝑖 and 𝜙𝑖), the infiltration of T-cells (𝜃) and
the activation of the T-cells (𝑚). Perturbations in these parameter values affected
the treatment outcome the most, and for example, the possibility of dormant disease
versus inevitable relapse depended on them.
It was assumed that an initial immune response had already occurred, and inves-
tigated virtual patients started with densities of 0.05 for all cell types. However, the
initial state affects the outcome and required treatment schedule. For example, if the
cancer cell density is high with nearly non-existent T-cell density, an intense initial
treatment schedule may be required to gain a treatment response.
In the case studies we included one or two cancer cell subpopulations (denoted
by 1 and 2). However, more subpopulations could be included along with normally
behaving cell types. One cancer cell type is considered to divide faster (𝑠𝐶1 > 𝑠𝐶2)
and with targeted therapy, one subpopulation is considered resistant (𝛽1 = 0) while
the other one remains sensitive (𝛽2 > 0).
2.3.1 No treatment
We investigated a few baseline cases to gain insights into possible outcomes when
treatment is not given. In some cases the treatment was not even needed as the
immune system is able to keep cancer in control [50; 51]. Depending on the qualities
of the immune system (e.g. the T-cell proliferation and activation rates), the cancer
cell density approached a fixed steady state, or the density fluctuated as the T-cells
kept defeating the cancer cells at each relapse attempt. However, in practice such
behavior is likely to happen in the early stages of cancer development and it is not a
guarantee, that cancer would not variate further and relapse.
Some patients are less fortunate, and despite a possible immune response, cancer
eventually dominates and potentially becomes fatal. In Figure 4 the T-cell density
quickly decreases without treatment, and the cancer cell densities increase. Eventu-
11
Anni S. Halkola
C2
C1
T
Ctot
Figure 4. Cancer cell densities for two subtypes
(𝐶1 in red, 𝐶2 in magenta) with the total cancer
cell density (𝐶𝑡𝑜𝑡, black dashed) and the active
killer T-cell density (𝑇 , green). (Adopted with per-
mission from Publication I: Fig. 2a.)
ally, the less aggressive subtype 𝐶2
starts to decrease as 𝐶1 starts to dom-
inate. Carrying capacity and resources
start to limit the cancer cell density,
however, a smaller cancer density may
have already proven fatal to the pa-
tient. For such patients, the treatment
is needed for a chance of recovery. Dif-
ferent treatment choices may lead to dif-
ferent outcomes and thus it is important
to analyze the options. Unfortunately,
even with careful planning, the treat-
ment may not always lead to full recov-
ery and relapse is inevitable.
2.3.2 Mono-immunotherapy
Immunotherapy has been used in the practice as a monotherapy [52; 46]. We inves-
tigated two different virtual patients who received mono-immunotherapy either as
a single-shot (case study 1) or with repeated treatment periods (case study 2). For
the case study 2, multiple different treatment schedules were considered. For both
of these case studies, the virtual patients would suffer from the failure of immune
response (see Figure 4).
Single-shot immunotherapy (case study 1)
For the case study 1, we considered one cancer cell subpopulation. We investigated
how the attractor landscape changed based on the treatment and how the cancer cell
density was affected by this change. The immunotherapy was initiated when the can-
cer cell density exceeded a pre-defined detection limit threshold (here set to 0.5). As
seen in Figure 5a, only one treatment period is enough to revive the T-cell density
and the immune system is able to regulate cancer to a dormant state [42; 43; 44; 45].
A single dose of immunotherapy (specifically anti-PD-1) has been observed to cause
major responses also in real patients [52; 53].
Here the patient’s cancer and immune system characteristics enabled the change
in the attractor (full blown cancer versus a dormant cancer). Without treatment, the
attractor landscape had two possible steady states. Due to the initial condition, the
cell densities started to reach the attractor with high cancer cell density and low killer
T-cell density (Fig. 6a top row). The immunotherapy changed the attractor landscape
for the duration of treatment period, allowing the cell densities to approach a differ-
ent steady state (Fig. 6a bottom row). After the treatment, the attractor landscape
12
Cancer cell population model with immune system
a) Treatment suitable for one-shot b) Treament unsuitable for one-shot
Figure 5. a) The one-shot immunotherapy is given for a period of 16 days (grey bar) and a dormant
state is reached. b) The immunotherapy is given for a longer period (24 days). The attractor is the
original with high cancer cell density. As the threshold of 0.5 is exceeded again, the treatment would
be repeated. (Panel a adopted with permission from Publication I: Fig. 3a.)
reverted back, however, the cell densities had changed and a different steady state
was reached as the original attractor had now become unattainable (Fig. 6a top row).
Thus, depending on the patient’s biological qualities, there may exist a manageable
dormant state, even if originally cancer is rapidly increasing. However, the timing
and dosage of the single-shot immunotherapy were crucial to achieve the change
from the original, unfavorable attractor to the manageable dormant state. With an
unfavorable treatment (e.g. too short or too long duration) cancer would have pre-
vailed after a small decrease in the cancer cell density (see Fig. 5b and Fig. 6b).
Repeated immunotherapy (case study 2)
Despite the success of one-shot immunotherapy in some cases, many cases require
repetitive treatment. For example, in practice nivolumab and pembrolizumab (immune-
checkpoint inhibitors) are administered every 2nd or 3rd week respectively. We in-
vestigated a virtual patient with two cancer cell subpopulations (𝐶1 and 𝐶2) with one
slightly more aggressive than the other (𝑠𝐶1 = 1, 𝑠𝐶2 = 0.95). The virtual patient
would suffer from a rapid increase in the cancer density if no treatment is used (Fig. 4
is for this patient).
First we investigated how different treatment initiation options affected the out-
come, and whether there was a superior scheduling technique. We compared a pre-set
treatment periods to threshold-based treatment initiation. Pre-set treatment periods
are widely used option, however, a more adaptive alternative is to initiate the treat-
ment when a total cancer cell count (i.e., total tumor burden) exceeds a threshold.
We investigated the options by changing the treatment period and, depending on the
scheduling, either interval between treatments (drug holiday) or the given threshold.
13
Anni S. Halkola
No treatment
Treatment
No treatment
a) Successful one-shot b) Unsuccessful one-shot
Treatment
Figure 6. Phase plane plots showing the isoclines (black lines) of the killer T-cell density T and
cancer cell density C, a) for a successful one-shot immunotherapy and b) for an unsuccessful one-
shot immunotherapy. The trajectories with respect to time are shown in blue. The line is solid when
the treatment status corresponds to the phase plane plot (top row: during no treatment, bottom row:
during treatment). The treatment period is a) 16 days or b) 24 days.
As expected, for the pre-set periods, short treatment periods required short drug
holidays. If the drug holiday was continued for too long, the treatment might fail
even if the treatment was repeated. Interestingly, for the threshold-based initiation,
the treatment period length had little effect on the outcome (total cancer cell mean
density) since short treatment periods were repeated depending on the cancer density
at the end of the previous treatment. If the threshold was still exceeded, the treatment
was repeated immediately without a drug holiday. However, the selection of thresh-
old affected the outcome, as for the higher thresholds the cancer cell density reached
a higher level before the treatment.
Both of the options had treatment regimen parameters that resulted in small can-
cer mean density (examples in Fig. 7). Choosing a specific treatment period along
with either a suitable drug holiday or a threshold might result in a better therapeu-
tic effect. However, a small cancer mean density did not necessarily imply a good
treatment, as it might require long treatment periods with short or non-existent drug
holidays. Such heavy treatment scheduling might be unbearable for the patient and
thus it is important to consider also the time spent under treatment. If the treatment
period and drug holiday were selected carefully, successful treatment result could be
reached without increasing the time in treatment in the chosen time interval (Fig. 7).
14
Cancer cell population model with immune system
0 20 40 60 80 100
0.2
0.4
0.6
0.8
Proportion of time in treatment (%)
C t
o
t m
e
a
n
 d
en
sit
y
ll
l
ll
l
l
ll
ll
l
ll
l Treatment period 4 days
Treatment period 10 days
Treatment period 20 days
Pre−set period
l Treatment period 4 days
Treatment period 10 days
Treatment period 20 days
Treatment threshold
Figure 7. Mean cancer cell density with respect to the proportion of time in treatment (% from the
total follow-up time interval). Either pre-set periods (blue) or threshold-based initiation (green) is
used with treatment periods of 4, 10 or 20 days. For the pre-set periods, the proportion of time in
treatment changes when the drug holiday is changed. For the threshold-based treatment initiation,
the threshold level determines the time in treatment. (Adopted with permission from Publication I:
Fig. 5.)
Overall, cancer was kept in control more likely with a threshold-based treatment ini-
tiation, since the treatment was repeated without drug holiday for as long as needed
to decrease the total cancer density. Despite the possibility of a long continuous treat-
ment period, the sufficient decrease in cancer density might have led to a longer drug
holiday before the density reached the given threshold again. In practice, however,
the use of pre-set treatment periods is more applicable since the threshold-based ini-
tiation requires heavy monitoring of the cancer cell density, and this would require
improvement in the diagnostic tests.
Adapting the treatment (case study 2)
To make the approach of pre-set periods more adaptive and reduce the time in treat-
ment (less side effect and costs) without losing the treatment success, we changed the
treatment schedule in the middle of the follow-up period. Since a long drug holiday
would lead to almost immediate treatment failure, the treatment was started with a
tighter schedule (treatment 12 days, drug holiday 6 days). After the third treatment
period, the drug holiday was increased (16 days). Compared to the non-adaptive
scheduling, the successful adapted treatment resulted in a similar total cancer cell
mean density while the time in treatment decreased (Fig. 8). Despite the successful
treatment adaptation in this case, the treatment outcome was sensitive for the selec-
tion of the new schedule. If the drug holiday was increased even more (to 20 days)
there was a notable increase in total cancer cell density, since the T-cell density would
decrease on an unrecoverable level during the drug holidays. Adapting the time of
the schedule change would not have helped, if the new drug holiday was too long.
15
Anni S. Halkola
Figure 8. The pre-set treatment period of 12 days, followed by 6 days of drug holiday, is changed
to schedule with 16 days of drug holiday after 48 days (black arrow). Right-side barplots: the
treatment success measures, with corresponding values for non-adaptive treatments with 6 days
of drug holiday (blue) and 16 days of drug holiday (cyan). Bottom bar: proportion of time spend in
treatment. (Adopted with permission from Publication I: Fig. 6a.)
2.3.3 Combination therapy
In practice, immunotherapy is also used in a combination with other treatments [54].
We investigated especially combinations with targeted treatment or chemotherapy.
We compared the effectiveness of combination treatments against monotherapies.
Immunotherapy with targeted treatment (case study 3)
The targeted treatment targets a specific molecular aberration in a cancer cell. Un-
fortunately, some cancer cells may be missing the aberration making them treatment
resistant against that treatment. With the virtual patient in case study 3, we inves-
tigated a situation with two cancer cell subtypes and assumed that only one of the
subtypes is sensitive (here 𝐶2) whereas the other is resistant (𝐶1). We tested multiple
different treatment strategies by simultaneous threshold-based treatment initiation,
and altering the durations of targeted and immunotherapy. The combined treatment
duration was set to a constant of 14 days, for simplicity.
Some duration combinations (e.g. both with 7 day duration), resulted in a chronic
disease with repeated treatments. Eventually the sensitive subpopulation diminished
due to the targeted treatment and the resistant subpopulation prevailed. After this,
the targeted treatment became useless and the immunotherapy was the main effec-
tive treatment. On the other hand, a stable disease was also reached with selected
combinations.
16
Cancer cell population model with immune system
We calculated the total cancer cell mean density and number of treatment initia-
tions for varied treatment thresholds and proportions of targeted treatment (Fig. 9).
As expected, smaller thresholds led to smaller cancer mean densities with multi-
ple treatment initiations. However, if the stable disease was reached, the number of
treatment initiations stayed low while the cancer mean density was still relatively
low. Such a stable disease was reached with mono-immunotherapy (0% of targeted
treatment), but also with carefully selected monotherapy of targeted treatment (100%
of targeted treatment). However, in general, the combination therapy worked better
than targeted therapy alone since after the sensitive subpopulation had disappeared,
the targeted treatment had no effect.
a) Total cancer cell mean deansity b) Number of treatment initiations
Figure 9. The effect of immunotherapy and targeted treatment combinations. a) The total cancer
cell mean density with respect to treatment threshold and the proportion of targeted treatment (here
0% equals 14 days of immunotherapy, 100% equals 14 days of targeted therapy, 50% equals 7 days
of both initiated simultaneously). b) The corresponding numbers of treatment initiations. (Adopted
with permission from Publication I: Fig. 7d and e.)
In Publication I we considered only a case where the two cancer subpopulations had
proliferation rates close to each other (1 vs. 0.95) and the less aggressive subpopu-
lation was sensitive to the targeted treatment. If the difference between the prolifer-
ation rates is higher (1 vs. 0.65) and the more aggressive subpopulation is sensitive,
there is another way to manage cancer with targeted therapy alone. Assuming that
there is no other targeted treatment to target the resistant subpopulation, the resistant
subpopulation could theoretically be controlled by the dominating, more aggressive
subpopulation. Then the sensitive subpopulation would be kept under control by
the targeted treatment. Figure 10 presents a case where the more rapidly dividing
subpopulation 𝐶1 is sensitive and the subpopulation 𝐶2 is resistant. The treatment
is initiated when the 𝐶1 cell density exceeds a given threshold. The resistant sub-
population is controlled by the sensitive subpopulation and T-cell increase during
treatments (caused by increase of cancer cell death and antigen delivery). However,
17
Anni S. Halkola
despite cancer being under control, this method would require monitoring the densi-
ties of specific subtypes instead of total cancer cell density, and thus the application
in practice may be challenging.
Figure 10. The dynamics in a case, where the treatment-sensitive subpopulation 𝐶1 is notably less
aggressive than the treatment-resistant subpopulation 𝐶2. Targeted treatment is initiated when the
density of the sensitive population exceeds the treatment threshold.
Immunotherapy with chemotherapy (case study 4)
For the case study 4, we considered a virtual patient who received chemotherapy
either alone or in combination with immunotherapy. To minimize the side effects, it
would be preferred to use the minimal effective treatment regimen (small dosages,
short treatment periods). Depending on the treatment choices, cancer was chronic
with repetitive treatment or, more rarely, a stable disease was reached with one or two
treatment periods. We observed that a combination treatment often resulted in shorter
time in treatment when compared to the chemotherapy alone. Simultaneously, there
was a possibility of smaller cancer cell mean density. The side effects were calculated
as the proportion of T-cell loss caused by cytostatic drug to the overall T-cell loss.
For short treatment periods, the combination therapy decreased the side effects when
compared to the chemotherapy. For longer treatment periods, the situation shifted
slightly to the opposite case, because the side effects were calculated as a proportion
and the extended immunotherapy decreased the overall loss of T-cells. However, the
differences were mild and other types of side effects and tolerability issues could
affect the comparison of chemotherapy and combination treatment.
18
3 Cancer metapopulation model with
angiogenesis
In Publication II, cancer was modelled as a metapopulation, in which cancer pop-
ulation patches along with the dispersal pool (here the circulatory system) form a
metapopulation. Cancer cells may inhibit patches, which are possible sites for clus-
ters of cancer cells. Cells may emigrate into the dispersal pool and immigrate from
there into patches (Fig. 11). Cancel cells have an emigration strategy, which tells
how eagerly the cell tries to emigrate and leave the patch. We included angiogenesis,
which is the formation of new vessels. Cancer cells may contribute to the angiogen-
esis trying to enhance the resource (e.g. glucose) availability. The elevated resource
inflow eases the cell growth when the tumor grows and the pre-existing vasculature
becomes insufficient. Anti-angiogenic treatment is designed to prevent the angiogen-
esis and thus to inhibit the cancer cell growth. In addition to angiogenesis mechanism
in our metapopulation model, we also included an anti-angiogenic treatment.
Figure 11. The metapopulation structure with cancer cells (C) inhabiting the habitat patches, i.e.,
possible sites for cluster of cells (black circles). Cell may emigrate into the dispersal pool, i.e., the
circulatory system (red rectangle) and immigrate from there into habitat patches. (Adopted with
permission from Publication II: Fig. 1.)
19
Anni S. Halkola
3.1 Literature review
As mentioned in Chapter 2, cancer has been widely modelled with different ap-
proaches. However, to our knowledge, our model in Publication II is the only one
considering the cancer ecosystem as a continuous-time metapopulation with varying
strategies. The metapopulation modelling approach has been previously used to in-
vestigate population dynamics in various other scenarios, that are not cancer specific
(see e.g. [55; 56; 57; 58; 59; 60; 61; 62]). The metapopulation dynamics along with
strategy evolution can be investigated using tools of adaptive dynamics [63; 64; 65].
Angiogenesis is the formation process of new blood vessels from pre-existing vascu-
lature. It is an important step in the cancer development, since the expanding tumor
requires new blood supply to maintain growth [66; 67; 68]. The cancer cells de-
velop signalling mechanisms, such as secreting VEGF (vascular endothelial growth
factor), that stimulate angiogenesis [69]. These mechanisms may be targeted by anti-
angiogenic treatments that prevent the angiogenesis (e.g., anti-VEGF) [33; 34; 35].
However, as with many cancer treatments, cancer cells may develop treatment re-
sistance to anti-angiogenic treatment [70; 33]. For example, the targeted angiogen-
esis mechanism might be bypassed with another signalling pathway (e.g., fibroblast
growth factors [33]).
3.1.1 Angiogenesis modelling
Angiogenesis has been modelled before in relation to cancer (see e.g. [11; 12; 13; 71;
72; 73]), and also in other contexts, such as wound-healing [74; 75]. For example,
Letellier et al. [11], Pinho et al. [12] and Yonucu et al. [13] modelled cancer with
ODEs for different cell types (e.g. cancer, normal/host) and vasculature. Angiogene-
sis mechanism was included through changes in the endothelial cells or blood vessel
distribution. Nagy and Armbruster [71] considered differential equations for tumor
mass, immature vascular endothelial cell mass and total microvessel length. They
also included cells’ energy status in form of ODEs for adenosine 5’ mono-, di- and
triphosphatases (AMP, ADP and ATP) to investigate the trade-offs between prolifer-
ation and angiogenesis contribution.
In our model in Publication II, we have included angiogenesis through changes
in resource inflow instead of including an additional cell type. The resources (e.g.
glucose or oxygen) flow into and out of the tumor microenvironment and higher an-
giogenesis contribution by the cells elevates the resource inflow, which further aids
the cell proliferation. However, secreting angiogenic factors (like VEGF) requires
energy, and we included a trade-off between birth rate and angiogenesis contribu-
tion. In our model, angiogenesis was also considered to increase the emigration rate
20
Cancer metapopulation model with angiogenesis
through increased vasculature [76]. However, there have been mixed results on the
relation of angiogenesis and emigration, since cancer angiogenesis usually leads to
distorted veins, through which the access is compromised [77].
3.1.2 Treatment modelling
When modelling cancer angiogenesis, an anti-angiogenic treatment could be in-
cluded to counter the effects. For example, Pinho et al. [12] and Yonucu et al. [13]
modelled the anti-angiogenic treatment with ODEs for therapy agents. The treat-
ment concentration then affected the amount of endothelial cells or vessels. In our
metapopulation model, the effect of anti-angiogenic treatment directly cancels the
effects of angiogenesis within the resource flow and the emigration rate. Since the
original (mature) vessels are assumed to be VEGF independent [78; 79], a baseline
resource inflow will remain even with high levels of anti-angiogenic treatment.
Thus, an anti-angiogenic treatment as a monotherapy does not typically lead to
a full response and a complete eradication requires a combination treatment [17; 80;
81]. Pinho et al. [12] and Yonucu et al. [13] included chemotherapy along with anti-
angiogenic therapy. Also other treatment options, like immunotherapy could be con-
sidered for a combination with anti-angiogenic treatment [82; 83]. Instead of regard-
ing another treatment as an adjuvant therapy for the anti-angiogenic treatment, the
situation could be considered other way around. Anti-angiogenic treatments could
be considered to reduce the need for more toxic treatments, since anti-angiogenic
treatments generally have manageable side effects (e.g., hypertension [34]), while
treatments like chemotherapy may have severe side effects. However, despite being
generally well tolerated, the anti-angiogenic treatment may also lead to increased risk
of bleeding complications and gastrointestinal perforation [84; 85]. In our model, we
included a cytotoxic treatment to investigate how the combination affected the out-
come compared to monotherapies.
A combination treatment may be required also since the cancer cells may develop
resistance to the anti-angiogenic treatment [70; 33]. We considered only a single
angiogenesis mechanisms (like VEGF secretion), but other proangiogenic factors
like fibroblast growth factors (FGFs) and platelet derived growth factors (PDGFs)
[33] could be used to promote angiogenesis. We included the resistance to anti-
angiogenic treatment by preventing the treatment effect on resistant cells.
3.2 Model
We constructed a metapopulation model for cancer with habitat patches and dispersal
pool (see Fig. 11). The state and living conditions in a patch are determined by the
resource inflow through vasculature and the cells living in the patch. The patches
have no direct interactions, and the cells immigrate only through the dispersal pool
21
Anni S. Halkola
(circulatory system). We omitted the spatial relation of patches and the distance is
not considered in the immigration of cells. However, the population size of the new
patch candidate affects the immigration probability, since the cells are more likely to
immigrate into a less-populated patch.
3.2.1 Within-patch dynamics
Dynamics within a patch depend on the resource availability and the cancer cells.
The resource dynamics are assumed to be fast compared to the cancer cell dynamics.
Resources
As in the Publication I, we included the resources in our metapopulation model in
Publication II. Resources (𝑅) (e.g. glucose) follow the chemostat dynamics [36]
when flowing into and out of the patch. The constant baseline inflow ?^? was in-
cluded with possibility of angiogenesis-based increase. Angiogenesis-contributing
cells enhance the vasculature, and the resource inflow by the factor
𝑓(𝑛,𝑎) =
(︃
1 +
𝐴max
∑︀𝑘
𝑖=1 𝑛𝑖𝑎𝑖
1 +
∑︀𝑘
𝑖=1 𝑛𝑖𝑎𝑖
)︃
, (7)
where 𝐴max is the upper bound, based on the limited space for the new vessels.
Vectors 𝑛 = (𝑛1, ..., 𝑛𝑘) and 𝑎 = (𝑎1, ..., 𝑎𝑘), where 𝑛𝑖 is the number of type 𝑖
cells and 𝑎𝑖 is the corresponding angiogenesis strategy (i.e., how much a type 𝑖 cell
contributes to angiogenesis). Cells collect resources by law of mass action with their
resource consumption rate 𝑠𝑖 [37]. Thus, including the outflow concentration 𝑅 and
the relative flow speed 𝜆, local resources follow the following differential equation:
𝑑𝑅
𝑑𝑡
= 𝜆
(︁
𝑓 (𝑛,𝑎) ?^?−𝑅
)︁
−𝑅
𝑘∑︁
𝑖=1
𝑛𝑖𝑠𝑖. (8)
However, if the anti-angiogenic treatment is administrated, it prevents the angiogen-
esis and the increase in resource concentration. Cancer cells may develop resistance
𝜌𝑖 to anti-angiogenic treatment [70; 33]. Resistant cells (0 < 𝜌𝑖 ≤ 1) decrease
the effect of anti-angiogenic treatment on their angiogenesis contribution. Modified
Equation (7) includes the anti-angiogenic treatment 𝜓𝑎 and the treatment resistance
𝜌𝑖:
𝑓(𝑛,𝑎,𝜌, 𝜓𝑎) =
(︃
1 + ℎ(𝜓𝑎)
𝐴max
∑︀𝑘
𝑖=1 𝑛𝑖𝑎𝑖𝑙(𝜌𝑖, 𝜓𝑎)
1 +
∑︀𝑘
𝑖=1 𝑛𝑖𝑎𝑖
)︃
, (9)
where 𝜌 = (𝜌1, ..., 𝜌𝑘) are the resistance strategies of cell types 1, ..., 𝑘. The treat-
ment effect is
ℎ(𝜓𝑎) =
1
1 + 𝜓𝑎
, (10)
22
Cancer metapopulation model with angiogenesis
where 𝜓𝑎 is the concentration of anti-angiogenic treatment. The resistance effect
𝑙(𝜌𝑖, 𝜓𝑎) = 1 + 𝜓𝑎𝜌𝑖 (11)
is determined by the treatment function ℎ(𝜓𝑎) and cancels the treatment effect only
on cell type 𝑖. Thus, if the cell type 𝑖 is resistant and the cell type 𝑗 is sensitive, only
the contribution of type 𝑖 transpires. It is assumed that the resource dynamics are
fast, and in the equilibrium (𝑑𝑅/𝑑𝑡 = 0) the resource concentration is
𝑅*(𝑛,𝑎) =
𝜆𝑓 (𝑛,𝑎) ?^?
𝜆+
∑︀𝑘
𝑖=1 𝑛𝑖𝑠𝑖
. (12)
Cancer cells
The within-patch cancer cell dynamics were modelled as a continuous-time Markov
chain (see Fig. 12 for a monomorphic population), in which the population’s state is
determined by the numbers of cell types and the state changes based on the cell pro-
liferation, death and migration. It was assumed that a patch could have a maximum
of 𝐾 cells and no division would happen in a full patch.
0 1 2 K
Immigration
Cell division or
immigration
Cell division or
immigration
Death or
emigration
Death or
emigration
Death or
emigration
Catastrophe
Figure 12. The states and transitions of the Markov chain describing monomorphic (i.e., all cells
are of the same type) resident population dynamics. (Adopted with permission from Publication II:
Fig. 2.)
Cells collect resources and proliferate based on the resource intake, thus the birth
rate of cell type 𝑖 is
𝑏𝑛,𝑎,𝑖 = 𝛾𝑠𝑖𝑅
*(𝑛,𝑎)𝑔(𝑎𝑖, 𝜌𝑖), (13)
where 𝑅* is the resource concentration in equilibrium, 𝑠𝑖 is the resource consump-
tion rate and 𝛾 is the conversion coefficient from resource intake to cell division.
Since angiogenesis contribution reserves resources from proliferation, the trade-off
23
Anni S. Halkola
𝑔(𝑎𝑖, 𝜌𝑖) was included. In addition, it was assumed that acquiring and maintaining
resistance (𝜌𝑖 > 0) requires energy and has a trade-off on proliferation. The trade-off
is
𝑔(𝑎𝑖, 𝜌𝑖) =
1
1 + 𝑎𝑖 + 𝛽𝜌𝑖
, (14)
where, 𝛽 is the trade-off factor for resistance. The trade-off factor of angiogenesis
was set to 1.
The cells die with the rate
𝑑𝑛𝑇 = 𝑑0 + 𝛿𝑛𝑇 + 𝜓𝑐, (15)
where 𝑑0 is the baseline death rate, 𝜓𝑐 is the effect of cytotoxic treatment and 𝑛𝑇 =∑︀𝑘
𝑖=1 𝑛𝑖 is the total number of cells in the patch. The factor 𝛿 = 1/𝐾 is positive, as
it was assumed that a higher population size leads to a higher death rate [86]. Cells
emigrate with the rate
𝑞𝑛,𝑎,𝑒𝑖 = 𝑒𝑖𝑓 (𝑛,𝑎) , (16)
where 𝑒𝑖 is the emigration strategy (i.e., how easily cell leaves the patch). The factor
𝑓 is as in Equation (9) since it was assumed that the changes in vasculature affect the
emigration probability [76; 87; 88].
3.2.2 Metapopulation dynamics
The state of a patch depends on the numbers of cells in the patch. A metapopula-
tion’s state is investigated by calculating the distribution of the cell numbers. In a
monomorphic population the probabilities 𝑝𝑛, i.e., the probabilities that a randomly
chosen patch has 𝑛 cells, satisfy the differential equations (forward Kolmogorov
equations):
𝑑𝑝0(𝑡)
𝑑𝑡
= −𝛼𝐷𝑆0𝑝0(𝑡) + [𝑑1 + 𝑞1,𝑎,𝑒] 𝑝1(𝑡) + 𝜇
𝐾∑︁
𝑛=1
𝑝𝑛(𝑡)
𝑑𝑝𝑛(𝑡)
𝑑𝑡
= − [𝛼𝐷𝑆𝑛 + 𝑛(𝑏𝑛,𝑎 + 𝑑𝑛 + 𝑞𝑛,𝑎,𝑒) + 𝜇] 𝑝𝑛(𝑡)
+ [𝛼𝐷𝑆𝑛−1 + (𝑛− 1)𝑏𝑛−1,𝑎] 𝑝𝑛−1(𝑡) + (𝑛+ 1) [𝑑𝑛+1 + 𝑞𝑛+1,𝑎,𝑒] 𝑝𝑛+1(𝑡)
𝑑𝑝𝐾(𝑡)
𝑑𝑡
= − [𝐾(𝑑𝐾 + 𝑞𝐾,𝑎,𝑒) + 𝜇] 𝑝𝐾(𝑡) + [𝛼𝐷𝑆𝐾−1 + (𝐾 − 1)𝑏𝐾−1,𝑎] 𝑝𝐾−1(𝑡),
(17)
where the dispersal pool size 𝐷 satisfies the differential equation:
𝑑𝐷(𝑡)
𝑑𝑡
= −
(︃
𝛼
𝐾∑︁
𝑛=0
𝑝𝑛(𝑡)𝑆𝑛 + 𝜈
)︃
𝐷(𝑡) +
𝐾∑︁
𝑛=1
𝑞𝑛,𝑎,𝑒𝑛𝑝𝑛(𝑡). (18)
24
Cancer metapopulation model with angiogenesis
Here 𝜈 is the death rate of cancer cells in the dispersal pool. Cells encounter patches
with the rate 𝛼 and immigrate with the probability 𝑆𝑛 = (𝐾−𝑛)/𝐾 upon encounter.
Thus, the cells are more likely to immigrate into less inhabited patches.
In a metapopulation-dynamical steady state
𝑑𝑝𝑛(𝑡)
𝑑𝑡
= 0 and
𝑑𝐷(𝑡)
𝑑𝑡
= 0. How-
ever, since 𝑝𝑛(𝑡) are probabilities, they must satisfy
∑︀𝐾
𝑛=0 𝑝𝑛(𝑡) = 1, and thus one
equation in (17) needs to be replaced by 𝑝0(𝑡)+𝑝1(𝑡)+ · · ·+𝑝𝐾(𝑡) = 1. In a steady
state, 𝑝𝑛(𝑡) and 𝐷(𝑡) do not depend on time, so that 𝑝𝑛(𝑡) = 𝑝𝑛 and 𝐷(𝑡) = ?¯?.
From the linear system of equations
0 = −𝛼𝐷𝑆0𝑝0 + [𝑑1 + 𝑞1,𝑎,𝑒] 𝑝1 + 𝜇(1− 𝑝0)
0 = − [𝛼𝐷𝑆𝑛 + 𝑛(𝑏𝑛,𝑎 + 𝑑𝑛 + 𝑞𝑛,𝑎,𝑒) + 𝜇] 𝑝𝑛
+ [𝛼𝐷𝑆𝑛−1 + (𝑛− 1)𝑏𝑛−1,𝑎] 𝑝𝑛−1 + (𝑛+ 1) [𝑑𝑛+1 + 𝑞𝑛+1,𝑎,𝑒] 𝑝𝑛+1
1 = 𝑝0 + 𝑝1 + · · ·+ 𝑝𝐾 .
(19)
The probabilities 𝑝𝑛 can be solved for a fixed value of 𝐷. Based on Equation (18)
the equilibrium value ?¯? must then satisfy
1 =
𝐾∑︀
𝑛=1
𝑞𝑛,𝑎,𝑒𝑛𝑝𝑛(?¯?)(︂
𝛼
𝐾∑︀
𝑛=0
𝑝𝑛(?¯?)𝑆𝑛 + 𝜈
)︂
?¯?
, (20)
where 𝑝𝑛(𝐷) denote the solution of (19). Solutions 𝑝𝑛 and ?¯? were searched numer-
ically.
3.2.3 Evolutionary dynamics
Cell’s strategies variate randomly, forming new variant cells that then experience the
environment set by the resident population. The variant’s invasion fitness is the long-
term exponential growth rate in this environment [89]. Only a variant with a positive
invasion fitness may invade and become a new resident. However, calculating the in-
vasion fitness is often complicated and, instead, the metapopulation fitness 𝑅metapop
[90; 59; 58] is calculated to determine the sign of the invasion fitness.
From an initially small variant population in the dispersal pool, a focal variant
may immigrate into a patch. The variant and all of its descendants form a variant
colony in the patch and, due to the variants rarity, no other variants are expected to
arrive during the lifetime of the variant colony. The variant colony sends emigrants
to the dispersal pool and the variant’s metapopulation fitness is the expected number
of variant emigrants produced by the variant colony during its lifetime (i.e. the time
before the variant colony goes extinct).
25
Anni S. Halkola
The state of the variant colony is determined by the number of variants 𝑛𝑚 and
residents 𝑛𝑟, which change according to a Markov chain. When the focal variant
immigrates, there is exactly one variant (𝑛𝑚 = 1). The number of residents depends
on chance, based on the resident population distribution (𝑝𝑛). The initial probability
distribution of the Markov chain is
𝛿0(𝑛𝑟, 1) =
𝑝𝑛𝑟𝑆𝑛𝑟∑︀𝐾−1
𝑛=0 𝑝𝑛𝑆𝑛
, for 0 ⩽ 𝑛𝑟 < 𝐾
𝛿0(𝑛𝑟, 𝑛𝑚) = 0, for 𝑛𝑚 > 1.
(21)
The transition intensities from the state (𝑛𝑟, 𝑛𝑚) to neighboring states are
𝑐0,+𝑛𝑟,𝑛𝑚 = 𝑛𝑚𝑏𝑛,𝑎,𝑚, if 𝑛𝑇 = 𝑛𝑟 + 𝑛𝑚 < 𝐾
𝑐0,−𝑛𝑟,𝑛𝑚 = 𝑛𝑚(𝑑𝑛𝑇 + 𝑞𝑛,𝑎,𝑚)
𝑐+,0𝑛𝑟,𝑛𝑚 = 𝛼𝐷𝑆𝑛𝑇 + 𝑛𝑟𝑏𝑛,𝑎,𝑟, if 𝑛𝑇 = 𝑛𝑟 + 𝑛𝑚 < 𝐾
𝑐−,0𝑛𝑟,𝑛𝑚 = 𝑛𝑟(𝑑𝑛𝑇 + 𝑞𝑛,𝑎,𝑟).
(22)
The subscript denotes the original state and the superscript indicates the change
in corresponding numbers of individuals. The transition intensity out of the state
(𝑛𝑟,𝑛𝑚) is
𝑐𝑛𝑟,𝑛𝑚 = 𝑐
0,+
𝑛𝑟,𝑛𝑚 + 𝑐
0,−
𝑛𝑟,𝑛𝑚 + 𝑐
+,0
𝑛𝑟,𝑛𝑚 + 𝑐
−,0
𝑛𝑟,𝑛𝑚 + 𝜇, (23)
where 𝜇 is the rate of catastrophes, which erase the local population entirely (e.g.
surgery). When a transition happens, the probability that the population changes
from the state (𝑛𝑟, 𝑛𝑚) to the state (𝑛𝑟, 𝑛𝑚 + 1) is 𝑐
0,+
𝑛𝑟,𝑛𝑚/𝑐𝑛𝑟,𝑛𝑚 . Other transi-
tion probabilities are formed correspondingly and they are collected to the transition
probability matrix 𝑃 .
The expected number of visits 𝑤 in each state are calculated from [91]:
𝑤 = 𝛿0(𝐼 − 𝑃 )−1, (24)
where 𝐼 is the identity matrix and the matrix inverse exists for transient states. The
fitness is calculated from
𝑅metapop =
𝛼𝑚
∑︀𝐾−1
𝑛=0 𝑝𝑛𝑆𝑛
𝜈 + 𝛼𝑚
∑︀𝐾−1
𝑛=0 𝑝𝑛𝑆𝑛
∑︁ 𝑛𝑚𝑞𝑚𝑤𝑛𝑟,𝑛𝑚
𝑐𝑛𝑟,𝑛𝑚
, (25)
where 𝑞𝑚 is the emigration rate of variants and 𝑤𝑛𝑟,𝑛𝑚 corresponds to the state
(𝑛𝑟, 𝑛𝑚) in the vector 𝑤. The invasion fitness is positive if and only if 𝑅metapop > 1.
Evolutionary dynamics are illustrated using pairwise invasibility plots (PIPs), in
which variant’s fitness is marked either negative or positive (𝑅metapop < 1 or
𝑅metapop > 1 correspondingly) with respect to resident’s and variant’s strategies.
Curves on which 𝑅metapop = 1, are called neutral contours. For example, on the
26
Cancer metapopulation model with angiogenesis
diagonal resident’s and variant’s strategies are equal and thus the fitness is neutral.
To investigate where the strategy (or strategies) would evolve, a variant’s fitness
gradient is calculated since a variant’s fitness increases the most in the direction given
by the fitness gradient. In reality, the strategies do not evolve straight along the gra-
dient, however, in that direction the invasion probability is higher and over time the
strategies will evolve to the gradient’s direction. For a multi-dimensional strategy
(strategy vector) the fitness gradient is
𝐷(𝑎𝑟, 𝑒𝑟, 𝜌𝑟) =(︂
𝜕
𝜕𝑎𝑚
𝑅metapop,
𝜕
𝜕𝑒𝑚
𝑅metapop,
𝜕
𝜕𝜌𝑚
𝑅metapop
)︂ ⃒⃒⃒⃒
⃒
𝑎𝑚=𝑎𝑟,𝑒𝑚=𝑒𝑟,𝜌𝑚=𝜌𝑟
.
(26)
In case of a one-dimensional strategy, the fitness function (25) is differentiated only
with respect to that strategy, after which the variant’s strategy is set to be equal to
resident’s strategy. A point where all the components of the fitness gradient become
zero is called evolutionarily singular strategy (or just singular strategy). In a PIP,
the singular strategies are positioned at the intersections of the diagonal and other
neutral contours. A singular strategy is uninvadable (evolutionarily stable) if there is
no strategy that could make an invasion. The opposite holds for invadable singular
strategy. A singular strategy is attracting, if near the point smaller values have a
positive gradient and higher values have a negative gradient. If the opposite holds,
the singular strategy is repelling.
In addition to analyzing the singular strategies, we also monitored their effect on
the population level by calculating the average population size
∑︀𝐾
𝑛=0 𝑛𝑝𝑛 and the
average emigration
∑︀𝐾
𝑛=0 𝑝𝑛𝑒𝑓(𝑛, 𝑎).
3.3 Results
In the main analyses, we focused on the strategies (angiogenesis, emigration and re-
sistance) and the treatment parameters’ effect on them. First, we investigated the
three strategies (angiogenesis, emigration and resistance) separately, while keep-
ing the remaining two as constants. Joint evolution of angiogenesis and emigration
strategies was analyzed with anti-angiogenic and cytotoxic monotherapies and with
the combination therapy. Finally, the joint evolution of all three strategies was inves-
tigated under anti-angiogenic treatment.
Parameter values were selected so that the metapopulation was viable. The se-
lected parameter values were further justified by comparing the observed angiogene-
sis dynamics qualitatively to real-life observations [17; 66; 80; 67; 81]. The treatment
parameters were increased from zero (no treatment) until the population became non-
viable, the angiogenesis strategy became zero or the main trend was detected. Other
27
Anni S. Halkola
parameters were mainly kept constant, however, additional analyses were conducted
also on their changes and effects on the strategy dynamics (Publication II Supple-
mentary Material). It was noted, that especially the joint dynamics of all three strate-
gies were sensitive to the changes in parameter values, indicating the importance of
patient-specific biology and need for personalized treatments.
3.3.1 Single strategies
Angiogenesis increased the birth rate indirectly through resources, however, there
was also the direct cost on the birth rate. The relation between benefit (resources) and
costs (trade-off) determined whether a positive singular strategy existed (Fig. 13a).
Despite the cost for an investing cell, the resource increase benefited the whole patch,
including the cell’s kin. Therefore, kin selection [92] also played a role in the angio-
genesis evolution. If the singular angiogenesis strategy was detected, it was attract-
ing and uninvadable. However, the attracting singular strategy did not necessarily
maximize the average population size because the best for individual might differ
from the best for the population. In the case anti-angiogenic treatment was added,
the angiogenesis strategy would decrease (Fig. 13b), eventually becoming zero if the
treatment was high enough.
Emigration evolution was also affected by the kin selection, since by emigrat-
ing, the cells left more resources for their relatives (and other cells) in the original
patch. The emigrating cells might also directly benefit from the emigration, if the
new patch was better than the original. Thus, a positive singular emigration strat-
egy was detected (Fig. 13c). The existence of a positive singular strategy was also
reasonable, because the positive catastrophe rate would cause non-viability if the em-
igration strategy was close to zero (narrow black bars in Figs. 13c and d). Similarly
to angiogenesis, the singular emigration strategy did not necessarily maximize the
average population size. Addition of anti-angiogenic treatment increased the singu-
lar emigration strategy 13d).
Without treatment, the treatment resistance had only a negative effect on pro-
liferation and the strategy evolved to zero (Fig. 13e). With treatment, the situation
shifted and a positive singular strategy might exist (Fig. 13f), however, this depended
on the level of treatment. We mainly analyzed the resistance strategy together with
angiogenesis and emigration strategies (see Section 3.3.3).
3.3.2 Joint evolution of angiogenesis and emigration strategies
We investigated the joint evolution of angiogenesis and emigration strategies, i.e.
the both were allowed to variate while the resistance strategy was kept zero. We
observed that increasing one strategy component made the other one less favorable
28
Cancer metapopulation model with angiogenesis
a) Angiogenesis PIP - no treatment b) Angiogenesis PIP - treatment 
c) Emigration PIP - no treatment d) Emigration PIP - treatment 
e) Resistance PIP - no treatment f) Resistance PIP - treatment 
Figure 13. Pairwise invasibility plots (PIP) for a) angiogenesis strategy without treatment and b)
with anti-angiogenic treatment. c) PIP for emigration strategy without treatment and d) with anti-
angiogenic treatment. e) and f) corresponding PIPs for resistance strategy. In each, the other
strategies are kept constant 𝑒 = 0.2, 𝑎 = 0.2 and 𝜌 = 0. The anti-angiogenic treatment 𝜓𝑎 = 0 for
left panel and 𝜓𝑎 = 7 for the right panel.
29
Anni S. Halkola
for the cancer cells. The joint singular strategy was numerically calculated from the
condition 𝐷(𝑎𝑟, 𝑒𝑟) = (0, 0). The singular strategy did not necessarily maximize
the average population size, as was observed in the case of a single strategy. Since
the parameters defined the conditions for the strategy dynamics, we investigated the
parameters’ effect on the singular strategies. Especially we focused on the treatment
parameters.
Monotherapies
As the anti-angiogenic treatment decreased the benefits of angiogenesis, an increase
in the treatment caused the singular angiogenesis strategy to decrease. Consequently,
the emigration was compromised and initially the average emigration rate decreased.
However, the singular emigration strategy increased with the anti-angiogenic treat-
ment and the decrease in the average emigration rate abated. The singular angio-
genesis strategy became zero, if the concentration of anti-angiogenic treatment was
increased enough. The increase in the emigration strategy also abated since the effect
of the treatment disappeared in the absence of angiogenesis.
The anti-angiogenic treatment did not necessarily lead to non-viability since the
baseline resource inflow ?^? kept the population viable. However, if the population
was originally non-viable (insufficient ?^?), increasing the anti-angiogenic treatment
enough would lead the population back to non-viability.
Cytotoxic treatment increased the death rate. Consequently the populations were
smaller and cells were more likely related. In such situation both angiogenesis and
emigration were favored based on kin selection. However, as noted in the begin-
ning of Section 3.3.2 the increase in one strategy component decreased the other.
In the joint evolution, the indirect effects could cause such phenomenon when in-
creasing the cytotoxic treatment. However, the direct effects of cytotoxic treatment
dominated, and both strategy components increased. The average emigration rate
initially increased along with emigration and angiogenesis strategies, but started to
decrease and approach the emigration strategy. This happened since the increase in
treatment decreased the average population size and eventually there were no cells
left to contribute to angiogenesis. Thus, regardless of the angiogenesis and emigra-
tion strategies, the population became non-viable when the cytotoxic treatment was
high enough.
Combination of anti-angiogenic and cytotoxic treatments
Both monotherapies have their downsides, as anti-angiogenic treatment does not
lead to non-viability and cytotoxic treatment would be preferred at minimal effec-
tive doses due to side effects. Therefore, we investigated whether the combination
30
Cancer metapopulation model with angiogenesis
treatment would overcome these downsides.
Similarly to the cytotoxic monotherapy, the population eventually became non-
viable when cytotoxic treatment was increased despite the increase in the angiogen-
esis and emigration strategies (Fig. 14). Compared with the cytotoxic monother-
apy (Fig. 14a) the combination therapy (Fig. 14b) lowered the need of cytotoxic
treatment as the non-viability was reached with smaller 𝜓𝑐. The concentration of
anti-angiogenic treatment affected how much the need for cytotoxic treatment de-
creased. However, increasing the anti-angiogenic treatment started to lose its benefit
and almost equal result could already be achieved with lower level of anti-angiogenic
treatment. A lower treatment level would also be preferred considering the additional
costs of combination therapy. See also Section 3.3.3 for the emergence of treatment
resistance.
Cytotoxic treatment, ψc
An
gi
og
en
es
is 
st
ra
te
gy
,
 
a
r
0 1 2 3 4 5
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
a) No anti−angiogenic treatment
Cytotoxic treatment, ψc
An
gi
og
en
es
is 
st
ra
te
gy
,
 
a
r
0 1 2 3 4 5
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
b) Anti−angiogenic treatment
Figure 14. Angiogenesis strategy (blue line) at the joint singular strategy of angiogenesis and em-
igration with respect to cytotoxic treatment a) without anti-angiogenic treatment and b) with anti-
angiogenic treatment. Black area denotes non-viability. (Adopted with permission from Publication
II: Fig. 7a and b.)
3.3.3 Joint evolution with treatment resistance
We considered the joint evolution of treatment resistance, angiogenesis and emi-
gration strategies. Since the resistance was modelled against the anti-angiogenic
treatment, we focused on how that treatment affected the joint singular strategy and
especially the treatment resistance.
Bistability scenario
For the selected parameter values, the number of singular strategies depended on the
concentration of anti-angiogenic treatment (Fig. 15a). For small concentrations, a
single attracting singular strategy existed with no treatment resistance. An additional
31
Anni S. Halkola
singular strategy with positive resistance was detected for higher treatment levels.
Between the two attracting singular strategies, there was a saddle point. The stable
manifold of this saddle point determined the areas in the strategy space from which
the strategies evolved to each attracting singular strategy. The maximum treatment
effect without a high risk of resistance would be before the positive resistance singu-
lar strategy appears (around 3.8 in Fig. 15a).
We further investigated a situation, where the resident population was not treat-
ment resistant despite the anti-angiogenic treatment (𝜓𝑎 = 7). The angiogenesis
and emigration strategies were at their joint singular strategy corresponding to the
selected treatment. For such resident, the resistance strategy component had a neg-
ative fitness gradient (black dot and arrow in Fig.15b). Since the resistance strategy
had to be non-negative, it was likely to remain in zero. However, some variants with
multiple radical variations could invade (area of positive fitness in Fig. 15b).
0 2 4 6 8 10
0.0
0.2
0.4
0.6
Anti−angiogenic treatment, ψa
R
es
is
ta
nc
e 
st
ra
te
gy
,
 
ρ r
a) Resistance strategy
R
es
is
ta
nc
e 
st
ra
te
gy
,
 
ρ m
Angiogenesis strategy, am
l0
0.12
0.25
0.38
0.5
0.03 0.07 0.1
0.95
1.00
0.85
1.05
_
b) Variant's fitness
Figure 15. a) Resistance strategy at the joint singular strategy of all three strategies. The solid
lines are attracting singular strategies, the dashed is repelling singular strategy. b) Slice of the 3D
fitness matrix calculated with respect to resident’s strategy (black dot). Here 𝑒𝑚 = 𝑒𝑟. The arrow
points the direction of fitness gradient from the resident’s strategy. The black lines are contour lines
for different fitness values. (Panel a adopted with permission from Publication II: Fig. 8a.)
The effect of cytotoxic treatment
As noted in the previous section, the anti-angiogenic treatment did not necessarily
lead to treatment resistance, even if the concentration was increased. If the com-
bination of anti-angiogenic and cytotoxic treatment was used, the bistability of sin-
gular strategies was present only for a narrow range of treatment regimens. The
singular strategy with zero resistance was repelling for a wide range of treatment
combinations leading to the emergence of treatment resistance. Compared to the
anti-angiogenic monotherapy, a lower average population was reached with the com-
bination, despite the higher risk of treatment resistance. However, as in the previous
32
Cancer metapopulation model with angiogenesis
sections, the combination treatment might come with additional burden for the pa-
tient.
In addition to treatment parameters, the bistability with anti-angiogenic monother-
apy shifted or disappeared when other parameters were changed. For example, when
decreasing the resistance trade-off on birth rate (𝛽), resistance strategy became more
profitable for cancer cells and, therefore, a positive resistance was attained even with
smaller treatment concentration.
33
4 Cancer survival prediction based on
clinical data
In Publication III, cancer survival was predicted with two methods: LASSO selection
and a greedy budget-constrained Cox regression (Greedy Cox). In Publication IV, a
new method OSCAR (Optimal Subset CArdinality Regression) was presented for
penalized (Cox) regression. In both Publications the aim was to provide new tools
for better prediction while simultaneously considering the model sparsity. The model
sparsity was also desired, since in real life applications the cost of measuring the
parameter values may considerably increase along with the number of predictors. In
both Publications III and IV, the cost aspect was considered with real life data from
advanced prostate cancer.
4.1 Literature review
Prognostic modelling is used to aid the prediction of cancer survival and assessing
who may benefit from a certain treatment regimen. Fitting the model into a data set
results in model coefficient values that can then be used to form predictions for new
patient data. An additional data can be used to further train the model to better fit a
target population. For survival prediction, Cox regression [9] has been widely used
as the base of modelling. In the Cox’s proportional hazards model, the hazard is
ℎ(𝑡,𝑥) = ℎ0(𝑡)𝑒
𝛽⊤𝑥, (27)
where ℎ0(𝑡) is a shared baseline hazard, 𝑥 is a feature vector and 𝛽 is the coefficient
vector.
Furthermore, to prevent over-fitting and enable model sparsity, penalized Cox
regression [93; 94] has been applied. For example, in the DREAM 9.5 Challenge for
survival prediction of metastatic castration-resistant prostate cancer (mCRPC), the
top model was based on an ensemble-based penalized Cox regression model (ePCR,
[95]). In addition to penalized methods, different approaches, such as other statisti-
cal, Bayesian, neural network or random forest methods have also been utilized to
ensure model sparsity (see e.g. [96; 97; 98; 99; 100; 101; 102]).
For penalized Cox regression, most commonly used penalization terms are 𝐿1-
norm (LASSO), 𝐿2-norm (ridge regression) or their sum (elastic net) [103]. Both 𝐿1
and 𝐿2 push coefficients (𝛽𝑖) towards zero, however, 𝐿1 selects only some nonzero
34
Cancer survival prediction based on clinical data
coefficients, while𝐿2 tends to set none to exactly zero. 𝐿1 has difficulties with highly
correlated parameters, selecting only one of them. 𝐿2 tends to give correlated param-
eters an equal coefficient value. Elastic net combines the two so that the penalization
to be minimized is
𝜆
(︁
𝛼
∑︁
|𝛽𝑖|+ (1− 𝛼)
∑︁
𝛽2
)︁
, (28)
where the first term is the effect of 𝐿1 and the latter is the effect of 𝐿2. Here 𝛼 de-
termines the trade-off between 𝐿1 and 𝐿2, and 𝜆 can be altered to tune the amount
of penalization. For example, Halabi et al. [99] and Goeman [104] utilized LASSO-
penalty, while Simon et al. [94] and Guinney et al. [95] considered elastic net with
varying 𝛼. Also Meier et al. [100] tested between the three penalties. In Publication
III, both 𝐿1 (LASSO selection) and 𝐿2 (Greedy Cox) were included.
However, another option for 𝐿1 and 𝐿2 is to restrict the number of predictors
with 𝐿0-pseudonorm, which calculates the number of nonzero coefficients. Unlike
𝐿1 and 𝐿2, 𝐿0 is not considered a proper norm since it is non-homogeneous [105].
In addition, it is discontinuous and noncovex, which leads to challenging, NP-hard
(nondeterministic polynomial hard) optimization problem [106; 107; 108; 109]. This
has led to under-presentation of 𝐿0 in the penalized regression. A few 𝐿0-approaches
are implemented into R-packages, that include generalized linear models with linear
and logistic regression [110; 111]. A model, augmented penalized minimization-𝐿0
(APM-𝐿0) [112] augments 𝐿0 for Cox’s proportional hazards model. However, the
𝐿0 is approached indirectly through minimizing the 𝐿0-pseudonorm of surrogate co-
efficients, which are compared to the actual coefficients used in the hazards model.
In Publication IV, the 𝐿0-pseudonorm was used directly in the optimization of model
coefficients. An exact, continuous representation of 𝐿0 [107] was included as a car-
dinality constraint that limits the nonzero coefficients.
In addition to differences in the penalty functions, the actual model fitting has been
done with different methods. One popular approach is to use coordinate descent
in the optimizing the selection of coefficient values. Coordinate descent has been
utilized by, for example, Simon et al. [94] and Li et al. [112]. Goeman [104] com-
bined gradient ascent and the Newton-Raphson algorithm. Our approach with 𝐿0-
pseudonorm resulted in a continuous exact representation of the penalization term,
which is also a DC function (difference of two convex functions). The DC structure
may be exploited in the optimization and thus we used DBDC (the double bundle
method) for DC optimization [113; 114]. DBDC suits the problem also since the
𝐿0 representation is still nonconvex. Another optimization algorithm LMBM (the
limited memory bundle method) [115; 116] was included as an option for larger
problems, however, LMBM does not utilize DC structure.
In Publication III and IV, the presented methods were applied on prostate cancer data.
35
Anni S. Halkola
As the prostate cancer is one of the most common cancers among men [1; 117], it is a
relevant subject for prognostic prediction. Furthermore, the need for cost-efficiency
of prediction was considered in both Publications, since the application of models
requires money. Some parameter values are also often measured as groups or kits
like basic blood count, which includes hemoglobin, white blood cells and red blood
cells, among others. In both Publications III and IV, the kit structure was included.
4.2 Materials and methods
Concordance index (C-index) [118] was used as the measure of model performance
in both Publications III and IV. In Publication III, also integrated AUC [119] was
calculated for supplementary analyses to enable comparison to previous research.
Cross-validation was used to further evaluate the model selection. For Publication
IV also bootstrap analysis was performed to interpret which features were the most
robust.
4.2.1 Data
We tested our survival prediction methods in four prostate cancer cohorts. Three data
cohorts, MAINSAIL, VENICE and ASCENT [120; 121; 122], were from random-
ized clinical trials and one data cohort was from real-world hospital registry data.
The clinical trial data were constructed in the DREAM competition (the Prostate
Cancer Challenge, PCC-DREAM), hosted by Project Data Sphere [123]. In Publica-
tion IV the trial cohorts were divided into validation (N=132, N=150 and N=119 for
MAINSAIL, VENICE and ASCENT respectively) and training sets (N=394, N=448
and N=357 respectively). Missing values were imputed separately for each cohort
using the median values calculated from the training sets.
The real-world hospital registry data were collected from patients with castration re-
sistant prostate cancer treated at the Turku University Hospital (TYKS). The data
was processed as in [124]. For Publication III, only variables with less than 40% of
missing values were selected. The missing values were imputed using medians as
the median imputation has previously given satisfactory results [124]. The k-nearest
neighbor imputation was also tested, however it did not improve model performance.
For Publication IV, the patients with diagnosis date before 2010 were discarded
due to high percentage of missing values. Unlike in Publication III, only laboratory
measurement data and the age group information were included. Furthermore, the
features with more than 50% of missing values were discarded. Similarly to the trials
data, TYKS cohort was divided into validation set (N=195) and training set (N=590).
The remaining missing values were imputed using medians calculated from the train-
ing data.
36
Cancer survival prediction based on clinical data
In both Publications III and IV the features were also considered as groups or
kits in which the features are measured together. Prices for the examinations were
provided by the Helsinki University Hospital. The prices were scaled in relation to
PSA, which was given a reference value of 100.
For both Publication III and IV the initially same four data cohorts were tested
in survival prediction, however, the data were processed slightly differently and dif-
ferent features were included. In addition, for Publication IV, the TYKS data was
updated to include new patient data that had been collected after the analyses per-
formed for Publication III.
4.2.2 LASSO selection and Greedy Cox methods
In Publication III the approaches were based on Cox regression [93] with either 𝐿1
(LASSO) or 𝐿2 (ridge regression) penalization. Penalization was varied to determine
the sparsity and select different sets of variables. In the Cox’s proportional hazards
model (27), the coefficients 𝛽 were selected by maximizing the penalized partial log
likelihood [103].
We denoted the set of available variables with 𝐹 , and formed subgroups 𝐺𝑖 that
included the variables measured in kit 𝑖. Variables measured on their own were
considered as groups containing one variable. A variable was allowed to belong to
multiple subgroups (e.g., hemoglobin individually or as a part of blood count). A
group was also allowed to be included partially, with only some of the group vari-
ables in the model. Selected model variables were denoted by 𝑆 =
⋃︀
𝑖∈𝐼 𝑠𝑖, where
𝐼 is the group indices included in the model and 𝑠𝑖 are the variables selected from
the group 𝐺𝑖. The price of the group 𝐺𝑖 was denoted by 𝑐𝑖 and the cost for a model
was calculated by 𝐶 =
∑︀
𝑖∈𝐼 𝑐𝑖. The cost was not reduced by exclusion of some
variables in a selected group since the group was measured as a whole regardless of
how many of its variables were included in the model.
Since the price of measuring the selected variables is not included in the LASSO
method, the price was computed afterwards using a heuristic that added variable
groups sequentially. Groups were added until all selected variables in 𝑆 were in-
cluded. At each step, the group with the smallest cost per new variable was selected.
However, we also proposed an alternative, a greedy budget-constrained Cox re-
gression (Greedy Cox) algorithm. The algorithm sequentially selects the variable
group that, with already selected variables, gives the best cross-validated (threefold)
model performance estimate. Since an added group may have also unfavorable vari-
ables, an inner selection is performed on the variables of the considered group. The
mechanism of inner selection of individual variables is similar to the group selection.
The selection of new groups is stopped if the budget does not fit any of the remain-
ing groups or the addition does not improve prediction performance. Unlike LASSO
with 𝐿1-penalization, Greedy Cox utilizes 𝐿2-penalization.
37
Anni S. Halkola
At the initialization, 𝑆 is initialized as an empty set and the remaining budget is
set to the total budget. The group selection step, followed by inner selection, require
the matrix of clinical variable measurements of the patients and the survival data, i.e.
time and type of event (death or right-censoring). In addition, the input should in-
clude the remaining budget, variable groups 𝐺𝑖 and their prices 𝐶 as well as already
selected variables 𝑆 and the selected group indices 𝐼 . The group and inner selections
return the index 𝑖 of the selected group and the variables 𝑠𝑖 selected from that group.
The remaining budget, 𝑆 and 𝐼 are then updated. The group and inner selections
are called repeatedly until the budget does not fit any new groups that improve the
performance.
4.2.3 OSCAR method
The OSCAR method presented in Publication IV implements the penalized Cox re-
gression [93]. However, unlike LASSO selection and Greedy Cox in Publication
III, the OSCAR method utilizes 𝐿0-pseudonorm in the penalization. In addition to
Cox regression and survival prediction, the OSCAR method includes also logistic
and linear regression for prediction of binomial and scalar responses. Similarly to
the methods in Publication III, the OSCAR method can consider the features either
separately or grouped into kits. The OSCAR method was implemented utilizing For-
tran, C and R and it is available as an R package OSCAR in the Central R Archive
Network (CRAN) at: https://CRAN.R-project.org/package=oscar.
In Cox’s proportional hazards model [9] the hazard for the patient 𝑖 at time 𝑡 is
ℎ𝑖(𝑡) = ℎ0(𝑡)𝑒
𝑥⊤𝑖 𝛽, where ℎ0(𝑡) is a shared baseline hazard, 𝑥𝑖 is the vector of ob-
served features and 𝛽 is the vector of predictor coefficients. In the OSCAR method,
the vector 𝛽 is estimated by maximizing the scaled log partial likelihood, which
leads to an equivalent solution as the Breslow approximation of partial likelihood
[125; 94]. The scaled log partial likelihood is
𝑙(𝛽) =
2
𝑛
𝑚∑︁
𝑖=1
⎧⎨⎩∑︁
𝑗∈𝐷𝑖
𝑥⊤𝑗 𝛽 − 𝑑𝑖 ln
⎛⎝∑︁
𝑗∈𝑅𝑖
𝑒𝑥
⊤
𝑗 𝛽
⎞⎠⎫⎬⎭ , (29)
where 𝑛 is the number of observations and 𝑚 is the number of unique event times.
Furthermore, 𝑅𝑖 is the set of patient indices at risk and 𝐷𝑖 is the set of indices with
the label equal to one at time 𝑡𝑖. As the label equal to one indicates an event (death)
and zero indicates right-censoring, 𝑑𝑖 = |𝐷𝑖| is the number of events at time 𝑡𝑖.
Maximizing the concave function 𝑙 is the same as minimizing the convex function
−𝑙 [126].
The minimization problem is modified further to include the regularization and
enable model sparsity. The 𝐿0-pseudonorm is included as a cardinality constraint
38
Cancer survival prediction based on clinical data
that restricts the number of nonzero predictors in 𝛽. Unlike in LASSO, the values of
nonzero coefficients in 𝛽 are not penalized and are allowed to vary freely.
The 𝐿0-pseudonorm ‖𝛽‖0 calculates the number of nonzero components, and
it gives the cardinality constraint ‖𝛽‖0 ≤ 𝐾. However, since the constraint is dis-
continuous, we used an approach [107] with the largest-𝑘 norm to form an exact
continuous representation. The constraint ‖𝛽‖0 ≤ 𝐾 is replaced by the equivalent
constraint ‖𝛽‖1−|||𝛽|||[𝐾] = 0 [107; 106]. Here |||𝛽|||[𝐾] is the largest-𝑘 norm, which
is the sum of the 𝑘 largest absolute value elements in 𝛽. Thus, the optimization prob-
lem is rewritten as ⎧⎨⎩min𝛽∈R𝑝 −𝑙(𝛽)s.t. ‖𝛽‖1 − |||𝛽|||[𝐾] = 0. (30)
The constraint is continuous but also noncovex due to its combinatorial structure.
The constrained optimization problem in (30) is rewritten as an unconstrained prob-
lem [127; 128]
min
𝛽∈R𝑝
𝑓(𝛽) = −𝑙(𝛽) + 𝜌
(︁
‖𝛽‖1 − |||𝛽|||[𝐾]
)︁
(31)
where 𝜌 is a positive penalization parameter. By suitably large 𝜌 the original con-
straint in (30) is forced to hold.
The objective in (31) is nonconvex and nonsmooth, but it is also a DC function
(Difference of two Convex functions) and can be represented as 𝑓 = 𝑓1− 𝑓2, where
𝑓1(𝛽) = −𝑙(𝛽) + 𝜌‖𝛽‖1 and 𝑓2(𝛽) = 𝜌|||𝛽|||[𝐾].
Here the functions 𝑓1 and 𝑓2 are convex.
The optimization problem (31) is solved with the new OSCAR algorithm presented
in Publication IV. Since the objective is DC, the algorithm utilizes the double bundle
method (DBDC) for DC optimization [114; 113]. The starting points for DBDC are
generated with and incremental approach. Solving the problem starts with allowing
only a single predictor (or kit) to be selected into the model. Then the number of
predictors (or kits) is increased until the maximal number of predictors (or kits) is
reached. The solution for the cardinality-constrained problem with 𝑘 − 1 predictors
is used to derive promising starting points for the problem with 𝑘 predictors. Since
such selection may fall in local optimum, we use several starting points and, from
the resulting solutions, the best is selected.
Similarly to Greedy Cox in Publication III, the OSCAR algorithm requires the
input matrix including the observed patient features and the survival data (time and
event/right-censoring) in case of the survival prediction. If logistic or linear regres-
sion is used, a vector of responses should be given instead. Additionally, the max-
imum number of predictors 𝐾 needs to be defined. If not given by the user, it is
39
Anni S. Halkola
calculated from the feature matrix in the R interface. A kit matrix may also be given
to fit the features as groups.
In the algorithm, the solution is initialized by solving the unconstrained problem
−𝑙 with DBDC. Starting from 𝑘 = 1 the starting points are constructed by varying
the previous solution 𝛽*𝑘−1 (or the initialized solution in case of 𝑘 = 1). Then for
all starting points the problem (31) is solved with DBDC method. The penalization
parameter is increased if the cardinality-constraint is not satisfied. The best solution
is selected from the solutions obtained from different starting points. The steps are
repeated up until the problem with cardinality 𝑘 = 𝐾 is solved. The algorithm re-
turns the solutions (i.e., the vectors 𝛽) for each cardinality-constrained problem with
𝑘 = 1, . . . ,𝐾 predictors.
In addition to DBDC, the method can use other solvers that are capable to handle
nonsmoothness and nonconvexity. In our R package, we included also the limited
memory bundle method (LMBM) [115; 116], which can be used in optimization
instead of DBDC. LMBM is for general nonconvex and nonsmooth minimization
problems and it can handle large-scale problems. However, unlike DBDC, LMBM
does not benefit from the DC structure.
4.3 Results
4.3.1 LASSO selection and Greedy Cox
The methods (LASSO selection and Greedy Cox) were tested in the four data co-
horts TYKS, VENICE, MAINSAIL and ASCENT (TYKS was named RWC, i.e.
real-world cohort). Different training and testing data combinations were tested. If
the training and testing data were from the same cohort, subsets of the cohort were
assigned to either training or testing set. C-index was calculated either for a single
fitting or cross-validation.
The budged was varied, and for low budgets Greedy Cox gave better results. In
cross-validation, the highest model performance was usually reached with LASSO
selection. However, if the budget was further increased, the performance of LASSO
selection deteriorated while the results of Greedy Cox stayed on the same level. This
may be because in LASSO selection the penalization is reduced to obtain sets with
higher budget while Greedy Cox has fixed penalization. However, with Greedy Cox,
a small increase in a low budget might result in lower model performance since with
the increased budget, an expensive variable group was selected first with no money
left for other non-free variables. If a lower budget was given, the first variable group
had to be less expensive and another low-cost, but non-free, group could be selected
as well. Peak model performance was achieved with a budget around 200-400 with
usually 10-15 variables.
In TYKS cohort with an unlimited budget, Greedy Cox selected prostate-specific
40
Cancer survival prediction based on clinical data
antigen (PSA) as the first predictor followed by the group of hemoglobin (HB), white
blood cells (WBC), ratio of lymphocytes and leukocytes (LYMperLEU) and ratio
of neutrophils and leukocytes (NEUperLEU). In the third step alkaline phosphatase
(ALP) was selected in addition to the previous variables. Also LASSO selection
started with PSA, HB, ALP and the age group (AGEGRP) (four steps in total). With
more steps, the order of selected variables started to differ between the two methods.
In Figure 16 the results are presented in relation to the costs and C-indices (until the
budget is 300). In the trial cohorts the selected variables varied, however, lactate
Greedy Cox
LASSOStep 1
Step 1
Step 2
Step 2
Steps 3 and 4
Figure 16. The C-index in TYKS cohort with respect to cost when the model was fitted with Greedy
Cox (red) or LASSO selection (red). The budget was unlimited for Greedy Cox and the cost was
calculated afterwards.
dehydrogenase (LDH) was selected in all trial cohorts with both Greedy Cox and
LASSO selection. However, LDH was not available in TYKS cohort. Most likely
variables to be selected in trial cohorts were LDH, HB, ALP, PSA and whether the
patient had suffered congestive heart failure (CHF). From these, PSA, HB and ALP
were also prominent predictors in TYKS cohort.
In addition, models for selected budgets were further investigated with the ePCR
package [124] after the variables were first selected by LASSO selection or Greedy
Cox. The further analysis with ePCR was conducted to optimize the amount of pe-
nalization (e.g. 𝜆 in Eq. (28)). Since the ePCR package can construct ensemble mod-
els combining results of multiple models, some ensembles were also investigated for
models fitted in different data cohorts. For example, when tested in TYKS cohort, the
best model was fitted in VENICE cohort with variables selected by LASSO selection
(Fig. 17). See Supplementary of Publication III for more details.
41
Anni S. Halkola
integrated AUC
0.50 0.60 0.70 0.80
VEN_LASSO_r_b380
ensemble_LASSO_r_b380
VEN_LASSO_r_b190
ensemble_LASSO_r_b190
VEN_greedy_r_b190
ensemble_greedy_r_b380
VEN_greedy_r_b380
ensemble_LASSO_f_b380
VEN_LASSO_f_b380
ensemble_greedy_f_b190
VEN_greedy_f_b190
Figure 17. Integrated AUC for models tested in TYKS cohort. Models were either fitted to VENICE
cohort or VENICE, MAINSAIL and VENICE + MAINSAIL for ensemble models. Two different bud-
gets are shown (190 and 380). Variables were first selected by LASSO selection or Greedy Cox in
the VENICE cohort. The letter ’f’ indicates that the full set of variables was allowed, and the letter
’r’ indicates that the set was reduced to those with less than 40% of missing values in TYKS cohort.
The colors indicate Bayes factors (BF) calculated to compare the models to the top model (yellow):
BF < 20 (blue), BF ≥ 20 (green). (Adopted with permission from Publication III: Suppl. Fig. 3.)
4.3.2 OSCAR
We tested the OSCAR method in the four data cohorts and analyzed the results ob-
tained from fitting into the training data (C-index), cross-validation and bootstrap
analysis. In addition, we investigated, which models were cost-efficient (i.e., maxi-
mum accuracy, minimal cost), thus making the underlying problem a multi-objective
optimization problem. We obtained an approximation of the Pareto-front from cross-
validation C-index and costs for corresponding cardinalities in the training data fit.
We further tested the models of the approximated Pareto-front in the validation sets
separated from the data before model fitting. We also compared the OSCAR method
to LASSO and APM-𝐿0 [112]. See Figure 18 for the schematic illustration of the
method and the post-fitting analyses.
In the TYKS cohort the PSA was the most robust predictor and in the bootstrap
analysis it was selected most often if only one predictor was allowed. HB, ALP and
age group were also promising predictors and the C-index improved when compared
to the PSA alone (0.654 vs 0.726). The cost also remained low with these predictors
(100 with only PSA vs. 160 with the four predictors). The model accuracy did not
dramatically improve with more predictors (maximum 0.735 in training data). In
42
Cancer survival prediction based on clinical data
1    1    0     0     0
0    0    1     0     0
0    1    0     1     1
0    0    0     0     1
Feature
Kit
Kit structure 
can be included
Multiple starting 
points and the best 
is selected for the 
fixed cardinality
oscar::oscar.bs
oscar::oscar.bs.visu
Solver in optimization: 
DBDC or LMBM
User-friendly R 
package in CRAN
OSCAR method
L0-pseudonorm to 
restrict the number
of predictors 
Bootstrap
Feature 
frequency (%)
PSA
HB
ALP
PULSE
CREAT
HEMAT
NA
HEIGHTBL
RBC
SYSTOLICBP
ALT
WBC
BMI
NEU
DIASTOLICBP
PLT
TESTO
POT
TBILI
WEIGHTBL
GLU
0           50          100
Cardinality (k)
Bootstrap for a better 
view on the most 
important predictive
parameters
Cross-validation to further 
validate the model 
performance and assist in 
the cardinality selection
Cardinality (k)
2 4 6 8 10 
0.60 
0.65 
0.70 
0.75 
0.80 
C
-in
de
x
Bootstrap Cross-validation
k=1
k=2
k=3
k=6
k=13
k=15
0.74
0.72
0.70
0.68
0.66
M
ea
n 
C
-in
de
x 
in
 C
V
0 10 20 30
Cost
Pareto-front
1 4 7 10 13 16 19
0.76
40
Cross-validation Pareto-front
oscar::oscar.cv
oscar::oscar.cv.visu
oscar::oscar.pareto
oscar::oscar.pareto.visu
oscar::oscar
Model family: 
- the Cox's proportional
  hazards model (cox)
- binomial model (logistic)
- linear regression model 
  (mean square error, mse)
Problem solved
for all cardinalities
Plot CV result with costs
to form a Pareto-front:
max C-index with min 
cost for clinical relevance
R interface
Fortran
Control the tuning
parameters
oscar::oscar.control
Result returned and 
saved as a S4-class
object in R
Post fitting
Fitting
Figure 18. Schematic illustration of the OSCAR method and post-fitting analyses. The bottom
panel results are from simulated real-world registry data included in the R package. (Figure adopted
with permission from Publication IV: Fig. 1.)
43
Anni S. Halkola
comparison to LASSO and APM-𝐿0, all three methods selected similar predictors,
e.g. PSA, HB and ALP. However, in the cross-validation, OSCAR resulted in higher
mean C-index than LASSO and APM-𝐿0. This might be due to OSCAR method not
penalizing the values of selected model coefficients, unlike LASSO which pushes
coefficients towards zero.
The cross-validation results were also investigated against the costs and, within
a selected cost, the OSCAR performance was the highest when compared to the
other two methods. The approximated Pareto-fronts (min cost, max accuracy) are
presented in Fig. 19. The models of the approximated Pareto-fronts were tested in
validation data and similar accuracy was reached when compared to the training data
(C-index around 0.70-0.73). If the kit structure was included, PSA was selected first,
followed by B-PVKT kit (basic blood count and trombocytes) consisting here of
HB, platelets (PLT), WBC, RBC and hematocrit (HEMAT). In this application, the
model performance was on a similar level than without the kit structures. However,
including a kit with highly prognostic feature (like HB in the B-PVKT) might have
less prognostic features that are included as well.
OSCAR
LASSO
APML0k=1
k=2
k=4
k=9
k=11 k=12
Figure 19. The cross-validation mean C-indices against the costs for the approximated Pareto-
front models obtained with OSCAR (red), LASSO (yellow) and APM-𝐿0 (red). The cardinalities
are marked for OSCAR. For LASSO the numbers of predictors are 2 and 4 starting from left, and
corresponding numbers for APM-𝐿0 are 4,6 and 14.
In trial cohorts (VENICE, MAINSAIL, ASCENT) the most robust predictors varied.
Unlike in TYKS, PSA was not as robust and, based on bootstrap analysis, it was
selected often only in ASCENT. If only one predictor was allowed, ALP was se-
lected most often in ASCENT and VENICE cohorts. In MAINSAIL, LDH was the
most robust predictor. HB was relevant predictor in VENICE and MAINSAIL. In
the VENICE cohort, potential additional predictors were aspartate aminotransferase
(AST), creatinine (CREAT), sodium (NA) and albumin (ALB). In the MAINSAIL
cohort magnesium (MG), body mass index (BMI), ALB, AST and weight were po-
44
Cancer survival prediction based on clinical data
tential predictor candidates. In the ASCENT cohort neutrophils (NEU), calcium
(CA), LDH and HB were selected often in the bootstrap analysis.
When compared to LASSO and APML-𝐿0, all three suggested similar predic-
tors. However, like with TYKS cohort, in the cross-validation analysis, the OS-
CAR method resulted in higher mean C-index than LASSO and APM-𝐿0. Similarly
to TYKS, the Pareto-front was approximated from cross-validation mean C-index
and costs and OSCAR resulted in higher accuracies when compared to LASSO and
APM-𝐿0 within the same cost. The models of the approximated Pareto-fronts were
tested in the validation data. Considering the objective of maximizing C-index and
minimizing the cost, OSCAR performed well in the validation data.
45
5 Discussion
Previous chapters summarize the approaches that were investigated to obtain predic-
tions of cancer behavior and treatment outcomes. Publications I and II have more
theoretical mechanistic approaches to achieve qualitative predictions of cancer be-
havior. Publications III and IV lean towards machine learning with training data and
patient outcomes.
The ODE model in Publication I (Chapter 2) incorporates the key aspects of cancer
and immune system communication: activation and infiltration of T-cells, elimina-
tion of cancer cells and cancer cells making T-cells inefficient. In addition, the im-
munotherapy was incorporated with targeted and chemotherapy. The model could be
easily altered for different patients to form an individualized model to make patient-
specific predictions and treatment decisions. For example, the cancer cells’ sensi-
tivity to either the T-cell attack (antigen missing) or different treatments could be
tailored. Also, it is possible to further develop the model and include other im-
munotherapy mechanisms, such as the lack and treatment of the initial T-cell activa-
tion (CTLA4 vs. anti-CTLA4).
In Publication I, different immunotherapy schedules were compared with each
other (pre-set treatment periods vs. threshold-based treatment initiation), and it was
noted that both have their advantages and disadvantages. However, the threshold-
based initiation led to a positive treatment outcome (low cancer burden) more easily,
as the treatment schedule responded to the need rather than sticking to a pre-set
schedule that might not have been suitable.
On the other hand, the pre-set periods may also be adapted to an extent, by start-
ing with a tighter schedule and lengthening the drug holidays (gaps between treat-
ments) later, when the initial treatment response is already achieved. In general, the
treatment combinations, with either immuno- and targeted therapy or immuno- and
chemotherapy led to a better treatment effect than monotherapies, however, despite
the additional possibility of fewer treatment initiations, the burden of two therapies
has to be considered. With careful treatment selection and based on the underlying
biological patient characteristics, a stable, possibly undetectable, disease was attain-
able [42; 43; 44; 45] or cancer remained chronic with repeated treatment periods
[129].
46
Discussion
The metapopulation model in Publication II (Chapter 3) adapted some aspects present
in the Publication I (e.g., the resources), however, instead of the immune response,
the process of angiogenesis was incorporated. Also, instead of investigating the treat-
ment effect on non-variating cancer cell subpopulations (Publication I), the metapop-
ulation model was utilized to investigate the change in the populations, i.e., how the
different treatment regimens affected the qualities of the population. We investigated
the treatments’ effect using tools of adaptive dynamics [63; 64; 65] and analyzed
how the cells’ angiogenesis contribution, emigration eagerness and resistance against
anti-angiogenic treatment would evolve under the treatments. A positive angiogene-
sis contribution was shown to be beneficial for cancer cells if the benefits outweigh
the costs [66; 67]. Similarly to real-life observations, the anti-angiogenic monother-
apy did not usually extinguish cancer completely and combination treatments were
required [17; 80; 81]. The main focus was on the anti-angiogenic therapy, however,
we modelled and analyzed also cytotoxic treatment. It was observed that the combi-
nation with anti-angiogenic treatment reduced the need for cytotoxic treatment when
the goals was to reach population non-viability. The treatments and their dosages
have to be carefully selected to avoid high risk of treatment resistance while main-
taining sufficient treatment effect.
Both immune response and angiogenesis play an important part in the cancer de-
velopment, the previous being the body’s own attempt to eradicate the anomaly, and
the latter being the key hallmark to achieve a higher level of cancer growth. Thus it
would be interesting to connect these key phenomena present separately in Publica-
tions I and II into a combined model. Since both of these models also incorporate
resources, it would be interesting to include the cancer cells’ ability to adapt to low
level of resources or hypoxia [130; 131]. The emergence of new (resistant) variants
often causes the relapse and treatment failure [132], and thus it would be interesting
to investigate it also in the ODE model in Publication I. Multiple dynamically adapt-
ing populations could be included, also to present normal or endothelial cells. Also
the treatment mechanisms could be expanded to include a combination of multiple
immune response mechanisms and corresponding therapies (e.g., anti-PD1 and anti-
CTLA4), or a combination of multiple angiogenesis mechanisms and therapies, such
as VEGF family and fibroblast growth factors (FGFs) [33]. In the case the angiogen-
esis and immune response mechanism were combined, the combination therapy of
anti-angiogenic treatment with immunotherapy could be considered [82; 83].
In Publications III and IV (Chapter 4) the cancer survival was predicted with meth-
ods based on penalized Cox regression [9; 93]. It was noted that the economic costs
could be reduced without a great loss in prediction accuracy. In Publication III,
Greedy Cox algorithm performed better than the LASSO selection, when the bud-
get was low, however, overall the LASSO selection gave the highest accuracy. Even
though not completely comparable due to slight differences in the training data, OS-
47
Anni S. Halkola
CAR method in Publication IV produced similar or slightly higher level of accuracy
as the methods in Publication III. For Greedy Cox and OSCAR, the grouped struc-
ture (kit structure) of features was included to accommodate the real-world practice
of measuring the features together (e.g. complete blood cell count).
The methods in Publication III used traditional 𝐿1 and 𝐿2 penalizations, how-
ever, for Publication IV we implemented 𝐿0-pseudonorm, which is discontinuous
and noncovex, making the optimization problem NP-hard [108]. The penalization
term was rewritten to an exact continuous presentation which was also in the form
of a DC composition. This enabled the use of DBDC algorithm [113; 114] in the
optimization. 𝐿0-pseudonorm also allowed the nonzero model coefficients to vary
freely away from zero (unlike, e.g., in LASSO).
We investigated the model performance in four data cohorts, one from hospital
registry data (TYKS) and three from randomized clinical trials (VENICE, MAIN-
SAIL, ASCENT). For TYKS data, all tested methods resulted in similar main pre-
dictors: PSA, HB, ALP and AGEGRP. If the number of predictors (or budget) was
increased, the selected predictors started to vary. For trial cohorts, methods in Pub-
lication III selected predictors such as LDH, HB, ALP, CHF and PSA. These were
distinguished also with OSCAR method: ALP, AST, CREAT, NA, HB and ALB
for VENICE; HB, ALB and LDH for MAINSAIL; PSA, ALP, LDH and HB for
ASCENT. Some of the detected differences between cohorts and methods may be
due that in Publication III, there were more features from which to select. In addi-
tion, some important predictors selected in trial cohorts were not available in TYKS
cohort. However, the selected parameters are reasonable and validate the method
functionality, when compared to previous research (e.g. [95]) and clinical practice.
For example, PSA reflects the disease severity and has been assessed as an important
predictor of disease state [133]. Furthermore, ALP is associated with metastases in
prostate cancer [134]. LDH is released by damaged tissues and has been linked to
cancer [135].
Since the cost was calculated after the model fitting for OSCAR and LASSO se-
lection, in future work it would be an interesting challenge to incorporate the cost as
a part of the objective function. This would require the inclusion of discrete coeffi-
cients in the optimization. The selection of starting points in OSCAR method affects
the outcome and also the running time, thus further development on the starting point
selection would be interesting and potentially beneficial for the overall applicability.
48
6 Conclusions
We modelled cancer dynamics with an ODE model incorporating the mechanism
of immune response and immunotherapy. We predicted the treatment outcome of
different treatment regimens, either as monotherapies or in combinations. Different
scheduling techniques were also compared to see whether the pre-set treatment peri-
ods usually used in practice would be outperformed by a threshold-based treatment
initiation. Generally, a more adaptive schedule, either with a threshold-initiation
or modifying the pre-set treatment schedule, resulted in smaller cancer burden along
with less time in treatment. We also considered cancer in a continuous-time metapop-
ulation model with migration through a dispersal pool. A mechanism of angiogenesis
and anti-angiogenic treatment was included and, using the tools provided by adaptive
dynamics, we investigated how the cancer cells affect the angiogenesis depending on
the treatment. We also analyzed the evolution on emigration and treatment resistance
against anti-angiogenic treatment.
Three variable selection methods were implemented and tested, one being a cost-
specified greedy algorithm (Publication III), one based on LASSO regularization
(Publication III) and one based on 𝐿0-pseudonorm regularization with DBDC op-
timization (Publication IV). All three methods (LASSO selection, Greedy Cox and
OSCAR) were designed to achieve prognostic models with varying levels of sparsity
to monitor, not only the prediction accuracy, but also the cost-efficiency. For Greedy
Cox, the budget was incorporated into the algorithm, whereas the other two methods
consider the cost afterwards. Group-wise feature selection was included in Greedy
Cox and OSCAR since group-wise measurement of features is general practice in
clinic. Greedy Cox and LASSO selection restrict the model parameters with well-
known 𝐿2 and 𝐿1 penalizations. The 𝐿0 penalization in OSCAR method has been
under-represented in regularized regression due to its discontinuous and nonconvex
nature. OSCAR method utilizes an exact DC-composition of the penalty term. OS-
CAR method was also made available as a user-friendly R package in CRAN. We
tested and validated the methods in prostate cancer patient data (real-world and trial
cohorts), which is clinically significant due to high occurrence among men.
We have presented highly applicable models and methods with varying approaches
and key aspects (i.e. angiogenesis and immune response) to gain insights into the
49
Anni S. Halkola
cancer development and response as well as the patient survival prediction. We have
incorporated multiple treatment options in the theoretical models and analyzed the
qualitative outcomes. We utilized a real-world and trial data cohorts from advanced
prostate cancer to predict the overall survival with three penalized regression meth-
ods. Hopefully these novel models and methods lead to more patient-specific treat-
ment choices and bring new tools into the quest of defeating cancer.
50
List of References
[1] Hyuna Sung, Jacques Ferlay, Rebecca L. Siegel, Mathieu Laversanne, Isabelle Soerjomataram,
Ahmedin Jemal, and Freddie Bray. Global Cancer Statistics 2020: GLOBOCAN Estimates of
Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA: A Cancer Journal for
Clinicians, 71(3):209–249, 2021. ISSN 0007-9235. doi: 10.3322/caac.21660.
[2] Ivana Bozic, Johannes G. Reiter, Benjamin Allen, Tibor Antal, Krishnendu Chatterjee, Preya
Shah, Yo Sup Moon, Amin Yaqubie, Nicole Kelly, Dung T. Le, Evan J. Lipson, Paul B. Chap-
man, Luis A. Diaz, Bert Vogelstein, and Martin A. Nowak. Evolutionary dynamics of cancer
in response to targeted combination therapy. eLife, 2013(2):1–16, 2013. ISSN 2050084X. doi:
10.7554/eLife.00747.
[3] Ivana Bozic and Martin A. Nowak. Timing and heterogeneity of mutations associated with drug
resistance in metastatic cancers. Proceedings of the National Academy of Sciences, 111(45):
15964–15968, 2014. ISSN 0027-8424. doi: 10.1073/pnas.1412075111.
[4] Andrej Fischer, Ignacio Va´zquez-Garcı´a, and Ville Mustonen. The value of monitoring to control
evolving populations. Proceedings of the National Academy of Sciences, 112(4):1007–1012,
2015. ISSN 0027-8424. doi: 10.1073/pnas.1409403112.
[5] Yoram Louzoun, Chuan Xue, Gregory B. Lesinski, and Avner Friedman. A mathematical model
for pancreatic cancer growth and treatments. Journal of Theoretical Biology, 351:74–82, 2014.
ISSN 00225193. doi: 10.1016/j.jtbi.2014.02.028.
[6] Franziska Michor and Kathryn Beal. Improving Cancer Treatment via Mathematical Model-
ing: Surmounting the Challenges Is Worth the Effort. Cell, 163(5):1059–1063, 2015. ISSN
10974172. doi: 10.1016/j.cell.2015.11.002.
[7] Jingsong Zhang, Jessica J. Cunningham, Joel S. Brown, and Robert A. Gatenby. Integrating
evolutionary dynamics into treatment of metastatic castrate-resistant prostate cancer. Nature
Communications, 8(1):1–10, 2017. ISSN 20411723. doi: 10.1038/s41467-017-01968-5.
[8] Boyang Zhao, Michael T. Hemann, and Douglas A. Lauffenburger. Modeling tumor clonal
evolution for drug combinations design. Trends in Cancer, 2(3):144–158, mar 2016. ISSN
24058033. doi: 10.1016/j.trecan.2016.02.001.
[9] D. R. Cox. Regression Models and Life-Tables. Journal of the Royal Statistical Society. Series
B, Methodological, 34(2):187–220, 1972. ISSN 0035-9246. doi: 10.1111/j.2517-6161.1972.
tb00899.x.
[10] Artur Ce´sar Fassoni, Ingo Roeder, and Ingmar Glauche. To Cure or Not to Cure : Consequences
of Immunological Interactions in CML Treatment. Bulletin of Mathematical Biology, (88881),
2019. ISSN 1522-9602. doi: 10.1007/s11538-019-00608-x.
[11] Christophe Letellier, Sourav Kumar Sasmal, Cle´ment Draghi, Fabrice Denis, and Dibakar
Ghosh. A chemotherapy combined with an anti-angiogenic drug applied to a cancer model
including angiogenesis. Chaos, Solitons and Fractals, 99:297–311, 2017. ISSN 09600779. doi:
10.1016/j.chaos.2017.04.013.
[12] S. T. R. Pinho, F. S. Bacelar, R. F. S. Andrade, and H. I. Freedman. A mathematical model
for the effect of anti-angiogenic therapy in the treatment of cancer tumours by chemotherapy.
Nonlinear Analysis: Real World Applications, 14(1):815–828, 2013. ISSN 14681218. doi:
10.1016/j.nonrwa.2012.07.034.
51
Anni S. Halkola
[13] Sirin Yonucu, Defne Y𝜄lmaz, Colin Phipps, Mehmet Burcin Unlu, and Mohammad Kohandel.
Quantifying the effects of antiangiogenic and chemotherapy drug combinations on drug delivery
and treatment efficacy. PLoS Computational Biology, 13(9):1–17, 2017. ISSN 15537358. doi:
10.1371/journal.pcbi.1005724.
[14] Stefano Avanzini and Tibor Antal. Cancer recurrence times from a branching process model.
PLoS Computational Biology, 15(11):1–30, 2019. ISSN 15537358. doi: 10.1371/journal.pcbi.
1007423.
[15] Emilia Kozłowska, Anniina Fa¨rkkila¨, Tuulia Vallius, Olli Carpe´n, Jukka Kemppainen, Seija
Gre´nman, Rainer Lehtonen, Johanna Hynninen, Sakari Hietanen, and Sampsa Hautaniemi.
Mathematical modeling predicts response to chemotherapy and drug combinations in ovar-
ian cancer. Cancer Research, 78(14):4036–4044, jul 2018. ISSN 15387445. doi: 10.1158/
0008-5472.CAN-17-3746.
[16] Xiaoran Lai, Oliver M. Geier, Thomas Fleischer, Øystein Garred, Elin Borgen, Simon W. Funke,
Surendra Kumar, Marie E. Rognes, Therese Seierstad, Anne-Lise Børresen-Dale, Vessela N.
Kristensen, Olav Engebra˚ten, Alvaro Ko¨hn-Luque, and Arnoldo Frigessi. Towards person-
alized computer simulation of breast cancer treatment : a multi-scale pharmacokinetic and
pharmacodynamic model informed by multi-type patient data. Cancer research, 2019. doi:
10.1158/0008-5472.CAN-18-1804.
[17] Toma´s Alarco´n, Markus R. Owen, Helen M. Byrne, and Philip K. Maini. Multiscale modelling
of tumour growth and therapy: The influence of vessel normalisation on chemotherapy. Compu-
tational and Mathematical Methods in Medicine, 7(2-3):85–119, 2006. ISSN 1748670X. doi:
10.1080/10273660600968994.
[18] N. L. Komarova, J. A. Burger, and D. Wodarz. Evolution of ibrutinib resistance in chronic lym-
phocytic leukemia (CLL). Proceedings of the National Academy of Sciences, 111(38):13906–
13911, 2014. ISSN 0027-8424. doi: 10.1073/pnas.1409362111.
[19] Daniel S. Chen and Ira Mellman. Oncology Meets Immunology: The Cancer-Immunity Cycle.
Immunity, 39(1):1–10, jul 2013. ISSN 10747613. doi: 10.1016/j.immuni.2013.07.012.
[20] Thomas A. Ferguson, Jayoung Choi, and Douglas R. Green. Armed response: How dying cells
influence T-cell functions. Immunological Reviews, 241(1):77–88, 2011. ISSN 01052896. doi:
10.1111/j.1600-065X.2011.01006.x.
[21] Jingyi Zhou, Gangyang Wang, Yinze Chen, Hongxia Wang, Yingqi Hua, and Zhengdong Cai.
Immunogenic cell death in cancer therapy: Present and emerging inducers. Journal of Cellular
and Molecular Medicine, 23(8):4854–4865, 2019. ISSN 1582-1838. doi: 10.1111/jcmm.14356.
[22] Daniel S. Chen, Bryan A. Irving, and F. Stephen Hodi. Molecular pathways: Next-generation
immunotherapy-inhibiting programmed death-ligand 1 and programmed death-1. Clinical
Cancer Research, 18(24):6580–6587, 2012. ISSN 10780432. doi: 10.1158/1078-0432.
CCR-12-1362.
[23] Drew M. Pardoll. The blockade of immune checkpoints in cancer immunotherapy. Nature
Reviews Cancer, 12(4):252–264, 2012. ISSN 1474175X. doi: 10.1038/nrc3239.
[24] F. Stephen Hodi, Steven J. O’Day, David F. McDermott, Robert W. Weber, Jeffrey A. Sosman,
John B. Haanen, Rene Gonzalez, Caroline Robert, Dirk Schadendorf, Jessica C. Hassel, Wallace
Akerley, Alfons J. M. van den Eertwegh, Jose Lutzky, Paul Lorigan, Julia M. Vaubel, Gerald P.
Linette, David Hogg, Christian H. Ottensmeier, Celeste Lebbe´, Christian Peschel, Ian Quirt,
Joseph I. Clark, Jedd D. Wolchok, Jeffrey S. Weber, Jason Tian, Michael J. Yellin, Geoffrey M.
Nichol, Axel Hoos, and Walter J. Urba. Improved Survival with Ipilimumab in Patients with
Metastatic Melanoma. New England Journal of Medicine, 363(8):711–723, aug 2010. ISSN
0028-4793. doi: 10.1056/nejmoa1003466.
[25] Greg T. Motz and George Coukos. Deciphering and Reversing Tumor Immune Suppression.
Immunity, 39(1):61–73, 2013. ISSN 10747613. doi: 10.1016/j.immuni.2013.07.005.
[26] Chengwen Liu, Weiyi Peng, Chunyu Xu, Yanyan Lou, Minying Zhang, Jennifer A. Wargo,
Jie Qing Chen, Haiyan S. Li, Stephanie S. Watowich, Yan Yang, Dennie Tompers Frederick,
52
LIST OF REFERENCES
Zachary A. Cooper, Rina M. Mbofung, Mayra Whittington, Keith T. Flaherty, Scott E. Wood-
man, Michael A. Davies, Laszlo G. Radvanyi, Willem W. Overwijk, Gregory Lize´e, and Patrick
Hwu. BRAF inhibition increases tumor infiltration by T cells and enhances the antitumor activ-
ity of adoptive immunotherapy in mice. Clinical Cancer Research, 19(2):393–403, 2013. ISSN
10780432. doi: 10.1158/1078-0432.CCR-12-1626.
[27] Carl H. June, Roddy S. O’Connor, Omkar U. Kawalekar, Saba Ghassemi, and Michael C.
Milone. CAR T cell immunotherapy for human cancer. Science, 359(6382):1361–1365, 2018.
ISSN 10959203. doi: 10.1126/science.aar6711.
[28] Young Kwang Chae, Ayush Arya, Wade Iams, Marcelo R. Cruz, Sunandana Chandra, Jaehyuk
Choi, and Francis Giles. Current landscape and future of dual anti-CTLA4 and PD-1/PD-L1
blockade immunotherapy in cancer; lessons learned from clinical trials with melanoma and non-
small cell lung cancer (NSCLC). Journal for ImmunoTherapy of Cancer, 6(1):39, dec 2018.
ISSN 2051-1426. doi: 10.1186/s40425-018-0349-3.
[29] Dev Karan and Peter Van Veldhuizen. Combination immunotherapy with prostate GVAX and
ipilimumab: safety and toxicity. Immunotherapy, 4(6):577–580, jun 2012. ISSN 1750-743X.
doi: 10.2217/imt.12.53.
[30] Ravi A. Madan, Mahsa Mohebtash, Philip M. Arlen, Matteo Vergati, Myrna Rauckhorst,
Seth M. Steinberg, Kwong Y. Tsang, Diane J. Poole, Howard L. Parnes, John J. Wright,
William L. Dahut, Jeffrey Schlom, and James L. Gulley. Ipilimumab and a poxviral vaccine
targeting prostate-specific antigen in metastatic castration-resistant prostate cancer: A phase
1 dose-escalation trial. The Lancet Oncology, 13(5):501–508, 2012. ISSN 14702045. doi:
10.1016/S1470-2045(12)70006-2.
[31] Jedd D. Wolchok, Harriet Kluger, Margaret K. Callahan, Michael A. Postow, Naiyer A. Rizvi,
Alexander M. Lesokhin, Neil H. Segal, Charlotte E. Ariyan, Ruth-Ann Gordon, Kathleen Reed,
Matthew M. Burke, Anne Caldwell, Stephanie A. Kronenberg, Blessing U. Agunwamba, Xi-
aoling Zhang, Israel Lowy, Hector David Inzunza, William Feely, Christine E. Horak, Quan
Hong, Alan J. Korman, Jon M. Wigginton, Ashok Gupta, and Mario Sznol. Nivolumab plus
Ipilimumab in Advanced Melanoma. New England Journal of Medicine, 369(2):122–133, 2013.
ISSN 0028-4793. doi: 10.1056/nejmoa1302369.
[32] Vikas Malhotra and Michael C. Perry. Classical chemotherapy: mechanisms, toxicities and the
therapeutic window. Cancer biology & therapy, 2(4 Suppl 1):4–6, 2003. ISSN 15384047. doi:
10.4161/cbt.199.
[33] Kabir A. Khan and Roy Bicknell. Anti-angiogenic alternatives to VEGF blockade. Clin-
ical and Experimental Metastasis, 33(2):197–210, 2016. ISSN 15737276. doi: 10.1007/
s10585-015-9769-3.
[34] Mehdi Rajabi and Shaker A. Mousa. The role of angiogenesis in cancer treatment. Biomedicines,
5(2), 2017. ISSN 22279059. doi: 10.3390/biomedicines5020034.
[35] Rakesh R. Ramjiawan, Arjan W. Griffioen, and Dan G. Duda. Anti-angiogenesis for cancer
revisited: Is there a role for combinations with immunotherapy? Angiogenesis, 20(2):185–204,
2017. ISSN 15737209. doi: 10.1007/s10456-017-9552-y.
[36] Hal L. Smith and Paul E. Waltman. The theory of the chemostat : dynamics of microbial compe-
tition. Cambridge University Press, 1995. ISBN 0521470277.
[37] Peter To´th and Ja´nos E´rdi. Mathematical Models of Chemical Reactions: Theory and Applica-
tions of Deterministic and Stochastic Models. Manchester University Press, 1989. ISBN ISBN
978-0-7190-2208-1.
[38] I. F. Tannock. The relation between cell proliferation and the vascular system in a transplanted
mouse mammary tumour. British Journal of Cancer, 22(2):258–273, 1968. ISSN 15321827.
doi: 10.1038/bjc.1968.34.
[39] Rudolf Gesztelyi, Judit Zsuga, Adam Kemeny-Beke, Balazs Varga, Bela Juhasz, and Arpad
Tosaki. The Hill equation and the origin of quantitative pharmacology. Archive for History of
Exact Sciences, 66(4):427–438, 2012. ISSN 00039519. doi: 10.1007/s00407-012-0098-5.
53
Anni S. Halkola
[40] Greg T. Motz and George Coukos. Deciphering and Reversing Tumor Immune Suppression.
Immunity, 39(1):61–73, 2013. ISSN 10747613. doi: 10.1016/j.immuni.2013.07.005.
[41] Mary L. Disis. Immune regulation of cancer. Journal of Clinical Oncology, 28(29):4531–4538,
2010. ISSN 0732183X. doi: 10.1200/JCO.2009.27.2146.
[42] Julio A. Aguirre-Ghiso. Models, mechanisms and clinical evidence for cancer dormancy. Nature
Reviews Cancer, 7(11):834–846, 2007. ISSN 1474175X. doi: 10.1038/nrc2256.
[43] Liliana Ossowski and Julio A. Aguirre-Ghiso. Dormancy of metastatic melanoma. Pigment Cell
and Melanoma Research, 23(1):41–56, 2010. ISSN 17551471. doi: 10.1111/j.1755-148X.2009.
00647.x.
[44] Robert D. Schreiber, Lloyd J. Old, and Mark J. Smyth. Cancer immunoediting: Integrating
immunity’s roles in cancer suppression and promotion. Science, 331(6024):1565–1570, 2011.
ISSN 00368075. doi: 10.1126/science.1203486.
[45] Daniela Senft and Ze’ev A. Ronai. Immunogenic, cellular, and angiogenic drivers of tumor
dormancy-a melanoma view. Pigment Cell and Melanoma Research, 29(1):27–42, 2016. ISSN
1755148X. doi: 10.1111/pcmr.12432.
[46] Suzanne L. Topalian, Mario Sznol, David F. McDermott, Harriet M. Kluger, Richard D. Carva-
jal, William H. Sharfman, Julie R. Brahmer, Donald P. Lawrence, Michael B. Atkins, John D.
Powderly, Philip D. Leming, Evan J. Lipson, Igor Puzanov, David C. Smith, Janis M. Taube,
Jon M. Wigginton, Georgia D. Kollia, Ashok Gupta, Drew M. Pardoll, Jeffrey A. Sosman, and
F. Stephen Hodi. Survival, durable tumor remission, and long-term safety in patients with ad-
vanced melanoma receiving nivolumab. Journal of Clinical Oncology, 32(10):1020–1030, 2014.
ISSN 15277755. doi: 10.1200/JCO.2013.53.0105.
[47] Z. J. Wolner, A. A. Marghoob, M. P. Pulitzer, M. A. Postow, and M. A. Marchetti. A case report
of disappearing pigmented skin lesions associated with pembrolizumab treatment for metastatic
melanoma. British Journal of Dermatology, 178(1):265–269, 2018. ISSN 13652133. doi:
10.1111/bjd.15354.
[48] J. Hamish E. Laing, George D. Wilson, and Christine A. Martindale. Proliferation rates in
human malignant melanoma. Melanoma Research, 13(3):271–277, 2003. ISSN 0960-8931. doi:
10.1097/00008390-200306000-00008.
[49] Ingunn Tufto and Einar K. Rofstad. Interstitial Fluid Pressure, Fraction of Necrotic Tumor
Tissue, and Tumor Cell Density in Human Melanoma Xenografts. Acta Oncologica, 37(3):
291–297, jan 1998. ISSN 0284-186X. doi: 10.1080/028418698429603.
[50] Gavin P. Dunn, Lloyd J. Old, and Robert D. Schreiber. The immunobiology of cancer immuno-
surveillance and immunoediting, 2004. ISSN 10747613.
[51] Jake S. O’Donnell, Michele W. L. Teng, and Mark J. Smyth. Cancer immunoediting and resis-
tance to T cell-based immunotherapy. Nature Reviews Clinical Oncology, 16(3):151–167, 2019.
ISSN 17594782. doi: 10.1038/s41571-018-0142-8.
[52] Alexander C. Huang, Robert J. Orlowski, Xiaowei Xu, Rosemarie Mick, Sangeeth M. George,
Patrick K. Yan, Sasikanth Manne, Adam A. Kraya, Bradley Wubbenhorst, Liza Dorfman, Kurt
D’Andrea, Brandon M. Wenz, Shujing Liu, Lakshmi Chilukuri, Andrew Kozlov, Mary Car-
berry, Lydia Giles, Melanie W. Kier, Felix Quagliarello, Suzanne McGettigan, Kristin Kreider,
Lakshmanan Annamalai, Qing Zhao, Robin Mogg, Wei Xu, Wendy M. Blumenschein, Jen-
nifer H. Yearley, Gerald P. Linette, Ravi K. Amaravadi, Lynn M. Schuchter, Ramin S. Herati,
Bertram Bengsch, Katherine L. Nathanson, Michael D. Farwell, Giorgos C. Karakousis, E. John
Wherry, and Tara C. Mitchell. A single dose of neoadjuvant PD-1 blockade predicts clinical
outcomes in resectable melanoma. Nature Medicine, 25(3):454–461, 2019. ISSN 1546170X.
doi: 10.1038/s41591-019-0357-y.
[53] Hirokazu Tokuyasu, Soichiro Ishikawa, Hiromitsu Sakai, Tomoyuki Ikeuchi, and Hiroshi Miura.
Single pembrolizumab treatment causing profound durable response in a patient with pulmonary
pleomorphic carcinoma. Respiratory Medicine Case Reports, 28(June):26–30, 2019. ISSN
22130071. doi: 10.1016/j.rmcr.2019.100879.
54
LIST OF REFERENCES
[54] Shaoming Zhu, Tian Zhang, Lei Zheng, Hongtao Liu, Wenru Song, Delong Liu, Zihai Li,
and Chong xian Pan. Combination strategies to maximize the benefits of cancer immunother-
apy. Journal of Hematology and Oncology, 14(1):1–33, 2021. ISSN 17568722. doi:
10.1186/s13045-021-01164-5.
[55] Mikko Heino and Ilkka Hanski. Evolution of migration rate in a spatially realistic metapopula-
tion model. American Naturalist, 157(5):495–511, 2001. ISSN 00030147. doi: 10.1086/319927.
[56] Tuomas Nurmi and Kalle Parvinen. On the evolution of specialization with a mechanistic under-
pinning in structured metapopulations. Theoretical Population Biology, 73(2):222–243, 2008.
ISSN 00405809. doi: 10.1016/j.tpb.2007.12.002.
[57] Kalle Parvinen. Evolution of dispersal in a structured metapopulation model in discrete time.
Bulletin of Mathematical Biology, 68(3):655–678, 2006. ISSN 00928240. doi: 10.1007/
s11538-005-9040-1.
[58] Kalle Parvinen and Johan A. J. Metz. A novel fitness proxy in structured locally finite metapop-
ulations with diploid genetics, with an application to dispersal evolution. Theoretical Population
Biology, 73(4):517–528, 2008. ISSN 00405809. doi: 10.1016/j.tpb.2008.01.002.
[59] Kalle Parvinen. Adaptive Dynamics of Altruistic Cooperation in a Metapopulation: Evolutionary
Emergence of Cooperators and Defectors or Evolutionary Suicide? Bulletin of Mathematical
Biology, 73(11):2605–2626, 2011. ISSN 00928240. doi: 10.1007/s11538-011-9638-4.
[60] Kalle Parvinen, Hisashi Ohtsuki, and Joe Yuichiro Wakano. Spatial heterogeneity and evolu-
tion of fecundity-affecting traits. Journal of Theoretical Biology, 454:190–204, 2018. ISSN
10958541. doi: 10.1016/j.jtbi.2018.06.005.
[61] Kalle Parvinen, Hisashi Ohtsuki, and Joe Yuichiro Wakano. Evolution of dispersal in a spa-
tially heterogeneous population with finite patch sizes. Proceedings of the National Academy
of Sciences of the United States of America, 117(13):7290–7295, 2020. ISSN 10916490. doi:
10.1073/pnas.1915881117.
[62] Helene C. Weigang and E´va Kisdi. Evolution of dispersal under a fecundity-dispersal trade-off.
Journal of Theoretical Biology, 371:145–153, 2015. ISSN 10958541. doi: 10.1016/j.jtbi.2015.
02.013.
[63] A˚ke Bra¨nnstro¨m, Jacob Johansson, and Niels von Festenberg. The Hitchhiker’s guide to adaptive
dynamics. Games, 4(3):304–328, 2013. ISSN 20734336. doi: 10.3390/g4030304.
[64] Odo Diekmann. A beginner’s guide to adaptive dynamics. Banach Center Publications, 63:
47–86, 2003. ISSN 0137-6934. doi: 10.4064/bc63-0-2.
[65] S. A. H. Geritz, E´ Kisdi, G. Mesze´na, and J. A. J. Metz. Evolutionarily singular strategies and
the adaptive growth and branching of the evolutionary tree. Evolutionary Ecology, 12(1):35–57,
1998. ISSN 02697653. doi: 10.1023/A:1006554906681.
[66] Judah Folkman. How Is Blood Vessel Growth Regulated in Normal and Neoplastic Tissue?- G.
H. A. Clowes Memorial Award Lecture. Cancer Research, 46:467–473, 1986.
[67] Lance Allen Liotta, Jerome Kleinerman, and Gerald M. Saidel. Quantitative Relationships of
Intravascular Tumor Cells, Tumor Vessels, and Pulmonary Metastases following Tumor Implan-
tation. Cancer Research, 34(5):997–1004, 1974. ISSN 15387445.
[68] Naoyo Nishida, Hirohisa Yano, Takashi Nishida, Toshiharu Kamura, and Masamichi Kojiro.
Angiogenesis in cancer. Vascular Health and Risk Management, 2(3):213–219, 2006. ISSN
11766344. doi: 10.2147/vhrm.2006.2.3.213.
[69] Yutaka Takahashi, Yasuhiko Kitadai, Corazon D. Bucana, Lee M. Ellis, Yutaka Takahashi,
Lee M. Ellis, and Karen R. Cleary. Expression of Vascular Endothelial Growth Factor and
Its Receptor, KDR, Correlates with Vascularity, Metastasis, and Proliferation of Human Colon
Cancer. Cancer Research, 55(18):3969–3972, 1995. ISSN 15387445.
[70] Gabriele Bergers and Douglas Hanahan. Modes of resistance to anti-angiogenic therapy. Nature
Reviews Cancer, 8(8):592–603, aug 2008. ISSN 1474-175X. doi: 10.1038/nrc2442.
[71] John D. Nagy and Dieter Armbruster. Evolution of uncontrolled proliferation and the angiogenic
switch in cancer. Mathematical Biosciences and Engineering, 9(4):843–876, oct 2012. ISSN
1551-0018. doi: 10.3934/mbe.2012.9.843.
55
Anni S. Halkola
[72] Rui D. M. Travasso, Eugenia Corvera Poire´, Mario Castro, Juan Carlos Rodrguez-Manzaneque,
and A. Herna´ndez-Machado. Tumor angiogenesis and vascular patterning: A mathematical
model. PLoS ONE, 6(5):1–10, 2011. ISSN 19326203. doi: 10.1371/journal.pone.0019989.
[73] Jiangping Xu, Guillermo Vilanova, and Hector Gomez. A mathematical model coupling tumor
growth and angiogenesis. PLoS ONE, 11(2):1–20, 2016. ISSN 19326203. doi: 10.1371/journal.
pone.0149422.
[74] Luke Olsen, Jonathan A. Sherratt, Philip K. Maini, and Frank Arnold. A mathematical model for
the capillary endothelial cell-extracellular matrix interactions in wound-healing angiogenesis.
IMA Journal of Mathemathics Applied in Medicine and Biology, 14(4):261–281, 1997. ISSN
02650746. doi: 10.1093/imammb/14.4.261.
[75] Richard C. Schugart, Avner Friedman, Rui Zhao, and Chandan K. Sen. Wound angiogenesis as a
function of tissue oxygen tension: A mathematical model. Proceedings of the National Academy
of Sciences of the United States of America, 105(7):2628–2633, 2008. ISSN 00278424. doi:
10.1073/pnas.0711642105.
[76] Diane R. Bielenberg and Bruce R. Zetter. The Contribution of Angiogenesis to the Process
of Metastasis. Cancer Journal (United States), 21(4):267–273, 2015. ISSN 1540336X. doi:
10.1097/PPO.0000000000000138.
[77] Dietmar W. Siemann. The Unique Characteristics of Tumor Vasculature and Preclinical Evidence
for its Selective Disruption by Tumor-Vascular Disrupting Agents. Cancer Treatment Reviews,
37(1):63–74, 2011. doi: 0.1016/j.ctrv.2010.05.001.
[78] Fabienne Baffert, Gavin Thurston, Michael Rochon-Duck, Tom Le, Rolf Brekken, and Don-
ald M. McDonald. Age-Related Changes in Vascular Endothelial Growth Factor Dependency
and Angiopoietin-1-Induced Plasticity of Adult Blood Vessels. Circulation Research, 94(7):
984–992, 2004. ISSN 00097330. doi: 10.1161/01.RES.0000125295.43813.1F.
[79] Michael S. Gee, William N. Procopio, Sosina Makonnen, Michael D. Feldman, Newman M.
Yeilding, and William M. F. Lee. Tumor vessel development and maturation impose limits
on the effectiveness of anti-vascular therapy. American Journal of Pathology, 162(1):183–193,
2003. ISSN 00029440. doi: 10.1016/S0002-9440(10)63809-6.
[80] Sajani S. Lakka and Jasti S. Rao. Antiangiogenic therapy in brain tumors. Expert Review Neu-
rotherapeutics, 8(10):1457–1473, 2008. doi: 10.1586/14737175.8.10.1457.
[81] Domenico Ribatti, Tiziana Annese, Simona Ruggieri, Roberto Tamma, and Enrico Crivellato.
Limitations of Anti-Angiogenic Treatment of Tumors. Translational Oncology, 12(7):981–986,
2019. ISSN 19365233. doi: 10.1016/j.tranon.2019.04.022.
[82] Kabir A. Khan and Robert S. Kerbel. Improving immunotherapy outcomes with anti-angiogenic
treatments and vice versa. Nature Reviews Clinical Oncology, 15(5):310–324, 2018. ISSN
17594782. doi: 10.1038/nrclinonc.2018.9.
[83] Won Suk Lee, Hannah Yang, Hong Jae Chon, and Chan Kim. Combination of anti-angiogenic
therapy and immune checkpoint blockade normalizes vascular-immune crosstalk to potenti-
ate cancer immunity. Experimental and Molecular Medicine, 52(9):1475–1485, 2020. ISSN
20926413. doi: 10.1038/s12276-020-00500-y.
[84] Francesca Elice and Francesco Rodeghiero. PL-09 side effects of anti-angiogenic drugs. Throm-
bosis Research, 129(SUPPL. 1):S50–S53, 2012. ISSN 00493848. doi: 10.1016/S0049-3848(12)
70016-6.
[85] Muhammad Wasif Saif, Aymen Elfiky, and Ronald R. Salem. Gastrointestinal perforation due to
bevacizumab in colorectal cancer. Annals of Surgical Oncology, 14(6):1860–1869, 2007. ISSN
10689265. doi: 10.1245/s10434-006-9337-9.
[86] Liang Qiao and Geofrrey C. Farrell. The effects of cell density, attachment substratum and
dexamethasone on spontaneous apoptosis of rat hepatocytes in primary culture. In Vitro Cellular
and Developmental Biology - Animal, 35(7):417–424, 1999. ISSN 10712690. doi: 10.1007/
s11626-999-0117-2.
[87] Dorothy H. Rowe, Jianzhong Huang, Mark L. Kayton, Richard Thompson, Andrea Troxel,
Kathleen M. O’Toole, Darrell Yamashiro, Charles J.H. Stolar, and Jessica J. Kandel. Anti-
56
LIST OF REFERENCES
VEGF antibody suppresses primary tumor growth and metastasis in an experimental model
of Wilms’ tumor. Journal of Pediatric Surgery, 35(1):30–33, 2000. ISSN 00223468. doi:
10.1016/S0022-3468(00)80008-1.
[88] Robert S. Warren, Hui Yuan, Mary R. Matli, Nancy A. Gillett, and Napoleone Ferrara. Regu-
lation by vascular endothelial growth factor of human colon cancer tumorigenesis in a mouse
model of experimental liver metastasis. Journal of Clinical Investigation, 95(4):1789–1797,
1995. ISSN 00219738. doi: 10.1172/JCI117857.
[89] J. A. J. Metz, R. M. Nisbet, and S. A. H. Geritz. How should we define ‘fitness’ for general
ecological scenarios? Trends in Ecology & Evolution, 7(6):198–202, jun 1992. ISSN 01695347.
doi: 10.1016/0169-5347(92)90073-K.
[90] J. A. J. Metz and M. Gyllenberg. How should we define fitness in structured metapopulation
models? Including an application to the calculation of evolutionarily stable dispersal strate-
gies. Proceedings of the Royal Society B: Biological Sciences, 268(1466):499–508, 2001. ISSN
14712970. doi: 10.1098/rspb.2000.1373.
[91] John G. Kemeny and J. Laurie Snell. Finite Markov chains. Princeton, 1960.
[92] W. D. Hamilton. The genetical evolution of social behaviour. I. Journal of Theoretical Biology,
7:1–16, 7 1964. ISSN 00225193. doi: 10.1016/0022-5193(64)90038-4.
[93] Robert Tibshirani. The lasso method for variable selection in the cox model. Statistics in
Medicine, 16(4):385–395, 1997. ISSN 02776715. doi: 10.1002/(SICI)1097-0258(19970228)16:
4⟨385::AID-SIM380⟩3.0.CO;2-3.
[94] Noah Simon, Jerome Friedman, Trevor Hastie, and Rob Tibshirani. Regularization paths for
Cox’s proportional hazards model via coordinate descent. Journal of Statistical Software, 39(5):
1–13, 2011. ISSN 15487660. doi: 10.18637/jss.v039.i05.
[95] Justin Guinney, Tao Wang, Teemu D. Laajala, Kimberly Kanigel Winner, J. Christopher Bare,
Elias Chaibub Neto, Suleiman A. Khan, Gopal Peddinti, Antti Airola, Tapio Pahikkala, Tuomas
Mirtti, Thomas Yu, Brian M. Bot, Liji Shen, Kald Abdallah, Thea Norman, Stephen Friend,
Gustavo Stolovitzky, Howard Soule, Christopher J. Sweeney, Charles J. Ryan, Howard I. Scher,
Oliver Sartor, Yang Xie, Tero Aittokallio, Fang Liz Zhou, and James C. Costello. Prediction of
overall survival for patients with metastatic castration-resistant prostate cancer: development of
a prognostic model through a crowdsourced challenge with open clinical trial data. The Lancet
Oncology, 18(1):132–142, 2017. ISSN 14745488. doi: 10.1016/S1470-2045(16)30560-5.
[96] Andrew J. Armstrong, Elizabeth S. Garrett-Mayer, Yi Chun Ou Yang, Ronald De Wit, Ian F.
Tannock, and Mario Eisenberger. A contemporary prognostic nomogram for men with hormone-
refractory metastatic prostate cancer: A TAX327 study analysis. Clinical Cancer Research, 13
(21):6396–6403, 2007. ISSN 10780432. doi: 10.1158/1078-0432.CCR-07-1036.
[97] Carine A Bellera, Gae¨tan MacGrogan, Marc Debled, Christine Tunon de Lara, Ve´ronique
Brouste, and Simone Mathoulin-Pe´lissier. Variables with time-varying effects and the cox model:
Some statistical concepts illustrated with a prognostic factor study in breast cancer. BMC Medi-
cal Research Methodology, 10:20, 12 2010. ISSN 1471-2288. doi: 10.1186/1471-2288-10-20.
[98] Jiadong Chu, Na Sun, Wei Hu, Xuanli Chen, Nengjun Yi, and Yueping Shen. The application of
bayesian methods in cancer prognosis and prediction. Cancer Genomics - Proteomics, 19:1–11,
12 2022. ISSN 1109-6535. doi: 10.21873/cgp.20298.
[99] Susan Halabi, Chen Yen Lin, W. Kevin Kelly, Karim Fizazi, Judd W. Moul, Ellen B. Kaplan,
Michael J. Morris, and Eric J. Small. Updated prognostic model for predicting overall survival in
first-line chemotherapy for patients with metastatic castration-resistant prostate cancer. Journal
of Clinical Oncology, 32(7):671–677, 2014. ISSN 15277755. doi: 10.1200/JCO.2013.52.3696.
[100] Richard Meier, Stefan Graw, Joseph Usset, Rama Raghavan, Junqiang Dai, Prabhakar Chalise,
Shellie Ellis, Brooke Fridley, and Devin Koestler. An ensemble-based Cox proportional haz-
ards regression framework for predicting survival in metastatic castration-resistant prostate can-
cer (mCRPC) patients. F1000Research, 5(0):1–16, 2016. ISSN 1759796X. doi: 10.12688/
f1000research.8226.1.
57
Anni S. Halkola
[101] Hideyuki Shimizu and Keiichi I. Nakayama. A 23 gene–based molecular prognostic score pre-
cisely predicts overall survival of breast cancer patients. EBioMedicine, 46:150–159, 8 2019.
ISSN 23523964. doi: 10.1016/j.ebiom.2019.07.046.
[102] Wan Zhu, Longxiang Xie, Jianye Han, and Xiangqian Guo. The application of deep learning
in cancer prognosis prediction. Cancers, 12:603, 3 2020. ISSN 2072-6694. doi: 10.3390/
cancers12030603.
[103] Mee Young Park and Trevor Hastie. L1-regularization path algorithm for generalized linear
models. Journal of the Royal Statistical Society. Series B: Statistical Methodology, 69(4):659–
677, 2007. ISSN 13697412. doi: 10.1111/j.1467-9868.2007.00607.x.
[104] Jelle J. Goeman. L1 penalized estimation in the Cox proportional hazards model. Biometrical
Journal, 52(1):70–84, 2010. ISSN 03233847. doi: 10.1002/bimj.200900028.
[105] Jean-Philippe Chancelier and Michel De Lara. Capra-Convexity, Convex Factorization and Vari-
ational Formulations for the ℓ0 Pseudonorm. Set-Valued and Variational Analysis, 30(2):597–
619, jun 2022. ISSN 1877-0533. doi: 10.1007/s11228-021-00606-z.
[106] Manlio Gaudioso, E. Gorgone, and J. B. Hiriart-Urruty. Feature selection in SVM via poly-
hedral k-norm. Optimization Letters, 14(1):19–36, 2020. ISSN 18624480. doi: 10.1007/
s11590-019-01482-1.
[107] Jun-ya Gotoh, Akiko Takeda, and Katsuya Tono. DC formulations and algorithms for sparse
optimization problems, volume 169. Springer Berlin Heidelberg, 2018. ISBN 1010701711. doi:
10.1007/s10107-017-1181-0.
[108] B. K. Natarajan. Sparse approximate solutions to linear systems. SIAM Journal on Computing,
24(2):227–234, 1995. ISSN 00975397. doi: 10.1137/S0097539792240406.
[109] Shai Shalev-Shwartz, Nathan Srebro, and Tong Zhang. Trading Accuracy for Sparsity in Opti-
mization Problems with Sparsity Constraints. SIAM Journal on Optimization, 20(6):2807–2832,
jan 2010. ISSN 1052-6234. doi: 10.1137/090759574.
[110] Wenchuan Guo, Shujie Ma, and Zhenqiu Liu. l0ara. URL https://cran.r-project.
org/package=l0ara.
[111] Hussein Hazimeh and Rahul Mazumder. Fast best subset selection: Coordinate descent and local
combinatorial optimization algorithms. Operations Research, 68(5):1517–1537, 2020. ISSN
15265463. doi: 10.1287/opre.2019.1919.
[112] Xiang Li, Shanghong Xie, Donglin Zeng, and Yuanjia Wang. Efficient ℓ0-norm feature selection
based on augmented and penalized minimization. Statistics in Medicine, 37(3):473–486, 2018.
ISSN 10970258. doi: 10.1002/sim.7526.
[113] Kaisa Joki and Adil M. Bagirov. Bundle methods for nonsmooth DC optimization. In Adil M.
Bagirov, Manlio Gaudioso, Napsu Karmitsa, Marko M. Ma¨kela¨, and Sona Taheri, editors, Nu-
merical Nonsmooth Optimization: State of the Art Algorithms, pages 263–296. Springer Inter-
national Publishing, Cham, 2020. ISBN 9783030349103. doi: 10.1007/978-3-030-34910-3 8.
[114] Kaisa Joki, Adil M. Bagirov, Napsu Karmitsa, Marko M. Ma¨kela¨, and Sona Taheri. Double Bun-
dle Method for finding Clarke Stationary Points in Nonsmooth DC Programming. SIAM Journal
on Optimization, 28(2):1892–1919, jan 2018. ISSN 1052-6234. doi: 10.1137/16M1115733.
[115] M. Haarala, K. Miettinen, and M. M. Ma¨kela¨. New limited memory bundle method for large-
scale nonsmooth optimization. Optimization Methods and Software, 19(6):673–692, 2004.
[116] Napsu Haarala, Kaisa Miettinen, and Marko M. Ma¨kela¨. Globally convergent limited memory
bundle method for large-scale nonsmooth optimization. Mathematical Programming, 109(1):
181–205, 2007.
[117] Nicolas Mottet, Joaquim Bellmunt, Michel Bolla, Erik Briers, Marcus G. Cumberbatch, Maria
De Santis, Nicola Fossati, Tobias Gross, Ann M. Henry, Steven Joniau, Thomas B. Lam, Mal-
colm D. Mason, Vsevolod B. Matveev, Paul C. Moldovan, Roderick C.N. van den Bergh, Thomas
Van den Broeck, Henk G. van der Poel, Theo H. van der Kwast, Olivier Rouvie`re, Ivo G. Schoots,
Thomas Wiegel, and Philip Cornford. EAU-ESTRO-SIOG Guidelines on Prostate Cancer. Part
1: Screening, Diagnosis, and Local Treatment with Curative Intent. European Urology, 71(4):
618–629, 2017. ISSN 18737560. doi: 10.1016/j.eururo.2016.08.003.
58
LIST OF REFERENCES
[118] Frank E. Harrell Jr, Robert M. Califf, David B. Pryor, Kerry L. Lee, and Robert A. Rosati.
Evaluating the yield of medical tests. JAMA, 247(18):2543–2546, 1982. doi: 10.1001/jama.
1982.03320430047030.
[119] Hung Hung and Chin-Tsang Chiang. Estimation methods for time-dependent AUC models with
survival data. Canadian Journal of Statistics, 38(1):8–26, 2010. ISSN 03195724. doi: 10.1002/
cjs.10046.
[120] Daniel P. Petrylak, Nicholas J. Vogelzang, Nikolay Budnik, Pawel Jan Wiechno, Cora N. Stern-
berg, Kevin Doner, Joaquim Bellmunt, John M. Burke, Maria Ochoa de Olza, Ananya Choud-
hury, Juergen E. Gschwend, Evgeny Kopyltsov, Aude Flechon, Nicolas Van As, Nadine Houede,
Debora Barton, Abderrahim Fandi, Ulf Jungnelius, Shaoyi Li, Ronald de Wit, and Karim
Fizazi. Docetaxel and prednisone with or without lenalidomide in chemotherapy-naive patients
with metastatic castration-resistant prostate cancer (MAINSAIL): A randomised, double-blind,
placebo-controlled phase 3 trial. The Lancet Oncology, 16(4):417–425, 2015. ISSN 14745488.
doi: 10.1016/S1470-2045(15)70025-2.
[121] Ian F. Tannock, Karim Fizazi, Sergey Ivanov, Camilla Thellenberg Karlsson, Aude Fle´chon,
Iwona Skoneczna, Francisco Orlandi, Gwenaelle Gravis, Vsevolod Matveev, Sevil Bavbek,
Thierry Gil, Luciano Viana, Osvaldo Are´n, Oleg Karyakin, Tony Elliott, Alison Birtle, Em-
manuelle Magherini, Laurence Hatteville, Daniel Petrylak, Bertrand Tombal, and Mark Rosen-
thal. Aflibercept versus placebo in combination with docetaxel and prednisone for treatment
of men with metastatic castration-resistant prostate cancer (VENICE): A phase 3, double-
blind randomised trial. The Lancet Oncology, 14(8):760–768, 2013. ISSN 14702045. doi:
10.1016/S1470-2045(13)70184-0.
[122] Howard I. Scher, Xiaoyu Jia, Kim Chi, Ronald De Wit, William R. Berry, Peter Albers, Brian
Henick, David Waterhouse, Dean J. Ruether, Peter J. Rosen, Anthony A. Meluch, Luke T.
Nordquist, Peter M. Venner, Axel Heidenreich, Luis Chu, and Glenn Heller. Randomized, open-
label phase III trial of docetaxel plus high-dose calcitriol versus docetaxel plus prednisone for
patients with castration-resistant prostate cancer. Journal of Clinical Oncology, 29(16):2191–
2198, 2011. ISSN 0732183X. doi: 10.1200/JCO.2010.32.8815.
[123] Project Data Sphere. URL https://www.projectdatasphere.org/.
[124] Teemu D. Laajala, Mika Murtoja¨rvi, Arho Virkki, and Tero Aittokallio. EPCR: An R-
package for survival and time-to-event prediction in advanced prostate cancer, applied to real-
world patient cohorts. Bioinformatics, 34(22):3957–3959, 2018. ISSN 14602059. doi:
10.1093/bioinformatics/bty477.
[125] N. E. Breslow. Contribution to the discussion of the paper by D.R. Cox. Journal of the Royal
Statistical Society B, 34:216–217, 1972.
[126] Stephen P. Boyd and Lieven Vandenberghe. Convex optimization. Cambridge University Press,
Cambridge ;, 2004. ISBN 0-521-83378-7.
[127] Jorge. Nocedal and Stephen J. Wright. Numerical optimization. Springer series in opera-
tions research and financial engineering. Springer, New York, 2nd ed. edition, 2006. ISBN
9780387400655. doi: 10.1201/b19115-11.
[128] Willard I. Zangwill. Non-Linear Programming Via Penalty Functions. Management Science, 13
(5):344–358, 1967.
[129] Evan J. Lipson, William H. Sharfman, Charles G. Drake, Ira Wollner, Janis M. Taube, Robert A.
Anders, Haiying Xu, Sheng Yao, Alice Pons, Lieping Chen, Drew M. Pardoll, Julie R. Brahmer,
and Suzanne L. Topalian. Durable cancer regression off-treatment and effective reinduction
therapy with an anti-PD-1 antibody. Clinical Cancer Research, 19(2):462–468, 2013. ISSN
10780432. doi: 10.1158/1078-0432.CCR-12-2625.
[130] K. L. Eales, K. E. R. Hollinshead, and D. A. Tennant. Hypoxia and metabolic adaptation of
cancer cells. Oncogenesis, 5(1):e190–e190, 2016. doi: 10.1038/oncsis.2015.50.
[131] Gabriele Grasmann, Ayusi Mondal, and Katharina Leithner. Flexibility and adaptation of can-
cer cells in a heterogenous metabolic microenvironment. International Journal of Molecular
Sciences, 22(3):1–19, 2021. ISSN 14220067. doi: 10.3390/ijms22031476.
59
[132] M. Gerlinger and C. Swanton. How Darwinian models inform therapeutic failure initiated by
clonal heterogeneity in cancer medicine. British Journal of Cancer, 103(8):1139–1143, 2010.
ISSN 15321827. doi: 10.1038/sj.bjc.6605912.
[133] Shahneen Sandhu, Caroline M. Moore, Edmund Chiong, Himisha Beltran, Robert G. Bristow,
and Scott G. Williams. Prostate cancer. The Lancet, 398(10305):1075–1090, 2021. ISSN
1474547X. doi: 10.1016/S0140-6736(21)00950-8.
[134] Daniel Heinrich, Oyvind Bruland, Theresa A. Guise, Hiroyoshi Suzuki, and Oliver Sar-
tor. Alkaline phosphatase in metastatic castration-resistant prostate cancer: Reassessment of
an older biomarker. Future Oncology, 14(24):2543–2556, 2018. ISSN 17448301. doi:
10.2217/fon-2018-0087.
[135] William R. Berry, John Laszlo, Edwin Cox, Ann Walker, and David Paulson. Prognos-
tic factors in metastatic and hormonally unresponsive carcinoma of the prostate. Cancer,
44(2):763–775, aug 1979. ISSN 0008-543X. doi: 10.1002/1097-0142(197908)44:2⟨763::
AID-CNCR2820440251⟩3.0.CO;2-5.
Anni S. H
alkola
A
I 674
A
N
N
A
LES U
N
IV
ERSITATIS TU
RK
U
EN
SIS
ISBN 978-951-29-9012-2 (PRINT)
ISBN 978-951-29-9013-9 (PDF)
ISSN 0082-7002 (PRINT)
ISSN 2343-3175 (ONLINE)
Pa
in
os
al
am
a 
O
y, 
Tu
rk
u,
 F
in
la
nd
 2
02
2
TURUN YLIOPISTON JULKAISUJA – ANNALES UNIVERSITATIS TURKUENSIS
SARJA – SER. AI OSA – TOM. 674 | ASTRONOMICA – CHEMICA – PHYSICA – MATHEMATICA | TURKU 2022
MATHEMATICAL 
MODELLING AND SURVIVAL 
PREDICTION IN CANCER
Anni S. Halkola