Critical parameters determining standard radiotherapy treatment outcome for glioblastoma multiforme: a computer simulation.

The aim of this paper is to investigate the most critical parameters determining radiotherapy treatment outcome in terms of tumor cell kill for glioblastoma multiforme tumors by using an already developed simulation model of in vivo tumor response to radiotherapy.

INTRODUCTION

Radiation therapy plays a central role in the management of cancer. Radiation is considered now a well understood cytotoxic agent [1: pg. 7]. Apart from advances in the so-called physical treatment planning, which refers to techniques used for improving radiation dosimetry and dose distribution, it is now widely recognized that further improvements in radiotherapy treatment will most probably require increasing grounding in fundamental radiobiological mechanisms [1].

Cancer biology plays a critical role in determining radiotherapy treatment outcome. In this context, efficient modeling of tumor response to radiotherapy can contribute to the elucidation of the involved biological mechanisms and to the emergence of patient-individualized treatments, on the presupposition of adequate clinical adaptation and testing.

During the last years researchers have enhanced the understanding of tumor behavior by means of various simulation models. Representative examples drawn from the extensive corresponding literature have been given in [2]. Analytical mathematical models are valuable tools to study some aspects of tumor radiobiology. Discrete time algorithmic descriptions (computer simulations) present the advantage of great adaptability in treating complex situations. They offer the possibility of accounting for a large number of mechanisms and interactions and they are therefore particularly suitable to describe in vivo tumor growth and response to irradiation. The stochastic nature of the involved biological phenomena favors the choice of stochastic techniques such as the generic Monte Carlo method for oncological simulations. An extensive presentation of computer simulation models of particular interest has been given in [3].

The goal of this study was an investigation of the relative impact of the most critical parameters determining radiotherapy treatment outcome in terms of tumor cell kill for glioblastoma multiforme (GBM) tumors by performing adequate simulations using an already developed simulation model of in vivo GBM response to radiotherapy. For each of the selected critical parameters a series of simulation runs covering the whole range of values that have appeared in the literature for GBM have been performed, while adjusting the remaining parameters at the most typical GBM values. GBM is a highly aggressive type of brain tumor. Prognosis for patients with GBM remains dismal despite efforts to improve current therapies and develop novel clinical approaches [4].

The comparison of the simulation results with clinical experience and experimental knowledge was selected as a means to reveal and substantiate the potential and flexibility of the particular simulation model in order to study biological phenomena related to cancer and in the long-term serve as a patient individualized treatment optimization tool, following a strict clinical validation procedure.

MATERIALS AND METHODOLOGY

Brief Outline of the Simulation Model

The In Silico Oncology Group (ISOG) simulation model of imageable GBM response to radiotherapy, which has been used for the purposes of this study, is based on the clinical, imaging, histopathologic, and molecular data of the patient and incorporates numerous fundamental biological mechanisms. The model is characterized by the combination of the following features: (i) possibility for simulation of both untreated in vivo tumor growth and in vivo tumor response to radiotherapy; (ii) possibility for consideration and use of the actual imaging, histopathologic and molecular data for each particular clinical case (clinical orientation); (iii) incorporation of numerous biological mechanisms by means of an explicit algorithmical description; (iv) introduction of the notion of the ‘‘geometrical cell’’ and its constituent compartments, called ‘‘equivalence classes’’, corresponding to discrete phases within or out of the cell cycle; (v) extensive use of random number generators to simulate the stochastic nature of various biological phenomena involved (Monte Carlo approach); (vi) discrete and modular character, which confers a high level of adaptability; (vii) possibility for 3D reconstruction and visualization of the results. A detailed description of the simulation model can be found in [2, 3, 5]. For the sake of clarity and completeness of information regarding the methods used in this work a brief outline of the model is presented below.

The clinician delineates the tumor and its inner regions of interest by using a dedicated computer tool. For GBM, high signal intensity in T1-weighted contrast-enhanced MRI with gadolinium corresponds to regions of actively proliferating tumor cells, while low intensity regions correspond to necrotic regions [4: pg. 109]. A 1mm-thick layer surrounding the necrotic region is assumed to correspond to the resting cells’ region. The distribution of the absorbed dose in the region of interest according to the treatment plan is also acquired.

A prototype system of quantizing cell clusters included within each elementary cubic volume of a discretizing mesh covering the anatomic area of interest lies at the heart of the simulation approach. During a simulation the geometrical mesh is scanned with a time step of an hour. The elementary cubic volume of the mesh is called “Geometrical Cell (GC)”. In each time step, the updated state of a given GC is determined on the basis of a number of algorithms describing the behavior of the cells constituting the tumor. More specifically, each GC of the mesh belonging to the tumor contains cells, which are distributed in a number of “classes” (compartments), each one characterized by the phase in which its cells are found (within or out of the cell cycle: G1, S, G2, M, G0, Necrosis, Apoptosis). Specially designed stochastic cellular automata describe tumor cell kinetics, by incorporating the following biological phenomena:

Cycling of proliferating cells through the subsequent phases of the cell cycle.

Spontaneous apoptosis.

Transition to a dormant (G0) phase due to inadequate supply with oxygen and nutrients.

Local reoxygenation and nutrient provision reestablishment.

Cell death through necrosis due to prolonged oxygen and nutrients deprivation.

Radiotherapy-induced cell death: possibility of distinction between two mechanisms of radiation-induced cell death, namely apoptotic and necrotic cell death.

Apoptotic cell death can be furthermore subdivided into radiation-induced interphase death (direct death through apoptosis) and radiation-induced mitotic apoptotic death. In most solid tumors, though, the majority of lethally damaged cells dies through a radiation-induced mitotic necrotic mechanism and is considered to undergo a few mitotic divisions prior to death and disappearance from the tumor (two mitotic divisions assumed in the current version of the simulation model) [1, 6].

The probability of cell death taking place through each specific mechanism depends primarily on cell type [7, 8].

Cell killing by irradiation is described by the Linear Quadratic or LQ Model, which is widely used in the pertinent literature [1]:

where S(D) is the surviving fraction after a (uniform) dose D (Gy) of radiation to a population of cells. The parameters α (alpha) (Gy-1) and β (beta) (Gy-2) are called the radiosensitivity parameters of the LQ model. Cell radiosensitivity varies considerably throughout the cell cycle [1, 9, 10]. The S phase is regarded as the most resistant. Cells in any proliferating phase are more radiosensitive than hypoxic cells residing in G0. Based on these observations the model currently uses different values for the radiosensitivity parameters of the LQ model for the S phase (αs,βs), the remaining proliferating phases G1, G2, M (αp,βp), and the G0 phase (αG0, βG0). The values of αs , βs and αG0, βG0 can be derived as perturbations of the (αp,βp) values, as for example in: αG0= αp /OER, βG0= βp /OER2, αs = 0.6 αp + 0.4αG0, βs = 0.6βp + 0.4βG0, where OER is the Oxygen Enhancement Ratio [9].

Tumor expansion or shrinkage is simulated on the basis of the following algorithms: If the number of tumor cells contained within a given GC drops below a given threshold, then a procedure which attempts to “unload” the remaining cells in the neighboring GCs takes place. Cells are preferentially placed within the neighboring GCs with the maximum available free space. If the given GC becomes empty, it is assumed to disappear from the tumor, in which case an appropriate shift of a chain of GCs, intended to fill the “vacuum”, leads to tumor shrinkage. This can happen after the killing of a number of cells by irradiation. On the other hand, if the number of alive and dead cells within a given GC exceeds a given threshold, then a similar procedure attempting to unload the excess cells in the surrounding GCs takes place. In case that the unloading procedure fails to sufficiently reduce the number of cells, then a new GC emerges. Its position relative to the “mother” GC is determined using a random number generator. An appropriate shifting of a chain of adjacent GCs leads to a differential expansion of the tumor.

A simplified flow diagram of the simulation model is depicted in Fig. ().

RESULTS

Literature Review

A thorough study of glioblastoma multiforme literature preceded the simulations, so as to define the real clinical range of the values of various model parameters. The parameters selected for the study of their relative impact in GBM’s response, based on accumulated clinical and experimental knowledge, were the α and β radiosensitivity parameters of the LQ model, the cell cycle duration (TC), the cell loss factor (CLF) and the growth fraction (GF) of the tumor. These constitute three parameters characterizing the kinetic behavior of tumors (TC, CLF, GF) and two parameters (α and β) widely used to describe tumor cell radiosensitivity. These parameters reflect characteristics of a tumor that are determined largely by the underlying genetic profile of the tumor cells.

More specifically, Table presents the values reported in the literature for the cell cycle duration of GBM cell lines, xenografts or primary tumor specimens. The corresponding references are [11, 12-17]. The distribution of the total cell cycle duration to the various phases constituting the cell cycle (G1, S, G2, M) has been assumed to follow the distribution in a typical malignant cell as described in [18]: G1≈TC40%, S≈ TC39%, G2≈TC19%, M≈ TC2%.

It should be noted that different distributions of cell cycle phase durations have been also encountered in the relevant literature. Many authors state that the varying length of G1 phase accounts for most of the differences reported in the total cell cycle duration between different cell types or between cells growing under different conditions, while S and M phases are normally quite fixed in duration [19]. A comparative study of the effect of the adoption of this alternative assumption has been also performed.

Table presents the results of the literature search for values of the α and β parameters of the LQ model for GBM, which are generally considered to reflect tumor cell intrinsic radiosensitivity. The corresponding references are [10, 20-24]. Cellular sensitivity to irradiation reflects a culmination of numerous molecular pathways including DNA repair and cell cycle checkpoint fidelity.

Table presents the values found in the literature for the cell loss factor of GBM, which is defined as the cell loss rate as a fraction of the cell birth rate [1, 13]. The corresponding references are [25-28]. The concept of the cell loss factor emphasizes that a tumor is a dynamic system and that its net growth is the cumulative effect of competing processes of cell proliferation, death and removal. The simulation model considers necrosis and apoptosis as the basic mechanisms of cell loss and assumes that the total cell loss factor is the sum of the cell loss factor due to necrosis and the cell loss factor due to apoptosis. Furthermore, in the particular case of GBM the contribution of necrosis to the total cell loss is considered far more important than that of apoptosis (necrosis 10-fold higher), as gliomas exhibit in general low levels of spontanteous or radiation-induced apoptotic cells [29-31]. It should also be noted that cell loss factor values are considered fixed throughout radiotherapy in the simulation model. This is also a reasonable first approximation for tumors that exhibit little or no regression in post-irradiation imaging studies, such as GBM [10]. This is thought to happen in cases where the blood flow to the tumor does not change during therapy and is usually observed for very slow growing tumors or those that have low cell loss factors at initiation of treatment, such as GBM [6, 32].

Finally, Table presents the values found in the literature for the growth fraction of GBM tumours, i.e. the fraction of actively cycling tumor cells. The corresponding references are [11, 13, 15-17, 33].

SIMULATION RESULTS

For each of the selected critical parameters a series of simulation runs covering the whole range of values that have been reported in the literature for GBM have been performed, while adjusting the remaining parameters at the most typical GBM values. For example, typical parameter values for GBM are: TC = 40h [11], α=0.36Gy-1, β=0.02Gy-2 [10, 34], CLF=0.3 [25], GF=25% [11]. The duration of the G0 phase has been adjusted to 25h [35] and the oxygen enhancement ratio has been taken equal to 3 [36, 37]. The initial tumor was a typical case of a 34.55cm3 GBM delineated based on MRI imaging data. The initial tumor volume corresponds to a total of 34550 GCs of 1mm3 each. Assuming a typical cell density of 106cells/mm3 [1:9-10] this corresponds to an initial number of cells equal to 3.455×1010 (or approximately 10gr). The outer boundary of the tumor and a central necrotic area were evident on the T1 imaging data, while an intermediate region with increased percentage of G0 cells has been assumed to surround the necrotic area in accordance with general biological knowledge. A central slice of the tumor is depicted in Fig. ().

Figs. () depict the simulation results in the form of the number of living tumor cells as a function of time. The time point t=0 corresponds to the start of the radiotherapy treatment. The standard fractionation scheme of 2Gy/day, 5 days/week, 60Gy in total, no irradiation during weekends, was used in all simulation cases.

According to the results obtained for the parameter ranges considered, the parameters with the major impact on the tumor’s response to radiotherapy were the alpha parameter of the LQ model and the cell cycle duration of tumor cells, while the effect of the other parameters was less pronounced. More specifically, starting with an otherwise identical tumor in all cases, by varying the values of the alpha parameter the resulting number of living tumor cells one week after the end of treatment (t=6 weeks) spans an interval of approximately 10 orders of magnitude (logs). The corresponding interval in the cases of TC, beta, CLF and GF parameters is approximately of the order of 7, 2, 1.5 and 0.5 logs correspondingly.

In Fig. () the effect of adopting the two alternative distributions for the cell cycle phase durations that have been discussed above, is depicted by means of a comparative simulation. More specifically, the values of TC=40h and TC=50h are considered for this exploratory simulations. The consideration of TC equal to 50h by attributing the extra duration compared to TC=40h in the G1 phase (G1:26h, S:15h, G2:8h, M:1h), instead of distributing the extra time in all phases according to [18] (G1:20h, S:19h, G2:10h, M:1h), seems to enhance the tumor cell kill effect, since G1 has been assumed as a relatively radiosensitive phase of the cell cycle.

DISCUSSION

The critical importance of the alpha parameter of the LQ model and of the cell cycle duration as major determinants of tumor response to radiotherapy is well documented in relevant literature [1:242, 38]. The cell cycle duration is directly linked to tumor cell proliferation during treatment. Rapidly proliferating tumors might be expected to repopulate more rapidly during treatment and thus be relatively resistant to some forms of treatment. Of course the proliferative potential of a tumor apart from the cell cycle time depends on the growth fraction and on the cell loss rate as well.

The linear component in the LQ model (α coefficient) is related to the initial slope of the survival curve and is thought to be due to single-track events. The quadratic component (β coefficient) is related to the slope of the survival curve at higher doses, and is considered to reflect sub-lethal damage repair. The ratio α/β is the dose for which the contribution of the linear and the quadratic terms are equal and is descriptive of the rate of the survival curve’s bending, with smaller values corresponding to higher rates. Clinically relevant radiation doses (around 2Gy, as in standard radiotherapy treatment for GBM) are generally dominated by the alpha component of the clonogenic survival curve, to which differences between radiosensitive and radioresistant cells are mainly attributed. For such dose fractions, the beta parameter seems to present little variation between different cells [1: pg. 190, 9: pg. 23-24, 91-95, 10, 38].

For example, in [39] the effect of radiosensitivity on tumor response to radiotherapy was examined by varying alpha from 0.2 to 0.35 to 0.5Gy-1 with α/β fixed at 10Gy. The model assumed that the bulk tumor was composed of viable cells and “doomed” cells and used a modified LQ model to compute “doomed” cells for a typical tumor of 5.8cm diameter (100gr), irradiated by a radiotherapy treatment of 30 fractions of 2Gy, given at a rate of 5 fractions per week (a fraction is delivered every 1.4days). The LQ model modification refers to the integration into the model of a description of tumor cell proliferation, whereby the tumor grows exponentially when its size is below a threshold value and in Gompertzian mode when its size is greater than the threshold. The tumor cell population is assumed homogeneous in terms of radiosensitivity.

The author of this study concludes that radiosensitivity has a very significant effect on the nadir of the surviving fraction. The depth of nadir of the surviving fraction increases with increasing radiosensitivity. The survival fraction interval at the end of treatment is approximately 9.5logs. This could be considered a rather good agreement with the results produced in our study taking into account the differences in the underlying simulation models. Furthermore after a parameter sensitivity analysis performed by considering 1% changes in the values of the model parameters, the author reaches the conclusion that radiosensitivity was, by far, the most dominant parameter, followed by proliferation rate and doomed cell loss rate, which is also in agreement with the results presented here.

On the other hand, the small contribution of the beta effect to the overall radiosensitivity of a tumour has been frequently reported in the literature [1: pg. 190, 9: pg. 24].

As far as the variation in the growth fraction of GBM tumors is concerned, it is observed from Fig. () that, as expected, as time increases the absolute difference between the numbers of surviving cells is constantly decreasing. Therefore, it is becoming progressively easier for the curves corresponding to different growth fractions to coincide at certain points if the stochasticity of the mathematical treatment is taken into account. It is obvious that after a point at which the number of surviving cells of two initially different tumors becomes the same their further relative time course is random given that the same treatment scheme is applied. It should be noted that different values of the growth fraction are computed in the model by appropriately adjusting the number of dead cells so as to have the same total number of cells (i.e. the hypoxic fraction remains the same).

To summarize, the relative impact on GBM response to standard radiotherapy treatment of three major parameters characterizing the kinetic behaviour of tumors (TC, GF, CLF) and two parameters widely used to describe radiosensitivity (α and β parameters of the LQ model) has been studied by means of comparative simulations with the help of a four-dimensional (in time and space) computer simulation model. A thorough literature review to investigate the real clinical ranges of these parameters for GBM tumors has been performed. According to the simulation results obtained, the parameters with the major impact on the tumor’s response to radiotherapy were the α parameter of the LQ model and the cell cycle duration of tumor cells, while the effect of the other parameters was less pronounced. In general, the simulation results are supported by clinical experience and reproduce the dismal prognosis for GBM in terms of tumor cell kill. Clinical reality has been reproduced in the sense that for almost all combinations of parameter values (apart from α=0.61Gy-1) reported in the literature for GBM, the standard radiotherapy treatment seems to be inadequate to kill all tumor cells that have been initially assumed for a 10gr-tumor.

It should be stressed that the values of these parameters that have been derived from literature reflect the specific molecular profiles of the cell lines or tumors for which they are computed, and thereby provide a means to reflect into the simulation model the effect of various genetic determinants of tumor response to radiotherapy. For example, using the values of the alpha and beta parameters found in [10], which refer to two isogenic cell lines differing only in their p53 status and performing comparative simulations, the influence of p53 gene status can be reflected in the simulations results, as has been performed in [5].

It should also be noted that as first approximation tumor cell heterogeneity in terms of cell cycle time or intrinsic radiosensitivity in the same cell cycle phase has not been taken into account. In addition, in the current version of the model all proliferating cells are assumed to be clonogenic; cells of limited proliferative potential or terminally differentiated cells are ignored. This is a logical first approximation for GBM which is a poorly differentiated tumor [25, 40]. Moreover, differentiated cells are generally not a danger to the patient, although their bulk contributes to the tumor volume. Nevertheless, versions of the model incorporating the distinction between stem cells, cells of limited proliferative potential and terminally differentiated cells are under preparation.

A further note of concern is the accuracy of the experimental values used in the simulation, which depends on the model system and the experimental procedure used for the measurement in the corresponding studies from which they were derived. Several causes of discrepancies between parameter values computed for established cell lines or tumor xenografts and real clinical tumors have been frequently reported in the literature [20, 26, 38, 41]. For example, established cancer cell lines some times exhibit different characteristics compared to in vivo tumors mainly as a result of the controlled environment of the experiments or the genetic drift observed over time.

CONCLUSIONS

Notwithstanding the uncertainties and simplifications mentioned in previous paragraphs, since the simulation results generally reproduce clinical reality, they thereby substantiate the clinical potential, the flexibility and the robustness of the simulation model, especially in view of its clinically oriented design. Its usefulness for gaining insight into critical biological mechanisms seems straightforward. From a wider perspective, the models’ translation into the clinical environment in order to serve as a patient individualized treatment optimization tool, following a strict validation procedure based on clinical data, can be envisaged. A similar process of clinical translation of computer simulation models of the ISOG group for the case of chemotherapy treatment of breast cancer and Wilms’ tumors is currently under way in the framework of the EC funded project “ACGT: Advancing Clinicogenomic Trials on Cancer” (FP6-2005-IST-026996).

A simplified case study which could characteristically explain the potential of the model to be used as a radiotherapy treatment optimization tool could be the following: 1.The imaging, histopathological and molecular data of the patient are collected. 2.The patient's imaging tumor data is used in order to three dimensionally reconstruct the macroscopic structure of the tumor and its metabolic regions of interest (well oxygenated, poorly oxygenated and badly oxygenated regions). 3.The patient's histopathological and molecular tumor data are preprocessed in order to estimate the patient individualized alpha and beta parameters of the linear quadratic model. 4.A number of candidate dose fractionations are simulated for the patient specific virtual tumor. 5.The level of normal tissue complications for each candidate schedule is calculated based on standard radiobiological techniques. 6.The fractionations which do not lead to normal tissue complications beyond the corresponding limits are selected. 7. The fractionation that computationally leads to the maximum tumor shrinkage (and minimum tumor viability) is the theoretically optimal scheme. 8.Provided that prospective and retrospective clinical studies agree with the theoretical outcome for the subpopulation in which the patient falls and that the clinician does not have any substantiated objections, the theoretically suggested treatment schedule may be applied to the patient.

Table 1

Values of Cell Cycle Duration for GBM Cell Lines/Tumors

TC (hours)	Reference	Source
36-48	[11]	Literature review
26.2±5.6	[12]	U87 cell line, 4 surgical specimens
24-48	[13]	13 surgical specimens
22	[14]	U251 cell line
75.6±45.7	[15]	24 glioma surgical specimens (of which 13 GBMs)
48-72	[16]	-
57±6	[17]	4 surgical specimens

Table 2

Values of the α and β Parameters of the Linear Quadratic Model for GBM Cell Lines

α (Gy-1)	β (Gy-2)	Reference	Cell line
0.05	0.05	[10]	U87MG
0.1	0.06	[10]	SF-126
0.31	0.04	[10]	SF-268
0.36	0.06	[10]	U-251MG
0.45	0.06	[10]	SF-188
0.61	0.02	[10]	U87-LUX.8
0.17	0.02	[10]	U87-175.4
0.167±0.052	0.021±0.003	[20]	GB A7(puro)
0.370±0.077	0.015±0.007	[20]	GB A7R
0.25(0.22-0.29)	-	[21]	HGL21
0.32(0.28-0.36)	-	[21]	D54MG
0.387±0.121	-	[22]	U87MG
0.226±0.023	-	[22]	U373MG
0.482±0.055	-	[22]	M059J
0.038±0.004	0.007±0.002	[23]	U87MG
0.05±0.007	0.018±0.003	[23]	T98G
0.031	0.0030	[24]	U251
0.017	0.019	[24]	U87
0.058	0.0061	[24]	T98G

Table 3

Values of the Cell Loss Factor for GBM Tumors

CLF	Reference	Source
0.278	[25]	5 U87 xenografts in nude mice
0.170.49	[26]	U87, HGL9 andU251-MG xenografts in nude mice
0.1710.450.2720.110.160.555	[27]	HGL9, HGL21, HGL29, U251-MG, MMC1, U87 and Hp56 xenografts in nude mice
0.750.320.130.440.040.37	[28]	15 ethylnitrosourea-induced rat gliomas

Table 4

Values of the Growth Fraction for GBM Tumors

GF (%)	Reference	Source
28 (14-4)	[11]	From literature review
8-46	[13]	13 surgical specimens
Approximately 15	[15]	24 glioma surgical specimens (of which 13 GBMs)
10-30	[16]	-
14-39 (31±10)	[17]	4 surgical specimens
10-40	[33]	22 surgical specimens