Introduction of hypermatrix and operator notation into a discrete mathematics simulation model of malignant tumour response to therapeutic schemes in vivo. Some operator properties.

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.


Introduction
An inelastic prerequisite for an effective treatment of cancer is understanding and modeling the corres ponding 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 3 must be investigated and treated as a selforganizing 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, biophy sical, engineering and statistical physics research. 4 This section provides a brief outline of several of the concepts and earlier research efforts to model tumour behaviour. Duechting et al 5 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 simula tion 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 equiva lent spherical tumour was considered in place of the generally arbitrarily shaped actual tumour. Addition ally, detailed cell cycle phase biology (phases G 1 , S, G 2 , 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, 10 have devel oped clinically oriented spatiotemporally models of tumour growth and invasion concerning glio blastoma 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 11 and Alarcon et al, [12][13][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 exten sion to clinical voluminous tumours is not an a priori manageable task. Wise et al 15 have developed a three dimensional multispecies 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. 16 An effort to overcome some of the above men tioned 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 radiothera peutic and chemotherapeutic schemes. [17][18][19][20][21][22][23][24][25][26][27][28][29] All param eters 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 integra tion of such knowledge. The obvious reasons for knowledge integration are both deeper understand ing 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). 22 Αs more and more complexity aspects of the natural phenomenon of cancer are incorporated into mathematical and computational models impor tant 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. 30 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 malig nant 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 repre sented by a hypermatrix α . A hypermatrix can be viewed as a matrix of [a matrix of [… of [matrices (or vectors) ]…]]. The hypermatrix α 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 rep resenting 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.

Discretization of the Biological problem
Collection of the appropriate imaging data (e.g. MRI T1 contrast enhanced, CT/PET, etc.), registration, inter polation and three dimensional reconstruction consti tute 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 charac terizes the subsets of each geometrical cell. 18,23,24

Operator and Hypermatrix notation
The following mathematical entities are considered in the proposed treatment: a stands for the hyper matrix corresponding to and dynamically describ ing the anatomical region of interest (including the tumour and possibly the surrounding normal tissue). Each vector element of the hypermatrix α is consid ered to have the following form that corresponds to a discretization mesh geometrical cell: a t a ( ) 0 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 mea sured in h) h p : number of therapy hit cells residing in phase p  h p number of non therapy hit cells residing in phase p where ξ min ,ξ max denote the minimum and maximum value respectively of the generic variable ξ during the simulation G 1 denotes the G 1 cell cycle phase, S denotes the DNA synthesis phase, G 2 denotes the G 2 cell cycle phase, M denotes mitosis, G 0 denotes the dormant G 0 phase, A denotes the apoptotic phase, N denotes the necrotic phase, D denotes the remnants of dead cells, s stands for suffcient oxygen and nutrient provision (for tumour cell proliferation),  s stands for insuff cient 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 0 is the set of non negative integers It should be noted that in general all physicobio logically 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 depen dent on the phase in which a cell resides when irradi ated 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 0 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 corre sponding 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 ) repre sents the cell phase within or out of the cell cycle in which a biological cell or a set of cells within a geo metrical 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 algo rithmic 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.
The various modules of algorithmic manipulations on the hypermatrix can be thought of as correspond ing to discrete operators acting on the hypermatrix in analogy to the action of continuous operators on a wave 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. 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 opera tors sequentially acting on the hypermatrix. There fore, the updated hypermatrix at the time point t n+1 is given by (14) The composite operator can be written as T stands for time update (i.e. just the increase of time by e.g. 1 h and not the updated state of the hyperma trix α 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 or in a more compact writing: 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 informa tion needed for the simulation to run. Their role is (at least at the present stage) rather to identify and decompose the major conceptual mathematical treat ment 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 pro posed 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).

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). 33 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 substan tially larger number of surviving cells than the oppo site 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 elementsgeometrical cells of the hypermatrix α (possibly already transformed by other partial operators) as there are in general interdepen dences among the various geometrical cells [e.g. dif ferential shifting of the overflowing biological cells that have emerged following mitoses].

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 opera tors" f J , J ∈{U, E, C, H, O, T} provides a concep tual 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 expres sions of a number of critical genes) can be used in order to perturb the population based average sur vival 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 treat ment intervention on the tumour. One could inter pret the procedure of perturbing the population based average values of the tumour biological parameters as a two step process. The first step refers to "sum marizing" 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.

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 8302 have been per formed. 34 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 radio therapy (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 param eters. In particular, the cases considered were the following: In all cases, we set (Kocher; 7 Perez and Brady, 45 p. 99): and The meaning of the symbols used is the following: α P , β p : the LQ Model parameters for all prolifera tive 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 0 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 adminis tration schedule was the same for all patients.

Construction of the hypermatrix
The first process that takes place before the appli cation 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 subregions and the anato mical 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 dedi cated 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 pur pose 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). 18 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 con vergence and CPU time demands points of view. 24 Typical clonogenic cell densities are 10 4 to 10 5 cells/mm 3 . 4 Since most GBM tumours are poorly differentiated and rapidly growing, we assume a clonogenic cell density of 2 × 10 5 cells/mm 3 in the proliferating cell region, 10 5 cells/mm 3 in the G 0 cell region and 0.2 × 10 5 cells/mm 3 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 biologi cal 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.

Operator based presentation of the simulation model basics
In the following the various processes constituting the entire simulation algorithm are briefly and sepa rately described with reference to the corresponding operators.

Operator f T
Time is discretized and incremented. One hour has been adopted as the unit of time since 1h is the approxi mate 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:

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 prolifera tive, 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 1 , S, G 2 , Mitosis) is considered based on experimental evidence (Katzung (Ed), 37 ).

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 distri bution 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 where α (Gy -1 ) and β (Gy -2 ) 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. 39 Note: Any eventual molecular perturbators of the cell surviving fraction are to be incorporated into this operator. The previously mentioned durations of the cell cycle phases TX, X ∈{G 1 , S, G 2 , M, G 0 } 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 stan dard deviations. However pseudorandom number generators are used in order to desynchronize the equivalence classes throughout the tumour.

Operator f C
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 37 the approximate percentage of the cell cycle time spent in each phase by a typical malignant cell is assumed as follows: TG 1 = 0.4T C , TS = 0.39T C , TG 2 = 0.19T C , TM = 0.02T C . The duration of the G 0 phase is taken to be TG 0 =25 h. 42 The cell loss factor (CLF) is considered equal to 0.3. 43 In 44 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 mor phological rules. These rules lead to tumour shrinkage or expansion conformal to the initial shape of the tumour (if the mechanical properties of the surround ing normal tissue are considered uniform around the tumour and the tumour is not in contact with practi cally 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 pos sible 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 ''6Neighbour'' 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 shift ing directions, according to the first version, results in a premature extensive fragmentation of the tumour region in case of radiotherapy, which is usually incom patible with clinical experience. The general trend is a conformal shrinking of most solid tumours (Perez and Brady, 45 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 consid ered 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.
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 previ ous metabolic potential field.

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, wildtype p53 and the tumour with intermediate radiosensitivity, correspondingly. Improved tumour control following highdose HF irradiation is evident in the diagrams and is in agreement with the conclu sions of the clinical study RTOG 8302. 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.
More specifically, the inspection of the simulation results reveals that AHF schedules, which employ a higher fraction dose compared to HF sched ules, 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 radio therapy 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 treat ment. 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.

Discussion
The introduction of operator notation into the pro cess of modeling malignant tumour response to therapeutic schemes has led to a brief and compre hensive description of the major steps of the simula tion 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 guid ance 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, 46 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" (FP62005IST026996). In particular the development of the software simulation tool named "Oncosimulator" 47 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" (FP7ICT20072223979).

conclusions
The introduction of the proposed hypermatrix and operator notation in order to denote, decompose and identify complex biological mechanisms that con tribute to the hypercomplex phenomenon of malig nant 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 com pact 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 8302 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.
publish with Libertas Academica and every scientist working in your field can read your article "I would like to say that this is the most author-friendly editing process I have experienced in over 150 publications. Thank you most sincerely." "The communication between your staff and me has been terrific. Whenever progress is made with the manuscript, I receive notice. Quite honestly, I've never had such complete communication with a journal." "LA is different, and hopefully represents a kind of scientific publication machinery that removes the hurdles from free flow of scientific thought." Your paper will be: • Available to your entire community free of charge • Fairly and quickly peer reviewed • Yours! You retain copyright http://www.la-press.com