E Ch

C Georgiadi

Desmedt

G S Stamatakos

E A Kolokotroni

D D Dionysiou

Georgiadi

C Desmedt

An advanced discrete state -discrete event multiscale simulation model of the response of a solid tumor to chemotherapy: Mimicking a clinical study

Keywords: Breast cancer, simulation, cancer modeling, chemotherapy, in silico oncology

published or not. The documents may come     An advanced discrete state -discrete event multiscale simulation model of the response of a solid tumor to chemotherapy: Mimicking a clinical study

Oncology Group, National Technical University of Athens, is characterized by several crucial new features, such as the explicit distinction of proliferating cells into stem cells of infinite mitotic potential and cells of limited proliferative capacity, an advanced generic cytokinetic model and an improved tumor constitution initialization technique.

A sensitivity analysis regarding critical parameters of the model has revealed their effect on the behavior of the biological system. The favorable outcome of an initial step towards the clinical adaptation and validation of the simulation model, based on the use of anonymized data from the TOP clinical trial, is presented and discussed. Two real clinical cases from the TOP trial with variable molecular profile have been simulated. A realistic time course of the tumor diameter and a reduction in tumor size in agreement with the clinical data has been achieved for both cases by selection of reasonable model

Introduction

Over the last years it has become clear that in order to understand both the disease and the complex natural phenomenon of cancer in a quantitative way sophisticated multiscale models are needed. Furthermore, multiscale cancer models are expected to substantially support treatment optimization in the patient individualized context. To this end several research groups have been working on this field on the global level. The major modeling approaches can be distinguished into the following three categories: continuous mathematics based models, discrete mathematics based models and hybrid models.

Continuous mathematics based models rely primarily on the diffusion equation, which is applied in order to describe the diffusion of molecules such as oxygen and glucose and/or the invasion of tumor cells into the surrounding tissue(s) in the case of diffusive tumors (e.g. glioblastoma) [START_REF] Murray | Mathematical Biology II. Spatial Models and Biomedical Applications[END_REF][START_REF] Swanson | Virtual brain tumours (gliomas) enhance the reality of medical imaging and highlight inadequacies of current therapy[END_REF][START_REF] Breward | A multiphase model describing vascular tumour growth[END_REF][START_REF] Clatz | Realistic simulation of the 3-D growth of brain tumors in MR images coupling diffusion with biomechanical deformation[END_REF][START_REF] Frieboes | An integrated computational/experimental model of tumor invasion[END_REF][START_REF] Guiot | The dynamic evolution of the power exponent in a universal growth model of tumors[END_REF][START_REF] Enderling | A mathematical model of breast cancer development, local treatment and recurrence[END_REF]. In many cases continuous equations are solved by finite difference methods such as the Crank-Nicholson method. Discrete mathematics based models rely primarily on the consideration of several discrete states in which cells may be found (e.g. cell cycle phases, stem cell state etc.) and the transitions between states (e.g. cell necrosis following a prolonged residence of a tumor cell in the G0 phase etc.) [START_REF] Duechting | Three-dimensional pattern generation applied to spheroidal tumor growth in a nutrient medium[END_REF][START_REF] Ginsberg | Computer simulation of fractionated radiotherapy: further results and their relevance to percutaneous irradiation and brachytherapy[END_REF][START_REF] Stamatakos | In silico radiation oncology: combining novel simulation algorithms with current visualization techniques[END_REF]Stamatakos et al., , 2006aStamatakos et al., , 2006b[START_REF] Stamatakos | 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[END_REF][START_REF] Dionysiou | A four-dimensional simulation model of tumour response to radiotherapy in vivo: parametric validation considering radiosensitivity, genetic profile and fractionation[END_REF], 2006, Dionysiou and Stamatakos, 2006;[START_REF] Deisboeck | In silico cancer modeling: is it ready for prime time?[END_REF]. State transitions can be decided by making use of several possible "decision calculators" such as cytokinetic diagrams, agent-based techniques etc.

Hybrid methods exploit the potential of both continuous and discrete methods in order to address in more detail several tumor dynamics mechanisms. A clinically oriented example of the latter is the combination of pharmacokinetics differential equations with discrete state/event simulations in order to simulate the response of large imageable tumors to chemotherapeutic treatment (Stamatakos et al., 2006b). A further example refers to the modeling of cancer cell invasion to tissue [START_REF] Ramis-Conde | Mathematical modelling of cancer cell invasion of tissue[END_REF]. As more and more biological mechanisms are being quantitatively understood most new models tend to fall in this category i.e. to make use of a combination of continuous and discrete mathematics. In that way the particular advantages of both approaches are exploited and therefore mechanisms of both predominantly continuous and predominantly discrete mechanisms can be addressed.

Hybrid models with a strong discrete character have been shown to be particularly adequate for the simulation of the response of large imageable clinical tumors to therapeutic interventions such as radiotherapy and chemotherapy [START_REF] Stamatakos | In silico radiation oncology: combining novel simulation algorithms with current visualization techniques[END_REF](Stamatakos et al., , 2006a(Stamatakos et al., , 2006b[START_REF] Stamatakos | 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[END_REF][START_REF] Dionysiou | A four-dimensional simulation model of tumour response to radiotherapy in vivo: parametric validation considering radiosensitivity, genetic profile and fractionation[END_REF]Dionysiou et al., , 2006)). Since clinical validation is a sine qua non necessity for clinically oriented multiscale cancer models, this paper presents an actual clinical trial driven model concerning the response of early breast cancer to epirubicin. A branch of the Trial of Principle (TOP: Topoisomerase II Alpha Gene Amplification and Protein Overexpression Predicting Efficacy of Epirubicin) trial (TOP trial) has been addressed as a paradigm of a chemotherapy optimization targeted clinicogenomic trial.

Carcinoma of the breast is the most commonly diagnosed cancer in women. Even though a variety of effective treatments exists, their benefits and adverse effects vary considerably. Early results indicate that gene expression profiles may contribute to predicting response to breast cancer treatment [START_REF] Potti | Genomic signatures to guide the use of chemotherapeutics[END_REF]. A number of clinical trials are running with the aim to identify gene expressions that correlate with chemotherapy results. The TOP trial, which is also part of the ACGT European project (ACGT), addresses a particular molecular subgroup of breast cancer patients and aims to determine the predictive rate of a single class of chemotherapeutic agent. More specifically, this study aims to evaluate the topoisomerase II alpha gene amplification and protein overexpression as markers predicting the efficacy of epirubicin monotherapy in the primary treatment of breast cancer patients. Inclusion criteria are non-metastatic, early breast cancer patients with estrogen receptor (ER) -negative tumors with a maximum dimension larger than 2 cm as defined by ultrasound. The patients are treated with single-agent epirubicin as neoadjuvant treatment, followed by surgery and adjuvant chemotherapy [START_REF] Durbecq | Topoisomerase-IIA expression as a predictive marker in a population of advanced breast cancer patients randomly treated either with single-agent doxorubicin or single-agent docetaxel[END_REF].

The model presented, although having its roots in previous work done and published by the In Silico Oncology Group, National Technical University of Athens (In Silico Oncology Group) is characterized by several crucial new features such as an advanced generic cytokinetic model for the response of a tumor cell to chemotherapeutic treatment, the explicit distinction of proliferating cells into stem cells and cells of limited mitotic potential and the adaptation of the relative populations of the various cell categories/phases to the cell category/phase transition rates (a novel tumor initialization strategy aiming at avoiding abnormal patterns of free tumor growth as described in subsequent paragraphs). Available imaging-based, histological, molecular and treatment data are exploited by the model in order to strengthen patient individualization modeling and in the future patient individualized treatment optimization. The expected use of the model following thorough clinical adaptation, optimization and validation is to simulate either several candidate treatment schemes for a particular patient and support the selection of the optimal one for him or her or to simulate the expected extent of tumor shrinkage for a given time instant and decide on the adequacy or not of the simulated scheme. The integrated simulation system has been termed "Oncosimulator" [START_REF] Stamatakos | 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[END_REF]. Two "oncosimulators" are currently being developed, The structure of the paper is as follows. First the cytokinetic model proposed and adopted is presented. Subsequently the basic algorithms of the spatiotemporal discretization of the tumor are provided, followed by the tumor initialization strategy.

Considerations on the pharmacokinetics and pharmacodynamics of epirubicin are presented. A description of the model parameters is provided. A sensitivity analysis of the model regarding the effect of critical parameters on the dynamics of the biological system is presented, followed by the description of an initial step towards the clinical adaptation and validation of the model. Simulations for two real case studies from the TOP trial are presented and discussed. The paper concludes with a discussion of the main points and the conclusions of the work.

The discrete statediscrete event cytokinetic model

The simulation model presented is based on the well documented assumption that tumor sustenance is due to the existence of stem cells i.e. cells with theoretically unlimited mitotic potential. Most cancers comprise a heterogeneous population of cells with differences in their proliferative potential. Cancer stem cells are generally thought to represent a small portion of the total tumor cell population [START_REF] Stingl | Molecular heterogeneity of breast carcinomas and the cancer stem cell hypothesis[END_REF].

They possess the ability of self renewal i.e. the ability to give birth to more stem cells. Fig. 1 depicts the generic cytokinetic model proposed for the case of tumor growth and response to chemotherapy. A cytokinetic model limited to free tumor growth has been presented in a previous publication [START_REF] Kolokotroni | 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[END_REF].

Proliferating cells, including stem and LIMP cells, go through the four phases of the cell cycle: gap 1 (G1) phase, DNA synthesis (S) phase, gap 2 (G2) phase and mitosis (M) phase. Depending on the nutrient and oxygen supply level of the local environment the daughter cells may re-enter the cell cycle at G1 phase or enter the dormant (resting) G0 phase following mitosis. Cells that are dormant due to insufficient nutrient supply and oxygenation can survive for a limited period of time. Subsequently, they die through necrosis unless in the meantime the local metabolic conditions regarding nutrition and oxygenation have become adequate. In the latter case dormant cells re-enter the G1 phase. Proliferating, dormant and differentiated cells may die due to aging and spontaneous apoptosis. Differentiated cells may also die through necrosis.

When a tumor is chemotherapeutically treated, a fraction of cancer cells are lethally hit by the drug or its metabolites. Lethally hit cycling tumor cells enter a rudimentary cell cycle that leads to apoptotic death via a specific phase dictated by the action mechanism of the chemotherapeutic agent used. Similarly, in the case of cell cycle non specific drugs, lethally hit dormant (G0) cells enter the G0hit phase. Marking of a cell as hit by the drug is assumed to take place at the instant of drug administration although its actual time of death is dictated by the specific pharmacokinetics and pharmacodynamics of the drug. For the special case of epirubicin treatment, the fraction of cells marked as hit by the chemotherapeutic agent is considered the same for all cell cycle phases and the G0 phase since epirubicin is a cell cycle non specific drug as well as a cell cycle phase non specific drug [START_REF] Fda | Food and Drug Administration[END_REF]. It is pointed out however that cell cycle phase specific drugs can be readily modeled by the cytokinetic model shown in Fig. 1 by appropriately selecting the "Chemo" induced exit from the normal cell cycle for both cases of stem and LIMP cells.

Spatiotemporal discretization of the tumor

Although macroscopically inhomogeneous tumors of generic shapes can be simulated by the full version of the model developed, in this paper the special case of a macroscopically homogeneous tumor of generic shape is presented. The application examples, however, concern the even simpler case of a macroscopically homogeneous spherical tumor since this has been the best approximation to the imaging data available from the TOP trial (TOP trial). Such a spatially simple tumor model is a reasonable first approximation to the modeling of a wide range of breast cancer tumors. The latter is supported by both the rather homogeneous normal soft tissue biomechanics of the breast and the accumulated clinical experience regarding breast tumors. In fact the tumor size as reported in the Case Report Forms (CRFs) of the TOP trial is defined as the maximum diameter based on ultrasound examination. It should also be noted that macroscopic spatial homogeneity implies that the model parameter values used refer actually to their spatial average throughout the tumor.

In contrast with its spatial simplicity the model is characterized by a rather high complexity regarding the number of mitotic potential cell categories and cycling phases considered, as well as the corresponding transitions. A three dimensional (3-D) cubic discretizing mesh is superimposed upon the anatomical region of interest. The elementary cube of the mesh is called geometrical cell (GC) and in the cases considered in this paper its size is 1 mm 3 . Each GC occupied by the tumor is assumed to initially contain a Number of Biological Cells (NBC) (see also section 6(i)). The biological cells residing within each geometrical cell of the mesh are distributed into the five categories mentioned above i.e. the stem, LIMP, differentiated, apoptotic and necrotic categories.

From the mathematical standpoint each cell category defines an equivalence class.

Distribution of the cells into the five equivalence classes creates one level of biological cell population partitioning within each GC. At each given instant each stem or LIMP cell can be either proliferating or dormant. Proliferation or dormancy creates another level of cell population partitioning. Cell cycle phases (G1, S, G2, M) introduce a finer partitioning of proliferating cells (stem and LIMP) into subclasses. A further partitioner in the case of therapeutic intervention is treatment hitting i.e. a boolean variable denoting whether a biological cell has been hit by treatment. The relative population (expressed as the fraction of the total tumor cell population) of each equivalence class and its equivalence subclasses is initialized based on cell category and cell phase transition rates as described in section 4 and Appendix A.

The initial distribution of the proliferating cells throughout the cell cycle phases is assumed analogous to the corresponding cell cycle phase durations (see also section 6(ii)). In order to tackle computing power limitations a number of carefully chosen quantizations are introduced. Time is discretized. The time unit in the practical cases considered in this paper is taken one hour since this is approximately the duration of mitosis, the shortest cell cycle phase [START_REF] Bast | Cancer Medicine 5th[END_REF]. For any given instant the biological cells belonging to the same cell category and cell cycle phase within a given GC are assumed synchronized. However, biological cells belonging to different GCs or to different categories and cell cycle phases within the same GC are not assumed synchronized. From the computational standpoint a sufficient number of registers are used to describe the state of each GC occupied by the tumor. They include i.a. the number of biological cells residing in each equivalence class and subclass and the time spent at each subclass. Mean time spent by the biological cells of a given equivalence subclass in the same subclass is initialized using a random number generator (Monte Carlo technique). The time under consideration can vary between 0 and the maximum time of the corresponding phase.

At each time step i.e. every hour the discretizing mesh covering the anatomical region of interest is virtually scanned in order to apply the basic metabolic, cytokinetic, pharmacokinetic/pharmacodynamic and mechanical rules that govern the spatiotemporal evolution of the tumor system. For practical reasons each complete virtual scan can be viewed as consisting of two mesh scans. A brief outline of the scanning procedure is given below.

(i) First Scan

The first scan aims at updating the state of each GC according to the proposed and adopted cytokinetic model of In any one of the cases of dormant (including stem and LIMP), differentiated, necrotic and apoptotic cells a fraction of the corresponding subclass(es) population may be transferred to another subclass or disappear from the tumor at each time step according to the cytokinetic model (Fig. 1). Therefore the following transitions may take place. I and sections 6(iii),6(iv)).

As far as chemotherapy is concerned, at any time instant corresponding to drug administration, the numbers of proliferating and dormant cells belonging to each one of the phases G1, S, G2, M or G0 and to each one of the stem or LIMP mitotic potential categories that are designated as hit by the drug are computed. The latter is achieved through the utilization of the cell kill ratio (CKR) parameter that corresponds to the drug and dose (per m 2 of the patient surface) considered. In terms of the simulation model's parameters the cell kill ratio is the percentage of LIMP and stem cells hit by the chemotherapeutic agent after each drug administration. The above mentioned cell numbers are added to the corresponding cell numbers of the drug affected equivalence subclasses, designated as "phase"hit e.g. G1hit in the cytokinetic diagram of Fig. 1. A general strategy for the calculation of the cell kill ratio for the case of epirubicin is described in section 5.

(ii) Second Scan

The second scan aims to preserve a roughly uniform cell density throughout the tumor volume. To this end, adequately shaped morphological rules are introduced, which may lead to tumor expansion, as is the case in free tumor growth, or no change in tumor volume or tumor shrinkage as is usually the case after treatment administration.

More specifically, at any given time point the total cell population that can be accommodated in each GC is allowed to fluctuate between a minimum (0.9*NBC) and a maximum (1.1*NBC) value. If the total population exceeds the maximum value of 1.1*NBC then a procedure is initiated that attempts to unload the total GC population minus NBC to neighboring GCs (26 GC neighborhood is considered) possessing empty space i.e. GCs with total cell population less than NBC. The procedure starts from the neighboring GC possessing the maximum free space. If two or more neighboring GCs possess the same free space then a random number generator is used so as to select the visiting order of the GCs. The procedure is repeated until all the excess cells have been transferred if possible. If the procedure fails to reduce the total population of the GC under consideration below the upper limit (maximum value) then an adjacent GC is freed from its contents which are moved outwards. The latter push the contents of a chain of geometrical cells outwards too. The excess contents of the GC under consideration are placed into the newly freed adjacent GC. The previous process leads to differential tumor expansion. The position of the GC to be freed from its contents relative to the GC with the excess contents is determined using a random number generator. The shifting of the chain of GCs mentioned above can take place along any randomly selected direction.

On the other hand if the GC's total cell population is below the minimum value then a similar procedure attempts to unload all cells to neighboring GCs possessing free space.

If the GC becomes empty then a chain of GC contents is shifted towards the GC under consideration so as to fill the vacuum generated. The latter leads to differential tumor shrinkage. Shifting of the GC content chain takes place in a way analogous to the previously described one.

The above procedure, however, may give rise to the following "side effects". chemotherapy is rather scarce in the relevant literature. It is possible to expect that in some cases some degree of fragmentation will result, although pressure from adjacent normal tissue is one possible counter-acting mechanism. The available postchemotherapy data from the TOP trial are in the form of tumors' maximum dimension, which led us to the assumption that at least the particular breast cancer tumors remained fairly cohesive.

The "second scan" algorithms introduced into the simulation model ensure uniform tumor cell density throughout the tumor volume (in the absence of any special information) and conformal tumor shrinkage, as depicted e.g. in (Perez and Brady, 1998). Our primary concern is to avoid artificial tumor fragmentation, namely an extensive premature fragmentation caused by the very procedure of shifting chains of adjacent geometrical cells along randomly selected directions. There is still the possibility to have geometrical cells with less or more than the typical cell content (currently used margin=0.1*typical number of cells in a geometrical cell) and this flexibility in cell number actually "incorporates" the phenomenon of tumor "fragmentation", in the sense that when chemotherapy is simulated geometrical cells become emptier over time. The above mentioned margin has been adjusted in order to achieve a rather uniform cell density throughout the tumor volume, but may be for example increased to permit a more extensive fragmentation in case that relevant information is available. These algorithms relate to geometrical aspects of the simulation and may be easily omitted from simulation runs in cases of macroscopically homogeneous tumors (such as those addressed in this paper) as they have no effect on (biological) cell numbers. The latter is the most critical aspect of the simulation both in case of untreated tumor growth and in tumor shrinkage induced by chemotherapy treatment.

A simplified flowchart of the entire simulation procedure pertaining to a macroscopically homogeneous solid tumor of arbitrary shape is provided in Fig. 2.

The core output of the simulation code is a number of matrices corresponding to the discretization mesh of the anatomic region of interest for a series of time points. These matrices contain all the necessary biological information to be used for the 3D metabolic, kinetic etc. reconstruction of the tumor at several time points or the creation of several diagrams describing tumor dynamics and response to therapy as a function of space and/or time.

It should be noted that the full model (not presented in this paper) supports the division of tumor area into different metabolic regions (e.g. necrotic and proliferative) and the handling of each region separately. In this case, GCs are initially characterized as necrotic or proliferative based on pertinent imaging data available and a process of image segmentation, interpolation and 3-D reconstruction. Different values of specific model parameters can be assigned to each region. A special procedure has been devised in order to conformally expand or shrink any internal necrotic region(s) during tumor evolution.

Tumor constitution initialization

If tumor GCs are initialized with arbitrary values of tumor cell populations belonging to the various mitotic potential categories and phases (stem, LIMP, differentiated, dead, proliferating, dormant etc.) whereas at the same time specific category / phase transition rates are used, it is very likely to observe an abnormal free tumor growth behavior. A decrease in tumor volume followed by a volume increase is a very common pattern [START_REF] Kolokotroni | 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[END_REF]. In order to avoid such artificial behavior the concept of a nomogram linking cell category / phase relative populations with cell category /phase transition rates was initially introduced [START_REF] Kolokotroni | 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[END_REF]. However, due to its shortcomings the nomogram method has been replaced by a new, more efficient and flexible technique, described in appendix A.

Condition for monotonic free growth

As far as untreated tumor growth simulations are concerned, depending on the values assigned to the code input parameters, a tumor that grows over time (which is normally the expected case) or a tumor that gradually diminishes (and therefore is usually biologically irrelevant) may result. The outcome strongly depends on the properties and the resulting behavior of the stem proliferating cells. In order to create a growing tumor the number of stem cells must increase over time and the following inequality must hold:

1 ) P 1 )( 1 ( ) 1 ( sleep t sym T A P R C (1)
where T c is the cell cycle duration, R A the apoptosis rate of cancer stem cells, P sym the fraction of stem cells that perform symmetric division and P sleep the fraction of newborn cells that enter G0 phase (Table I). The rationale behind the above inequality, which holds true for the majority of cases, can be found in a previous publication [START_REF] Kolokotroni | 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[END_REF]. Further analysis concerning free tumor growth preconditions is underway.

Pharmacokinetics and pharmacodynamics of the chemotherapeutic agent

Epirubicin is an anthracycline chemotherapeutic agent, derivative of doxorubicin. It exerts its cytotoxic action through various mechanisms; the most established one is intercalation between bases of double stranded DNA, thereby inhibiting nucleic acid (DNA -RNA) and protein synthesis. It also interferes with DNA replication and transcription, blocks helicase activity, inhibits DNA cleavage by topoisomerase II, and forms toxic oxygen-free radicals causing DNA and cell damage (FDA 1999). In specific cases epirubicin is favored over other anthracycline drugs (doxorubicin), as it appears to cause fewer side-effects due to its less toxic nature at equivalent therapeutic doses.

Epirubicin is considered a cell cycle non-specific as well as a cell cycle phase non specific drug [START_REF] Fda | Food and Drug Administration[END_REF].

In the simulation model tumor cells are assumed to absorb the drug (to be treatment hit) at all cycling phases as well as at G0, whereas apoptotic death of hit cells is assumed to occur at the end of S phase (i.e. between S and G2 phases). The latter is supported by the observation that the maximal cytotoxic effect of epirubicin has been observed on S and G2 phases (FDA, 1999, p.40).

Following intravenous administration, epirubicin is rapidly and widely distributed into tissues. Its higher concentration is observed in liver, spleen, kidney and small intestine.

Epirubicin undergoes extensive hepatic metabolism [START_REF] Danesi | Pharmacokinetics and pharmacodynamics of combination chemotherapy with paclitaxel and epirubicin in breast cancer patients[END_REF] and is also metabolized by other organs and cells, including red blood cells. It is eliminated mainly through biliary excretion and, to a lesser extent, by urinary excretion [START_REF] Fda | Food and Drug Administration[END_REF].

Epirubicin pharmacokinetics has been described in various studies by an open threecompartment model with elimination from the central compartment [START_REF] Danesi | Pharmacokinetics and pharmacodynamics of combination chemotherapy with paclitaxel and epirubicin in breast cancer patients[END_REF][START_REF] Fda | Food and Drug Administration[END_REF]. If a specific value of CKR was to be used in a simulation run, then based on the parameter values of such a model and making use of e.g. the SAAM pharmacokinetics software (SAAM II, 2009) various pharmacokinetic quantities of interest could be estimated such as the Area Under Curve (AUC) for a given drug dose.

Subsequently, a reasonably good initial approximation would be to compute the survival fraction for epirubicin treated tumor cells on the basis of experimental FDA data concerning the pharmacodynamics of epirubicin, and more specifically the in vitro cytotoxicity of epirubicin on HeLa cells [START_REF] Fda | Food and Drug Administration[END_REF]. The survival fraction for any realistic dose for which no experimental data are available could be approximately calculated through linear interpolation. A lower CKR might be expected for the in vivo case due to imperfect vascularization, unless the particular genetic profile of a specific tumor somehow increases responsiveness to epirubicin treatment. It is noted that in cases of human breast cancer steep concentration gradients have been shown for the highly relevant drug doxorubicin in pertinent studies [START_REF] Lankelma | Doxorubicin gradients in human breast cancer[END_REF].

In the simulations presented in the following sections of the paper, the suggested value of the CKR for each case study (the "apparent" CKR) is the CKR that produces good agreement between the evolution of the simulated tumor and that of the real tumor according to the clinical data. The individual patient's genetic profile data, if available, could be used to perturb the population based mean cell survival fraction. In this way an appropriate molecular signature of the patient will be exploited and therefore further individualization of the treatment plan will be achieved.

Description of the model parameters

Table I presents the simulation model's input parameters and their range of values according to pertinent literature or based on logical assumptions supported by basic science or clinical experience in case of lack of literature data. In the following a summarizing description of the model's parameters and their adopted values is provided for reasons of clarity:

(i) NBC As has been described previously, this parameter refers to the number of biological cells initially occupying each GC belonging to the tumor. In order to preserve a roughly constant mean cell density throughout the tumor volume, the NBC parameter is allowed to fluctuate between a minimum and a maximum value. Therefore, as described in previous paragraphs, accordingly shaped morphological rules governing the special evolution of the tumor have been introduced. A typical value of tumor cell density found in literature is 10 6 biological cells/mm 3 (Steel, 1997, p.9), which corresponds to 10 6 biological cells per occupied GC since the considered GCs' dimensions for the simulations that have been performed so far is 1mm x 1mm x1mm.

(ii) Cell cycle phases duration

According to literature the duration of the various cell cycle phases follows a normal distribution. Although this observation has been explicitly taken into account in an in vitro tumor growth and treatment response model by our group [START_REF] Zacharaki | Simulating growth dynamics and radiation response of avascular tumour spheroid model validation in the case of an EMT6/Ro multicellular spheroid[END_REF], in the present model cell cycle phase durations are considered constant and equal to their literature based mean values in order to accelerate executions. It is pointed out however that algorithmically constructing a normal distribution of phase durations is a pretty simple task.

Based on literature, breast cancer cell cycle duration may vary from approximately 20h to 96h [START_REF] Sterin | Levels of phospholipid metabolites in breast cancer cells treated with antimitotic drugs[END_REF][START_REF] Cos | Modulation of the lenght of cell cycle time of MCF-7 human breast cancer cells by melatonin[END_REF][START_REF] Descamps | Nerve Growth Factor Is Mitogenic for Cancerous but Not Normal Human Breast Epithelial Cells[END_REF][START_REF] Barnes | Expression of inducible Hsp70 enhances the proliferation of MCF-7 breast cancer cells and protects against the cytotoxic effects of hyperthermia[END_REF][START_REF] Meyer | Practical breast carcinoma cell kinetics: review and update[END_REF]. The duration of mitosis is considered constant and equal to 1h [START_REF] Bast | Cancer Medicine 5th[END_REF]. The rest of the cell cycle phases durations are computed based on [START_REF] Salmon | Cancer Chemotherapy[END_REF], after having taken into consideration the above assumption regarding the constant duration of mitosis. More specifically the following equations are used:

T G1 = T S = 0.41(T c -T M ), T G2 = 0.18(T c -T M ), T M = 1h.

(iii) Duration of G0, duration of apoptosis and necrosis

According to literature dormant cells resting in G0 phase can survive under hypoxic conditions for 4-10 days [START_REF] Maseide | Mathematical modeling of chronical hypoxia in tumors considering potential doubling time and hypoxic cell lifetime[END_REF]. Tumor apoptotic cells are generally considered to be rapidly phagocytosed in vivo [START_REF] Dewey | Radiation-induced apoptosis: relevance to radiotherapy[END_REF], contrary to the time-consuming process of necrosis products removal.

(iv) Cell category/phase transition rates/fractions

The values of these parameters have been selected based on both qualitative or semiquantitative information and the dictates of the accumulated basic science and clinical experience. Systematic use of TOP clinical trial data is expected to permit a quantitative refinement of the initial assumptions.

(v) Number of mitotic divisions that LIMP (committed progenitor) cells undergo before they become terminally differentiated (also called LIMP mitotic stages)

By varying the number of LIMP mitotic stages (assumed range 1 to 10) tumors characterized by different differentiation degrees can be simulated. To the best of our knowledge, there exists no specific information regarding possible values for this parameter in the pertinent literature. Nevertheless, the choice of its value is governed by the need to achieve specific relative percentages of stem and LIMP cells (which is also a reflection of a tumor's degree of differentiation). A larger number in this interval leads to lower percentages of stem cells.

Sensitivity analysis

The presented model is characterized by a considerable number of input parameters as presented in the previous section and summarized in Table I. A number of parametric studies have been performed in order to study the model's behavior in relation to the value of each input parameter for tumor free growth and chemotherapeutic response simulations and thereby detect the most critical parameters of the model. An extensive sensitivity analysis has been performed and reported in the ACGT project deliverable D8.3 'Report on the Refinement and Optimization of the Algorithms and Codes, and the Initial Clinical Validation and Adaptation of the "Oncosimulator". The following section presents indicative results. A detailed sensitivity analysis involving the role of several model parameters on tumor free growth and response to therapy will be the topic of a separate paper (to appear). A characteristic approach that will be presented in this separate paper is to adopt a ±10% variation around a reference value of each model parameter and the inspection of the variation in the output. Such analyses permit for example a sorting of the model's parameters according to the magnitude of their effect on the output.

Table II presents the values assigned to the model's input parameters for the indicative sensitivity analyses presented in this paper.

Tumor free growth parametric studies

Doubling time

Different types of tumors with respect to the degree of aggressiveness can be simulated by assigning suitable values to the model's parameters. In the following, an indicative number of exploratory simulations are presented. These have been performed in order to study the effect of various input parameters on tumor growth rate. The model's parameters considered in these simulation studies are the fraction of cells entering the G0 phase following mitosis (sleep fraction, P sleep ), the fraction of stem cells that perform symmetric division (symmetric division fraction, P sym ) and the necrosis rate of differentiated tumor cells (R NDiff, ) (see Table I).

For the analyses presented below, homogeneous tumors have been considered, consisting of a single proliferating region with no distinguishable necrotic or other areas. This assumption has been dictated by the availability of only volumetric and not necrosis distribution data in the TOP trial. This implies that microscopic necrotic regions which are primarily due to imperfect neovascularization of the tumor are postulated, rather homogeneously distributed over the entire tumor.

The cell category transition rates are assumed to be constant. Such an approximation is considered applicable for a relatively short time interval compared to the tumor lifetime (for example during the chemotherapeutic treatment course), and the constant values are assumed to reflect the means of the actual cell category transition rates over the interval.

The use of the means of several tumor dynamics parameters over a substantial time interval is quite customary in radiobiology (Begg, 1997, pp.14-16). Such a tumor would be characterized by a grossly exponential growth pattern which in fact approximates a segment of the Gompertzian curve (Steel, 1997, p.10;[START_REF] Zoubek | Natural course of a Wilms' tumour[END_REF][START_REF] Graf | What are the expectations of a Clinician from In Silico Oncology[END_REF]. Therefore for a short time interval the population of the various cell categories evolves over time according to the equation:

N(t)=N o e at

(2)

where N o is the population at time 0 and a is the growth rate constant.

Only for theoretical purposes in certain simulation cases is the considered time interval deliberately extended beyond what would seem a short segment of the Gompertzian curve.

The doubling time, T d , defined as the period of time required for the total tumor cell population to double, is a widely used parameter characteristic of the growth rate of a tumor. Under the previous restrictions the following equation holds:

T d = (ln2)/a (3)
The doubling time is the same for all cell category populations and can be computed based on the following relationship:

T d = (ln2)*(t 2 -t 1 )/ln(N 2 /N 1 ) (4)
where N 1 and N 2 are the populations at times t 1 and t 2 respectively.

The effect of the parameters considered on tumor doubling time is revealed in Fig. 3.

The following remarks can be made:

(a) Fraction of cells entering the G0 phase after mitosis (Sleep fraction, P sleep ): For low values of the sleep fraction, the simulated tumor is characterized by a small doubling time of a few days implying a rather aggressive tumor. As P sleep increases, the doubling time increases as well, initially slowly for a wide range of values, and then rapidly towards infinity. There exists an upper limit of the sleep fraction that corresponds to the equality in ( 1). This limit depends on the cell cycle duration and corresponds to a tumor that remains constant over time. As the cell cycle duration increases, the permitted value range of the sleep fraction narrows, since the upper limit decreases. Values of P sleep larger than the upper limit result in a tumor that diminishes over time (Fig. 3a).

(b) Fraction of stem cells that perform symmetric division (Symmetric division fraction, P sym ): In contrast to the sleep fraction, in the case of the symmetric division fraction there exists a lower limit that corresponds to the equality in (1) and increases as the cell cycle duration increases. As P sym decreases, the doubling time approaches this limit asymptotically (Fig. 3b).

(c) Necrosis rate of differentiated cells (R NDiff ): According to the simulations that have been performed, this parameter appears to have no noticeable effect on tumor doubling time. There exists no upper or lower limit for its values, as the necrosis rate of differentiated cells does not appear in (1) (Fig. 3c). (vi) The sum of LIMP and stem cells relative populations decreases slightly as P sleep increases (in the order of 10 -2 ) to account for the increase of necrotic cells relative population (results not shown).

Response to Chemotherapy Parametric Studies

The effect of the cell kill ratio (i.e. the percentage of LIMP and stem cells hit by the chemotherapeutic agent after each drug administration) on the tumor's volume and diameter is shown in Fig. 5 for various cell cycle durations. The simulation algorithms address the case of primary chemotherapy ("neo-adjuvant" chemotherapy) with singleagent epirubicin (administration once every 3 weeks for 4 consecutive cycles) for early breast cancer patients, in accordance with the TOP trial. The simulated tumor grows freely for 3 weeks after the end of chemotherapy. The percentage decrease in tumor volume/diameter is calculated based on the following equation:

Percentage decrease = ((X initial -X final ) ×100%)/X initial (5) where X refers to tumor volume or diameter.

The following observations can be noted:

(i) Obviously the outcome of a chemotherapeutic scheme strongly depends on the pharmacodynamics of the chemotherapeutic agent used. The model successfully simulates the shrinkage of a tumor to a higher degree as the cell kill ratio increases. The value of the cell kill ratio could be thought of as summarizing important genetic determinants influencing a particular patient's tumor response to chemotherapy.

(ii) The free growth volume doubling time of a tumor is expected to influence the effectiveness of a chemotherapeutic fractionation scheme. In general, the time for tumor volume to double (volume doubling time) is determined by three main factors: the cell cycle duration, T c , the growth fraction, GF, which is the percentage of actively proliferating cells, and the rate of cell loss (e.g. through necrosis or apoptosis) [START_REF] Steel | Clonogenic cells and the concept of cell survival[END_REF]. Our simulation studies involved the cell cycle duration, as a determinant of tumor doubling time. By increasing the value of the cell cycle duration, while keeping the rest of the model parameters constant, a slower evolving tumor can be simulated.

The effectiveness of the chemotherapeutic scheme is pronounced for tumors with higher doubling times, due to the restrained repopulation of the tumor between two therapeutic sessions.

(iii) There exists a limit in the cell kill ratio values, below which the specific therapeutic scheme fails to shrink the tumor. This limit depends on the tumor's doubling time. As doubling time increases (due to an increase in the cell cycle duration) this limit decreases.

Mimicking a branch of the Trial of Principle (TOP) trial. Towards clinical adaptation of the simulation model.

As presented in previous paragraphs, the model's behavior substantiates its potential to serve as the basis of a treatment optimization system following a successful completion of the clinical adaptation and validation process which will rely on the use of anonymized real data before and after chemotherapeutic treatment for the case of the TOP breast cancer clinical trial. The TOP trial aims to evaluate the topoisomerase II alpha gene amplification and protein overexpression as markers predicting the efficacy of epirubicin in the primary treatment of breast cancer patients. Regarding the correlation of topoisomerase II alpha (TOP2A) amplification/ expression with response to anthracyclines, and epirubicin in particular, controversial results have been reported by several investigators during the last decade (e.g. [START_REF] Knoop | Retrospective Analysis of Topoisomerase IIa Amplifications and Deletions As Predictive Markers in Primary Breast Cancer Patients Randomly Assigned to Cyclophosphamide, Methotrexate, and Fluorouracil or Cyclophosphamide, Epirubicin, and Fluorouracil: Danish Breast Cancer Cooperative Group[END_REF][START_REF] Gennari | HER2 Status and Efficacy of Adjuvant Anthracyclines in Early Breast Cancer: A Pooled Analysis of Randomized Trials[END_REF][START_REF] Pritchard | HER-2 and Topoisomerase II As Predictors of Response to Chemotherapy[END_REF].

The presented simulation model is characterized by the ability to incorporate actual molecular data in order to predict the outcome of a chemotherapeutic treatment in a patient-individualized context. The effect of a patient's molecular profile on its tumor's response to chemotherapy can be modeled through adequate adaptation of the parameters related to the pharmacodynamics of the chemotherapeutic agent. Depending on the status and expression of critical genes such as topo II α, p53 HER-2/neu, etc., the drug pharmacodynamics is expected to vary around its population based mean value.

In the present work, two real clinical cases from the TOP trial with variable molecular profile (including topo II α amplification status) have been simulated. More specifically, the patient specific data that have been provided and utilized by the model are the following:

(i) Maximum dimension of the initial tumor before the beginning of the therapy as measured by ultrasound examination and the date that the above examination was performed.

(ii) Maximum dimension of the tumor during or after the completion of chemotherapy (depending on the available data) as measured by ultrasound examination and the date that the above examination was performed.

(iii) Chemotherapeutic scheme (number of therapeutic sessions, time and dose of epirubicin administration)

(iv) Molecular profile information, including topo II α gene amplification status as measured by FISH (fluorescent in situ hybridization). Amplification of the gene was defined as a relative copy number ratio ≥2.

(v) Tumor histological grade (determining the degree of tumor's differentiation). A tumor is described in TOP trial's CRF forms as poorly, moderately or well differentiated according to the the Elston-Ellis grade system [START_REF] Elston | Pathologic prognostic factors in breast caner. I. The valued of histological grade in breast cancer: experience from a large study with longterm follow-up[END_REF]. The values assigned to the code input parameters and the resulting tumor characteristics are presented in Table III. Critical properties of the resulting tumor, such as doubling time (232days) and percentage of proliferating cells (2.6%) are in accordance with breast cancer literature data (e.g. [START_REF] Spratt | Association between the actual doubling times of primary breast cancer with histopathologic characteristics and Wolfe's parenchymal mammographic patterns[END_REF][START_REF] Al-Hajj | Prospective indentification of tumorigenic breast cancer cells[END_REF][START_REF] Liu | Mammary stem cells, self-renewal pathways and carcinogenesis[END_REF][START_REF] Stingl | Molecular heterogeneity of breast carcinomas and the cancer stem cell hypothesis[END_REF].

For each case study one possible combination of the model input parameter values that successfully simulates the clinical case i.e. predicts the actual tumor volume shrinkage and satisfies several biological boundary conditions is reported below. Fig. 7 shows the time course of the following quantities of interest for the two clinical cases: population of proliferating cells (stem and LIMP) and total cell population (Fig. 7a), and tumor diameter (Fig. 7b). The above populations include both intact cells and cells affected by IV).

The values assigned to the cell kill ratio are shown in Table V. The value of the CKR that produced good agreement with the clinical data is reported as the suggested "apparent CKR" for each particular clinical case, which can be thought of as summarizing important genetic determinants influencing the tumor's response to epirubicin monotherapy. Table V includes topo II a amplification status for demonstration purposes (since topo II a amplification status has a particular role in the TOP trial), but the apparent CKR value incorporates the effect of all critical genetic factors constituting the tumor's molecular profile. It is noted that no statistical conclusions can be drawn concerning specifically the effect of topo II α amplification status on the therapeutic outcome, a study that is out of the scope of the present work. In conclusion, by properly adjusting the pharmacodynamics of epirubicin as reflected in the cell kill ratio value, the patient's molecular profile data (including topo II a amplification status) can be taken into account. The relevant methodology is to be refined according to the details provided by the TOP trial clinical cases.

Discussion

The model presented in this paper, having its roots in previous work of the In Silico [START_REF] Dean | Tumour stem cells and drug resistance[END_REF]. Cancer stem cells may be naturally resistant to chemotherapy; indeed considerable evidence exists for many types of tumors, including breast cancer [START_REF] Dean | Tumour stem cells and drug resistance[END_REF][START_REF] Diehn | Analyzing the Sensitivity of Breast Cancer Stem Cells to Ionizing Radiation and Chemotherapy[END_REF][START_REF] Fillmore | Human breast cancer cell lines contain stem-like cells that self-renew, give rise to phenotypically diverse progeny and survive chemotherapy[END_REF]. A major advantage of the present model is the fact that such a diversification of chemotherapeutic resistance is facilitated by its discrete nature and the characteristics of the particular cytokinetic model adopted.

A the GCs are initialized using the cell category /phase relative population values corresponding to this equilibrium state, a mathematically monotonous and biologically normal free growth behavior will be achieved. The challenge is to successfully locate the point beyond which equilibrium has been achieved and use the relative populations (or "fractions of populations") after that point for the correct initialization of the tumor.

Certain combinations of category / phase transition rates cannot sustain tumor growth (see section 4 for the condition for monotonic free growth). In such cases the method will correctly fail to create the initial tumor and a relevant warning message will be issued by the simulation system.

More specifically, the technique consists in the following: A limited number of geometrical cells N GCs are considered. Each GC initially contains a small number of stem cells, e.g. 100, residing in the various cell cycle phases (G1, S, G2, M) and the G0 phase. Time initialization, i.e. the time already spent by clustered stem cells in the phase they reside, is assigned using a pseudo-random generator. Different random number sets are assigned to different GCs. The aim is to avoid artificial synchronizations which would result in the group of GCs considered behaving as one big GC. The system is left to evolve and produce all cell category populations (distributed to the various cell phases). The code execution has to continue until equilibrium is reached and the various cell categories population percentages have been stabilized. The equilibrium condition applied is described by the following inequality:

for the rest of the mitotic potential cell categories for various proliferation statuses.

Taking into consideration the dependence of the execution time on the number of GCs II. II. Higher values of cell cycle duration correspond to more slowly evolving tumors, hence tumors with larger doubling time.

  clinically adapted using real clinical trial multiscale data and clinically validated within the framework of the EC funded projects ACGT (ACGT) and ContraCancrum (ContraCancrum).

  Two types of stem cell division are possible: symmetric and asymmetric. An asymmetric stem cell division gives rise to one daughter cell with stem cell fate and one daughter cell (limited mitotic potential or committed progenitor cell) that follows the path towards differentiation. A symmetric stem cell division gives rise to two daughter cells both with a stem cell fate (Morisson and Kimble, 2006). Specifically, the following five categories of cancer cells are considered in the model: i. Stem cells ii. Limited mitotic potential (LIMP) or restricted/committed progenitor cells: cells able to perform a limited number of divisions before becoming terminally differentiated iii. Differentiated cells: terminally differentiated cells iv. Apoptotic cells: cells that have died through apoptosis v. Necrotic cells: cells that have died through necrosis Stem, LIMP and differentiated cells constitute three categories with distinct mitotic potential.

Fig. 1 .

 1 The time registers of the various cell subclasses within each geometrical cell are updated and the cytokinetic diagram is applied within each GC as follows. Spontaneous apoptosis induced cell loss from each non treatment perturbed cell cycle phase and the G0 phase is calculated for each cell category based on the spontaneous apoptotic rates assumed. Any necessary transitions between equivalence subclasses (G1oS, SoG2, G2oM, MoG1 or MoG0) take place for biological cells clustered in the same subclass. If the mean time that the clustered cells have spent in the corresponding phase has become equal to or larger than the phase duration then the cells enter a new phase and equivalence subclass.

For

  stem and LIMP cells: G0oG1 or G0oNecrosis or G0oApoptosis. For differentiated cells: DifferentiatedoNecrosis or DifferentiatedoApoptosis. For dead cells of any mitotic potential category: ApoptosisoCell disappearance, NecrosisoCell disappearance. Most of the corresponding rates are parameters of the model (see Table

  (a) Premature extensive tumor fragmentation: some GCs belonging to the tumor become artificially separated from the main tumor mass. (b) Vacuum enclosures: holes that correspond to empty GCs are created inside the tumor. The algorithm developed in order to avoid the occurrence of the above side effects is given in appendix B. It should be noted that explicit information regarding tumor fragmentation induced by

( d )

 d Cell cycle duration (T c ): Higher values of T c are associated with slowly evolving tumors and higher values of doubling time, on the condition that the rest of the model parameters are kept constant. However the same doubling time can result from different values of cell cycle duration after proper adjustment of the rest of the model parameters.7.1.2 Relative cell category populationsAs far as the effect of the model's parameters on the relative populations of cells of various categories is concerned, P sleep and T c have been chosen as characteristic parameters for such a study. A presentation of a more detailed analysis covering all the model's parameters is out of the scope of the present paper. Based on Fig.4the following observations can be made:(i) It is expected that variations in the value of the fraction of cells that enter G0 phase after mitosis would mainly result in a redistribution of LIMP and stem cells in the dormant and the proliferative cell cycle phases. Indeed, according to our simulation results, as P sleep increases the percentage of proliferating cells (proliferating LIMP and stem cells residing in the various proliferative cell cycle phases) decreases whereas the percentage of cells resting in G0 phase (dormant LIMP and stem cells) increases.(ii) Increasing cell cycle duration T c results in an increase of proliferating cells (the sum of LIMP proliferating and stem proliferating cells) at the expense of mainly the G0 and dead cell populations.(iii) The effect of P sleep on differentiated cells relative population appears to be insignificant. For higher values of the cell cycle duration there is a slight increase in the population, whereas for low values of the cell cycle duration the differentiated cells relative population slightly decreases with the increase of P sleep (in the order of 10 -3 for the whole range of P sleep ).(iv) For high values of the sleep fraction the population of dead cells increases. This is due to the increase of necrotic cells, a class that is populated by cells formerly residing in the G0 phase.(v) The relative population (fraction) of stem to LIMP cells remains unaffected by the parameter P sleep (results not shown).

Fig. 6

 6 Fig. 6 describes each case study in terms of the above presented data. The date of the

  the drug (and therefore destined to die) but not yet actually dead. Qualitatively a fairly expected and reasonable tumor dynamics behavior can be easily noticed. Shrinkage of the tumor after each chemotherapeutic session and tumor repopulation are successfully demonstrated. Most importantly though, in quantitative terms, a reduction in tumor size in agreement with the clinical data has been achieved in both cases. According to the clinical data the tumors in case A and B have a maximum diameter of 12mm in day 62 and 15mm in day 111 respectively, which is clearly reproduced by the simulation results (Fig 7b and Table

Oncology

  Group, National Technical University of Athens (In Silico Oncology Group) is an actual clinical trial driven model concerning the response of early breast cancer to epirubicin monotherapy. Available imaging-based, histological, molecular and treatment data are exploited in order to strengthen patient individualization modeling and in the future patient individualized treatment optimization. Clinical validation constitutes an undeviating prerequisite for clinically oriented multiscale cancer models. A branch of the Trial of Principle (TOP) trial has been addressed as a paradigm of a chemotherapy optimization targeted clinicogenomic trial. A basic feature of the simulation model is the explicit consideration of cells of distinct mitotic potential, i.e. stem cells of infinite proliferative potential, LIMP (or committed progenitor) cells of limited proliferative activity and terminally differentiated cells. It is well documented in the relevant literature that tumor sustenance is due to the existence of stem cells, which generally represent a small but potentially critical portion of the total tumor cell population

  further prominent feature of the model is the introduction of an improved tumor initialization technique, which dispenses shortcomings of previously used methods and is now incorporated into the hard code, thereby offering high flexibility to the user of the model. This is of particular importance since the clinically oriented nature of the model implies that the evolution of already fully developed clinical tumors is to be simulated. Proper tumor initialization excludes the possibility of artificial tumor behavior which could interfere with the interpretation of the simulation results.A number of parametric studies have been performed in order to study the model's behavior in relation to the value of each chosen input parameter with the primary aim to deepen and advance quantification of our understanding of tumor response to chemotherapeutic treatment in the breast cancer and more specifically the TOP trial context. Similar studies are to be performed for all code parameters in order to thoroughly analyze the sensitivity of the model's behavior to its parameters variations.An extensive use of the clinical trial data is expected to crucially support this procedure, while pertinent optimization techniques such as artificial neural networks, genetic algorithms etc. have been planned to be used in this context. By using real medical data in conjunction with plausible values for the model parameters a reasonable prediction of the actual tumor volume shrinkage for two clinical cases of the TOP trial has been made possible. The value of the cell kill ratio parameter that produced good agreement with the clinical data is suggested as the "apparent CKR" for each particular clinical case, which can be thought of as summarizing important genetic determinants influencing the tumor's response to epirubicin monotherapy. The whole system can be seen as a paradigmal case of hybridizing clinicogenomics trials with multiscale computer modeling primarily focused on the cell and higher levels. For each clinical case one possible combination of the model input parameter values that successfully predicts the actual tumor volume shrinkage has been reported. It should be pointed out that the large number of inherent biological boundary conditions (e.g. monotonic increase of all tumor cell class populations for an imageable freely and smoothly growing tumor) dramatically limits the number of possible solutions (i.e. sets of parameter values that are able to predict real tumor shrinkage following treatment). This is particularly important since such limitations drastically facilitate the approach to the solution best representing clinical reality for each given medical data case. The model has been developed to support and incorporate individualized clinical data such as imaging data (e.g., CT, MRI, PET slices, possibly fused) (e.g. Marias et al 2007), including the definition of the tumor contour and internal tumor regions (proliferating, necrotic), histopathologic (e.g., type of tumor) and the genetic data (e.g., p53 status, if available). The use of anonymized real data before and after chemotherapeutic treatment for the case of the TOP breast cancer clinical trials constitute the basis of the clinical adaptation and validation process. Obviously as more and more sets of medical data are exploited the reliability of the model "tuning" is expected to increase. Finally, a detailed description of various technical issues involved in the present work (such as image processing, grid execution and visualization techniques) will be the topic of a dedicated paper. 10. Conclusions An advanced multiscale model of clinical tumor response to chemotherapy has been presented. The paradigm of early breast cancer treated by epirubicin according to a branch of the TOP (Topoisomerase II Alpha Gene Amplification and Protein Overexpression Predicting Efficacy of Epirubicin) trial has been addressed. The model, stemming from previous work of our research group, is characterized by several crucial new features such as an advanced generic cytokinetic model for the response of a tumor cell to chemotherapeutic treatment, the explicit distinction of proliferating cells into stem cells and committed progenitor cells of limited proliferative capacity and the adaptation of the relative populations of the various cell categories/phases to the cell category/phase transition rates for free tumor growth, which constitutes an important novel tumor initialization strategy. A sensitivity analysis regarding critical parameters of the model has revealed their effect on the behavior of the biological system. In particular, the joint effect of the cell cycle duration and the fraction of newborn tumor cells that enter the G0 phase on the relative cell category populations, the joint effect of the cell cycle duration and the fraction of cells that perform symmetric division on the doubling time and the joint effect of the cell cycle duration and the necrosis rate of differentiated cells on the tumor doubling time have been studied. Furthermore, the joint effect of the cell kill ratio and the cell cycle duration on the therapeutic outcome (volume and diameter decrease percentages) has been also presented. Finally, a realistic time course of the tumor diameter for two clinical cases of the TOP trial has been achieved by selection of reasonable model parameter values, thus demonstrating a possible adaptation process of the model to real clinical trial data. The above studies and their results support the potential of the model to serve as both a theoretical investigation and a patient specific treatment optimization tool following completion of an ongoing clinical adaptation, optimization and validation process.

  FIGURE CAPTIONS

Fig. 2 .

 2 Fig.2. Simplified flow chart of the simulation procedure for a macroscopically

Fig. 3 .

 3 Fig. 3. The effects of (a) the fraction of cells entering the G0 phase following mitosis

Fig. 4 .

 4 Fig. 4. The effect of the fraction of cells entering the G0 phase following mitosis (P sleep )

Fig. 5

 5 Fig.5The effect of cell kill ratio (CKR) on tumor volume and tumor diameter. Three

Fig. 6 .

 6 Fig.6. The epirubicin administration schedules considered in case studies A and B.

Fig. 7 .

 7 Fig.7. Simulated time course of the population of proliferating cells and total cell

Fig

  Fig. A.3. Study of the model convergence. Relative variation of the average of N = Tc =

Acknowledgements

This work has been 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). Fruitful discussions with Norbert Graf, University Hospital of Saarland, Manolis Tsiknakis, Foundation for Research and Technology Hellas and Nikolaos Uzunoglu, National Technical University of Athens as well as constructive feedback from the anonymous reviewers are duly acknowledged.

Appendix A

The principle of the tumor constitution initialization technique is to start with a small number of stem cells and with specific cell category transition rates that are assumed to hold true for a relatively small time interval about the treatment baseline. Specific values are assigned to the phase durations and transition rates. Gradually, all cell categories and phases become populated and after sufficient time the relative cell categories populations tend to reach an equilibrium state. If in subsequent simulations Here again the fluctuation is still not significant.

(d) The small positive real number ε Ideally, equilibrium of the relative cell populations has been reached when the relative variation of the average of N consecutive values of the relative population of each cell equivalence subclass over time has become zero. Since this may happen at infinity due to small quantization errors, the small real positive value of ε = 0.00001 = 10 -5 has been adopted. The latter has been shown to ensure convergence in the breast cancer parameter subspace considered.

(e) The number k of consecutive average values considered

Satisfaction of the equilibrium condition for k = 10 consecutive averages has been proved to ensure convergence in the parameter subspace considered.

In summary, the optimal parameter values for achieving convergence in the tumor constitution initialization procedure, for the specific parameter subspace considered in this paper, are N GCs = 64, N = T c , M = T c , ε = 10 -5 and k =10. It is noted that the initialization code execution time ranges from a few fractions of a second to a few seconds on a standard laptop depending on the particular characteristics of the tumor (e.g. cell cycle duration, fraction of cells entering the G0 phase following mitosis etc.).

Appendix B

The algorithm developed for the avoidance of premature extensive fragmentation and vacuum enclosures artifacts:

(I) detects tumor occupied GCs that are surrounded by empty GCs in a "6-GCs Neighborhood" and moves their contents (by one GC at each time step) towards the tumor's center of mass. The direction of movement is chosen based on the minimum distance of the GC under consideration from the center of mass along the x, y, z coordinates. The corresponding quantity to be calculated each time is the following: min{abs(GC.x-center.x), abs(GC.y-center.y), abs(GC.z-center.z)} where abs() denotes absolute value, GC.x , GC.y and GC.z are the x , y and z coordinates of the GC respectively and center.x, center.y and center.z are the x, y and z coordinates of the tumor's center of mass respectively.

If more than one direction is characterized by the same minimum distance then a random number generator is used for the selection of the movement direction.

(II) detects empty GCs that are surrounded by occupied GCs in a "6-GCs Neighborhood" and fills them with the contents of adjacent GCs by applying the tumor shrinkage procedure described above.

M time steps (hours) over time. Triangles correspond to M equal to the numerical value of the cell cycle duration (96 h). Circles correspond to M = 50. Both curves correspond to the same tumor. The convergence of this parameter to zero is essential for the equilibrium condition described in appendix A to be applied.

4. Figure 1 4. Figure 2 4. Figure 3a 4. Figure 3b 4. Figure 3c 4. Figure 4 4. Figure 5 4. Figure 6 4. Figure 7a 4. Figure 7b 4.