Table of Contents Author Guidelines Submit a Manuscript
Computational and Mathematical Methods in Medicine
Volume 2014, Article ID 437094, 14 pages
http://dx.doi.org/10.1155/2014/437094
Research Article

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

1Institute of Medical Engineering, University of Luebeck, Ratzeburger Allee 160, 23562 Luebeck, Germany
2Graduate School for Computing in Medicine and Life Sciences, University of Luebeck, Ratzeburger Allee 160, 23562 Luebeck, Germany
3Institute for Computational Engineering and Sciences, University of Texas at Austin, 201 East 24th Street, Stop C0200, Austin, TX 78712-1229, USA
4Centre of Excellence for Technology and Engineering in Medicine (TANDEM), Ratzeburger Allee 160, 23562 Luebeck, Germany

Received 16 January 2014; Accepted 22 March 2014; Published 8 May 2014

Academic Editor: Volkhard Helms

Copyright © 2014 Tina A. Schuetz et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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.

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 [68]. 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 the Go 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 [1517]. 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 [1921], or the response to radiotherapy [2224] 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) [2528]. 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 our in silico simulation results to in 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 (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.

437094.fig.001
Figure 1: The molecular interaction network. Different arrow types pointing from one molecular species (given by the variable notation 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 to in 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 () 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 where 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 Section 2.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 () and mTORC1 (). 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 Figure 2. Representative simulation results are shown in Figure 3. In these simulations, a constant initial glucose concentration 3 × 10−1 gL−1, 1.125 gL−1, 2.25 gL−1, 4.5 gL−1} is assumed (a setting comparable to in 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 from in vitro experiments.

437094.fig.002
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).
437094.fig.003
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 output is recorded. Next, we multiplied each reaction parameter () with different factors . For each of these combinations, of a parameter and a factor separate simulations are performed and the model output 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 are calculated as follows (see [37]): with and . 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 coefficients .(1)If (i.e., the parameter is decreased) and (i.e., the model output is decreased as well), then the respective sensitivity coefficient is positive (i.e., ). The same holds true if and . If, on the other hand, and or and , this results in a negative sensitivity coefficient (i.e., ).(2)If the same parameter is separately multiplied with two different factors and for which and if additionally the resulting outcomes and do not differ (i.e., ), then the resulting sensitivity coefficients are the same except for the sign (i.e., ).(3)If the same parameter is separately multiplied with two different factors and that result in the same simulation outcome , then the smaller factor results in the larger respective sensitivity coefficient (i.e., and vice versa).

For the model output , we choose four different endpoints 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 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 ) and vice versa. The final tumour volume is estimated by the number of cells in the last time step . 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 and the number of proliferating cells in the last time step were documented. For all these four model outputs , the corresponding sensitivity coefficients were calculated and were named accordingly , and .

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 Figure 3, 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 . 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 steps , the final number of all tumour cells , and the number of migrating and proliferating cells in the last time step were recorded, and the respective sensitivity coefficients were calculated for all four endpoints . For the ease of notation, the reference to the parameter , the factor , and the initial glucose concentration will be omitted in the following, if the sensitivity coefficients are discussed in general; that is, they will be referred to as .

4.1. Sensitivity Coefficients

Here, only the results of the sensitivity coefficient calculations for one parameter () 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).

tab1
Table 1: The maximal absolute sensitivity coefficients. This table summarizes for each parameter the maximum of the absolute values () of the four different sensitivity coefficients (, and ) among all factors and all tested initial glucose values (). The values are dimensionless (d.u.) and the initial glucose concentration values are given in g L−1.
fig4
Figure 4: The sensitivity coefficients for parameter for the tested factors in the interval 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.

In Figure 4, the sensitivity coefficients for the four glucose settings are shown for the parameter . The range of the factors shown in Figure 4 is limited to factors . By this means, an equal number of factors are shown that decrease and increase the original parameter 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 and the observations made in Section 3).

Figure 4(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 coefficients from zero. However, if the parameter is multiplied with factors close to zero or two, this results in low sensitivity coefficients. The absolute maximum value of all sensitivity coefficients 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 is taken. Therefore, Table 1 provides the maximal absolute sensitivity coefficients together with the respective glucose values 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 coefficients , the maxima of the absolute values are primarily reached for the lowest glucose level of  gL−1 for all parameters. Only a few parameters exist for the absolute maximum 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 and , respectively. For , 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 take maximum values between and . The maximum over all parameters of all maximal sensitivity coefficients () is taken for the parameter . Except for the endpoint , the parameter results in the lowest .

4.2. Correlation between Changes in and

Besides investigating the sensitivity coefficients , and , we also analysed the correlation between changes in and . 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 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 is a direct measure for the tumour volume.

We define with and as in Section 3. If the total number of cells decreases for a certain parameter variation, we can observe that this corresponds to . Furthermore, if the variation of parameter by multiplication with a factor causes an increase in the number of time steps necessary for a simulation to terminate, this comes along with . In general, can be translated into .

Figure 5 shows the correlation of the number of time steps and the number of tumour cells with on the -axis and on the -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.

437094.fig.005
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 of corresponds to a decrease of . 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 ) and a decrease of the number of tumour cells (along with a slight decrease of ). 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 Figure 5, it would be desirable to end up in the upper left quadrant). Therefore, we identified the parameter-factor combinations whose corresponding simulations resulted in 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 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 Table 2. The column Factor+ includes those factors whose multiplication with a given parameter resulted in (for 3 × 10−1 gL−1 glucose concentration), and Factor− refers to those factors whose multiplication with a given parameter resulted in (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 (for 3 × 10−1 gL−1 glucose concentration) and (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.

tab2
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−).

Next, we explore all possible combinations of two parameter-factor pairs (that are listed in Table 2) at a time; that is, simulations are run in which two parameters are modified: one parameter () with an associated factor () from the Factor+ group and another parameter () with a corresponding factor () 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 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 and are calculated based on the mean values. The results are shown in Figure 6 in the same manner as the correlation between and was previously presented in Figure 5 for just one parameter change.

437094.fig.006
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 () and of the number of tumour cells (). 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 and at the same time a negative . 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 ( combined with ) that resulted in an increase of the number of time steps () along with a decrease of the total number of tumour cells (). 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.

tab3
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 () and factors () whose corresponding variations () combined with the second parameter scalings () (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 and .

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) combined with or ,(ii) combined with .

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 (, and ) 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 ( as a measure for the expansion velocity, as a measure for the tumour volume, and and 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 (, , and ).

In general, it can be observed (see Figure 4 and Figures S2–S13) that, measured in relative terms (cf. the definition of , (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 Table 1, 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. Figure 5), the first general observation is that an increase of the number of time steps () goes along with an increase of the total number of tumour cells () 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, and 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 ( and ) 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 Table 2) 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 coincides with an increase of . 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 ( and ). 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.

Appendix

Equations and Parameters Underlying the Molecular Model

In the following, the ODEs are given which determine the molecular interaction model that is shortly described in Section 2.1. Furthermore, in Table 4, the mapping of the variables to the corresponding molecular species is summarized. Finally, Table 5 lists all involved reaction parameters with their original value and unit. More details on the model can be found in [41]:

tab4
Table 4: Molecular species, their variable names, and units.
tab5
Table 5: The reaction parameters (DC: dimensionless constant).

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

Tina A. Schuetz was supported by the Graduate School for Computing in Medicine and Life Sciences funded by Germanys Excellence Initiative [DFG GSC 235/1] and Stefan Becker and Alina Toma were supported by the State Schleswig-Holstein, Germany (Programme for the Future-Economy: 122-09-024).

References

  1. A. M. Dubuc, S. Mack, A. Unterberger, P. A. Northcott, and M. D. Taylor, “The epigenetics of brain tumors,” in Cancer Epigenetics, vol. 863 of Methods in Molecular Biology, pp. 139–153, Humana Press, New York, NY, USA, 2012. View at Google Scholar
  2. A. Quick, D. Patel, M. Hadziahmetovic, A. Chakravarti, and M. Mehta, “Current therapeutic paradigms in glioblastoma,” Reviews on Recent Clinical Trials, vol. 5, no. 1, pp. 14–27, 2010. View at Publisher · View at Google Scholar · View at Scopus
  3. D. Krex, B. Klink, C. Hartmann et al., “Long-term survival with glioblastoma multiforme,” Brain, vol. 130, no. 10, pp. 2596–2606, 2007. View at Publisher · View at Google Scholar · View at Scopus
  4. T. A. Dolecek, J. M. Propp, N. E. Stroup, and C. Kruchko, “CBTRUS statistical report: primary brain and central nervous system tumors diagnosed in the United States in 2005–2009,” Neuro-Oncology, vol. 14, supplement 5, pp. v1–v49, 2012. View at Google Scholar
  5. F. B. Furnari, T. Fenton, R. M. Bachoo et al., “Malignant astrocytic glioma: genetics, biology, and paths to treatment,” Genes and Development, vol. 21, no. 21, pp. 2683–2710, 2007. View at Publisher · View at Google Scholar · View at Scopus
  6. A. Esquela-Kerscher, “MicroRNAs function as tumor suppressor genes and oncogenes,” in MicroRNAs in Development and Cancer, vol. 1 of Molecular Medicine and Medicinal Chemistry, 6, pp. 149–184, Imperial College Press, London, UK, 2011. View at Google Scholar
  7. B. Kefas, J. Godlewski, L. Comeau et al., “microRNA-7 inhibits the epidermal growth factor receptor and the akt pathway and is down-regulated in glioblastoma,” Cancer Research, vol. 68, no. 10, pp. 3566–3572, 2008. View at Publisher · View at Google Scholar · View at Scopus
  8. O. Kovalchuk, J. Filkowski, J. Meservy et al., “Involvement of microRNA-451 in resistance of the MCF-7 breast cancer cells to chemotherapeutic drug doxorubicin,” Molecular Cancer Therapeutics, vol. 7, no. 7, pp. 2152–2159, 2008. View at Publisher · View at Google Scholar · View at Scopus
  9. J. Godlewski, M. O. Nowicki, A. Bronisz et al., “MicroRNA-451 regulates LKB1/AMPK signaling and allows adaptation to metabolic stress in glioma cells,” Molecular Cell, vol. 37, no. 5, pp. 620–632, 2010. View at Publisher · View at Google Scholar · View at Scopus
  10. R. J. Shaw, “LKB1 and AMP-activated protein kinase control of mTOR signalling and growth,” Acta Physiologica, vol. 196, no. 1, pp. 65–80, 2009. View at Publisher · View at Google Scholar · View at Scopus
  11. A. Alexander and C. L. Walker, “The role of LKB1 and AMPK in cellular responses to stress and damage,” FEBS Letters, vol. 585, no. 7, pp. 952–957, 2011. View at Publisher · View at Google Scholar · View at Scopus
  12. A. Giese, M. A. Loo, N. Tran, D. Haskett, S. W. Coons, and M. E. Berens, “Dichotomy of astrocytoma migration and proliferation,” International Journal of Cancer, vol. 67, pp. 275–282, 1996. View at Google Scholar
  13. L. Edelstein-Keshet, Mathematical Models in Biology, Society for Industrial and Applied Mathematics, Philadelphia, Pa, USA, 2005.
  14. N. Bellomo, N. K. Li, and P. K. Maini, “On the foundations of cancer modelling: selected topics, speculations, and perspectives,” Mathematical Models and Methods in Applied Sciences, vol. 18, no. 4, pp. 593–646, 2008. View at Publisher · View at Google Scholar · View at Scopus
  15. L. H. Abbott and F. Michor, “Mathematical models of targeted cancer therapy,” British Journal of Cancer, vol. 95, no. 9, pp. 1136–1141, 2006. View at Publisher · View at Google Scholar · View at Scopus
  16. N. Bellomo, M. Chaplain, and E. D. Angelis, Eds., Selected Topics in Cancer Modeling—Genesis, Evolution, Immune Competition and Therapy, Birkhäauser, Boston, Mass, USA, 2008.
  17. J. Clairambault, “Can theorems help treat cancer?” Journal of Mathematical Biology, vol. 66, no. 7, pp. 1555–1558, 2013. View at Publisher · View at Google Scholar · View at Scopus
  18. E. Konukoglu, O. Clatz, B. H. Menze et al., “Image guided personalization of reaction-diffusion type tumor growth models using modified anisotropic eikonal equations,” IEEE Transactions on Medical Imaging, vol. 29, no. 1, pp. 77–95, 2010. View at Publisher · View at Google Scholar · View at Scopus
  19. C. Hogea, G. Biros, F. Abraham, and C. Davatzikos, “A robust framework for soft tissue simulations with application to modeling brain tumor mass effect in 3D MR images,” Physics in Medicine and Biology, vol. 52, no. 23, pp. 6893–6908, 2007. View at Publisher · View at Google Scholar · View at Scopus
  20. S. Becker, A. Mang, A. Toma, and T. M. Buzug, “In-silico oncology: an approximate model of brain tumor mass effect based on directly manipulated free form deformation,” International Journal of Computer Assisted Radiology and Surgery, vol. 5, no. 6, pp. 607–622, 2010. View at Publisher · View at Google Scholar · View at Scopus
  21. A. Mang, A. Toma, T. A. Schuetz, S. Becker, and T. M. Buzug, “A generic framework for modeling brain deformation as a constrained parametric optimization problem to aid non-diffeomorphic image registration in brain tumor imaging,” Methods of Information in Medicine, vol. 51, pp. 429–440, 2012. View at Google Scholar
  22. D. D. Dionysiou, G. S. Stamatakos, N. K. Uzunoglu, K. S. Nikita, and A. Marioli, “A four-dimensional simulation model of tumour response to radiotherapy in vivo: parametric validation considering radiosensitivity, genetic profile and fractionation,” Journal of Theoretical Biology, vol. 230, no. 1, pp. 1–20, 2004. View at Publisher · View at Google Scholar · View at Scopus
  23. R. Rockne, E. C. Alvord Jr., M. Szeto, S. Gu, G. Chakraborty, and K. R. Swanson, “Modeling diffusely invading brain tumors: an individualized approach to quantifying glioma evolution and response to therapy,” in Selected Topics in Cancer Modeling, pp. 207–221, Birkhäauser, Boston, Mass, USA, 2008. View at Google Scholar
  24. A. Bertuzzi, C. Bruni, F. Papa, and C. Sinisgalli, “Optimal solution for a cancer radiotherapy problem,” Journal of Mathematical Biology, vol. 66, pp. 311–349, 2012. View at Publisher · View at Google Scholar · View at Scopus
  25. A. R. A. Anderson, A. M. Weaver, P. T. Cummings, and V. Quaranta, “Tumor morphology and phenotypic evolution driven by selective pressure from the microenvironment,” Cell, vol. 127, no. 5, pp. 905–915, 2006. View at Publisher · View at Google Scholar · View at Scopus
  26. P. Macklin and J. Lowengrub, “Nonlinear simulation of the effect of microenvironment on tumor growth,” Journal of Theoretical Biology, vol. 245, pp. 677–704, 2007. View at Google Scholar
  27. M. Robertson-Tessi, A. El-Kareh, and A. Goriely, “A mathematical model of tumor-immune interactions,” Journal of Theoretical Biology, vol. 294, pp. 56–73, 2012. View at Publisher · View at Google Scholar · View at Scopus
  28. A. Toma, A. Mang, T. A. Schuetz, S. Becker, and T. M. Buzug, “A novel method for simulating the extracellular matrix in models of tumour growth,” Computational and Mathematical Methods in Medicine, vol. 2012, Article ID 109019, 11 pages, 2012. View at Publisher · View at Google Scholar
  29. W. W. Chen, B. Schoeberl, P. J. Jasper et al., “Input-output behavior of ErbB signaling pathways as revealed by a mass action model trained against dynamic data,” Molecular Systems Biology, vol. 5, article 239, 2009. View at Publisher · View at Google Scholar · View at Scopus
  30. J. E. Purvis, A. J. Shih, Y. Liu, and R. Radhakrishnan, “Cancer cell: linking oncogenic signaling to molecular structure,” in Multiscale Cancer Modeling, pp. 31–44, CRC Press, Taylor & Francis Group, London, UK, 2011. View at Google Scholar
  31. T. S. Deisboeck, Z. Wang, P. MacKlin, and V. Cristini, “Multiscale cancer modeling,” Annual Review of Biomedical Engineering, vol. 13, pp. 127–155, 2011. View at Publisher · View at Google Scholar · View at Scopus
  32. C. Athale, Y. Mansury, and T. S. Deisboeck, “Simulating the impact of a molecular “decision-process” on cellular phenotype and multicellular patterns in brain tumors,” Journal of Theoretical Biology, vol. 233, no. 4, pp. 469–481, 2005. View at Publisher · View at Google Scholar · View at Scopus
  33. L. Zhang, C. A. Athale, and T. S. Deisboeck, “Development of a three-dimensional multiscale agent-based tumor model: simulating gene-protein interaction profiles, cell phenotypes and multicellular patterns in brain cancer,” Journal of Theoretical Biology, vol. 244, no. 1, pp. 96–107, 2007. View at Publisher · View at Google Scholar · View at Scopus
  34. Z. Wang, L. Zhang, J. Sagotsky, and T. S. Deisboeck, “Simulating non-small cell lung cancer with a multiscale agent-based model,” Theoretical Biology and Medical Modelling, vol. 4, article 50, 2007. View at Publisher · View at Google Scholar · View at Scopus
  35. I. Ramis-Conde, D. Drasdo, A. R. A. Anderson, and M. A. J. Chaplain, “Modeling the influence of the E-cadherin-β-catenin pathway in cancer cell invasion: a multiscale approach,” Biophysical Journal, vol. 95, no. 1, pp. 155–165, 2008. View at Publisher · View at Google Scholar · View at Scopus
  36. G. Goel, I.-C. Chou, and E. O. Voit, “System estimation from metabolic time-series data,” Bioinformatics, vol. 24, no. 21, pp. 2505–2511, 2008. View at Publisher · View at Google Scholar · View at Scopus
  37. H. Rabitz, M. Kramer, and D. Dacol, “Sensitivity analysis in chemical kinetics,” Annual Review of Physical Chemistry, vol. 34, pp. 419–461, 1983. View at Google Scholar
  38. Y. Zhang and A. Rundell, “Comparative study of parameter sensitivity analyses of the TCR-activated Erk-MAPK signalling pathway,” IEE Proceedings: Systems Biology, vol. 153, no. 4, pp. 201–211, 2006. View at Publisher · View at Google Scholar · View at Scopus
  39. Z. Wang, V. Bordas, J. Sagotsky, and T. S. Deisboeck, “Identifying therapeutic targets in a combined EGFR-TGFβR signalling cascade using a multiscale agent-based cancer model,” Mathematical Medicine and Biology, vol. 29, no. 1, pp. 95–108, 2012. View at Publisher · View at Google Scholar · View at Scopus
  40. Z. Wang, C. M. Birch, and T. S. Deisboeck, “Cross-scale sensitivity analysis of a non-small cell lung cancer model: linking molecular signaling properties to cellular behavior,” BioSystems, vol. 92, no. 3, pp. 249–258, 2008. View at Publisher · View at Google Scholar · View at Scopus
  41. T. A. Schuetz, S. Becker, A. Mang, A. Toma, and T. M. Buzug, “Modelling of glioblastoma growth by linking a molecular interaction network with an agent based model,” Mathematical and Computer Modelling of Dynamical Systems, vol. 19, pp. 417–433, 2013. View at Google Scholar
  42. T. A. Schuetz, S. Becker, A. Mang, A. Toma, and T. M. Buzug, “A computational multiscale model of glioblastoma growth: regulation of cell migration and proliferation via microRNA-451, LKB1 and AMPK,” in Proceedings of the Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC '12), pp. 6620–6623, 2012.
  43. J. H. Lee, H. Koh, M. Kim et al., “Energy-dependent regulation of cell structure by AMP-activated protein kinase,” Nature, vol. 447, no. 7147, pp. 1017–1020, 2007. View at Publisher · View at Google Scholar · View at Scopus
  44. J. Sabari, D. Lax, D. Connors et al., “Fibronectin matrix assembly suppresses dispersal of glioblastoma cells,” PLoS ONE, vol. 6, no. 9, Article ID e24810, 2011. View at Publisher · View at Google Scholar · View at Scopus