# Identification of Crucial Parameters in a Mathematical Multiscale Model of Glioblastoma Growth

Authors: Tina A. Schuetz; Andreas Mang; Stefan Becker; Alina Toma; Thorsten M. Buzug
Journal: Computational and Mathematical Methods in Medicine
 (2014)

---

## Abstract

Glioblastomas are highly malignant brain tumours. Mathematical models and their analysis
provide a tool to support the understanding of the development of these tumours as well as
the design of more effective treatment strategies. We have previously developed a multiscale
model of glioblastoma progression that covers processes on the cellular and molecular scale.
Here, we present a novel nutrient-dependent multiscale sensitivity analysis of this model
that helps to identify those reaction parameters of the molecular interaction network that influence
the tumour progression on the cellular scale the most. In particular, those parameters
are identified that essentially determine tumour expansion and could be therefore used as potential
therapy targets. As indicators for the success of a potential therapy target, a deceleration
of the tumour expansion and a reduction of the tumour volume are employed.
From the results, it can be concluded that no single parameter variation results in a less
aggressive tumour. However, it can be shown that a few combined perturbations of two
systematically selected parameters cause a slow-down of the tumour expansion velocity accompanied
with a decrease of the tumour volume. Those parameters are primarily linked to
the reactions that involve the microRNA-451 and the thereof regulated protein MO25.


---

## Body

## 1. Introduction

The glioblastoma (GB) still can be considered to be the most aggressive primary brain tumour [1] in humans. The current standard therapy is composed of a resection (if possible) and a combined radio- and chemotherapy [2]. However, the current median survival is only about 12 months [3]. Only 4.7% of all patients are still alive five years after diagnosis [4]. Therefore, current research focuses on gaining a better understanding of the disease and developing new treatment options. Among others, several relevant mutations and signalling pathways for GB emergence and progression have been identified and are one of the primary fields of research [5]. Consequently, the efficacy of molecular targeted drugs for GB treatment has been the focus of recent investigations. Yet, in practice, none of the currently available molecular targeted drugs resulted in a significant improvement of the survival expectancy [2].In recent years, the significance of microRNAs and the thereof dependent signalling pathways for cancer has been discovered [6–8]. MicroRNAs are noncoding RNAs that influence the expression of proteins on the posttranscriptional level [6]. In [9], the role of microRNA-451 (miR-451) in glioblastoma is discussed. Godlewski et al. discovered that, on the one hand, in glioblastoma, the extracellular glucose concentration influences the level of miR-451. On the other hand, miR-451 regulates the signalling of the (well-known) LKB1-AMPK-mTOR pathway (LKB1: liver kinase B1, AMPK: AMP activated kinase, and mTOR: mammalian target of rapamycin) [10, 11]. This intracellular pathway is considerably involved in a cell’s decision for either a migrating or proliferating phenotype [9]. Thus, Godlewski et al. introduced a switch for the dichotomy between a migrating and proliferating phenotype in glioblastoma that is also known as theGo or Grow principle [12].Mathematical modelling has evolved as a useful tool for biology and medicine in gaining a better understanding of diseases and biological processes [13, 14]. Consequently, one application of mathematical models is the investigation of existing therapies and the development of new treatment options. This holds in particular true for models of cancer and cancer treatment [15–17]. In general, these models can be distinguished by the scale (macroscopic, microscopic or molecular) that they describe [14]. While on the macroscopic scale mainly general growth processes [18], deformation effects [19–21], or the response to radiotherapy [22–24] are examined, models on the microscopic scale deal among others with the interaction of tumour cells with their microenvironment (e.g., the immune system, the extracellular matrix, or nutrients) [25–28]. Molecular scale models on the contrary focus on effects of mutations and on molecular interaction and signalling pathways [29, 30]. By definition, all models constitute a simplification of reality. However, to obtain a more accurate model of reality, the coupling of models of different scales is unavoidable. In the case of cancer modelling, a few multiscale approaches exist [31]. In [32, 33], a model for glioblastoma growth was introduced that combines an agent based model with an EGFR (epidermal growth factor receptor) signalling network and focuses on the determination of the cell phenotypes “migrating” and “proliferating.” Later, also a model for the progression of lung cancer was developed [34] in that essentially only the molecular interaction network was interchanged. In [35], a model of cell-cell adhesion is presented that links intracellular E-cadherin signalling to a lattice-free model on the cellular scale.All mathematical models have in common that they contain parameters (e.g., diffusion constants, radiosensitivity parameters, absorption, secretion, mutation, or kinetic reaction rates). These parameters might vary from patient to patient or even from cell to cell and need to be determined/estimated from experimental or patient data in the model development process [36]. One important problem is to examine the dependence of the modelled system on the model parameters, that is, to investigate how the system behaviour changes if one or more parameters are varied. This is realized by means of a sensitivity analysis [37, 38]. The goal of such an analysis is to identify those parameters that considerably influence the model, that is parameters are identified to that the system reacts sensitively. On the one hand it is important to estimate these parameters accurately, on the other hand these parameters indicate potential targets for new therapies [39]. So far, sensitivity analyses have been conducted nearly exclusively on a single scale; that is, the parameters that are varied and the system output that is measured belong to the same modelling scale [40]. This might be sufficient for a single scale model; however, for models covering more than one scale, such an analysis would cover only parts of the model and is therefore nonsatisfying. In [39, 40], an approach for a multiscale sensitivity analysis of a multiscale growth model for lung cancer is introduced that couples the molecular and cellular scale. For this analysis, parameters of the molecular scale part of the model are modified (in particular, they variegated the initial conditions for a system of ordinary differential equations (ODEs)). The dependence of the model on these parameters is evaluated at the cellular scale. Thus, the sensitivity analysis lives up to the multiscale modelling concept.In [41], we presented a multiscale model of glioblastoma growth. On the cellular scale, an agent based model (ABM) represents individual cell actions and a partial differential equation (PDE) models the diffusion of a nutrient (glucose). Each cell is further equipped with a molecular interaction network that implements the regulation of cell migration and cell proliferation via the above described glucose-miR-451 signalling pathway. This interaction network is represented by a system of nine ODEs that describe the dependence of the nine relevant molecular species. In total, these nine equations contain 31 reaction parameters that were either taken from the literature or estimated to fit experimental data. We validated our multiscale model by comparing ourin silico simulation results toin vitro experiments.In the work at hand we analyse how the cellular scale of our model depends on the reaction parameters of the molecular interaction network. Motivated by the idea presented in [39, 40], parameters on the molecular scale are modified and the effect on the cellular scale is measured. In particular, we vary each of the 31 reaction parameters by multiplication with factors in the range [0.01,100] (i.e., also very extreme changes are considered) and consider different nutrient conditions. As a model output, we measure the total number of tumour cells and the number of migrating and of proliferating cells as a degree for the tumour volume and tumour build-up, respectively. Furthermore, the number of time steps necessary for a simulation to finish is evaluated as a measure for the expansion velocity. This allows for investigating which parameter changes have a significant impact on the model and, in particular, identifying the parameters and their respective changes that have a “positive” effect from the therapy development view.In the following, we will first briefly recapitulate the major aspects of our multiscale model, before we introduce the method for the multiscale sensitivity analysis. Finally, we will present the results of the analysis and draw some conclusions for future research.


## 2. A Cross-Scale Model of Tumour Growth

To model the progression of GB in a region of a few square millimetres (to be exact 3 mm × 3 mm), a cross-scale model bridging the molecular and cellular scale has been previously developed [41, 42]. In this model, an intracellular molecular interaction network (represented by a system of ODEs) is coupled with an agent based model. Each agent represents a single tumour cell as is described below.


### 2.1. Molecular Interaction Model

The molecular interaction model describes the influence of the extracellular glucose level on the proliferating versus migrating phenotype emergence of GB cells. The glucose level controls the concentration of miR-451 that in turn regulates the concentration of the active LKB1-MO25-STRAD complex [9]. This complex catalyses the phosphorylation of AMPK that influences cell proliferation via the mTORC1 pathway [10] and cell polarity through further signalling cascades [43]. The interaction model is depicted in Figure 1.Figure 1
The molecular interaction network. Different arrow types pointing from one molecular species (given by the variable notationxi and the biological term) to another indicate different kinds of reactions (see legend in figure for more details). The labels next to the reaction arrows denote the reaction parameters that are involved in the respective reactions. Adapted from [41]. Adapted with permission.The molecular interactions are modelled by mass-balance reactions and Michaelis-Menten equations for the enzyme kinetics. In total, the concentration development of the nine molecular species is mathematically represented by a nondimensionalised system of nine nonlinear ODEs. These ODEs involve 31 reaction parameters which are given together with the details on the system of ODEs in [41]. The ODEs and the reaction parameters are also summarized in the appendix.


### 2.2.  Cellular and Microenvironment Model

To model cellular processes (in terms of migration and proliferation) and the influence of environmental factors (in terms of the available nutrient concentration) on a tumour cell, an agent based modelling (ABM) approach is employed. In this model, each agent represents a single tumour cell.Comparable toin vitro experiments in a petri dish, in this computational model, tumour growth is investigated in a two-dimensional region. The region of interest (3 mm × 3 mm) is represented by a regular grid of 200 × 200 grid points. Each of these grid points either can hold a tumour cell (the grid spacing of 1.5 × 101 μm × 1.5 × 101 μm corresponds to the average size of a GB cell [44]) or is empty.Furthermore, each grid point holds information on the local nutrient level, in particular on the concentration of glucose. Initially, the glucose concentration (x10(0)) is assumed to be either constant or randomly distributed across the whole 3 mm × 3 mm region. Throughout the simulation, glucose is consumed by the tumour cells, cannot pass through the boundaries of the grid, and diffuses across the grid. The last is modelled by the use of a PDE of the form
(1)∂x10∂t=αD∇2x10,
where D is the diffusion coefficient for glucose and α a factor that realizes slower diffusion in case a grid point is occupied by a tumour cell.Tumour cells can migrate on the grid by moving to an empty neighbouring site or proliferate by placing a daughter cell on an empty neighbouring site. Migration and proliferation are governed by chemotaxis; that is, the movement is directed along a chemotactic gradient (here, glucose acts as the chemotactic agent).


### 2.3.  Bridging the Modelling Scales

To incorporate information from the gene and protein scale into the cellular level, each cell/agent on the grid (as described in Section2.2) is provided with a molecular interaction network (as given in Section 2.1). In each time step, each cell is provided with information on the local glucose level. Based on this information, the system of ODEs is evaluated independently for each individual cell. Then, the new phenotype of each cell (migrating, proliferating, or quiescent) is determined on basis of the concentration of AMPK (x6) and mTORC1 (x9). According to this phenotype, the cells are moved to a new location, a daughter cell is placed at a neighbouring site (after a certain delay to incorporate different time scales of migration and proliferation), or no interaction takes place.The general model setup is depicted in Figure2. Representative simulation results are shown in Figure 3. In these simulations, a constant initial glucose concentration x10(0)∈{3 × 10−1 gL−1, 1.125 gL−1, 2.25 gL−1, 4.5 gL−1} is assumed (a setting comparable toin vitro experiments). Initially, 797 cells are placed in a circular shape in the center of the grid. All simulations were run until the first cell reached the boundary. As could be seen in [41], these simulation results are in good agreement with results fromin vitro experiments.Figure 2
General model setup. The scheme shows how the different scales are combined. From the cellular grid on which the model is described by an ABM and a PDE, information about the local glucose concentration is transferred to each individual cell. Each cell is equipped with a system of ODEs describing the molecular interaction network. After simulating this system of ODEs, the cell’s phenotype (migrating, proliferating, or quiescent) is determined. This phenotype then decides on the action of each individual cell on the cellular scale (i.e., the phenotype is realized on the grid).Figure 3
The spatiotemporal tumour development for four different initial glucose concentrations (from top to bottom: 3 × 10−1 gL−1, 1.125 gL−1, 2.25 gL−1, and 4.5 gL−1). Black denotes quiescent, dark gray migrating, and light gray proliferating tumour cells. Adapted from [41]. Adapted with permission.


## 2.1. Molecular Interaction Model

The molecular interaction model describes the influence of the extracellular glucose level on the proliferating versus migrating phenotype emergence of GB cells. The glucose level controls the concentration of miR-451 that in turn regulates the concentration of the active LKB1-MO25-STRAD complex [9]. This complex catalyses the phosphorylation of AMPK that influences cell proliferation via the mTORC1 pathway [10] and cell polarity through further signalling cascades [43]. The interaction model is depicted in Figure 1.Figure 1
The molecular interaction network. Different arrow types pointing from one molecular species (given by the variable notationxi and the biological term) to another indicate different kinds of reactions (see legend in figure for more details). The labels next to the reaction arrows denote the reaction parameters that are involved in the respective reactions. Adapted from [41]. Adapted with permission.The molecular interactions are modelled by mass-balance reactions and Michaelis-Menten equations for the enzyme kinetics. In total, the concentration development of the nine molecular species is mathematically represented by a nondimensionalised system of nine nonlinear ODEs. These ODEs involve 31 reaction parameters which are given together with the details on the system of ODEs in [41]. The ODEs and the reaction parameters are also summarized in the appendix.


## 2.2.  Cellular and Microenvironment Model

To model cellular processes (in terms of migration and proliferation) and the influence of environmental factors (in terms of the available nutrient concentration) on a tumour cell, an agent based modelling (ABM) approach is employed. In this model, each agent represents a single tumour cell.Comparable toin vitro experiments in a petri dish, in this computational model, tumour growth is investigated in a two-dimensional region. The region of interest (3 mm × 3 mm) is represented by a regular grid of 200 × 200 grid points. Each of these grid points either can hold a tumour cell (the grid spacing of 1.5 × 101 μm × 1.5 × 101 μm corresponds to the average size of a GB cell [44]) or is empty.Furthermore, each grid point holds information on the local nutrient level, in particular on the concentration of glucose. Initially, the glucose concentration (x10(0)) is assumed to be either constant or randomly distributed across the whole 3 mm × 3 mm region. Throughout the simulation, glucose is consumed by the tumour cells, cannot pass through the boundaries of the grid, and diffuses across the grid. The last is modelled by the use of a PDE of the form
(1)∂x10∂t=αD∇2x10,
where D is the diffusion coefficient for glucose and α a factor that realizes slower diffusion in case a grid point is occupied by a tumour cell.Tumour cells can migrate on the grid by moving to an empty neighbouring site or proliferate by placing a daughter cell on an empty neighbouring site. Migration and proliferation are governed by chemotaxis; that is, the movement is directed along a chemotactic gradient (here, glucose acts as the chemotactic agent).


## 2.3.  Bridging the Modelling Scales

To incorporate information from the gene and protein scale into the cellular level, each cell/agent on the grid (as described in Section2.2) is provided with a molecular interaction network (as given in Section 2.1). In each time step, each cell is provided with information on the local glucose level. Based on this information, the system of ODEs is evaluated independently for each individual cell. Then, the new phenotype of each cell (migrating, proliferating, or quiescent) is determined on basis of the concentration of AMPK (x6) and mTORC1 (x9). According to this phenotype, the cells are moved to a new location, a daughter cell is placed at a neighbouring site (after a certain delay to incorporate different time scales of migration and proliferation), or no interaction takes place.The general model setup is depicted in Figure2. Representative simulation results are shown in Figure 3. In these simulations, a constant initial glucose concentration x10(0)∈{3 × 10−1 gL−1, 1.125 gL−1, 2.25 gL−1, 4.5 gL−1} is assumed (a setting comparable toin vitro experiments). Initially, 797 cells are placed in a circular shape in the center of the grid. All simulations were run until the first cell reached the boundary. As could be seen in [41], these simulation results are in good agreement with results fromin vitro experiments.Figure 2
General model setup. The scheme shows how the different scales are combined. From the cellular grid on which the model is described by an ABM and a PDE, information about the local glucose concentration is transferred to each individual cell. Each cell is equipped with a system of ODEs describing the molecular interaction network. After simulating this system of ODEs, the cell’s phenotype (migrating, proliferating, or quiescent) is determined. This phenotype then decides on the action of each individual cell on the cellular scale (i.e., the phenotype is realized on the grid).Figure 3
The spatiotemporal tumour development for four different initial glucose concentrations (from top to bottom: 3 × 10−1 gL−1, 1.125 gL−1, 2.25 gL−1, and 4.5 gL−1). Black denotes quiescent, dark gray migrating, and light gray proliferating tumour cells. Adapted from [41]. Adapted with permission.


## 3. Sensitivity Analysis

In total, 31 parameters determine the behaviour of the tumour growth model by regulating nine ODEs. One important question one should consider is whether there exist one or more parameters that have a crucial influence on the whole system. In particular, when considering novel targets for therapeutic interventions, one searches for parameters that slow down the tumour expansion and decrease the tumour volume.Sensitivity analyses are a tool to investigate the dependence of a system (e.g., a specific model outcome) on individual model parameters. The general procedure is to vary one or more parameters and measure the impact on the system. This allows identifying parameters that significantly influence the whole system. Generally, local sensitivity analyses (LSAs) and global sensitivity analyses (GSAs) can be distinguished [37]. In an LSA, all parameters are varied independently and, for each single variation, the impact on the system is measured. This allows for a rather simple and fast realization of the whole analysis. In contrast, a GSA implies that several parameter variations are combined (up to the exhaustion of all possible combinations) and the impact of these combined simultaneous perturbations is measured. This involves rather complex calculations and is computationally expensive. For the purpose of this study, it was sufficient to conduct an LSA and to further analyse the impact of the combined variation of two parameters at a time.The general idea of this sensitivity analysis was to investigate the effects of changes in the intracellular setting of the tumour cells on the cell behaviour and overall tumour expansion. If, for example, a certain parameter significantly slowed down the tumour expansion, the respective reaction could be an indicator for future therapy targets. To accommodate for such an analysis across different spatial scales, we followed the idea introduced in [40]. The parameters were varied on the subcellular scale; that is, the reaction parameters of the system of ODEs that represents the molecular interaction network were altered. However, the system output was measured on the cellular level; that is, values that embody the tumour expansion were recorded. These values were then compared to the original standard setting to allow for an assessment of the influence of the parameter.Initially, the simulation is performed with the original parameter setting and the model outputM is recorded. Next, we multiplied each reaction parameter kl (kl∈{k1,…,k18m}) with different factors bn∈ℝ. For each of these combinations, klbn:=kl·bn of a parameter kl and a factor bn separate simulations are performed and the model output Mkl,bn is recorded. For each parameter-factor setting, the simulation is run three times and the mean of the model output is calculated to account for the random elements of the model. Finally based on the model outputs and the parameter changes sensitivity coefficient Skl,bn are calculated as follows (see [37]):
(2)Skl,bn:=ΔMkl,bn/MΔklbn/kl=(Mkl,bn-M)/M(kl·bn-kl)/kl=Mkl,bn/M-1bn-1
with ΔMkl,bn:=Mkl,bn-M and Δklbn=kl·bn-kl. These sensitivity coefficients facilitate the evaluation of the model outcome relative to the responsible parameter change.Following this definition, the undermentioned observations can be made regarding the properties of the sensitivity coefficientsSkl,bn.(1)
Ifbn<1 (i.e., the parameter kl is decreased) and Mkl,bn<M (i.e., the model output is decreased as well), then the respective sensitivity coefficient is positive (i.e., Skl,bn>0). The same holds true if bn>1 and Mkl,bn>M. If, on the other hand, bn<1 and Mkl,bn>M or bn>1 and Mkl,bn<M, this results in a negative sensitivity coefficient (i.e., Skl,bn<0).(2)
If the same parameterkl is separately multiplied with two different factors b1 and b2 for which b1-1=1-b2 and if additionally the resulting outcomes Mkl,b1 and Mkl,b2 do not differ (i.e., Mkl,b1=Mkl,b2), then the resulting sensitivity coefficients are the same except for the sign (i.e., Skl,b1=-Skl,b2).(3)
If the same parameterkl is separately multiplied with two different factors b1 and b2 that result in the same simulation outcome Mkl,b1=Mkl,b2, then the smaller factor results in the larger respective sensitivity coefficient (i.e., b1<b2⇒Skl,b1>Skl,b2 and vice versa).For the model outputMkl,bne, we choose four different endpoints e that are representative for the respective tumour expansions and tumour characteristics. As a measure for the tumour expansion velocity, we recorded the number of simulation time steps Mkl,bntime necessary for the first tumour cell to reach the boundary of the simulation region. A high expansion velocity corresponds to a short time until the first cell reaches the boundary (i.e., a small Mkl,bntime) and vice versa. The final tumour volume is estimated by the number of cells in the last time step Mkl,bntotal. A large tumour volume is associated with a large number of tumour cells and a small volume with a small number of tumour cells. To further get a better understanding of the tumour build-up, additionally, the number of migrating cells in the last time step Mkl,bnmig and the number of proliferating cells in the last time step Mkl,bnprolif were documented. For all these four model outputs Mkl,bne, the corresponding sensitivity coefficients were calculated and were named accordingly Skl,bntime,Skl,bntotal,Skl,bnmig, and Skl,bnprolif.This setup allows for the assessment of the effects of changes in the reaction parameters on the overall behaviour of the tumour (tumour expansion and tumour build-up) as well in absolute as in relative terms.


## 4. Results

The model is implemented in C++. For the sensitivity analysis, each of the 31 reaction parameters of the molecular interaction model was separately multiplied with the factors 0.01, 0.1, 0.5, 0.8, 0.9, 0.95, 0.98, 0.99, 1.01, 1.02, 1.05, 1.1, 1.2, 1.5, 1.9, 1.99, 5, 10, 50, and 100.As could be seen in Figure3, the tumour growth behaves very different in different glucose environments. Therefore, the sensitivity analysis was also carried out for tumours growing under different initial glucose conditions x10(0). In particular, for each parameter-factor combination, the tumour growth simulation was run three times for each of the four initial glucose settings (3 × 10−1 gL−1, 1.125 gL−1, 2.25 gL−1, and 4.5 gL−1). Thus, in total, simulations for 2480 different settings were carried out.Each simulation was run until the first cell reached the boundary, the necessary number of time stepsMkl,bntime(x10(0)), the final number of all tumour cells Mkl,bntotal(x10(0)), and the number of migrating Mkl,bnmig(x10(0)) and proliferating Mkl,bnprolif(x10(0)) cells in the last time step were recorded, and the respective sensitivity coefficients Skl,bne(x10(0)) were calculated for all four endpoints e. For the ease of notation, the reference to the parameter kl, the factor bn, and the initial glucose concentration x10(0) will be omitted in the following, if the sensitivity coefficients are discussed in general; that is, they will be referred to as Se.


### 4.1. Sensitivity Coefficients

Here, only the results of the sensitivity coefficient calculations for one parameter (k1) are shown in detail (see Figure 4). Results for the other parameters are summarized in Table 1 and in the electronic supplementary material (see Figures S2–S13 in Supplementary Material available online at http://dx.doi.org/10.1155/2014/437094).Table 1
The maximal absolute sensitivity coefficients. This table summarizes for each parameterkl the maximum of the absolute values (𝔪kle) of the four different sensitivity coefficients (Stime,Stotal,Smig, and Sprolif) among all factors bn and all tested initial glucose values (x10(0)). The max⁡(∣·∣) values are dimensionless (d.u.) and the initial glucose concentration values x10(0) are given in g L−1.
k
l
S
time
S
total
S
mig
S
prolif
𝔪
k
l
time
x
10
(
0
)
𝔪
k
l
total
x
10
(
0
)
𝔪
k
l
mig
x
10
(
0
)
𝔪
k
l
prolif
x
10
(
0
)
d.u.
[g L−1]
d.u.
[g L−1]
d.u.
[g L−1]
d.u.
[g L−1]
k
1
5.98
0.3
11.91
0.3
12.28
0.3
11.13
0.3
k
1
i
6.60
0.3
14.88
0.3
17.38
0.3
13.06
0.3
k
2
13.20
0.3
21.69
0.3
22.99
0.3
20.78
0.3
k
3
11.75
0.3
21.53
0.3
22.49
0.3
21.61
0.3
k
4
17.94
0.3
38.42
0.3
43.66
0.3
34.89
0.3
k
5
6.80
0.3
13.55
0.3
16.12
0.3
11.37
0.3
k
6
6.80
0.3
9.88
0.3
10.55
0.3
9.56
0.3
k
7
11.18
0.3
15.95
1.125
17.10
1.125
15.08
1.125
k
8
5.57
0.3
13.63
2.25
15.01
2.25
12.92
2.25
k
9
11.14
0.3
11.14
2.25
12.16
2.25
10.59
2.25
k
10
8.25
0.3
13.76
0.3
14.77
0.3
13.16
0.3
k
11
m
1
5.58
1.125
12.54
1.125
13.45
1.125
11.71
1.125
k
11
c
1
3.71
0.3
11.33
2.25
11.95
2.25
10.95
2.25
k
11
m
2
7.73
0.3
14.87
0.3
16.39
0.3
14.91
0.3
k
11
c
2
12.43
0.3
14.21
0.3
16.55
0.3
12.31
0.3
k
12
i
8.66
0.3
18.44
0.3
20.67
0.3
16.93
0.3
k
12
8.25
0.3
15.22
0.3
17.52
0.3
13.79
0.3
k
13
m
1
10.93
0.3
21.62
0.3
23.88
0.3
19.88
0.3
k
13
c
1
5.36
0.3
10.57
0.3
10.65
0.3
10.56
0.3
k
13
m
2
11.96
0.3
19.50
0.3
19.63
0.3
19.31
0.3
k
13
c
2
7.01
1.125
15.41
0.3
15.51
0.3
15.46
0.3
k
14
8.66
0.3
18.14
0.3
21.16
0.3
14.95
0.3
k
15
i
6.49
1.125
12.11
1.125
10.14
1.125
14.23
1.125
k
15
14.85
0.3
32.23
0.3
36.54
0.3
29.02
0.3
k
16
9.04
0.3
15.40
0.3
17.48
0.3
13.18
0.3
k
17
m
1
7.01
0.3
13.86
0.3
15.18
0.3
12.85
0.3
k
17
c
1
6.19
0.3
12.30
1.125
13.24
1.125
10.97
1.125
k
17
m
2
3.92
0.3
8.52
0.3
9.44
0.3
8.54
1.125
k
17
c
2
9.95
0.3
20.23
1.125
22.69
1.125
18.10
1.125
k
18
m
9.90
0.3
11.81
0.3
11.73
0.3
11.71
0.3
k
18
c
6.80
0.3
14.18
0.3
15.80
0.3
12.30
0.3The sensitivity coefficients for parameterk1 for the tested factors bn in the interval ]0;2[ for all four tested initial glucose levels 3 × 10−1 gL−1, 1.125 gL−1, 2.25 gL−1, and 4.5 gL−1. (a) The sensitivity coefficients for the number of time steps necessary for the first cell to reach the boundary, (b) the sensitivity coefficients for the total number of cells in this last time step, (c) the sensitivity coefficients for the number of migrating cells, and (d) the sensitivity coefficient for the number of proliferating cells, respectively. Each marker represents the mean result of three simulations with the same parameter scaling with the same initial level of glucose. In all these figures, the black crosses (×) represent the data for an initial glucose level of 3 × 10−1 gL−1, the blue circles (●) for an initial glucose concentration of 1.125 gL−1, and the red diamond (◆) and the yellow plus sign (+) for an initial glucose level of 2.25 gL−1 and 4.5 gL−1, respectively.
(a)
(b)
(c)
(d)In Figure4, the sensitivity coefficients for the four glucose settings are shown for the parameter k1. The range of the factors shown in Figure 4 is limited to factors bn<2. By this means, an equal number of factors are shown that decrease and increase the original parameter k1 by the same intervals, respectively. Furthermore, multiplication of the original parameter with a high factor (5, 10, 50, and 100) tends to result in rather low sensitivity coefficients (compare the definition of Skl,bn and the observations made in Section 3).Figure4(a) shows the relative change of the number of time steps necessary to finish a simulation. In Figure 4(b), the sensitivity coefficients for the total number of tumour cells in the last time step can be seen. Figures 4(c) and 4(d) display the sensitivity coefficients for the number of migrating and proliferating cells in the last time step of the simulation, respectively. In all these figures, the crosses (×) represent the data for an initial glucose level of 3 × 10−1 gL−1, the circles (●) for an initial glucose concentration of 1.125 gL−1, and the diamond (◆) and the plus sign (+) for an initial glucose level of 2.25 gL−1 and 4.5 gL−1, respectively.For all of the four initial glucose settings, factors close to one cause the largest deviation of the four sensitivity coefficientsSk1,bne from zero. However, if the parameter k1 is multiplied with factors bn close to zero or two, this results in low sensitivity coefficients. The absolute maximum value of all sensitivity coefficients Sk1,bne(3)𝔪k1e=max⁡bn,x10(0)|Sk1,bne(x10(0))|
is reached for a low initial glucose concentration of 3 × 10−1 gL−1. With an increasing initial glucose level, the amplitudes of the sensitivity coefficients tend to decrease.The graphs of the sensitivity coefficients for the other parameters look similar to the ones explained above. The main differences affect the range of the sensitivity coefficients and the initial glucose concentration for which the maximum𝔪kle is taken. Therefore, Table 1 provides the maximal absolute sensitivity coefficients 𝔪kle together with the respective glucose values x10(0) that result in these maximal values for all parameters. In contrast to Figure 4, all factors from 0.01 to 100 were explored. The complete results for all reaction parameters and all glucose settings are given in the electronic supplementary material (Figures S2–S13).From the table, it can be seen that, for all sensitivity coefficientsSkl,bne, the maxima of the absolute values 𝔪kle are primarily reached for the lowest glucose level of x10(0)=3×10-1 gL−1 for all parameters. Only a few parameters exist for the absolute maximum 𝔪kle for a medium glucose level (1.125 gL−1 or 2.25 gL−1). The absolute maxima vary between 3.71–17.94 and 8.52–38.42 for Stime and Stotal, respectively. For Smig, the absolute maxima are in the range 9.44–43.66 and the sensitivity coefficients for the number of proliferating cells in the last time step Sprolif take maximum values between 8.54 and 34.89. The maximum over all parameters of all maximal sensitivity coefficients 𝔪kle (max⁡kl(𝔪kle)) is taken for the parameter k4. Except for the endpoint e=time, the parameter k17m2 results in the lowest 𝔪kle.


### 4.2. Correlation between Changes inMtime and Mtotal

Besides investigating the sensitivity coefficientsStime,Stotal,Smig, and Sprolif, we also analysed the correlation between changes in Mtime and Mtotal. Thus, we examined the correlation between tumour expansion velocity and tumour volume for different parameter variations. This is motivated by the fact that the number of time steps necessary for the first tumour cell to reach the boundary Mtime can be considered as a measure for the velocity of the tumour expansion and the total number of tumour cells in the last time step Mtotal is a direct measure for the tumour volume.We define(4)δMkl,bnj:=ΔMkl,bnjMj=Mkl,bnjMj-1,j∈{time,total}
with ΔMkl,bnj and M as in Section 3. If the total number of cells decreases for a certain parameter variation, we can observe that this corresponds to δMkl,bntotal<0. Furthermore, if the variation of parameter kl by multiplication with a factor bn causes an increase in the number of time steps necessary for a simulation to terminate, this comes along with δMkl,bntime>0. In general, δMkl,bnj=a can be translated into Mkl,bnj=(a+1)·M.Figure5 shows the correlation of the number of time steps and the number of tumour cells with δMkl,bntotal on the x-axis and δMkl,bntime on the y-axis. Each marker corresponds to a simulation with a different parameter scaling. As before, the simulations were run for four different glucose settings and the crosses (×, in black) represent the data for an initial glucose level of 3 × 10−1 gL−1, the circles (●, in blue) for an initial glucose concentration of 1.125 gL−1, and the diamond (◆, in red) and the plus sign (+, in yellow) for an initial glucose level of 2.25 gL−1 and 4.5 gL−1, respectively.Figure 5
The correlation between the changes in the number of time steps necessary for a simulation to complete and the number of tumour cells in this last time step for all tested parameter scalings. Each marker represents the mean result of three simulations with the same parameter scaling and with the same initial level of glucose. The black crosses (×) represent the data for an initial glucose level of 3 × 10−1 gL−1, the blue circles (●) for an initial glucose concentration of 1.125 gL−1, and the red diamond (◆) and the yellow plus sign (+) for an initial glucose level of 2.25 gL−1 and 4.5 gL−1, respectively.For each glucose setting, an increase in the number of time steps comes along with an increase in the number of tumour cells. A decrease ofMtime corresponds to a decrease of Mtotal. The increases are most prominent for a low initial glucose level of 3 × 10−1 gL−1. For medium glucose levels (1.125 gL−1 and 2.25 gL−1), one mainly observes an increase in the number of time steps (for a slight increase of Mtotal) and a decrease of the number of tumour cells (along with a slight decrease of Mtime). Assuming a high glucose level (4.5 gL−1), primarily a decrease of the number of time steps and a decrease of the final total number of tumour cells become apparent. None of the parameter scalings results in a significant decrease of the final total number of tumour cells at the same time with an increase in the number of time steps.


### 4.3.  Combined Parameter Changes

The goal of therapies should be to reduce the tumour expansion velocity and to reduce the tumour volume at the same time (i.e., in Figure5, it would be desirable to end up in the upper left quadrant). Therefore, we identified the parameter-factor combinations whose corresponding simulations resulted in δMtime>9 for an initial glucose value of 3 × 10−1 gL−1 (i.e., the marker group in the upper right of Figure 5). Furthermore, we determined the parameter-factor combinations that result in simulations with δMtime<-0.8 assuming an initial glucose level of 4.5 gL−1. The goal is to combine two of these different variations at a time to allow for a semiglobal sensitivity analysis.The respective parameters and factors are given in Table2. The column Factor+ includes those factors whose multiplication with a given parameter resulted in δMtime>9 (for 3 × 10−1 gL−1 glucose concentration), and Factor− refers to those factors whose multiplication with a given parameter resulted in δMtotal<-0.8 (for 4.5 gL−1 glucose concentration). It can be noted that Table 2 lists 14 out of the 31 reaction parameters. The variation of 9 out of these 14 parameters resulted in δMtime>9 (for 3 × 10−1 gL−1 glucose concentration) and δMtotal<-0.8 (for 4.5 gL−1 glucose concentration) by multiplication with different factors. Out of all 20 factors (from 0.01 to 100), the highest and lowest scaling factors (0.01, 0.1, 10, 50, and 100) predominantly resulted in a high number of time steps for a low glucose value or a low total number of cells for a high glucose value. No factors close to one appear in Table 2.Table 2
Parameters and corresponding scaling factors that result in a high number of time steps for a low glucose value (Factor+) or a low number of total cells for a high glucose value (Factor−).
Parameter
Factor+
Factor−
k
1
10, 50, 100
0.01
k
1
i
—
0.01
k
2
10, 50, 100
0.01
k
3
0.01, 0.1
50, 100
k
4
0.01, 0.1
50, 100
k
7
0.01, 0.1
50, 100
k
8
10, 50, 100
—
k
9
0.01, 0.1
50, 100
k
10
10, 50, 100
0.01
k
11
c
1
—
10, 50, 100
k
11
m
2
50, 100
—
k
11
c
2
0.01, 0.1
5, 10, 50, 100
k
12
i
—
0.01, 0.1
k
12
5, 10, 50, 100
0.01, 0.1Next, we explore all possible combinations of two parameter-factor pairs (that are listed in Table2) at a time; that is, simulations are run in which two parameters kl1,kl2 are modified: one parameter (kl1) with an associated factor (bn1) from the Factor+ group and another parameter (kl2) with a corresponding factor (bn2) out of the Factor− group. Obviously, there are some parameter-factor pairs that cannot be combined, since, for example, it is impossible to vary parameter k1 at the same time by the factors 10 and 0.01. In total, 603 combinations are possible and, for each of these combinations, three simulations are carried out for each of the four glucose settings. For each of the simulations, the necessary number of time steps for the first cell to reach the boundary and the total number of cells in the last time step are evaluated and δMtime and δMtotal are calculated based on the mean values. The results are shown in Figure 6 in the same manner as the correlation between Mtime and Mtotal was previously presented in Figure 5 for just one parameter change.Figure 6
The correlation between the changes in the number of time steps necessary for a simulation to complete and the number of tumour cells in this last time step for all tested combined parameter scalings. Each marker represents the mean result of three simulations with the same combination of parameter scalings with the same initial level of glucose. The black crosses (×) represent the data for an initial glucose level of 3 × 10−1 gL−1, the blue circles (●) for an initial glucose concentration of 1.125 gL−1, and the red diamond (◆) and the yellow plus sign (+) for an initial glucose level of 2.25 gL−1 and 4.5 gL−1, respectively.For a low initial glucose concentration (3 × 10−1 gL−1), most variations resulted in a significant increase of the number of time steps (δMtime>0) and of the number of tumour cells (δMtotal>0). For the medium and high glucose levels (1.125 gL−1, 2.25 gL−1, and 4.5 gL−1), some combined parameter variations exist that yield a positive δMtime>0 and at the same time a negative δMtotal<0. These combined parameter variations, therefore, cause a decrease of the tumour expansion velocity and of the tumour volume. For the low initial glucose concentration (3 × 10−1 gL−1), no such combination of parameter changes exists.For initial glucose levels of 1.125 gL−1, 2.25 gL−1, and 4.5 gL−1, Table 3 lists the combinations of parameter scalings (kl1·bn1 combined with kl2·bn2) that resulted in an increase of the number of time steps (δMtime>0.01) along with a decrease of the total number of tumour cells (δMtotal<-0.01). Instead of using exactly zero as the threshold, an additional margin of 0.01 is applied. Thus, combinations are excluded that influence the system output only very slightly.Table 3
The parameter scalings whose combinations result in a decrease of the tumour expansion velocity and the tumour volume. The first and second column list the parameters (kl1) and factors (bn1) whose corresponding variations (kl1·bn1) combined with the second parameter scalings (kl2·bn2) (given for different initial glucose concentrations in the third (1.125 gL−1), fourth (2.25 gL−1), and fifth (4.5 gL−1) column) resulted in δMtime>0.01 and δMtotal<-0.01.
k
l
1
b
n
1
k
l
2
*
b
n
2
1.125 g L−1
2.25 g L−1
4.5 g L−1
k
1
0.01
—
—
k
7
·
100
10
—
k
11
c
2
·
5
—
50
k
10
·
0.01
k
9
·
100; k11c2·10
k
1
i
·
0.01; k9·{50,100}; k10·0.01
k
1
i
0.01
—
k
8
·
50
k
8
·
50
k
2
0.01
—
—
k
4
·
0.01; k7·0.01
50
k
1
i
·
0.01; k9·100; k10·0.01
k
1
i
·
0.01; k9·100; k10·0.01
k
9
·
100; k10·0.01; k11c2·10
k
4
50
k
2
·
100; k3·0.01; k10·100
k
9
·
0.01
k
8
·
100
k
7
50
—
—
k
8
·
100; k9·0.01
k
11
c
2
5
k
8
·
10; k9·0.1;
k
3
·
0.01; k4·0.1;
k
8
·
10; k9·0.1
k
10
·
10; k12·5
k
7
·
0.1; k10·50
10
k
12
·
10
—
—
50
k
8
·
100; k9·0.01
k
8
·
100; k12·50
k
8
·
100; k9·0.01; k10·100; k12·50
100
—
k
12
·
100
—A few combined parameter scalings result in a decrease of the tumour expansion velocity and the tumour volume for all medium and high glucose levels (1.125 gL−1, 2.25 gL−1, and 4.5 gL−1), in particular, as follows:(i)
k2·50 combined with k9·100 or k10·0.01,(ii)
k11c2·50 combined with k8·100.On the one hand, some of the parameter scalings caused a decrease of the tumour expansion velocity and the tumour volume only in combination with a few (1–3) other parameter modifications across all initial glucose concentrations. On the other hand, four parameter scalings (k1·50,k2·50,k4·50, and k11c2·50) could be combined with a quite large selection of other parameter variations (≥5 combinations across all initial glucose settings).


## 4.1. Sensitivity Coefficients

Here, only the results of the sensitivity coefficient calculations for one parameter (k1) are shown in detail (see Figure 4). Results for the other parameters are summarized in Table 1 and in the electronic supplementary material (see Figures S2–S13 in Supplementary Material available online at http://dx.doi.org/10.1155/2014/437094).Table 1
The maximal absolute sensitivity coefficients. This table summarizes for each parameterkl the maximum of the absolute values (𝔪kle) of the four different sensitivity coefficients (Stime,Stotal,Smig, and Sprolif) among all factors bn and all tested initial glucose values (x10(0)). The max⁡(∣·∣) values are dimensionless (d.u.) and the initial glucose concentration values x10(0) are given in g L−1.
k
l
S
time
S
total
S
mig
S
prolif
𝔪
k
l
time
x
10
(
0
)
𝔪
k
l
total
x
10
(
0
)
𝔪
k
l
mig
x
10
(
0
)
𝔪
k
l
prolif
x
10
(
0
)
d.u.
[g L−1]
d.u.
[g L−1]
d.u.
[g L−1]
d.u.
[g L−1]
k
1
5.98
0.3
11.91
0.3
12.28
0.3
11.13
0.3
k
1
i
6.60
0.3
14.88
0.3
17.38
0.3
13.06
0.3
k
2
13.20
0.3
21.69
0.3
22.99
0.3
20.78
0.3
k
3
11.75
0.3
21.53
0.3
22.49
0.3
21.61
0.3
k
4
17.94
0.3
38.42
0.3
43.66
0.3
34.89
0.3
k
5
6.80
0.3
13.55
0.3
16.12
0.3
11.37
0.3
k
6
6.80
0.3
9.88
0.3
10.55
0.3
9.56
0.3
k
7
11.18
0.3
15.95
1.125
17.10
1.125
15.08
1.125
k
8
5.57
0.3
13.63
2.25
15.01
2.25
12.92
2.25
k
9
11.14
0.3
11.14
2.25
12.16
2.25
10.59
2.25
k
10
8.25
0.3
13.76
0.3
14.77
0.3
13.16
0.3
k
11
m
1
5.58
1.125
12.54
1.125
13.45
1.125
11.71
1.125
k
11
c
1
3.71
0.3
11.33
2.25
11.95
2.25
10.95
2.25
k
11
m
2
7.73
0.3
14.87
0.3
16.39
0.3
14.91
0.3
k
11
c
2
12.43
0.3
14.21
0.3
16.55
0.3
12.31
0.3
k
12
i
8.66
0.3
18.44
0.3
20.67
0.3
16.93
0.3
k
12
8.25
0.3
15.22
0.3
17.52
0.3
13.79
0.3
k
13
m
1
10.93
0.3
21.62
0.3
23.88
0.3
19.88
0.3
k
13
c
1
5.36
0.3
10.57
0.3
10.65
0.3
10.56
0.3
k
13
m
2
11.96
0.3
19.50
0.3
19.63
0.3
19.31
0.3
k
13
c
2
7.01
1.125
15.41
0.3
15.51
0.3
15.46
0.3
k
14
8.66
0.3
18.14
0.3
21.16
0.3
14.95
0.3
k
15
i
6.49
1.125
12.11
1.125
10.14
1.125
14.23
1.125
k
15
14.85
0.3
32.23
0.3
36.54
0.3
29.02
0.3
k
16
9.04
0.3
15.40
0.3
17.48
0.3
13.18
0.3
k
17
m
1
7.01
0.3
13.86
0.3
15.18
0.3
12.85
0.3
k
17
c
1
6.19
0.3
12.30
1.125
13.24
1.125
10.97
1.125
k
17
m
2
3.92
0.3
8.52
0.3
9.44
0.3
8.54
1.125
k
17
c
2
9.95
0.3
20.23
1.125
22.69
1.125
18.10
1.125
k
18
m
9.90
0.3
11.81
0.3
11.73
0.3
11.71
0.3
k
18
c
6.80
0.3
14.18
0.3
15.80
0.3
12.30
0.3The sensitivity coefficients for parameterk1 for the tested factors bn in the interval ]0;2[ for all four tested initial glucose levels 3 × 10−1 gL−1, 1.125 gL−1, 2.25 gL−1, and 4.5 gL−1. (a) The sensitivity coefficients for the number of time steps necessary for the first cell to reach the boundary, (b) the sensitivity coefficients for the total number of cells in this last time step, (c) the sensitivity coefficients for the number of migrating cells, and (d) the sensitivity coefficient for the number of proliferating cells, respectively. Each marker represents the mean result of three simulations with the same parameter scaling with the same initial level of glucose. In all these figures, the black crosses (×) represent the data for an initial glucose level of 3 × 10−1 gL−1, the blue circles (●) for an initial glucose concentration of 1.125 gL−1, and the red diamond (◆) and the yellow plus sign (+) for an initial glucose level of 2.25 gL−1 and 4.5 gL−1, respectively.
(a)
(b)
(c)
(d)In Figure4, the sensitivity coefficients for the four glucose settings are shown for the parameter k1. The range of the factors shown in Figure 4 is limited to factors bn<2. By this means, an equal number of factors are shown that decrease and increase the original parameter k1 by the same intervals, respectively. Furthermore, multiplication of the original parameter with a high factor (5, 10, 50, and 100) tends to result in rather low sensitivity coefficients (compare the definition of Skl,bn and the observations made in Section 3).Figure4(a) shows the relative change of the number of time steps necessary to finish a simulation. In Figure 4(b), the sensitivity coefficients for the total number of tumour cells in the last time step can be seen. Figures 4(c) and 4(d) display the sensitivity coefficients for the number of migrating and proliferating cells in the last time step of the simulation, respectively. In all these figures, the crosses (×) represent the data for an initial glucose level of 3 × 10−1 gL−1, the circles (●) for an initial glucose concentration of 1.125 gL−1, and the diamond (◆) and the plus sign (+) for an initial glucose level of 2.25 gL−1 and 4.5 gL−1, respectively.For all of the four initial glucose settings, factors close to one cause the largest deviation of the four sensitivity coefficientsSk1,bne from zero. However, if the parameter k1 is multiplied with factors bn close to zero or two, this results in low sensitivity coefficients. The absolute maximum value of all sensitivity coefficients Sk1,bne(3)𝔪k1e=max⁡bn,x10(0)|Sk1,bne(x10(0))|
is reached for a low initial glucose concentration of 3 × 10−1 gL−1. With an increasing initial glucose level, the amplitudes of the sensitivity coefficients tend to decrease.The graphs of the sensitivity coefficients for the other parameters look similar to the ones explained above. The main differences affect the range of the sensitivity coefficients and the initial glucose concentration for which the maximum𝔪kle is taken. Therefore, Table 1 provides the maximal absolute sensitivity coefficients 𝔪kle together with the respective glucose values x10(0) that result in these maximal values for all parameters. In contrast to Figure 4, all factors from 0.01 to 100 were explored. The complete results for all reaction parameters and all glucose settings are given in the electronic supplementary material (Figures S2–S13).From the table, it can be seen that, for all sensitivity coefficientsSkl,bne, the maxima of the absolute values 𝔪kle are primarily reached for the lowest glucose level of x10(0)=3×10-1 gL−1 for all parameters. Only a few parameters exist for the absolute maximum 𝔪kle for a medium glucose level (1.125 gL−1 or 2.25 gL−1). The absolute maxima vary between 3.71–17.94 and 8.52–38.42 for Stime and Stotal, respectively. For Smig, the absolute maxima are in the range 9.44–43.66 and the sensitivity coefficients for the number of proliferating cells in the last time step Sprolif take maximum values between 8.54 and 34.89. The maximum over all parameters of all maximal sensitivity coefficients 𝔪kle (max⁡kl(𝔪kle)) is taken for the parameter k4. Except for the endpoint e=time, the parameter k17m2 results in the lowest 𝔪kle.


## 4.2. Correlation between Changes inMtime and Mtotal

Besides investigating the sensitivity coefficientsStime,Stotal,Smig, and Sprolif, we also analysed the correlation between changes in Mtime and Mtotal. Thus, we examined the correlation between tumour expansion velocity and tumour volume for different parameter variations. This is motivated by the fact that the number of time steps necessary for the first tumour cell to reach the boundary Mtime can be considered as a measure for the velocity of the tumour expansion and the total number of tumour cells in the last time step Mtotal is a direct measure for the tumour volume.We define(4)δMkl,bnj:=ΔMkl,bnjMj=Mkl,bnjMj-1,j∈{time,total}
with ΔMkl,bnj and M as in Section 3. If the total number of cells decreases for a certain parameter variation, we can observe that this corresponds to δMkl,bntotal<0. Furthermore, if the variation of parameter kl by multiplication with a factor bn causes an increase in the number of time steps necessary for a simulation to terminate, this comes along with δMkl,bntime>0. In general, δMkl,bnj=a can be translated into Mkl,bnj=(a+1)·M.Figure5 shows the correlation of the number of time steps and the number of tumour cells with δMkl,bntotal on the x-axis and δMkl,bntime on the y-axis. Each marker corresponds to a simulation with a different parameter scaling. As before, the simulations were run for four different glucose settings and the crosses (×, in black) represent the data for an initial glucose level of 3 × 10−1 gL−1, the circles (●, in blue) for an initial glucose concentration of 1.125 gL−1, and the diamond (◆, in red) and the plus sign (+, in yellow) for an initial glucose level of 2.25 gL−1 and 4.5 gL−1, respectively.Figure 5
The correlation between the changes in the number of time steps necessary for a simulation to complete and the number of tumour cells in this last time step for all tested parameter scalings. Each marker represents the mean result of three simulations with the same parameter scaling and with the same initial level of glucose. The black crosses (×) represent the data for an initial glucose level of 3 × 10−1 gL−1, the blue circles (●) for an initial glucose concentration of 1.125 gL−1, and the red diamond (◆) and the yellow plus sign (+) for an initial glucose level of 2.25 gL−1 and 4.5 gL−1, respectively.For each glucose setting, an increase in the number of time steps comes along with an increase in the number of tumour cells. A decrease ofMtime corresponds to a decrease of Mtotal. The increases are most prominent for a low initial glucose level of 3 × 10−1 gL−1. For medium glucose levels (1.125 gL−1 and 2.25 gL−1), one mainly observes an increase in the number of time steps (for a slight increase of Mtotal) and a decrease of the number of tumour cells (along with a slight decrease of Mtime). Assuming a high glucose level (4.5 gL−1), primarily a decrease of the number of time steps and a decrease of the final total number of tumour cells become apparent. None of the parameter scalings results in a significant decrease of the final total number of tumour cells at the same time with an increase in the number of time steps.


## 4.3.  Combined Parameter Changes

The goal of therapies should be to reduce the tumour expansion velocity and to reduce the tumour volume at the same time (i.e., in Figure5, it would be desirable to end up in the upper left quadrant). Therefore, we identified the parameter-factor combinations whose corresponding simulations resulted in δMtime>9 for an initial glucose value of 3 × 10−1 gL−1 (i.e., the marker group in the upper right of Figure 5). Furthermore, we determined the parameter-factor combinations that result in simulations with δMtime<-0.8 assuming an initial glucose level of 4.5 gL−1. The goal is to combine two of these different variations at a time to allow for a semiglobal sensitivity analysis.The respective parameters and factors are given in Table2. The column Factor+ includes those factors whose multiplication with a given parameter resulted in δMtime>9 (for 3 × 10−1 gL−1 glucose concentration), and Factor− refers to those factors whose multiplication with a given parameter resulted in δMtotal<-0.8 (for 4.5 gL−1 glucose concentration). It can be noted that Table 2 lists 14 out of the 31 reaction parameters. The variation of 9 out of these 14 parameters resulted in δMtime>9 (for 3 × 10−1 gL−1 glucose concentration) and δMtotal<-0.8 (for 4.5 gL−1 glucose concentration) by multiplication with different factors. Out of all 20 factors (from 0.01 to 100), the highest and lowest scaling factors (0.01, 0.1, 10, 50, and 100) predominantly resulted in a high number of time steps for a low glucose value or a low total number of cells for a high glucose value. No factors close to one appear in Table 2.Table 2
Parameters and corresponding scaling factors that result in a high number of time steps for a low glucose value (Factor+) or a low number of total cells for a high glucose value (Factor−).
Parameter
Factor+
Factor−
k
1
10, 50, 100
0.01
k
1
i
—
0.01
k
2
10, 50, 100
0.01
k
3
0.01, 0.1
50, 100
k
4
0.01, 0.1
50, 100
k
7
0.01, 0.1
50, 100
k
8
10, 50, 100
—
k
9
0.01, 0.1
50, 100
k
10
10, 50, 100
0.01
k
11
c
1
—
10, 50, 100
k
11
m
2
50, 100
—
k
11
c
2
0.01, 0.1
5, 10, 50, 100
k
12
i
—
0.01, 0.1
k
12
5, 10, 50, 100
0.01, 0.1Next, we explore all possible combinations of two parameter-factor pairs (that are listed in Table2) at a time; that is, simulations are run in which two parameters kl1,kl2 are modified: one parameter (kl1) with an associated factor (bn1) from the Factor+ group and another parameter (kl2) with a corresponding factor (bn2) out of the Factor− group. Obviously, there are some parameter-factor pairs that cannot be combined, since, for example, it is impossible to vary parameter k1 at the same time by the factors 10 and 0.01. In total, 603 combinations are possible and, for each of these combinations, three simulations are carried out for each of the four glucose settings. For each of the simulations, the necessary number of time steps for the first cell to reach the boundary and the total number of cells in the last time step are evaluated and δMtime and δMtotal are calculated based on the mean values. The results are shown in Figure 6 in the same manner as the correlation between Mtime and Mtotal was previously presented in Figure 5 for just one parameter change.Figure 6
The correlation between the changes in the number of time steps necessary for a simulation to complete and the number of tumour cells in this last time step for all tested combined parameter scalings. Each marker represents the mean result of three simulations with the same combination of parameter scalings with the same initial level of glucose. The black crosses (×) represent the data for an initial glucose level of 3 × 10−1 gL−1, the blue circles (●) for an initial glucose concentration of 1.125 gL−1, and the red diamond (◆) and the yellow plus sign (+) for an initial glucose level of 2.25 gL−1 and 4.5 gL−1, respectively.For a low initial glucose concentration (3 × 10−1 gL−1), most variations resulted in a significant increase of the number of time steps (δMtime>0) and of the number of tumour cells (δMtotal>0). For the medium and high glucose levels (1.125 gL−1, 2.25 gL−1, and 4.5 gL−1), some combined parameter variations exist that yield a positive δMtime>0 and at the same time a negative δMtotal<0. These combined parameter variations, therefore, cause a decrease of the tumour expansion velocity and of the tumour volume. For the low initial glucose concentration (3 × 10−1 gL−1), no such combination of parameter changes exists.For initial glucose levels of 1.125 gL−1, 2.25 gL−1, and 4.5 gL−1, Table 3 lists the combinations of parameter scalings (kl1·bn1 combined with kl2·bn2) that resulted in an increase of the number of time steps (δMtime>0.01) along with a decrease of the total number of tumour cells (δMtotal<-0.01). Instead of using exactly zero as the threshold, an additional margin of 0.01 is applied. Thus, combinations are excluded that influence the system output only very slightly.Table 3
The parameter scalings whose combinations result in a decrease of the tumour expansion velocity and the tumour volume. The first and second column list the parameters (kl1) and factors (bn1) whose corresponding variations (kl1·bn1) combined with the second parameter scalings (kl2·bn2) (given for different initial glucose concentrations in the third (1.125 gL−1), fourth (2.25 gL−1), and fifth (4.5 gL−1) column) resulted in δMtime>0.01 and δMtotal<-0.01.
k
l
1
b
n
1
k
l
2
*
b
n
2
1.125 g L−1
2.25 g L−1
4.5 g L−1
k
1
0.01
—
—
k
7
·
100
10
—
k
11
c
2
·
5
—
50
k
10
·
0.01
k
9
·
100; k11c2·10
k
1
i
·
0.01; k9·{50,100}; k10·0.01
k
1
i
0.01
—
k
8
·
50
k
8
·
50
k
2
0.01
—
—
k
4
·
0.01; k7·0.01
50
k
1
i
·
0.01; k9·100; k10·0.01
k
1
i
·
0.01; k9·100; k10·0.01
k
9
·
100; k10·0.01; k11c2·10
k
4
50
k
2
·
100; k3·0.01; k10·100
k
9
·
0.01
k
8
·
100
k
7
50
—
—
k
8
·
100; k9·0.01
k
11
c
2
5
k
8
·
10; k9·0.1;
k
3
·
0.01; k4·0.1;
k
8
·
10; k9·0.1
k
10
·
10; k12·5
k
7
·
0.1; k10·50
10
k
12
·
10
—
—
50
k
8
·
100; k9·0.01
k
8
·
100; k12·50
k
8
·
100; k9·0.01; k10·100; k12·50
100
—
k
12
·
100
—A few combined parameter scalings result in a decrease of the tumour expansion velocity and the tumour volume for all medium and high glucose levels (1.125 gL−1, 2.25 gL−1, and 4.5 gL−1), in particular, as follows:(i)
k2·50 combined with k9·100 or k10·0.01,(ii)
k11c2·50 combined with k8·100.On the one hand, some of the parameter scalings caused a decrease of the tumour expansion velocity and the tumour volume only in combination with a few (1–3) other parameter modifications across all initial glucose concentrations. On the other hand, four parameter scalings (k1·50,k2·50,k4·50, and k11c2·50) could be combined with a quite large selection of other parameter variations (≥5 combinations across all initial glucose settings).


## 5. Discussion and Conclusion

Molecular targeted therapies are theoretically a promising approach for glioblastoma treatment. However, so far, none of the existing agents could significantly improve the overall survival [2]. Therefore, further research is inevitable: already well-known signaling networks as well as newly discovered ones need to be examined further for their therapeutic potential. For this purpose, mathematical modelling is a useful tool. In the work at hand, we have presented an approach for identifying influential parameters in a cross-scale model of glioblastoma growth. This previously developed model [41] describes the relevant cellular processes (migration and proliferation) in combination with an intracellular molecular interaction network. A multiscale sensitivity analysis allows us to measure the effect on the tumour growth if reaction parameters on the molecular scale are perturbed. Parameter variations that “positively” influence the tumour progression (i.e., slower expansion velocity and smaller tumour volume) indicate potential targets for new therapeutic interventions. So far, only sensitivity analyses that cover processes on a single scale or that solely variegate the initial conditions for a lung cancer model have been presented [39, 40]. To our knowledge, here, we have presented for the first time a multiscale sensitivity analysis for a model of glioblastoma growth that not only accounts for a systematic variation of all reaction parameters but also considers different nutrient conditions. We could restrict the analysis to an LSA since for therapeutic purposes it is impossible to modify too many parts of a signaling network.Godlewski et al. [9] discovered that microRNA-451 plays an essential role for the nutrient level depending decision of glioma cells to either migrate or proliferate. This miR-451 directly regulates the expression level of the protein MO25 that is relevant for further downstream signaling. Therefore, it was to be expected that those parameters that control the miR-451 and MO25 concentration will have a significant influence on the whole system. Indeed, Table 1 supports this presumption. For all four sensitivity coefficients under examination (Stime as a measure for the expansion velocity, Stotal as a measure for the tumour volume, and Smig and Sprolif as indicators for the build-up of the tumour), the maximum values among all parameter variations are taken for parameters that are relevant for the control of the level of miR-451 and MO25 (k1, k2, and k4).In general, it can be observed (see Figure4 and Figures S2–S13) that, measured in relative terms (cf. the definition of Skl,bn, (2)), small parameter perturbations cause the largest changes for all parameters. For the same parameter-factor combination as well positive as negative sensitivity coefficients can be taken, depending on the initial glucose concentration. Furthermore, since the graphs for all parameters have a very similar structure (no parameter sets itself noticeably apart from the others), further conclusions on the effects of a parameter variation by means of these plots are hardly possible.Nevertheless, on the basis of Table1, it can be observed that the sensitivity coefficients primarily take their absolute maximum values (maximum over all 20 factors and 4 initial glucose values) for a low initial glucose concentration (3 × 10−1 gL−1). Since tumours growing in a low glucose environment expand very fast (which corresponds to a small number of time steps until the first cell reaches the boundary), any changes that only slightly slow down or speed up the expansion are considerably noticeable. The same holds true for the dependence of the tumour volume (corresponding to the number of tumour cells) on the parameter variations.If one does not consider the relative but the absolute changes in the model output (cf. Figure5), the first general observation is that an increase of the number of time steps (δMtime>0) goes along with an increase of the total number of tumour cells (δMtotal>0) and vice versa. An increase of the number of tumour cells by a parameter variation is ultimately caused by an increase of the number of proliferating cells. Proliferation, in turn, is a process that happens much slower than migration and therefore the tumour expansion is slowed down. Thus, this relationship was to be expected in a more or less explicit form. The extreme changes for a low initial glucose condition (3 × 10−1 gL−1) are due to the fact that, in the original setting, the tumour expands very fast and is only very loosely packed, that is, consisting of rather few cells. Therefore, the potential for an increase of the number of time steps and the total number of tumour cells due to an increase of the number of proliferating cells is very high. Consequently, δMtime and δMtotal reach high values. On the contrary, under a high glucose level (4.5 gL−1), the number of time steps and tumour cells is already rather high and it is hardly possible to increase these. However, also, the decrease due to the parameter variations turns out to be only modest. Unfortunately, none of the parameter perturbations resulted in a slower growing tumour that consists of fewer cells (no markers in the upper left quadrant of Figure 5) for any of the glucose settings. Thus, from a single parameter variation, no indicator for a potential therapy could be identified.The identification of those parameters that resulted in the most extreme absolute changes (δMtime<-0.8 and δMtime>9) of the model output (Table 2) yielded mainly parameters that control the “early” parts of the network. That is, primarily, reactions that involve miR-451 and MO25 cause these high variations. Furthermore, it can be noticed that only the multiplication of these parameters with very large or very small factors did change the tumour volume and expansion velocity considerably. Thus, in relative terms, small perturbations of the parameters cause large changes. However, in absolute terms, large parameter modifications are responsible for the large changes.The combined scaling of two selected parameters (see Table2) resulted in an overall similar relationship of the number of time steps and the total number of tumour cells (see Figure 6) as for a single parameter variation: an increase of Mtime coincides with an increase of Mtotal. However, a few pairs of parameter scalings could be identified that, for a medium to high glucose level, cause a less aggressive tumour expansion in the sense that the expansion velocity could be slowed down and the tumour volume could be decreased (δMtime>0.01 and δMtotal<-0.01). For a low initial glucose condition (3 × 10−1 gl−1), none of the variations resulted in the desirable effect. Table 3 further supports the key role of the reactions that regulate the miR-451 and MO25 level for a less aggressive tumour expansion. Most combinations listed in this table contain at least one parameter that controls the reactions that involve either miR-451 or MO25.


### 5.1. Conclusion

For a previously developed multiscale model of glioblastoma growth, we could identify a few combinations of reaction parameter modifications that result in a less aggressive tumour progression. Therefore, these parameters can also be thought of as potential targets for therapies. The parameters are primarily part of reactions that involve the microRNA-451 and the thereof dependent protein MO25 that are responsible for the further downstream signalling. Thus, the analysis supports the relevance of these two molecules for future therapy research.


## 5.1. Conclusion

For a previously developed multiscale model of glioblastoma growth, we could identify a few combinations of reaction parameter modifications that result in a less aggressive tumour progression. Therefore, these parameters can also be thought of as potential targets for therapies. The parameters are primarily part of reactions that involve the microRNA-451 and the thereof dependent protein MO25 that are responsible for the further downstream signalling. Thus, the analysis supports the relevance of these two molecules for future therapy research.



---

*Source: 437094-2014-05-08.xml*