Introduction of Hypermatrix and Operator Notation into a Discrete Mathematics Simulation Model of Malignant Tumour Response to Therapeutic Schemes In Vivo . Some Operator Properties

Abstract

The tremendous rate of accumulation of experimental and clinical knowledge pertaining to cancer dictates the development of a theoretical framework for the meaningful integration of such knowledge at all levels of biocomplexity. In this context our research group has developed and partly validated a number of spatiotemporal simulation models of in vivo tumour growth and in particular tumour response to several therapeutic schemes. Most of the modeling modules have been based on discrete mathematics and therefore have been formulated in terms of rather complex algorithms (e.g. in pseudocode and actual computer code). However, such lengthy algorithmic descriptions, although sufficient from the mathematical point of view, may render it difficult for an interested reader to readily identify the sequence of the very basic simulation operations that lie at the heart of the entire model. In order to both alleviate this problem and at the same time provide a bridge to symbolic mathematics, we propose the introduction of the notion of hypermatrix in conjunction with that of a discrete operator into the already developed models. Using a radiotherapy response simulation example we demonstrate how the entire model can be considered as the sequential application of a number of discrete operators to a hypermatrix corresponding to the dynamics of the anatomic area of interest. Subsequently, we investigate the operators' commutativity and outline the “summarize and jump” strategy aiming at efficiently and realistically address multilevel biological problems such as cancer. In order to clarify the actual effect of the composite discrete operator we present further simulation results which are in agreement with the outcome of the clinical study RTOG 83–02, thus strengthening the reliability of the model developed.

Keywords

cancer multiscale modeling computer models tumour growth tumour response to treatment oncosimulator radiobiology operator notation hypermatrix

1. Introduction

An inelastic prerequisite for an effective treatment of cancer is understanding and modeling the corresponding spatiotemporal natural phenomenon of tumour growth and response to therapeutic schemes concurrently on several biocomplexity levels. The usually fast growth and resilience of tumours suggest that they are emerging, opportunistic systems rather than random, disorganized and diffuse cell masses.^1,2 Therefore, the entire in vivo growing tumour rather than only a single cell³ must be investigated and treated as a self-organizing complex dynamic system. In this context there is need for advanced computational models to simulate the complexity of solid tumour growth, invasion and metastasis combining a range of disciplines including medical, biological, biophysical, engineering and statistical physics research.⁴

This section provides a brief outline of several of the concepts and earlier research efforts to model tumour behaviour. Duechting et al⁵ introduced a simulation model which concerns the in vitro case or the early avascular stages of small in vivo tumours and is based on a consideration of the distinct phases of the cell cycle. Kocher et al^6,7 presented a simulation model of the development of a tumour spheroid and its response to radiosurgery. The detailed imaging based geometry of the clinical tumour which might facilitate the direct clinical validation of the model had not been considered however. Instead an equivalent spherical tumour was considered in place of the generally arbitrarily shaped actual tumour. Additionally, detailed cell cycle phase biology (phases G₁, S, G₂, M) had not been taken into account, grouping of the cells into only proliferating and dormant classes being considered instead. It is noted that none of the above mentioned models has been applied to large clinical tumours (of varied geometrical shapes) and none of them simulates shrinkage for an arbitrarily shaped clinical tumour undergoing treatment. In the pure tumour growth models presented by Kansal et al^1,2 a discretising grid is used in which each geometrical cell is able to contain a large number of biological cells, although the grid has not been used to discretise clinical tumours of arbitrary shape. Swanson et al,^8,9 Mandonnet et al,¹⁰ have developed clinically oriented spatiotemporally models of tumour growth and invasion concerning glioblastoma multiforme (GBM). Although growth and invasion constitute fundamental phenomena related to GBM treatment optimization, the investigators have not focused neither on the radiobiological nor on the pharmacodynamic mechanisms that determine the cell survival probabilities and the subsequent shrinkage. Byrne et al¹¹ and Alarcon et al,^12–14 have developed mathematical models of avascular tumour growth and angiogenesis evolution pertinent mainly to the initial stages of tumour development. Valuable insight can be gained using such models, but extension to clinical voluminous tumours is not an a priori manageable task. Wise et al¹⁵ have developed a three dimensional multi-species non linear tumour growth model. A review of significant efforts to model cancer in silico (=on the computer) has also been presented by Deisboeck et al.¹⁶

An effort to overcome some of the above mentioned limitations has been previously made by our research group through the development of four-dimensional patient-specific in vivo simulation models of imageable tumour response to radiotherapeutic and chemotherapeutic schemes.^17–29 All parameters used in the models have already been defined and can be determined (in principle) experimentally or clinically. Therefore, use of new mathematically dictated parameters of ambiguous physical meaning has been avoided.

Looking now at the problem of cancer from a broader theoretical perspective we realize that the impressive rate of production of experimental and clinical knowledge pertaining to the disease dictates the development of a generic and desirably universal theoretical framework for the meaningful integration of such knowledge. The obvious reasons for knowledge integration are both deeper understanding of cancer and optimization of therapeutic schemes on the patient individualized context by performing in silico experiments. At the same time as cancer is an excellent paradigm of multilevel biological phenomena any rigorous theoretical treatment of cancer dynamics could provide important hints for the modeling and simulation of other biological phenomena including other homeostatic imbalances (diseases).²² As more and more complexity aspects of the natural phenomenon of cancer are incorporated into mathematical and computational models important discrete mathematics modeling treatments tend to become difficult to understand by the wider cancer modeling community. Therefore, the need for a symbolic mathematical notation has become obvious. Such an approach could be viewed as following the formalist school founded by D. Hilbert.³⁰ According to the formalist thesis mathematics is concerned with formal symbolic systems.

Stimulated by these remarks we present the introduction of an operator notation into the malignant tumour response to therapeutic schemes in vivo making use of several versions of the simulation model developed by our group (see citations above). In order to proceed to the introduction of a discrete mathematics operator notation the anatomic region of interest and its biological dynamics are represented by a hypermatrix $\bar{\underline{α}}$ . A hypermatrix can be viewed as a matrix of [a matrix of [… of [matrices (or vectors)]…]]. The hypermatrix $\bar{\underline{α}}$ is created by the superposition of a discretization mesh on the anatomic region of interest and the consideration of equivalence classes within each geometrical cell of the mesh representing the various phases within or out of the cell cycle that a biological cell can be found. Discrete time represents a further dimension of the hypermatrix.

2. Discretization of the Biological Problem

Collection of the appropriate imaging data (e.g. MRI T1 contrast enhanced, CT/PET, etc.), registration, interpolation and three dimensional reconstruction constitute the initial steps of an in vivo discrete mathematics treatment of the tumour growth and therapy response phenomenon. A discretizing mesh is superimposed on the anatomic region of interest and the contents of each geometrical cell of the mesh are distributed into equivalence classes corresponding to the various phases within or out of cell cycle. The mean time spent within each phase is another parameter that characterizes the subsets of each geometrical cell.^18,23,24

3. Operator and Hypermatrix Notation

The following mathematical entities are considered in the proposed treatment: $\bar{\underline{α}}$ stands for the hypermatrix corresponding to and dynamically describing the anatomical region of interest (including the tumour and possibly the surrounding normal tissue). Each vector element of the hypermatrix $\bar{\underline{α}}$ is considered to have the following form that corresponds to a discretization mesh geometrical cell:

\bar{\underline{α}} (x_{i}, y_{j}, z_{k}, p_{l}, t_{n}) = (g^{i j k l n}, N_{p}^{i j k l n}, t_{p}^{i j k l n}, h_{p}^{i j k l n}, {\tilde{h}}_{p}^{i j k l n})

(1)

\bar{\underline{a}} (t_{0}) = {\bar{\underline{a}}}_{0} initial state of the tumour (just before the start of the treatment course to be simulated)

(2)

where the following symbols have been introduced:

p: phase within or out of the cell cycle

g: oxygen and nutrient provision

N_p: number of biological cells in phase p

t_p: mean time spent in phase p (time is usually measured in h)

h_p: number of therapy hit cells residing in phase p

${\tilde{h}}_{p}$ number of non therapy hit cells residing in phase p

x_{i} \in [x_{\min}, x_{\max}]

(3)

y_{j} \in [y_{\min}, y_{\max}]

(4)

z_{k} \in [z_{\min}, z_{\max}]

(5)

t_{n} \in [0, t_{\max}]

(6)

p_{l} \in [G_{1}, S, G_{2}, M, G_{0}, A, N, D]

(7)

where

ξ_min, ξ_max denote the minimum and maximum value respectively of the generic variable ξ during the simulation

G₁ denotes the G₁ cell cycle phase,

S denotes the DNA synthesis phase,

G₂ denotes the G₂ cell cycle phase,

M denotes mitosis,

G₀ denotes the dormant G₀ phase,

A denotes the apoptotic phase,

N denotes the necrotic phase,

D denotes the remnants of dead cells,

g \in {s, \tilde{s}}

(8)

s stands for sufficient oxygen and nutrient provision (for tumour cell proliferation), $\tilde{s}$ stands for insufficient oxygen and nutrient provision (for tumour cell proliferation).

Obviously this binary character of the oxygen and nutrient provision is to be considered only a first simplifying approximation.

N_{p} \in N_{0}

(9)

N₀ is the set of non negative integers

t_{p} \in [0, t_{p \max}]

(11)

h_{p} \in [0, N_{p}]

(12)

{\tilde{h}}_{p} \in [0, N_{P}]

(13)

It should be noted that in general all physicobiologically different components of each hypermatrix element can be considered dependent on all five dimensions of the proposed abstract space of tumour dynamics. For example oxygen and nutrient provision can change dramatically in space and time within the tumour. The treatment outcome is generally dependent on the phase in which a cell resides when irradiated or treated with chemotherapy. Cell cycle phases have generally different durations and therefore the mean time spent within each phase equivalence class of a given geometrical cell is dependent on the cell phase. Even the oxygen and nutrient provision (to the biological cells belonging to the same phase within a geometrical cell) may be microscopically related to the phase under consideration. A relatively large number of tumour cells within the G₀ phase located around a given point may imply inadequate oxygen and nutrient provision, although in this particular case dormancy is normally the outcome rather than the cause of inadequate oxygen and nutrient provision

Figure 1 provides a schematic representation of the proposed five dimensional discrete abstract space of tumour dynamics. Three dimensions (those corresponding to the variables x_i, y_j, z_k) represent space, another one (corresponding to the variable t_n) time and the fifth one (corresponding to the variable p_l) represents the cell phase within or out of the cell cycle in which a biological cell or a set of cells within a geometrical cell of the discretization mesh is found at a given instant. The entire simulation can be viewed as the periodic application of a number of discrete algorithmic operators on the hypermatrix of the anatomic region of interest. The period of application is equal to the time separating two consecutive discretization mesh scans. This has been taken equal to 1h in all applications referred to in this particular paper.

Figure 1.

A schematic representation of equation 1 showing the location of an indicative hypermatrix element $\bar{\underline{a}} (x_{i}, y_{j}, z_{k}, p_{l}, t_{n})$ (x_i, y_j, z_k, p_l, t_n) and its physically inhomogeneous and multidimensional content $(g^{i j k l n}, N_{p}^{j k l n}, t_{p}^{j k l n}, h_{p}^{j k l n}, {\tilde{h}}_{p}^{j k l n})$ [see text for symbols; ∼ h_G2 ≡ $\tilde{h}$ _G ]. The proposed five dimensional discrete abstract space of tumour dynamics (corresponding to the localization of each hypermatrix element) is shown on the bottom right of the diagram. Three dimensions (corresponding to variables x_i, y_j, z_k) represent space, another one (corresponding to variable t_n) represents time and the fifth one (corresponding to variable p_l) represents the cell phase within or out of the cell cycle in which a biological cell or a set of biological cells within a geometrical cell of the discretization mesh is found at a given instant.

The various modules of algorithmic manipulations on the hypermatrix can be thought of as corresponding to discrete operators acting on the hypermatrix in analogy to the action of continuous operators on a wave function in quantum mechanics (Schiff³¹ pp. 148–186). The critical importance of considering abstract (vector) spaces and operators has been made clear by practically all fields of physics (Morse and Feshbach,³² Part 1, pp. 76–92).

In order to proceed to a symbolic formulation of the operator application we make use of the following symbols:

f stands for the composite discrete operator i.e. the operator formed by the synthesis of all partial operators sequentially acting on the hypermatrix. Therefore, the updated hypermatrix at the time point t_n+1 is given by

f (\bar{\underline{a}} (t_{n})) = \bar{\underline{a}} (t_{n + 1})

(14)

The composite operator can be written as

f = f^{U} f^{E} f^{C} f^{H} f^{O} f^{T}

(15)

where

f^{J}, J \in {U, E, C, H, O, T} stands for a "partial operator"

(16)

T stands for time update (i.e. just the increase of time by e.g. 1 h and not the updated state of the hypermatrix $\bar{\underline{α}}$ at any time t_n).

O stands for the oxygen and nutrient provision status

H stands for the effect of therapy (mainly cell survival)

C stands for the eventually perturbed cell cycling due to therapy

E stands for differential expansion or shrinkage

U stands for oxygen and nutrients field update

\begin{matrix} f^{U} (f^{H} (f^{C} (f^{H} (f^{O} (f^{T} \bar{\underline{a}} (t_{n})))))) \\ = \underline{\bar{a}} (t_{n + 1}), n = 0, 1, ….. \end{matrix}

(17)

or in a more compact writing:

\begin{matrix} f^{U} f^{E} f^{C} f^{H} f^{O} f^{T} \bar{\underline{a}} (t_{n}) \\ = \bar{\underline{a}} (t_{n + 1}), n = 0, 1, ….. \end{matrix}

(18)

where the application of the operators takes place from the right to the left.

It is noted that the term partial operator as used in this work essentially denotes the application of complex algorithmic manipulations. The entire model has been constructed with a number of algorithmic manipulations (partial operators) applied repeatedly in a given order. Therefore, the whole model (entire algorithm) can be considered as the application of a composite operator. Subsequently this composite operator can always be decomposed into partial operators since the former is nothing more than a conceptual clustering of the latter.

It is quite obvious that the above mentioned concepts and symbols cannot include all the information needed for the simulation to run. Their role is (at least at the present stage) rather to identify and decompose the major conceptual mathematical treatment steps than to represent any assumption details. The proposed approach is to be seen as a continually evolving and optimized process. Examples of such evolutionary stages could be the following. In order to address the non imageable components of highly invasive tumours such as glioblastoma multiforme at a large time scale further operators could be proposed so as to explicitly handle diffusion phenomena at the cellular/tissue level. Additional operators could handle the refined biomechanics of the tumour and adjacent normal tissues as well complex molecular networks which largely determine the response of a single tumour cell to treatment (see Section 5).

4. Non Commutativity of the Operators

A careful study of the behaviour of any two of the “partial operators” f^J, f^K where J, K ∈ {U, E, C, H, O, T} reveals that they are non commutative. It is noted that the terms commutative and permutable are used interchangeably when applied to operators (but not always when applied to subgroups).³³ For example applying the cell cycle clock (f^C) to the proliferating cells first and subsequently calculating the effect of a radiotherapy fraction (f^H) may lead to a substantially larger number of surviving cells than the opposite sequence. The reason for such a difference could be the successful completion of mitosis for far more tumor cells in the first case which would lead to a larger number of tumour cells (tumour burden) before treatment is applied. Furthermore, it must be stressed that in general each partial operator has to be applied concurrently on all elements-geometrical cells of the hypermatrix $\bar{\underline{α}}$ (possibly already transformed by other partial operators) as there are in general interdependences among the various geometrical cells [e.g. differential shifting of the overflowing biological cells that have emerged following mitoses],

5. Multilevel Biology Considerations: The “Summarize and Jump” Strategy

In order to achieve a quite realistic prediction of the response of a tumour to therapeutic interventions several levels of biocomplexity have to be addressed at the same time. The decomposition of the composite discrete operator f to its constituent “partial operators” f^J, J ∈ {U, E, C, H, O, T} provides a conceptual tool useful for the analysis and superposition of several critical mechanisms that may take place on different biocomplexity levels. For example the molecular profile of a given tumour (e.g. the expressions of a number of critical genes) can be used in order to perturb the population based average survival fraction following irradiation with dose D so that the genetically produced radiosensitivity or radioresistance of the particular tumour of a given patient is taken into account during the simulation. This “summary” of the molecular level phenomena is currently incorporated into the partial operator f^H which represents the microscopic effect of the treatment intervention on the tumour. One could interpret the procedure of perturbing the population based average values of the tumour biological parameters as a two step process. The first step refers to “summarizing” what is happening on one biocomplexity level (here the molecular level) and providing the amount of perturbation whereas the second one refers to the “jumping” to another level (here the cellular level). In this case the “summary” refers to the percentage by which the survival fraction will have to be perturbed whereas the “jumping” refers to the fact that the individualized survival fraction is related to the cellular level on which the main bulk of the simulation process takes place. Furthermore, it is worth noting that it is at the cellular level that complete success or failure of tumour treatment is defined since complete success of tumour treatment implies that not even a single proliferative or dormant tumour cell has been left alive.

6. Some Indicative Points of the Model and their Correspondence to Specific “Partial Operators”

The introduction of the suggested notation is demonstrated through a version of a simulation model of (T1 gadolinium enhanced) imageable glioblastoma response to radiotherapeutic schemes. As a thorough account of the assumptions made would be beyond the scope of the present paper we only provide an indication of the operator notation introduction by referring to selected modeling points. Further modeling details can be found i.a. in^{18,19,23,24,26,27}

In order to clinically validate the simulation model, a series of simulation executions corresponding to the various arms of the RTOG study 83–02 have been performed.³⁴ This was a randomized Phase I/II study of escalating doses for Hyperfractionated radiotherapy (HF, 1.2 Gy twice daily to doses of 64.8, 72, 76.8, or 81.6 Gy) and Accelerated Hyperfractionated radiotherapy (AHF, 1.6 Gy twice daily to doses of 48 or 54.4 Gy) with carmustine (BCNU) for adults with supratentorial glioblastoma multiforme (GBM) or anaplastic astrocytoma. The study has revealed that GBM patients who received the higher HF doses had survival superior to the patients in the AHF arms or lower HF doses.

The in silico experiments performed involve three hypothetical imageable GBM tumours, otherwise identical except for their radiosensitivity parameters. In particular, the cases considered were the following:

a.
A GBM tumour with intact wild-type (wt) p53 function, and accordingly adjusted LQ Model parameters.³⁵
$α_{p} = 0.61 {Gy}^{- 1}, β_{p} = 0.02 {Gy}^{- 2}$

b.
A GBM tumour with mutant (mt) p53 gene³⁵:
$α_{p} = 0.17 {Gy}^{- 1}, β_{p} = 0.02 {Gy}^{- 2}$

c.
A GBM tumour with intermediately adjusted radiosensitivity:
$α_{p} = 0.36 {Gy}^{- 1}, β_{p} = 0.02 {Gy}^{- 2}$

In all cases, we set (Kocher;⁷ Perez and Brady,⁴⁵ p. 99):
$α_{G0} = α_{p} / OER, β_{G0} = β_{p} / {OER}^{2}, OER=3$
(19)

and
$α_{S} = 0.6 α_{p} + 0.4 α_{G0}, β_{S} = 0.6 β_{p} + 0.4 β_{G0} .$
(20)

The meaning of the symbols used is the following:
α_P, β_p: the LQ Model parameters for all proliferative cell cycle phases except for the DNA synthesis phase (S phase).

α_S, β _S : the LQ Model parameters for the S phase.

α_G0, β_G0. the LQ Model parameters for the resting G₀ phase.

The delivery of irradiation takes place at 08:00 and 16:00 every day, 5 days per week (no irradiation during weekends). The distribution of the absorbed dose in the tumour region is assumed to be uniform. It should also be noted that carmustin, which was administered to all patients enrolled in the RTOG—83–02 study, is assumed not to significantly modify the relative effectiveness of the radiation therapy schedules considered, as the chemotherapy administration schedule was the same for all patients.
6.1. Construction of the Hypermatrix

The first process that takes place before the application of the operators is the discretization of the anatomic region of interest and the construction of the corresponding hypermatrix. The imaging data (e.g. T1 gadolinium enhanced MRI slices, PET slices etc.) including the definition of the tumour contour, its metabolically active sub-regions and the anatomical structures of interest, the histopathological (e.g. type of tumour) and the genetic data (e.g. p53 status and other molecular data) of the patient are collected. The clinician delineates the tumour and the anatomical structures of interest by using a dedicated computer tool. In the case of radiotherapy, the planned distribution of the absorbed dose (e.g. in Gy) in the region of interest is also acquired. For the purpose of the 3D reconstruction and visualization of both the initial tumour and the simulation outcome, the 3D visualization package AVS/Express 4.2 has been used, (for details concerning the use of AVS/Express in the simulation model refer to Stamatakos et al).¹⁸ The description of the biological activity of the tumour is implemented by introducing the notion of the “geometrical cell (GC)”, the elementary cubic volume of the 3D discretizing mesh covering the region of interest as previously mentioned.

We assume that each GC of the mesh initially and normally accommodates a Number of Biological Cells (NBC). However, the maximum number of biological cells that is allowed to be accommodated within a GC is assumed to be NBC + [a fraction of NBC]. NBC apparently depends on the selected size of the GC and determines the quantization error of the model. The fraction of NBC considered in the code executions in this paper was 1/10 as this was shown to be the optimal one from both the convergence and CPU time demands points of view.²⁴

Typical clonogenic cell densities are 10⁴ to 10⁵ cells/mm^3,4 Since most GBM tumours are poorly differentiated and rapidly growing, we assume a clonogenic cell density of 2 × 10⁵ cells/mm³ in the proliferating cell region, 10⁵ cells/mm³ in the G₀ cell region and 0.2 × 10⁵ cells/mm³ in the dead cell region of the tumour.

It is noted that each multidimensional element of the proposed hypermatrix has a clear physical and/or biological meaning (e.g. number of biological cells within a given phase, number of therapy hit cells residing in a given phase etc.). Therefore, the hypermatrix (such as the exemplary one outlined above) is used to describe the distribution of several critical physical and biological quantities over space and time.

6.2. Operator Based Presentation of the Simulation Model Basics

In the following the various processes constituting the entire simulation algorithm are briefly and separately described with reference to the corresponding operators.

6.2.1. Operator f^T

Time is discretized and incremented. One hour has been adopted as the unit of time since 1h is the approximate duration of mitosis, the shortest cell cycle phase. According to the prescribed radiotherapy scheme the specific instants corresponding to the delivery of a radiation dose to the tumour region are defined. In each time step the geometrical mesh is scanned and the updated state of a given GC is determined on the basis of a number of behaviour algorithms:

6.2.2. Operator f^O

During each scan of the discretization mesh the effect of the oxygen and nutrient provision on the cells of each geometrical cell is taken into account. This provision, determining the metabolic potential of a region, is based on the imaging data and determines the distribution of the tumour cells in the proliferative, dormant and necrotic states within the regions without taking into account the eventual therapeutic interventions effects. Furthermore, distribution of the cells over the cell cycle phases (G₁, S, G₂, Mitosis) is considered based on experimental evidence (Katzung (Ed),³⁷).

6.2.3. Operator f^H

(I) At the time instants corresponding to the delivery of a specific radiation dose to the tumour (according to the prescribed radiotherapy scheme and the acquired distribution of the absorbed dose in the region of interest) the number of cells killed in a particular GC is calculated based on the Linear Quadratic (LQ) Model, which is widely used in the pertinent literature.^5,7,38,39 The fraction of cells surviving from a radiation dose D is given by

{S(D)=exp[-(αD+βD}^{2})]

(21)

where α (Gy⁻¹) and β (Gy⁻²) characterize the initial slope and the curvature, respectively, of the survival curve.

In an untreated tumour simulation case, the dose D would be set to zero.

(II) Lethally damaged cells following exposure to radiation undergo two mitotic divisions prior to death and disappearance from the tumour.³⁹

Note: Any eventual molecular perturbators of the cell surviving fraction are to be incorporated into this operator.

6.2.4. Operator f^C

At each time step the time registers of all GCs are incremented by one hour. Cell loss due to apoptosis and necrosis is computed. According to the cytokinetic model appearing in²³ possible transitions of the cells within a GC include: G₁ →S (if time spent in the G₁ cell cycle phase ≥ TG₁; TG₁ = duration of the G₁ phase), S→G₂ (if time spent in the S cell cycle phase ≥ TS; TS = duration of the S phase), G₂→M (if time spent in the G₂ cell cycle phase ≥ TG₂; TG₂ = duration of the G₂ phase), M→G₁ or M→G₀ (if time spent in the M cell cycle phase ≥ TM; TM = duration of mitosis), G₀→G₁ (if adequate oxygen and nutrient provision has been re-established) or G₀→N (if time spent in the G₀ cell cycle phase ≥ TG₀; TG₀ = maximum duration of the G₀ phase before a cell enters necrosis). The previously mentioned durations of the cell cycle phases TX, X ∈ {G₁, S, G₂, M, G₀} seem to follow the normal distribution according to pertinent literature. As a first approximation, we use the mean values of the duration of each cell cycle phase and neglect standard deviations. However pseudo-random number generators are used in order to de-synchronize the equivalence classes throughout the tumour.

The cell cycle duration T_C has been taken equal to 40 h. This is the average of the cell cycle durations we have found in the literature for GBM cell lines.^40,41 In Katzung³⁷ the approximate percentage of the cell cycle time spent in each phase by a typical malignant cell is assumed as follows: TG₁ = 0.4T_C, TS = 0.39T_C, TG₂ = 0.19T_C, TM = 0.02T_C. The duration of the G₀ phase is taken to be TG₀ = 25 h.⁴²

The cell loss factor (CLF) is considered equal to 0.3.⁴³ In⁴⁴ the authors note that cell loss is mainly due to necrosis (CLFN) and apoptosis (CLFA) and that gliomas have a low CLF in general. We assume that the total CLF (0.3) is the sum of the CLFN (0.27) and CLFA (0.03). We hypothesize low levels of apoptotic cells for GBM, as we have found that this is in general the case for gliomas.^35,36,44

6.2.5. Operator f^E

The differential tumour expansion and shrinkage algorithms are based on the use of random number generators in conjunction with adequately formed morphological rules. These rules lead to tumour shrinkage or expansion conformal to the initial shape of the tumour (if the mechanical properties of the surrounding normal tissue are considered uniform around the tumour and the tumour is not in contact with practically undeformable tissues such as bone). Two versions of the expansion and shrinkage algorithms have been tested. First version (a): For each GC, one out of the six possible directions of shrinkage or expansion is randomly chosen (Cartesian coordinate system XYZ centered at the current GC. Each axis defines two possible directions of movement). Second version (b1) Shrinkage: The outermost tumour GC is detected along each one of the six possible directions of shrinkage (Cartesian coordinate system XYZ centered at the current GC. Each axis defines two possible directions of movement). Its “6-Neighbour” GCs belonging to the Tumour (NGCT) are counted. The direction corresponding to the maximum NGCT is finally selected out of the six possible directions as the direction along which the shifting of the GCs will take place (shifting direction). In case that more than one shifting directions have the same maximum NGCT, then the selection is based on the use of a random number generator. Second version (b2) Expansion: A similar, though inverse, morphological—mechanical rule can be applied in the case of tumour expansion. The need for the formulation of the second improved version of the tumour shrinkage and expansion algorithm has arisen from the inspection of the macroscopic results of the simulation algorithms. Specifically, the completely random selection of one out of the six possible shifting directions, according to the first version, results in a premature extensive fragmentation of the tumour region in case of radiotherapy, which is usually incompatible with clinical experience. The general trend is a conformal shrinking of most solid tumours (Perez and Brady,⁴⁵ Figs. 1—4, p. 10). Using the second version of the algorithms this problem is solved. The mechanical properties of the surrounding normal tissue are considered uniform around the tumour, with the exception of an absolute lack of deformability of the bone. As a first approximation immunological reactions, invasion and formation of metastases have been ignored.

Figure 2.

The proposed ten levels of biocomplexity.

Figure 3.

number of surviving tumour cells as a function of time for the glioblastoma tumour with mutant p53 gene (see section 6). The radiotherapeutic schemes correspond to schemes considered by the RTOG 83–02 clinical study.

Figure 4.

Number of surviving tumour cells as a function of time for the glioblastoma tumour with wild-type p53 gene (see section 6). The radiotherapeutic schemes correspond to schemes considered by the RTOG 83–02 clinical study. The curves on the right correspond to the HF schedules whereas the curves on the left correspond to the AHF schedules.

6.2.6. Operator f^U

After having completed a scan of the discretizing mesh the oxygen and nutrient field is updated based on the criterion determining the relative position of the proliferative, dormant and necrotic regions of the tumour. The reason for this process is to take into account any eventual expansion or shrinkage of the tumour that would lead to a perturbation of the previous metabolic potential field.

7. Results

In order to clarify the actual effect of the composite discrete operator f introduced above the following indicative simulation predictions are included in this paper. Figures 3, 4 and 5 present the results of the in silico experiments in the form of the number of surviving tumour cells (proliferating and dormant) as a function of time, for the tumours with mutant p53, wild-type p53 and the tumour with intermediate radiosensitivity, correspondingly. Improved tumour control following high-dose HF irradiation is evident in the diagrams and is in agreement with the conclusions of the clinical study RTOG 83–02. In fact, the higher the total dose in an HF schedule, the better the result in terms of tumour cell kill. It should be noted that the clinical study compared the treatments based on patient survival whereas the simulation results presented refer to tumour responsiveness. Although a correspondence of patient survival with tumour responsiveness may not always hold on a single patient basis, when taking into account relatively large clinical trial populations such a correspondence becomes much more reasonable.

Figure 5.

Number of surviving tumour cells as a function of time for the glioblastoma tumour with intermediately adjusted radiosensitivity parameters (see section 6). The radiotherapeutic schemes correspond to schemes considered by the RTOG 83–02 clinical study.

More specifically, the inspection of the simulation results reveals that AHF schedules, which employ a higher fraction dose compared to HF schedules, seem at first to be beneficial as they achieve the maximum tumour cell kill at some instant. Nevertheless, the duration of the AHF schedules is smaller; as a result, if they fail in eradicating “all” tumour cells, tumour repopulation begins earlier (as in the cases of the tumour with mt p53 and the tumour with intermediate radiosensitivity). Only in the case of the tumour with wt p53, do all radiotherapy schemes kill all the clonogenic cells we have initially assumed (although this may not be the case in reality due to considerable instabilities of the simulation when it comes to the last few living tumour cells (chaotic behaviour limits)), so the tumour does not regrow after the end of the treatment. Of course, regions of potential microscopic disease have not been considered, and the accuracy of the simulation model decreases as the number of tumour cells is reduced.

8. Discussion

The introduction of operator notation into the process of modeling malignant tumour response to therapeutic schemes has led to a brief and comprehensive description of the major steps of the simulation process. In this way highly complex algorithmic treatments are decomposed into simpler procedures which are readily identifiable by the wider research community. The use of mathematical symbols to denote complex algorithmic processes is expected to function as a stimulant for the advancement of multilevel biological modeling through symbolic mathematical expressions. Such expressions could by themselves provide hints for further questions, investigations and optimizations due to their inherent logical and quantitative associations. Especially the analogy of the operator notation in biomedicine with their use in several fields of physics such as classical and quantum mechanics can act as a source of guidance and quantitative insight into the hypercomplex biological phenomena.

Obviously the treatment presented should be viewed only as an initial step of a rather long term modeling process as more and more experimental and clinical knowledge could be incorporated into the models which are under continuous development and optimization. Specific aspects that are currently addressed in parallel include i.a. the mitotic potential categories of the cancer stem, progenitor and differentiated cells (leading to a new dimension in the discrete abstract space),^28,29 adjacent normal tissue response,⁴⁶ molecular networks (adapting the cell survival probability to the patient individualized context), chemotherapy optimization^20,21,28,29 etc. A possible application of the approach presented in this paper is on the development of the in silico oncology action of the European Commission (EC) funded project “ACGT: Advancing Clinicogenomic Trials on Cancer” (FP6–2005–IST-026996). In particular the development of the software simulation tool named “Oncosimulator”⁴⁷ can be described based on the notation proposed in this paper. An analogous simulation tool is being developed within the framework of the EC funded project “ContraCancrum: Clinically Oriented Translational Cancer Multilevel Modelling” (FP7-ICT–2007–2–223979).

9. Conclusions

The introduction of the proposed hypermatrix and operator notation in order to denote, decompose and identify complex biological mechanisms that contribute to the hypercomplex phenomenon of malignant tumour growth and response to therapeutic schemes has been presented through the use of a discrete mathematics simulation model concerning radiotherapy response. Several aspects of the model had been developed by our research group in the past. Symbolic operator notation has provided a compact way of describing the most crucial simulation steps thus offering a possible basis for quantitative multilevel biology based on discrete mathematics. Furthermore, simulation results mimicking branches of the RTOG 83–02 clinical study have provided both clarification of the actual content of the composite discrete operator and additional evidence for the potential of the simulation approach presented.

Disclosures

The authors report no conflicts of interest.

Footnotes

Acknowledgements

This work was supported in part by the European Commission under the projects “ACGT: Advancing Clinicogenomic Trials on Cancer” (FP6–2005–IST-026996) and ContraCancrum: Clinically Oriented Translational Cancer Multilevel Modelling” (FP7-ICT–2007–2–223979). The authors acknowledge N. Graf, University Hospital of Saarland and N. Uzunoglu, National Technical University of Athens for stimulating discussions as well as the anonymous reviewers for their highly constructive comments and suggestions.

References

Kansal

A.R.

, Torquato

, Harsh

G.R.

, Chiocca

E.A.

, Deisboeck

T.S.

Simulated brain tumour growth dynamics using a three-dimensional cellular automaton.

J Theor Biol. 2000a; 203: 367–82.

Kansal

A.R.

, Torquato

, Harsh

I.V.

, Chiocca

E.A.

, Deisboeck

T.S.

Cellular automaton of idealized brain tumour growth dynamics.

BioSystems. 2000b; 55: 119–27.

Kraus

, Wolf

Emergence of self-organization in tumour cells: relevance for diagnosis and therapy.

Tumour Biol. 1993; 14: 338–53.

Jones

, Dale

R.G.

Mathematical models of tumour and normal tissue response.

Acta Oncol. 1999; 38: 883–93.

Duechting

, Ginsberg

, Ulmer

Modelling of radiogenic responses induced by fractionated irradiation in malignant and normal tissue.

Stem Cells. 1995; 13(S1): 301–6.

Kocher

, Treuer

Reoxygenation of hypoxic cells by tumour shrinkage during irradiation. A computer simulation.

Strahlentherapie und Onkologie. 1995; 171: 219–30.

Kocher

, Treuer

, Voges

, Hoevels

, Sturm

, Mueller

R.P.

Computer simulation of cytotoxic and vascular effects of radiosurgery in solid and necrotic brain metastases, Radiother Oncol. 2000; 54: 149–56.

Swanson

K.R.

, Alvord

E.C.

Jr , Murray

J.D.

Virtual brain tumours (gliomas) enhance the reality of medical imaging and highlight inadequacies of current therapy.

Br J Cancer. 2002; 86: 14–8.

Swanson

K.R.

, Bridge

, Murray

J.D.

, Alvord

E.C.

Virtual and real brain tumors: using mathematical modelling to quantify glioma growth and invasion.

J Neurol Sci. 2003; 216: 1–10.

10.

Mandonnet

, Delattre

J.Y.

, Tanguy

M.L.

, Swanson

K.R.

, Carpentier

A.F.

, Duffau

Continuous growth of mean tumor diameter in a subset of grade II gliomas.

Ann Neurol. 2003; 53: 524–8.

11.

Byrne

, Matthews

Assymetric growth of avascular solid tumours: exploiting symmetries.

IMA J Mathematics Applied in Medicine and Biology. 2002; 19: 1–29.

12.

Alarcon

, Byrne

H.M.

, Maini

P.K.

A cellular automaton model for tumour growth in inhomogeneous environment.

J Theor Biol. 2003; 225: 257–74.

13.

Alarcon

, Byrne

H.M.

, Maini

P.K.

Towards whole-organ modelling of tumour growth.

Prog Biophys Mol Biol. 2004a; 85: 451–72.

14.

Alarcon

, Byrne

H.M.

, Maini

P.K.

A mathematical model of the effects of hypoxia on the cell-cycle of normal and cancer cells.

J Theor Biol. 2004b; 229: 395–411.

15.

Wise

S.M.

, Lowengrub

J.S.

, Frieboes

H.B.

, Cristini

Three-dimensional multispecies nonlinear tumor growth—I Model and numerical method.

J Theor Biol. 2008; 253: 524–43.

16.

Deisboeck

T.S.

, Zhang

, Yoon

, Costa

In silico cancer modeling: is it ready for prime time?

Nat Clin Pract Oncol. 2009; 6: 34–42.

17.

Stamatakos

, Dionysiou

, Nikita

, Zamboglou

, Baltas

, Pissakas

In vivo tumour growth and response to radiation therapy: a novel algorithmic description.

Int J Radiat Oncol Biol Phys. 2001; 51 Suppl 1: 240.

18.

Stamatakos

G.S.

, Dionysiou

D.D.

, Zacharaki

E.I.

, Mouravliansky

N.A.

, Nikita

, Uzunoglu

In silico radiation oncology: combining novel simulation algorithms with current visualization techniques.

Proc. IEEE 2002; 90(11): 1764–77.

19.

Stamatakos

G.S.

, Antipas

V.P.

, Uzunoglu

N.K.

, Dale

R.G.

A four dimensional computer simulation model of the in vivo response to radiotherapy of glioblastoma multiforme: studies on the effect of clonogenic cell density.

Brit J Radiology. 2006a; 79: 389–400.

20.

Stamatakos

G.S.

, Antipas

V.P.

, Uzunoglu

N.K.

Simulating chemotherapeutic schemes in the individualized treatment context: The paradigm of glioblastoma multiforme treated by temozolomide in vivo.

Comput Biol Med. 2006b; 36: 1216–34.

21.

Stamatakos

G.S.

, Antipas

V.P.

, Uzunoglu

N.K.

A spatiotemporal, patient individualized simulation model of solid tumor response to chemotherapy in vivo: the paradigm of glioblastoma multiforme treated by temozolomide.

IEEE Trans Biomed Engineering. 2006c; 53: (8): 1467–77.

22.

Stamatakos

G.S.

Towards a collaborative formulation of the Mathematical Principles of Natural Philosophy: Living Matter. The paradigm of In Silico Oncology. DIMACS Workshop on Computational Tumor Modeling, DIMACS Center, Rutgers University, Piscataway, NJ, USA. 2006d. http://dimacs.rutgers.edu/Workshops/TumorModeling/abstracts.html.

23.

Dionysiou

D.D.

, Stamatakos

G.S.

, Uzunoglu

N.K.

, Nikita

K.S.

, Marioli

A four-dimensional simulation model of tumour response to radiotherapy in vivo: parametric validation considering radiosensitivity, genetic profile and fractionation.

J Theor Biol. 2004; 230: 1–20.

24.

Dionysiou

D.D.

, Stamatakos

G.S.

, Uzunoglu

N.K.

, Nikita

K.S.

A computer simulation of in vivo tumour growth and response to radiotherapy: New algorithms and parametric results.

Computers in Biology and Medicine. 2006; 36: 448–64

25.

Dionysiou

D.D.

, Stamatakos

G.S.

, Gintides

, Uzunoglu

, Kyriaki

Critical Parameters Determining Standard Radiotherapy Treatment Outcome for Glioblastoma Multiforme: A Computer Simulation.

The Open Biomedical Engineering Journal. 2008; 2: 43–51.

26.

Dionysiou

D.D.

, Stamatakos

G.S.

Applying a 4D multiscale in vivo tumor growth model to the exploration of radiotherapy scheduling: the effects of weekend treatment gaps and p53 gene status on the response of fast growing solid tumors.

Cancer Informatics. 2006; 2: 113–21.

27.

Stamatakos

G.S.

, Uzunoglu

Computer Simulation of Tumour Response to Therapy. Cancer Bioinformatics: from therapy design to treatment. Edited by Sylvia Nagl. John Wiley & Sons, Ltd., Chichester, UK. 2006.

28.

Kolokotroni

E.A.

, Stamatakos

G.S.

, Dionysiou

D.D.

, Georgiadi

E Ch

, Desmedt

, Graf

N.M.

Translating Multiscale Cancer Models into Clinical Trials: Simulating Breast Cancer Tumor Dynamics within the Framework of the “Trial of Principle” Clinical Trial and the ACGT Project. Proc. 8th IEEE International Conference on Bioinformatics and Bioengineering (BIBE 2008), Athens, Greece, 2008 Oct 8–10. IEEE Catalog Number: CFP08266, ISBN: 978-1-4244-2845-8, Library of Congress: 2008907441, Paper No. BE-2.1.1, length: 8 pages (in electronic format). 2008.

29.

Georgiadi

E Ch

, Stamatakos

G.S.

, Graf

N.M.

Multilevel Cancer Modeling in the Clinical Environment: Simulating the Behavior of Wilms Tumor in the Context of the SIOP 2001/GPOH Clinical Trial and the ACGT Project. Proc. 8th IEEE International Conference on Bioinformatics and Bioengineering (BIBE 2008), Athens, Greece, 2008 Oct 8–10. IEEE Catalog Number: CFP08266, ISBN: 978-1-4244-2845-8, Library of Congress: 2008907441, Paper No. BE-2.1.2, length: 8 pages (in electronic format). 2008.

30.

Eves

Foundations and Fundamental Concepts of Mathematics. Third Edition. Dover, Mineola, US. 1990.

31.

Schiff

L.I.

Quantum Mechanics. Third Edition. McGraw Hill Kogakusha, Tokyo. 1981. pp. 148–86.

32.

Morse

P.M.

, Feshbach

Methods of Theoretical Physics. Mc Graw-Hill, New York. 1953.

33.

Miller

G.A.

Groups which are the products of two permutable proper subgroups.

PNAS. 1935; 21: 469–72.

34.

Werner-Wasik

, Scott

C.B.

, Nelson

D.F.

Final report of a phase I/II trial of hyperfractionated and accelerated hyperfractionated radiation therapy with carmustine for adults with supratentorial malignant gliomas.

Cancer. 1996; 77: 1535–43.

35.

Haas-Kogan

D.A.

, Yount

, Haas

; p53-dependent G₁ arrest and p53 independent apoptosis influence the radiobiologic response of glioblastoma. Int J Radiat Oncol Biol Phys. 1996; 36: 95–103.

36.

Tribius

, Pidel

, Casper

ATM protein expression correlates with radioresistance in primary glioblastoma cells in culture.

Int J Radiat Oncol Biol Phys. 2001; 50: 511–23.

37.

Katzung

B.G.

, (Ed.). Basic and Clinical Pharmacology, 8th ed, New York: McGraw-Hill/Lange Medical Books. 2001.

38.

Fowler

J.F.

The linear-quadratic formula and progress in fractionated radiotherapy.

Br J Radiol. 1989; 62: 679–94.

39.

Steel

(Ed.). Basic Clinical Radiobiology. Arnold, London, UK. 2002.

40.

Dillehay

L.E.

A model of cell killing by low-dose-rate radiation including repair of sublethal damage, G₂ block, and cell division.

Rad Res. 1990; 124: 201–7.

41.

Hegedues

, Czirok

, Fazekas

, Babel

, Madarasz

, Viscsek

Locomotion and proliferation of glioblastoma cells in vitro: statistical evaluation of videomicroscopic observations.

J Neurosurgery. 2000; 92: 428–34.

42.

Duechting

, Ulmer

, Lehrig

, Ginsberg

, Dedeleit

Computer simulation and modelling of tumour spheroid growth and their relevance for optimization of fractionated radiotherapy.

Strahlenther Onkol. 1992; 168: 354–60.

43.

Huang

, Allam

, Perez

, Taghian

, Freeman

, Suit

The effect of combining recombinant human tumor necrosis factor-alpha with local radiation on tumor control probability of a human glioblastoma multiforme xenograft in nude mice.

Int J Radiat Oncol Biol Phys. 1995; 32(1): 93–8.

44.

Nakajima

, Nakasu

, Morikawa

, Inubushi

Estimation of volume doubling time and cell loss in an experimental rat glioma model in vivo.

Acta Neurochir. 1998; 140: 607–13.

45.

Perez

, Brady

(Eds.). Principles and Practice of Radiation Oncology. Philadelphia: Lippincott-Raven. 1998.

46.

Antipas

V.P.

, Stamatakos

G.S.

, Uzunoglu

N.K.

A patient-specific in vivo tumor and normal tissue model for prediction of the response to radiotherapy: a computer simulation approach.

Methods Inf Med. 2007; 46: 367–75.

47.

Stamatakos

G.S.

, Dionysiou

D.D.

, Graf

N.M.

The Oncosimulator: a multilevel, clinically oriented simulation system of tumor growth and organism response to therapeutic schemes. Towards the clinical evaluation of in silico oncology. Proceedings of the 29th Annual International Conference of the IEEE EMBS Cite Internationale, August 23–26, SuB07.1: 6628–6631, Lyon, France. 2007.

Introduction of Hypermatrix and Operator Notation into a Discrete Mathematics Simulation Model of Malignant Tumour Response to Therapeutic Schemes In Vivo . Some Operator Properties

Abstract

Keywords

1. Introduction

2. Discretization of the Biological Problem

3. Operator and Hypermatrix Notation

5. Multilevel Biology Considerations: The “Summarize and Jump” Strategy

6. Some Indicative Points of the Model and their Correspondence to Specific “Partial Operators”

6.2. Operator Based Presentation of the Simulation Model Basics

6.2.1. Operator fT

6.2.2. Operator fO

6.2.3. Operator fH

6.2.5. Operator fE

7. Results

9. Conclusions

Disclosures

Footnotes

Acknowledgements

References

6.2.1. Operator f^T

6.2.2. Operator f^O

6.2.3. Operator f^H

6.2.5. Operator f^E