Heike Meiselbach

Heinrich Sticht
email: h.sticht@biochem.uni-erlangen.de

Effect of the SH3-SH2 domain linker sequence on the structure of Hck kinase

Keywords:

The coordination of activity in biological systems requires the existence of different signal transduction pathways that interact with one another and must be precisely regulated. The Src-family tyrosine kinases, which are found in many signaling pathways, differ in their physiological function despite their high overall structural similarity. In this context, the differences in the SH3-SH2 domain linkers might play a role for differential regulation, but the structural consequences of linker sequence are yet poorly understood. We have therefore performed comparative molecular dynamics simulations of wildtype Hck and of a mutant Hck, in which the SH3-SH2 domain linker is replaced by the sequence from the homologous kinase Lck. The simulations reveal that linker replacement does not only affect the orientation of the SH3 domain itself, but also leads to an alternative conformation of the activation segment in the Hck kinase domain. The sequence of the SH3-SH2 domain linker thus exerts a remote effect on the active site geometry and might therefore play a role in modulating the structure of the inactive kinase or for the fine-tuning of the activation process itself.

Introduction

The Src family of non-receptor tyrosine kinases to date comprises nine members (Src, Fyn, Lck, Hck, Yes, Fgr, Yrk, Blk, and Lyn) which play critical roles in eukaryotic signal pathways that control a diverse array of processes such as cell growth, differentiation, activation and transformation [START_REF] Thomas | Cellular functions regulated by src family kinases[END_REF][START_REF] Blume-Jensen | Oncogenic kinase signalling[END_REF][START_REF] Brown | Regulation, substrates and functions of src[END_REF]. These proteins have a common molecular architecture that includes five distinct regions [START_REF] Tatosyan | Kinases of the src family: Structure and functions[END_REF][START_REF] Engen | Structure and dynamic regulation of src-family kinases[END_REF]: a unique N-terminal region with sequences for lipid attachment, the regulatory SH3 and SH2 domains and a kinase domain, followed by a negative regulatory C-terminal tail (Fig. 1b). The SH2 and SH3 domains bind specific phosphotyrosine and proline-rich motifs, respectively. The catalytic domain consists of two lobes, which form a cleft that serves as a docking site for ATP and other substrates representing the active site of the kinase.

The crystal structures of the inactive ("closed") forms of c-Src [START_REF] Williams | The 2.35 a crystal structure of the inactivated form of chicken src: A dynamic molecule with multiple regulatory interactions[END_REF][START_REF] Xu | Three-dimensional structure of the tyrosine kinase csrc[END_REF][START_REF] Xu | Crystal structures of c-src reveal features of its autoinhibitory mechanism[END_REF] and Hck [START_REF] Sicheri | Crystal structure of the src family tyrosine kinase hck[END_REF][START_REF] Schindler | Crystal structure of hck in complex with a src family-selective tyrosine kinase inhibitor[END_REF] reveal that intramolecular interactions involving the SH2 and SH3 domains are essential for suppression of kinase activity despite being far away from the active site. In fact, these regulatory domains are bound to the distal side of the catalytic domain, restricting its conformational flexibility to hinder productive ATP binding. The SH3 domain is docked to the N-terminal lobe of the kinase domain via a proline-containing motif in the SH2-kinase linker (Fig. 1b).

Activation of Src-family kinases (SFKs) requires the disruption of the intramolecular interactions, the dephosphorylation of the C-terminal tail by phosphatases, or the binding of high-affinity ligands to the SH2 or the SH3 domain [START_REF] Sicheri | Crystal structure of the src family tyrosine kinase hck[END_REF][START_REF] Schindler | Crystal structure of hck in complex with a src family-selective tyrosine kinase inhibitor[END_REF][START_REF] Lerner | Sh3-dependent stimulation of src-family kinase autophosphorylation without tail release from the sh2 domain in vivo[END_REF][START_REF] Moarefi | Activation of the src-family tyrosine kinase hck by sh3 domain displacement[END_REF]. Additionally, the active conformation becomes stabilized upon autophosphorylation of Tyr416 (in Src) in the activation segment and is required for maximal kinase activity [START_REF] Sw | Structural biology of protein tyrosine kinases[END_REF][START_REF] Hubbard | Protein tyrosine kinase structure and function[END_REF]. In this case, the phosphorylated activation segment, which contains the highly conserved aspartatephenylalanine-glycine (DFG) motif, becomes less ordered and moves out of the active site, thereby allowing substrate access [START_REF] Yamaguchi | Structural basis for activation of human lymphocyte kinase lck upon tyrosine phosphorylation[END_REF].

An additional feature, which turned out to be important for the activation process, is the sequence of the linker that connects the SH3 and SH2 domain. In wildtype Src and Hck, a short rigid linker joins the regulatory SH3 and SH2 domains and leads to dynamic coupling of the two domains in the closed conformation, thereby forming a stiff clamp [START_REF] Harrison | Variation on an src-like theme[END_REF][START_REF] Young | Dynamic coupling between the sh2 and sh3 domains of c-src and hck underlies their inactivation by cterminal tyrosine phosphorylation[END_REF].

Computational studies and site-directed mutagenesis experiments have shown that the disruption of the structural properties of the rigid linker by introduction of glycine mutations reduces the strength of the clamp and facilitates activation of the kinase [START_REF] Young | Dynamic coupling between the sh2 and sh3 domains of c-src and hck underlies their inactivation by cterminal tyrosine phosphorylation[END_REF].

In Hck and most other SFKs the connecting linker stabilizes the SH3-SH2 domain orientation by a salt-bridge between a conserved glutamate of the linker (E147 in Hck) and a lysine (K104 in Hck) of the SH3 domain. In the homologous kinase Lck, E147 is replaced by proline and therefore the salt-bridge cannot be formed. Accordingly, the isolated Lck SH3-SH2 domain pair shows an ill-defined domain orientation in NMR-spectroscopic studies [START_REF] Arold | The role of the src homology 3-src homology 2 interface in the regulation of src kinases[END_REF][START_REF] Hofmann | Binding, domain orientation, and dynamics of the lck sh3-sh2 domain pair and comparison with other src-family kinases[END_REF].

The effect of the Lck linker sequence on the intact kinase, however, is yet poorly understood due to the lack of an intact Lck crystal structure containing both regulatory and kinase domains.

To assess whether the differences in linker sequence also affect the structure of inactive SFKs, we have investigated the effect of linker replacement using Hck as a model system.

Comparative molecular dynamics simulations were performed for wildtype Hck and for a mutant Hck (Hck*), in which the SH3-SH2 domain linker was replaced by the respective sequence from Lck. To address the effect of the inhibitor PP1 on kinase dynamics, simulations were performed in the presence and absence of PP1 totalling up to four systems studied. Our results show that an altered linker sequence does not only result in local changes, but also has a remote effect on the geometry of the active site.

Materials and methods

Preparation of starting structures

The crystal structures of inactive Hck complexed with the nucleoside analog PP1 (Protein Data Bank ID 1QCF, [START_REF] Schindler | Crystal structure of hck in complex with a src family-selective tyrosine kinase inhibitor[END_REF]), was used to generate the starting structures. For consistency, the sequence numbering from entry 1QCF was also adapted in the present work. To obtain the model of mutant Hck (termed Hck*), the wildtype SH3-SH2 linker sequence ( 140 V-D-S-L-E-T-E-E 147 ) was replaced with the linker sequence form Lck (A-N-S-L-E-P-E-P) (Fig. 1a) using the Sybyl 7 [20]. Inspection of the model revealed no steric clashed indicating that the Lcklinker sequence can be readily accommodated in the Hck scaffold. The systems Hck and Hck* were simulated each with and without inhibitor PP1.

The force field parameters for PP1, which are not available in the AMBER package, were obtained as follows: The initial coordinates of PP1 were extracted from the 1QCF pdb file.

ArgusLab [21] was utilized to add missing hydrogen atoms. Then, the structure was subjected to two consecutive geometry optimizations with Gaussian03 [START_REF] Frisch | Gaussian03[END_REF] using the ab initio methods HF/MIDI!, and HF/6-31G(d). For all quantum mechanical geometry optimizations the stationary points found were ensured to be true minima by the calculation of the vibrational frequencies. The atomic charges for PP1 were then obtained following the established procedure [23,[START_REF] Cornell | Application of resp charges to calculate conformational energies, hydrogen bond energies and free energies of solvation[END_REF] by fitting the charges to the HF/6-31G(d) computed electrostatic potential using the Antechamber tool from the AMBER program suite.

Molecular dynamics (MD) simulations

All MD simulations were performed using the PMEMD module of AMBER 8.0 [START_REF] Case | Amber 8. University of California[END_REF][START_REF] Pearlman | Amber, a package of computer programs for applying molecular mechanics, normal mode analysis, molecular dynamics and free energy calculations to simulate the structural and energetic properties of molecules[END_REF] with the parm99 force field [START_REF] Cheatham | A modified version of the cornell et al. Force field with improved sugar pucker phases and helical repeat[END_REF][START_REF] Cornell | A second generation force fiel for the simulaton of proteins, nucleic acids and organic molecules[END_REF]. For the organic compound PP1 the general AMBER force field (gaff) [START_REF] Wang | Development and testing of a general amber force field[END_REF] was used. The parameters for the phosphotyrosine present in the C-terminal kinase tail were assigned according to Homeyer et al. [START_REF] Homeyer | Amber force-field parameters for phosphorylated amino acids in different protonation states: Phosphoserine, phosphothreonine, phosphotyrosine, and phosphohistidine[END_REF].

An appropriate number of Na+ counter ions was added to neutralize the system and afterwards the molecules were solvated in a box of water using the TIP3P water model [START_REF] Jorgensen | Comparison of simple potential functions for simulating liquid water[END_REF] with at least 10 Å of water around every atom of the solute. All structures were minimized in a three-step procedure by using the SANDER module of AMBER following a previously established protocol [START_REF] Wartha | Molecular dynamics simulations of hiv-1 protease suggest different mechanisms contributing to drug resistance[END_REF].

After energy minimization, the system was equilibrated for 0.1 ns, raising the temperature gradually to 298 K by coupling to a temperature bath with a time constant of 0.2 ps, as described by Berendsen [START_REF] Berendsen | Molecular dynamics with coupling to an external bath[END_REF]. Subsequently, 30-ns MD simulations with standard NPT conditions were performed for data collection. Classical equations of motion were propagated numerically with a time step of 1 fs. A cutoff distance of 10 Å was used for the non-bonded interactions, and the particle mesh Ewald method [START_REF] Darden | Particle mesh ewald. An n.Log(n) method for ewald sums in large systems[END_REF] was employed to calculate the longrange electrostatic interactions. All bonds involving hydrogen atoms were constrained using the SHAKE procedure [START_REF] Ryckaert | Numerical integration of the cartesian equations of motion of a system with constraints: Molecular dynamics of n-alkanes[END_REF].

The convergence of temperature, pressure, and energy of the systems, as well as the atomic root mean square deviations of the structure, were used to verify the stability of the systems.

For the visualization and structural analysis, the programs Sybyl [20], AMBER [START_REF] Case | Amber 8. University of California[END_REF], VMD [START_REF] Humphrey | Vmd --visual molecular dynamics[END_REF], and DS ViewerPro 6 [37] were used.

Results and discussion

Influence of the SH3-SH2 linker sequence on the global structure and dynamics of PP1inhibited Hck

To determine the effect of the SH3-SH2 linker on the global structure of Src kinases, 30-ns molecular dynamics simulations of wildtype Hck and of Hck with a mutated SH3-SH2 linker (Hck*) were carried out. To address the effect of the inhibitor PP1 bound to the active site, comparative MD simulations were performed for Hck and Hck* with and without PP1 totalling up to four systems studied.

The mobility of the distinct domains and the conformational changes of the kinase structure during the molecular dynamics simulations were first analyzed by monitoring the distances between the centers-of-mass of the different domains. In this analysis, the two subdomains of kinase domain (N-lobe and C-lobe) were considered separately (Fig. 2, right panel).

For the PP1-bound kinases (Fig. 2) the most prominent difference is detected for the distance between the SH3-domain and the N-lobe of the kinase domain. The distance of approximately 29 Å in the start structure of the inhibited Hck* increases up to 32 Å during the simulation (Fig. 2a). Since the only difference between the Hck and Hck* structure is the composition of the SH3-SH2 kinase linker, it is most likely that the linker exchange is responsible for the different behaviour of the Hck* system. Inspection of the distance between the SH3-SH2 linker and the active site reveals that the linker itself also becomes significantly displaced during the MD simulation (Fig. 3).

The rearrangement of the SH3-SH2 linker and of the SH3 domain can also be seen as the most prominent motion in the overlay of the initial and final structure of Hck* (Fig. 4b). In contrast, only minor changes are observed for wildtype Hck during the simulation (Fig. 4a).

The changes of the linker orientation in Hck* also affect the adjacent parts of the SH2 domain (in particular helix B), while the SH2 interface with the C-lobe of the kinase domain is almost unchanged (Fig. 4b). The kinase domain itself is also only marginally affected by the linker replacement and shows only minor fluctuations both in the Hck and in the Hck* simulation compared to the starting structure (Fig. 4a,b).

Effect of linker replacement on the uninhibited Hck kinase

We addressed the question whether linker replacement also has similar structural effects in the absence of the PP1 inhibitor. The MD simulations performed in the absence of PP1 generally reveal fluctuations for both systems indicating that Hck and Hck* need to adapt to the absence of the inhibitor (Fig. 5). Similarly to the PP1-bound simulation, the distance between SH3 domain and N-lobe becomes also larger in the unliganded Hck* compared to Hck (see Fig. 5a vs. 2a). This difference, however, is only in part caused by an increased distance in unliganded Hck* ─ a second effect is a decrease of the respective distance in wildtype Hck (Fig. 5a). The latter effect might be a consequence of the tight SH3-SH2 clamp in the wildtype, which results in a tighter compression of the active site after removal of the inhibitor (Fig. 6a). Such an effect is not observed for Hck* (Fig. 6b) suggesting that the SH3-SH2 clamp on the kinase domain is less strong than in Hck.

The simulations with and without PP1 revealed that the mutation of the linker sequence influences the position of the SH3 domain with respect to the N-lobe. Thus, the identity of the linker sequence does not only affect the structure of isolated SH3-SH2 domain pair as shown in previous studies [START_REF] Arold | The role of the src homology 3-src homology 2 interface in the regulation of src kinases[END_REF][START_REF] Hofmann | Binding, domain orientation, and dynamics of the lck sh3-sh2 domain pair and comparison with other src-family kinases[END_REF], but also the structure of the full-length kinase in its closed conformation. In the next step, we investigated whether linker replacement also has an effect on the geometry of the active site itself.

Effects of the SH3-SH2 linker sequence on the activation segment

The active site, which is located between the N-and C-lobe of the kinase domain, contains an activation segment with a highly conserved aspartate-phenylalanine-glycine (DFG) motif. This segment does not only undergo significant structural changes during the activation process itself, but can also adopt at least two distinct conformations in the inactive kinase. This structural feature also plays a critical role for the drug binding properties of tyrosine kinases [START_REF] Levinson | A src-like inactive conformation in the abl tyrosine kinase domain[END_REF][START_REF] Nagar | Crystal structures of the kinase domain of c-abl in complex with the small molecule inhibitors pd173955 and imatinib (sti-571)[END_REF][START_REF] Nagar | Structural basis for the autoinhibition of c-abl tyrosine kinase[END_REF][START_REF] Dar | Small molecule recognition of c-src via the imatinib-binding conformation[END_REF][START_REF] Jacobs | Classifying protein kinase structures guides use of ligand-selectivity profiles to predict inactive conformations: Structure of lck/imatinib complex[END_REF][START_REF] Seeliger | C-src binds to the cancer drug imatinib with an inactive abl/c-kit conformation and a distributed thermodynamic penalty[END_REF].

Due to this important structural role of the DFG-motif, the subsequent analysis focused on the conformation of the respective residues in Hck and Hck*. Figure 7 shows the analysis of the reflects a local energy minimum which is formed during a transition from the DFG-in to the DFG-out conformation.

Since our results indicate that mutations of the SH3-SH2 linker affect the active site geometry of Hck, we have checked the literature, whether remote effects of mutations on the active site conformation have already been reported for other SFKs in the past. An interesting example comes from c-Src, which can adopt either a DFG-in or a DFG-out conformation. For wildtype c-Src, the DFG-out conformation is energetically quite unfavourable [START_REF] Aleksandrov | Molecular dynamics simulations show that conformational selection governs the binding preferences of imatinib for several tyrosine kinases[END_REF][START_REF] Seeliger | C-src binds to the cancer drug imatinib with an inactive abl/c-kit conformation and a distributed thermodynamic penalty[END_REF]. This is reflected in the low affinity of c-Src for the drug imatinib, which exclusively binds the DFGout conformation. The structural origin for the low imatinib binding affinity of c-Src was also investigated by mutational studies [START_REF] Seeliger | C-src binds to the cancer drug imatinib with an inactive abl/c-kit conformation and a distributed thermodynamic penalty[END_REF]. The strongest increase in drug sensitivity was observed for three mutations (W260A, F405A, L407G), which destabilize the DFG-in conformation of c-Src. While F405 and L407 are located close to the imatinib-binding site, W260 is located in the N-lobe approximately 15 Å apart from the ligand indicating that longrange effects actually affect the drug binding properties.

Our present results further support the importance of long-range effects by demonstrating that the inactive conformation of the kinase domain cannot only be destabilized by mutations in the kinase domain itself, but also by mutations within the regulatory domains. The SH3-SH2 domain linker mutations studied here are more than 30 Å apart from the active site (Fig. 1b, Fig. 4). The sequence of the SH3-SH2 domain linker thus exerts a remote effect on the active site geometry and might therefore play a role in modulating the structure of the inactive kinase or for the fine-tuning of the activation process itself.

Figure captions

Fig. 1 (Fig. 2 Fig. 3 Fig. 4 Fig. 5 Fig. 6 Fig. 7

 1234567 Fig. 1 (a) Structures of the linker between regulatory domains from the Hck kinase (left), and Lck kinase (right). (b) Overlay of the crystal structure of the Hck kinase (orange, 1QCF) in complex with the inhibitor PP1 (blue), and the Hck* kinase (grey)

Acknowledgments

The authors thank Pia Rücker and Dr. Holger Dinkel for fruitful discussions, and A. Jens Meiselbach-Wilke for critically reading the manuscript. This work was supported by grants from the Deutsche Forschungsgemeinschaft.

backbone angles of the DFG motif and the N-terminally adjacent amino acid A403 during the MD simulations of the inhibited (Fig. 7 a-d) and uninhibited (Fig. 7 e-h) kinase.

No significant changes of the torsion angles φ and ψ were detected in the PP1-bound kinases Hck and Hck* (Fig. 7a-d). This indicates that the inhibitor efficiently fixes the active site conformation regardless of the SH3-SH2 linker sequence present. This rigidity is also consistent with the structural overlay showing that the geometry of the kinase domain is almost unaffected by the SH3-SH2 linker sequence (Fig. 4).

The analysis of the torsion angles of uninhibited wildtype Hck generally shows no major fluctuations (Fig. 7 e-h; black curves) and the angles are very similar to the inhibited kinases (Fig. 7 a-d). In contrast, the respective segment of Hck* adopts an alternative backbone geometry during the simulation (Fig. 7 e-h; red curves). After 0.3 ns a correlated change is detected for the neighboring amino acids A403(ψ) and D404(φ). Correlated changes of the backbone angles are observed for the amino acids D404(ψ), F405(φ ,ψ), and G406(φ) after approximately 4 ns. For the remaining simulation time these novel φ/ψ combinations remain remarkably stable and only minor fluctuations are detected.

Comparison of the active site conformation to that of other inactive SFKs

The conformation of the activation loop detected in our simulations was compared to that of other SFKs. At least two distinct stable conformations of this segment are observed in inactive kinases, which differ in the backbone conformation of the DFG-motif. These conformations were termed "DFG-in" and "DFG-out", respectively [START_REF] Levinson | A src-like inactive conformation in the abl tyrosine kinase domain[END_REF].

All starting structures of the present simulations correspond to the "DFG-in" conformation, which is also retained in the simulations of PP1-bound Hck/Hck* as well as in the simulation of uninhibited wildtype Hck. For the uninhibited Hck*, a distinct conformation forms during the MD simulation and the backbone conformation formed by Asp404 is somewhat reminiscent of the "DFG-out" conformation, which is also characterized by a negative Asp404 φ angle of approximately -130 to -150 degree. For residues F405 and G406 the conformation of Hck*, however, is clearly distinct from the DFG-out conformation previously observed in inactive c-Src and Lck [START_REF] Jacobs | Classifying protein kinase structures guides use of ligand-selectivity profiles to predict inactive conformations: Structure of lck/imatinib complex[END_REF][START_REF] Seeliger | C-src binds to the cancer drug imatinib with an inactive abl/c-kit conformation and a distributed thermodynamic penalty[END_REF]. One might speculate that the Hck* backbone conformation either represents a novel stable topology or that the 30-ns structure 20


Anna Fyrberg
email: anna.fyrberg@liu.se

Curt Peterson

Bertil Kågedal

Kourosh Lotfi

Induction of fetal hemoglobin and ABCB1 gene expression in 9-β-D-arabinofuranosylguanineresistant MOLT-4 cells

Keywords: 9-β-D-arabinofuranosylguanine, fetal hemoglobin, P-glycoprotein, microarray, hypomethylation

Induction of fetal hemoglobin and ABCB1 gene expression in 9-β-arabinofuranosylguanine-resistant MOLT-4 cells

Introduction

The anti-leukemic drug 9-β-D-arabinofuranosylguanine (AraG, Atriance ® ) was synthesized in the 1950s but only recently brought into clinical use for the treatment of relapsed or refractory T-cell acute lymphoblastic leukemia (T-ALL) or lymphoblastic lymphoma [START_REF] Kurtzberg | The long and winding road of the clinical development of Nelarabine[END_REF][START_REF] Deangelo | Nelarabine induces complete remissions in adults with relapsed or refractory T-lineage acute lymphoblastic leukemia or lymphoblastic lymphoma: Cancer and Leukemia Group B study 19801[END_REF][START_REF] Cohen | Approval summary: nelarabine for the treatment of T-cell lymphoblastic leukemia/lymphoma[END_REF]. The drug is used clinically as a prodrug, nelarabine, which is more water-soluble than AraG and is demethylated in blood to AraG [START_REF] Deangelo | Nelarabine induces complete remissions in adults with relapsed or refractory T-lineage acute lymphoblastic leukemia or lymphoblastic lymphoma: Cancer and Leukemia Group B study 19801[END_REF]. AraG is then activated intracellularly by the cytosolic enzyme deoxycytidine kinase (dCK) and the mitochondrial enzyme deoxyguanosine kinase [START_REF] Lotfi | Low level of mitochondrial deoxyguanosine kinase is the dominant factor in acquired resistance to 9-beta-Darabinofuranosylguanine cytotoxicity[END_REF][START_REF] Zhu | Incorporation of nucleoside analogs into nuclear or mitochondrial DNA is determined by the intracellular phosphorylation site[END_REF]. The high content of mitochondria in neurons can partially explain why severe side effects of AraG treatment include life-threatening neurotoxicity [START_REF] Kurtzberg | The long and winding road of the clinical development of Nelarabine[END_REF]. Inherited or acquired resistance to nucleoside analog treatment is a common phenomenon leading to treatment failure. In order to study possible resistance mechanisms to AraG, we performed a gene expression microarray analysis of the human leukemic T-ALL cell line MOLT-4 and its AraG-resistant variant.

In this study we show that traditional nucleoside analog resistance genes were only moderately affected in the resistant MOLT-4/AraG cells compared to parental cells. However, we could see a massive induction of fetal hemoglobin genes and the multidrug resistance gene ABCB1 coding for the drug efflux pump P-glycoprotein (P-gp). We suggest that AraG induces the expression of hemoglobin, ABCB1, and other genes by hypomethylation of DNA and that P-gp may be induced to help excrete fetal hemoglobin waste products from the cells since AraG is not considered a substrate of this pump.

Materials and methods

Reagents and drugs

9-β-D-arabinofuranosylguanine (AraG, Atriance ® ) was obtained from R.I. Chemical (Orange, CA, USA), 6-mercaptopurine (6-MP, Puri-Nethol ® ) was from GlaxoSmithKline (Solna, Sweden) 2-chloro-2´-arabino-fluoro-2´-deoxyadenosine (clofarabine, CAFdA, Evoltra ® ) was a kind gift from Dr. Howard Cottam, University of California (San Diego, CA. USA), 2fluoroarabinosyladenine (FaraA, Fludara ® ) was a gift from Dr. Zéev Shaked, Berlex (Alameda, CA, USA), gemcitabine (dFdC, Gemzar ® ) was from Eli Lilly (Solna, Sweden), daunorubicin (Dnr, Cerubidin ® ) was from Rhône-Poulenc Rorer (Bristol, UK), cladribine (CdA, Leustatin ® ) and laniquidar were from Johnson & Johnson (Schaffhausen, Switzerland), cytarabine (AraC, Cytosar ® ) was from Pharmacia & Upjohn (Stockholm, Sweden), and cyclosporin A (Sandiummun ® ) was from Novartis (Stockholm, Sweden). Radioactive substances were from Moravek Biochemicals (Brea, CA. USA), Western blot equipment and reagents were from Bio-Rad Laboratories (Hercules, CA, USA), cell-culturing reagents were from Gibco, Life Technologies (Paisley, UK), real-time PCR equipment and reagents were from Applied Biosystems (Foster City, CA, USA), and additional chemicals were purchased from Sigma-Aldrich (Stockholm, Sweden).

Cell lines

Cell lines used in this study were the T-cell acute lymphoblastic leukemia cell line MOLT-4 (American Type Culture Collection, ATCC) and variants made resistant to AraG (0.9 M) or 6-MP (5 M) as earlier described [START_REF] Lotfi | Low level of mitochondrial deoxyguanosine kinase is the dominant factor in acquired resistance to 9-beta-Darabinofuranosylguanine cytotoxicity[END_REF], and the human chronic myelogenous leukemia cell line K562 (ATCC). The human lung adenocarcinoma epithelial cell line A549 was a gift from Dr. S A Coulthard (Newcastle University, Northern Institute for Cancer Research, Newcastle upon Tyne, UK). Cells were cultured in RPMI-1640 media supplemented with 10% fetal calf serum, 100 IU/ml penicillin, 100 g/ml streptomycin, and 2 mM L-glutamine. Cells were kept at 37 o C with 5% CO 2 in a humidified atmosphere and were subcultured twice weekly.

Cytotoxicity assay

Sensitivity of the cells to various drugs was determined by the methyl thiazole tetrazolium (MTT) assay as previously described [START_REF] Fyrberg | Cell cycle effect on the activity of deoxynucleoside analogue metabolising enzymes[END_REF]. Triplicate analyses of 20,000 cells per well for MOLT-4 cells and 10,000 cells/well for A549 in 100 l aliquots in 96-well microtiter plates with 5 l of drug dilutions were performed. Sensitivity to the different substances, expressed as an IC 50 value, was calculated as the percent survival compared to control cells.

RNA extraction, reverse transcriptase PCR, and real-time PCR

RNA was extracted using the RNeasy Mini Kit (Qiagen, Solna, Sweden) and complementary DNA was produced using the High Capacity Archive kit (Applied Biosystems). Details of PCR primers and Taqman probes for dCK and dGK have been published previously [START_REF] Chandra | Resistance of leukemic cells to 2-chlorodeoxyadenosine is due to a lack of calciumdependent cytochrome c release[END_REF][START_REF] Fyrberg | Cell cycle dependent regulation of deoxycytidine kinase, deoxyguanosine kinase, and cytosolic 5'-nucleotidase I activity in MOLT-4 cells[END_REF]. To measure hemoglobin gamma (Hs00361131_g1) and ABCB1 (Hs00184500_m1) expression, Taqman® Gene Expression Assays (20X) were used, and the glyceraldehyde-3 phosphate dehydrogenase (GAPDH) gene (Hs00266705_g1) was used as a housekeeping control by dividing the mean quantity of the gene studied with the mean quantity of GAPDH expression after analysis with the Sequence Detection Software 1.

Protein extraction and enzyme activity measurements

Cells were centrifuged and washed in cold PBS, and the pellets were resuspended in a protein extraction buffer containing 50 mM Tris-HCl (pH 7.6), 4 mM dithiothreitol, 0.5 mM phenylmethylsulfonyl fluoride, 20% glycerol, and 0.5% Nonidet P40 and then disrupted using liquid nitrogen. After a final 15 min centrifugation at 15,000 × g, the cell supernatants were collected and assayed for dCK and dGK activities.

Enzyme activities were measured by a method published previously based on a modified radiochemical method [START_REF] Fyrberg | Cell cycle effect on the activity of deoxynucleoside analogue metabolising enzymes[END_REF]. Concentrations of the substrates were 10 M deoxycytidine ([8- Ci/mmol) for dGK. To inhibit the activity of thymidine kinase s, 1 mM thymidine was added to the dCK reaction and to inhibit the activity of dCK, 50 mM deoxycytidine was added to the dGK reaction. Total radioactivity was measured in a liquid scintillation counter (RackBeta, LKB Wallac, Turku, Finland).

Enzyme activities were expressed as pmol of formed product/min/mg of protein, and the method of Lowry was used to determine protein content in the individual samples using the DC protein assay (Bio-Rad Laboratories).

Western blot

Hemoglobin gamma protein was visualized using western blot. Criterion Precast 4-15% Tris-HCl gels were loaded with 40 g of protein and run in a sodium dodecyl sulfate (SDS) running buffer containing 0.4 M glycine, 50 mM Tris-HCl, and 14 mM SDS, pH 8.5 at 100 V for 2 h. The proteins were transferred onto an Immune Blot PVDF membrane in a transfer buffer (25 mM Tris-HCl, 192 mM glycine, 20% methanol, pH 8.3) at 50 V for 1 h. Mouse antihemoglobin gamma (1:200) antibody (Santa Cruz Biotechnology Inc, Santa Cruz, CA, USA) and rabbit anti--actin antibody (1:1000, Abcam, Cambridge, UK) were used. All secondary peroxidase-labeled antibodies (1:1000) were from Amersham Biosciences (Little Chalfont, UK). The blots were developed using Supersignal ® West Pico Chemiluminescent Substrate (Pierce Biotechnology, Rockford, IL, USA).

Microarray analysis

Microarray analyses were performed with RNA from four MOLT-4 samples and four MOLT-4/AraG900 samples using Human Genome Focus Array (Affymetrix, Buckinhamshire, UK) containing more than 8,500 of the most well-characterized human genes. Samples were prepared according to the manufacturer's instructions. Chip data were analyzed using Genespring GX 11.0 (Agilent Technologies, Inc, Santa Clara, CA, USA). Samples were subjected to baseline transformation and thereafter background correction, normalization, and probe summarization were done using the robust multichip averaging (RMA) algorithm [START_REF] Irizarry | Exploration, normalization, and summaries of high density oligonucleotide array probe level data[END_REF].

Significantly differentially expressed genes were found using Welch´s t test, followed by the Benjamin-Hochberg post hoc test. Samples were filtered on expression, p-value computation was asymptotic, and the cut-off limit was drawn at 0.05. Data were also analyzed by ALMAC Diagnostics (ALMAC Group Ltd, Graigavon, UK) with essentially the same genes found significantly affected.

Global methylation assay

Cells were incubated with AraG (5 or 10 M) or 6-MP (5 or 10 M) for 4 weeks (eight passages) and then pulse-labeled to measure the global methylation status. They were maintained at 0.5 10 6 cells/ml for MOLT-4 cells and approximately 70% confluency for A549 cells. Cells were washed in pre-heated PBS (37°C) in order to wash out all traces of drugs and then resuspended in pulse-labeling media containing 1.5 mmol/l L-[methyl-3 H]methionine (37 MBq, specific activity 75.6 Ci/mmol, Perkin Elmer, Waltham, MA, USA) as methyl donor, and 0.2 mol/l [methyl-14 C]thymidine (1.85 MBq, specific activity 61 mCi/mmol, Amersham Biosciences, Little Chalfont, UK) to be incorporated into newly synthesized DNA. After 4 h of incubation, the cells were harvested and DNA extracted using the QIAamp ® DNA Blood Mini Kit (Qiagen, Solna, Sweden) according to the manufacturer's instructions. A ratio between 3 H and 14 C was then used to estimate methylation and the results were compared to the results from cells cultured with no drug.

Statistical analysis

Additional statistical analyses were performed using the GraphPad Prism, version 4.00, for Windows (GraphPad Software, San Diego, CA, USA). Data were expressed as means and standard deviations and analyzed using Student´s t test with a p-value of 0.05 being considered significant.

Results

Cytotoxicity of nucleoside analogs and expression and activity of activating enzymes

We have previously described the MOLT-4/AraG900-resistant subclone and compared it to the parental MOLT-4 cell line [START_REF] Lotfi | Low level of mitochondrial deoxyguanosine kinase is the dominant factor in acquired resistance to 9-beta-Darabinofuranosylguanine cytotoxicity[END_REF]. We reported then on resistance to AraG but also crossresistance against other nucleoside analogs and to anthracyclines. Once again we investigated the pattern of nucleoside-analog cytotoxicity and expression and activity of nucleoside analog activating enzymes in these two cell lines. The AraG-resistant cells showed the highest resistance to AraG, Dnr, and AraC (350-, 1400-and 1500-fold less sensitivity respectively) compared to the parental cells. They also showed resistance to CdA, CAFdA, and FaraA, although moderately (10-70-fold less sensitivity), while the effect of dFdC was slightly increased in the resistant cells as was also shown before by us [START_REF] Lotfi | Low level of mitochondrial deoxyguanosine kinase is the dominant factor in acquired resistance to 9-beta-Darabinofuranosylguanine cytotoxicity[END_REF].

The resistant AraG900 subclone had reduced dCK expression and activity (less than half of control cells). However, the expression of dGK was elevated two-fold compared to parental cells while the activity of dGK measured with AraG as substrate was approximately 10% of the activity in parental cells.

Microarray expression data analysis

The data analysis revealed 1223 genes with an at least two-fold differential expression between the MOLT-4 and MOLT-4/AraG900 cells. On considering the most up-regulated genes, several fetal hemoglobin genes were represented, as well as the ABCB1 gene (Table 1) encoding for the membrane transporter P-gp involved in multidrug resistance [START_REF] Lofgren | Mechanisms of cross-resistance between nucleoside analogues and vincristine or daunorubicin in leukemic cells[END_REF]. Upregulated genes also included KLF1, involved in hemoglobin gamma expression [START_REF] Perrine | Erythroid Kruppel-like factor (EKLF) is recruited to the gamma-globin gene promoter as a co-activator and is required for gamma-globin gene induction by short-chain fatty acid derivatives[END_REF], and BLVRB, encoding biliverdin reductase B, which is involved in hemoglobin degradation. The genes BCL11A and BCL11B, probably involved in the regulation of fetal hemoglobin gene expression [START_REF] Sankaran | Human fetal hemoglobin expression is regulated by the developmental stage-specific repressor BCL11A[END_REF], were 10-and 12-fold down-regulated, respectively (Table 1). The expression of the nucleoside analog metabolizing enzyme, dCK, was moderately down-regulated, and the expression of dGK was 1.6-fold elevated and was not included in the list of significantly altered genes. The list of significantly altered genes did not reveal any changes in other genes traditionally involved in the resistance mechanisms against AraG or any other nucleoside analogs.

Hemoglobin induction and hypomethylation

In order to determine if AraG could induce fetal hemoglobin and ABCB1 expression in the short-term perspective, and if that was associated with global DNA hypomethylation, MOLT-4 cells and A549 cells were incubated for eight passages with AraG and 6-MP, the latter being known to hypomethylate DNA [START_REF] Lambooy | Inhibition of DNA methylation in malignant MOLT F4 lymphoblasts by 6mercaptopurine[END_REF]. Pulse-labeling of DNA was then performed to assess global methylation, and RNA was prepared to measure hemoglobin gamma expression. Since the MOLT-4 and MOLT-4/AraG900 cell lines had been cultured for a long period of time (up to 6 months during resistance development), new MOLT-4 cells (cultured less than 10 passages since received from ATCC) were used for these experiments. As shown in Table 2, hemoglobin gamma was induced in both MOLT-4 and A549 cells treated with AraG or 6-MP, and so was the ABCB1 gene in both cell lines treated with AraG. The same was seen in the MOLT-4/AraG900 cell line used in the microarray experiment and in a 6-MP-resistant variant generated previously by us (Table 2).

Global methylation was significantly decreased in both cell lines incubated with AraG, but only significantly decreased in the A549 cells when incubated with 6-MP (Fig. 1).

Western blot analysis of hemoglobin gamma

In the MOLT-4/AraG900 samples used in the gene expression microarray analysis, hemoglobin gamma protein expression was induced as well as in the 6-MP-resistant MOLT-4 cells (Fig. 2). We have previously shown that the P-gp protein expression is induced in these cells as well [START_REF] Lotfi | Low level of mitochondrial deoxyguanosine kinase is the dominant factor in acquired resistance to 9-beta-Darabinofuranosylguanine cytotoxicity[END_REF]. These alterations could not be detected in the MOLT-4 and A549 cells that had been incubated with drugs for a short period of time. This is probably due to the shorter induction time, and the fact that these proteins are not normally expressed in these cells, and so it may take a higher amount of protein in order for them to be detectable using these antibodies.

Altered ABCB1 and hemoglobin gamma expression after P-gp inhibition

Several cytotoxic drugs can induce P-gp expression even though they are not substrates of this efflux pump, and many drugs also induce fetal hemoglobin. To investigate whether there is a connection between P-pg and fetal hemoglobin induction, AraG-resistant MOLT-4 cells and K562 cells, the latter known to constitutively express both fetal hemoglobin [START_REF] Bianchi | Accumulation of gamma-globin mRNA and induction of erythroid differentiation after treatment of human leukaemic K562 cells with tallimustine[END_REF] and low amounts of P-gp [START_REF] Arora | Reversal of P-glycoprotein-mediated multidrug resistance by diallyl sulfide in K562 leukemic cells and in mouse liver[END_REF], were incubated with two inhibitors of this pump. The cells were incubated with cyclosporin A or laniquidar at different concentrations for 96 h in order to see if the P-gp efflux pump might be involved in the excretion of hemoglobin or its degradation products. Cells were then counted in an automated Beckman Coulter Counter (Beckman Coulter, Inc., Brea, CA, USA) and the expression of hemoglobin gamma and ABCB1 was measured.

All three concentrations tested reduced cell growth of the two cell lines (Fig. 3). RNA was extracted from the cells incubated with the highest concentration of the drugs (10 M) and analyzed for hemoglobin gamma and ABCB1 expression. We speculated that inhibiting P-gp would lead to toxic concentrations of hemoglobin or waste products of hemoglobin and slow cell proliferation and possibly affect the expression of these genes. On inhibiting P-gp, hemoglobin gamma expression decreased to 5 percent of control in MOLT/AraG900 cells and to 15-21 percent in K562 cells (Table 3), while the expression of ABCB1 was markedly increased in the MOLT-4/AraG900 cells, possibly as a compensation for the inhibition.

Discussion

In the present study, using gene expression microarray analysis, we have demonstrated that fetal hemoglobin genes as well as the drug efflux pump, ABCB1, involved in multidrug resistance are up-regulated in the AraG-resistant variant of the human leukemic MOLT-4 cell line. We expected major changes in genes involved in nucleoside analog resistance, such as nucleoside analog-metabolizing enzymes, and hoped to obtain information about additional genes involved in the mechanisms of resistance. However, genes known to be involved in resistance development were only moderately affected together with changes in more than 1000 additional genes. The reason for these unexpected results it not entirely clear. It may be due to the fact that the cells were selected for resistance during several months and that the mechanisms for resistance change over time. It may also be because the concentrations used for resistance induction are of such a nature that low concentrations induce certain resistance mechanisms and higher concentrations induce other mechanisms. Also, the resistance mechanisms may not be reflected by the mRNA expression of the genes. The fact that the mRNA expression of dGK, which is known to phosphorylate and activate AraG, was, in fact, slightly induced while the enzymatic activity was 1/10 of that in control cells suggests posttranslational changes not seen at the mRNA level.

Fetal hemoglobin [START_REF] Enver | Developmental programs of human erythroleukemia cells: globin gene expression and methylation[END_REF] and ABCB1 gene induction [START_REF] Kantharidis | Altered methylation of the human MDR1 promoter is associated with acquired multidrug resistance[END_REF][START_REF] Nakayama | Hypomethylation status of CpG sites at the promoter region and overexpression of the human MDR1 gene in acute myeloid leukemias[END_REF] is correlated with hypomethylation of DNA, and we speculated that AraG has hypomethylating properties since these genes were induced in AraG-resistant cells. This was investigated in MOLT-4 cells and the lung adenocarcinoma cell line A549 by means of a four-week drug induction with AraG and the anti-metabolite 6-MP, a drug known to hypomethylate DNA [START_REF] Lambooy | Inhibition of DNA methylation in malignant MOLT F4 lymphoblasts by 6mercaptopurine[END_REF] and, like AraG, is incorporated in DNA. We used the A549 as a control cell line since many nucleoside analogs, and particularly AraG, have a very modest effect on solid tumors, and we wanted to see if the effects of AraG were applicable to other cells than leukemic. Fetal hemoglobin gamma was induced in both cell lines, while ABCB1 was induced in cells treated with AraG only. The global methylation of DNA was also significantly decreased in both cell lines treated with AraG. Drug concentrations that inhibited approximately 50% cell growth in the cultured cells were used since we wanted to mimic the development of resistance when some cells more tolerant to the treatment survive and expand into a less drug-sensitive population. The effect on DNA methylation may depend on drug concentrations, exposure times, and cell lines used.

The gamma-globin gene is silenced through methylation shortly after birth, a process known as the fetal switch, but can easily be re-expressed after treatment with demethylating agents [START_REF] Saunthararajah | Effects of 5-aza-2'-deoxycytidine on fetal hemoglobin levels, red cell adhesion, and hematopoietic differentiation in patients with sickle cell disease[END_REF][START_REF] Koshy | 2-deoxy 5azacytidine and fetal hemoglobin induction in sickle cell anemia[END_REF]. Substances shown to induce fetal hemoglobin in cell lines include sodium buturate, cisplatin, imatinib, cytosine arabinoside, and hydroxyurea [START_REF] Bianchi | Accumulation of gamma-globin mRNA and induction of erythroid differentiation after treatment of human leukaemic K562 cells with tallimustine[END_REF][START_REF] Jacquel | Imatinib induces mitochondria-dependent apoptosis of the Bcr-Abl-positive K562 cell line and its differentiation toward the erythroid lineage[END_REF][START_REF] Huang | CPEC induces erythroid differentiation of human myeloid leukemia K562 cells through CTP depletion and p38 MAP kinase[END_REF][START_REF] Bianchi | Induction of erythroid differentiation of human K562 cells by cisplatin analogs[END_REF][START_REF] Cortesi | In vitro effect on human leukemic K562 cells of co-administration of liposome-associated retinoids and cytosine arabinoside (Ara-C)[END_REF][START_REF] Park | Involvement of p38 kinase in hydroxyurea-induced differentiation of K562 cells[END_REF][START_REF] Witt | Induction of fetal hemoglobin expression by the histone deacetylase inhibitor apicidin[END_REF][START_REF] Witt | Butyrate-induced erythroid differentiation of human K562 leukemia cells involves inhibition of ERK and activation of p38 MAP kinase pathways[END_REF], 5-azacytidine, and decitabine [START_REF] Charache | Treatment of sickle cell anemia with 5-azacytidine results in increased fetal hemoglobin production and is associated with nonrandom hypomethylation of DNA around the gamma-deltabeta-globin gene complex[END_REF][START_REF] Lowrey | Brief report: treatment with azacitidine of patients with end-stage beta-thalassemia[END_REF][START_REF] Ley | 5-Azacytidine increases gamma-globin synthesis and reduces the proportion of dense cells in patients with sickle cell anemia[END_REF][START_REF] Desimone | 5-Azacytidine stimulates fetal hemoglobin synthesis in anemic baboons[END_REF]. Two genes, BCL11A and BCL11B, working as transcriptional repressors of the fetal hemoglobin genes [START_REF] Sankaran | Human fetal hemoglobin expression is regulated by the developmental stage-specific repressor BCL11A[END_REF][START_REF] Miccio | The role of the GATA-1/FOG-1/NuRD pathway in the expression of human {beta}-like globin genes[END_REF] were down-regulated in the resistant cells. Therefore, transcriptional repression of these two genes may also be one explanation for the massive fetal hemoglobin induction.

The expression of the nucleoside analog-activating enzymes dCK and dGK was only moderately altered according to the microarray analysis and real-time PCR measurements, while enzymatic activities were reduced. However, we have previously demonstrated a lack of correlation between mRNA expression and the enzymatic activity of these enzymes [START_REF] Lotfi | The pattern of deoxycytidine-and deoxyguanosine kinase activity in relation to messenger RNA expression in blood cells from untreated patients with B-cell chronic lymphocytic leukemia[END_REF].

The multidrug resistance gene, ABCB1, was one of the most up-regulated genes. We have shown before, that this gene is induced in MOLT-4 cells selected for resistance to AraG [START_REF] Lotfi | Low level of mitochondrial deoxyguanosine kinase is the dominant factor in acquired resistance to 9-beta-Darabinofuranosylguanine cytotoxicity[END_REF] and is probably responsible for the cross-resistance to daunorubicin seen in these cells. We have also studied the K562 cell line and a variant resistant to vincristine, a substrate of P-gp [START_REF] Green | Spontaneous reversal of p-glycoprotein expression in multidrug resistant cell lines[END_REF]. The presence of ABCB1 and P-gp expression after drug withdrawal was evaluated over 11 months. These levels remained unchanged for 5 months without vincristine, but then rapidly decreased, at which time cells also became more sensitive to P-gp substrates [START_REF] Green | Spontaneous reversal of p-glycoprotein expression in multidrug resistant cell lines[END_REF].

However, these cells remained highly resistant to AraG. It was shown that after 11 months of drug withdrawal, some cells still contained low levels of ABCB1/P-gp which were easily reinduced [START_REF] Green | Spontaneous reversal of p-glycoprotein expression in multidrug resistant cell lines[END_REF] and we showed that AraG had the ability to re-induce P-gp during short-term incubations. Like other nucleoside analogs, AraG is not considered to be a substrate of P-gp [START_REF] Litman | From MDR to MXR: new understanding of multidrug resistance systems, their properties and clinical significance[END_REF], but can yet induce P-gp expression. Prenkert et al. showed that a 10 min incubation with 0.5 M of the nucleoside analog cytarabine, induced a 1.7-fold increase in ABCB1 expression in human HL60 cells, while, for daunorubicin, which is a P-gp substrate [START_REF] Litman | From MDR to MXR: new understanding of multidrug resistance systems, their properties and clinical significance[END_REF], a 1.6-fold increased ABCB1 induction was not seen until after 24 h [START_REF] Prenkert | Rapid induction of P-glycoprotein mRNA and protein expression by cytarabine in HL-60 cells[END_REF]. It has previously been shown that HL60 cells resistant to the anthracycline doxorubicin had increased P-pg expression and showed cross-resistance to AraC, and upon inhibition of P-gp with cyclosporine A the intracellular retention of AraC was increased [START_REF] Mansson | Cross-resistance to cytosine arabinoside in a multidrug-resistant human promyelocytic cell line selected for resistance to doxorubicin: implications for combination chemotherapy[END_REF]. The authors therefore speculated that resistance to AraC in these cells was due both to decreased AraC activation by dCK which was slightly reduced, and efflux of AraC by P-gp. This may be the reason for AraG, and nucleoside analog-resistance in general, in the present study as well. Studies have also shown that AraC can activate the protein kinase C (PKC) signaling pathway [START_REF] Kharbanda | Regulation of c-jun gene expression in HL-60 leukemia cells by 1-beta-D-arabinofuranosylcytosine. Potential involvement of a protein kinase C dependent mechanism[END_REF][START_REF] Kharbanda | 1-beta-Darabinofuranosylcytosine activates serine/threonine protein kinases and c-jun gene expression in phorbol ester-resistant myeloid leukemia cells[END_REF] and the mitogen-activated protein (MAP) kinase p38, [START_REF] Kharbanda | 1-beta-Darabinofuranosylcytosine activates serine/threonine protein kinases and c-jun gene expression in phorbol ester-resistant myeloid leukemia cells[END_REF] eventually leading to activation of AP-1 which is the transcription factor for the ABCB1 gene and the hemoglobin gene locus [START_REF] Ikuta | Globin gene switching. In vivo protein-DNA interactions of the human beta-globin locus in erythroid cells expressing the fetal or the adult globin gene program[END_REF][START_REF] Guo | Increased p38-MAPK is responsible for chemotherapy resistance in human gastric cancer cells[END_REF].

Activation of p38 kinase is also activated by the hypomethylating cytosine analog decitabine, [START_REF] Lavelle | Decitabine induces cell cycle arrest at the G1 phase via p21(WAF1) and the G2/M phase via the p38 MAP kinase pathway[END_REF] and this intracellular pathway seems to be a common mechanism for hypomethylating drugs, and possibly shared by AraG as well.

It is tempting to speculate that ABCB1 expression is induced in response to AraG and other hypomethylating drugs due to the increased hemoglobin concentration in the cells since when hemoglobin is degraded and excreted, different P-gp transporters are expected to be involved in its transportation [START_REF] Beukeveld | Relationship between biliary lipid and protoporphyrin secretion; potential role of mdr2 P-glycoprotein in hepatobiliary organic anion transport[END_REF] and also since inhibitors of P-gp increase bilirubin entry into the rat brain [START_REF] Hanko | Administration of drugs known to inhibit P-glycoprotein increases brain bilirubin and alters the regional distribution of bilirubin in rat brain[END_REF]. In that case, this would explain why nucleoside analog-resistant and other resistant cells show a multidrug-resistant phenotype over-expressing P-gp. On analyzing the gene expression microarray analysis data, we also found up-regulated genes involved in hemoglobin synthesis and degradation. The ABCB6 gene encodes a mitochondrial transporter required for mitochondrial porphyrin uptake necessary for hemoglobin synthesis [START_REF] Krishnamurthy | Identification of a mammalian mitochondrial porphyrin transporter[END_REF], and this gene was slightly up-regulated, as well as ABCG2 (not shown) also involved in this transport [START_REF] Tamura | Functional validation of the genetic polymorphisms of human ATP-binding cassette (ABC) transporter ABCG2: identification of alleles that are defective in porphyrin transport[END_REF]. On inhibiting P-gp in the MOLT-4/AraG900 and the K562 cells, the expression of hemoglobin gamma was decreased, and in the MOLT-4/AraG900 cells the expression of ABCB1 was greatly induced. K562 cells constitutively express low levels of P-gp [START_REF] Arora | Reversal of P-glycoprotein-mediated multidrug resistance by diallyl sulfide in K562 leukemic cells and in mouse liver[END_REF][START_REF] Yague | Pglycoprotein (MDR1) expression in leukemic cells is regulated at two distinct steps, mRNA stabilization and translational initiation[END_REF] and fetal hemoglobin, [START_REF] Arora | Reversal of P-glycoprotein-mediated multidrug resistance by diallyl sulfide in K562 leukemic cells and in mouse liver[END_REF] and this is probably the reason for the modest effect on ABCB1 and hemoglobin gamma expressions after P-gp inhibition in this cell line. We speculate that this may occur due to increased retention of hemoglobin and its degradation products in the cells when inhibiting P-gp, with decreased hemoglobin production and increased P-gp expression in order to try to compensate for the decreased efflux of hemoglobin or hemoglobin waste products. The P-gp inhibitors slightly affected the cell growth in the parental MOLT-4 cells at the highest concentrations used, but since these cells do not express P-gp this effect must be due to some other actions of the drugs or possibly off-target effects due to the high concentration used.

One study showed that short-term treatment (96 h) with different cytotoxic drugs induced ABCB1 mRNA expression in K562 cells, but that this was due to mRNA stabilization which did not lead to P-gp expression at the membrane surface, but that long-term selection for drug resistance did lead to mRNA association with polyribosomes and translation of P-pg [START_REF] Yague | Pglycoprotein (MDR1) expression in leukemic cells is regulated at two distinct steps, mRNA stabilization and translational initiation[END_REF]. We have previously characterized the MOLT-4 and the MOLT-4/AraG900 cells with regard to Pgp expression. These experiments demonstrated an absence of P-gp in MOLT-4 cells, but a high expression in the AraG-resistant cells measured by western blot, [START_REF] Lotfi | Low level of mitochondrial deoxyguanosine kinase is the dominant factor in acquired resistance to 9-beta-Darabinofuranosylguanine cytotoxicity[END_REF] and we therefore assume that the massive up-regulation of ABCB1 seen in this study is accompanied by an increased P-pg expression. However, the fact that AraG induced ABCB1 expression in both MOLT-4 and A549 cells during short-term treatment, along with hemoglobin gamma expression suggests that AraG and other DNA-hypomethylating drugs should be used with caution together with P-gp substrate drugs due to the ABCB1 induction.

Leukemia arises partly due to inactivation of genes controlling differentiation and growth [START_REF] Lotem | Different blocks in the differentiation of myeloid leukemic cells[END_REF], and this may be a consequence of epigenetic silencing of the genes. These genes may become reactivated by such inhibitors of DNA methylation as 5-aza-2´-deoxycytidine or, possibly, AraG. It would be interesting to find out what genes are hypomethylated and reexpressed in patients after short-term AraG treatment and whether these genes could be used for targeted therapies, or if they only render the patients less susceptible to therapy by inducing genes involved in nucleoside analog resistance.

3 . 1 (

 31 Applied Biosystems). The reaction mixture contained 1X Taqman Universal Mastermix, forward and reverse primers (800 nM), and a probe (200 nM) or 1X Gene Expression Assay mix. Reaction conditions were: 2 min at 50 o C, 10 min at 95 o C, and 40 cycles with 15 s at 95 o C and 1 min at 60 o C in an ABI Prism 7500 Sequence Detection System.

3 H

 3 ]-2´-deoxycytidine, specific activity 8.3 Ci/mmol) for dCK and 20 M 9-β-Darabinofuranosylguanine ([8-3 H]-guanine-β-D-arabinofuranoside, specific activity 6.7

Figure 1

 1 Figure 1 Global methylation status was investigated in MOLT-4 cells and A549 cells cultured without drug and with AraG or 6-MP. Methylation was assessed using a DNA pulse-labeling technique with thymidine and methionine measuring methylated cytosine in newly synthesized DNA.

Figure 2

 2 Figure 2 Protein expression of hemoglobin gamma (Hb ) was investigated in the MOLT-4 and MOLT-4 AraG and 6-MP-resistant cells using Western blot. Hemoglobin gamma protein was highly induced in the resistant cells and compared to β-actin expression.

Figure 3

 3 Figure 3 The resistant MOLT-4/AraG900 and K562 cells were incubated during 96 h with the P-gp inhibitors cyclosporin A (cyc A) and laniquidar (laniq) at different concentrations (1, 5, 10 M). Cells were then counted and a graph was plotted showing growth inhibition of both drugs at all concentrations.Y-axis showing amount of million cells/ml.

Acknowledgments

This work was supported by grants from the Swedish Cancer Foundation, the Swedish Childhood Cancer Foundation, the Signe and Olof Wallentin Foundation, Capio's Research Foundation, the County Council of Östergötland, and the Swedish Fund for Research without Animal Experiments.


T Bastogne

R Keinj

P Vallois

A stochastic model of the cell and tumor lifespans in radiotherapy

Keywords: dynamic models, stochastic models, multinomial distribution, cell lifespan, target theory, cancer

, a multinomial model based on a discrete-time Markov chain able to take into account both cell repair and cell damage heterogeneity due to radiations has been proposed. In this paper, we introduce the notion of lifespan for a single cell, a tumor and a normal tissue. We also determine the cumulative distribution functions of these random variables. These results lead to original formulations of the tumor control probability (TCP) and the normal tissue complication probability (NTCP). Finally, we propose a new characteristic: the Efficiency-Complication Trade-off diagram which provides to the radiotherapeutist all the information needed to choose the most appropriate treatment plan, i.e. the suited number of dose fractions to be applied to control the tumor while preserving the adjacent normal tissue.

INTRODUCTION

Cancer is a disease that affects millions of people worldwide. One of the common therapies used to treat cancer is external beam radiotherapy. In this treatment, ionizing radiations are used to eliminate tumors. The ionizations induced by radiations cause a variety of possible lesions in cells, [START_REF] Curtis | Lethal and potentially lethal lesions induced by radiation -a unified repair model[END_REF]. The most harmful damages are the lesions which affect the DNA structure, [START_REF] Feinendegen | Systems biology and its potential role in radiobiology[END_REF]. Mathematical modelling may help to quantify the effects of the radioactive treatments on cell populations. It can be used to predict tumor growth and cancer spread, but also allows to determine the effectiveness of a specific treatment. Two probabilities are generally involved in the design of a treatment plan in radiotherapy : (i) the tumor control probability (TCP) and (ii) the normal tissue complication probability (NTCP). The TCP is the probability that all cancer cells are dead in the irradiated region, see [START_REF] Zaider | Tumour control probability: a formulation applicable to any temporal protocol of dose delivery[END_REF]; [START_REF] Dawson | Derivation of the tumour control probability (tcp) from a cell cycle model[END_REF]. The NTCP is another probability that measures the sensitivity of normal tissue to radiations. The optimum choice of radiation dose delivery technique in the treatment of a given tumor is those that maximizes the TCP while keeping the NTCP lower than a given tolerance threshold. A classical method to obtain TCP, and probably the most used one, is to evaluate the probability that the surviving cell number, in the tumor, equals to zero. The expression of TCP and NTCP requires stochastic models and we are interested here in the simplest class of models used in radiotherapy: the hit models associated with the target theory. Note that other model structures, e.g. survival curves in [START_REF] Fowler | The linear-quadratic formula and progress in fractionated radiotherapy[END_REF], stochastic cell population-dynamic models in [START_REF] Sachs | Simple ODE models of tumor growth and anti-angiogenic or radiation treatment[END_REF]or the cell-cycle models, in [START_REF] Kirkby | Mathematical modelling of the response of tumour cells to radiotherapy[END_REF], are also used to determine TCP and NTCP.

The target theory and hit-modeling paradigm were introduced in the 1920s when biologists were beginning to apply quantum physics to biology, see [START_REF] Dessauer | Uber einige Wirkungen von Strahlen[END_REF]; [START_REF] Blau | Uber einige Wirkungen von Strahlen[END_REF]. The modeling of radiation effects on living cells was initiated theoretically and experimentally by [START_REF] Lea | The mechanism of induction by radiation of chromosome aberrations in tradescantia[END_REF], and continued by [START_REF] Atwood | On the interpretation of multi-hit survival curves[END_REF][START_REF] Atwood | On the interpretation of multi-hit survival curves[END_REF][START_REF] Pollard | The direct action of ionizing radiation on enzymes and antigen[END_REF]. In this theoretical setting, it is assumed that specific regions of the DNA called targets must be all inactivated to kill the cell. Each target is deactivated when it is hit by a number of radiation particules. Several classes of hit models classified by the number of targets and the number of hits were proposed, see [START_REF] Satow | Hit and target models for dna damage with indirect action[END_REF]; [START_REF] Chapman | Target theory revisited: Why physicists are essential for radiobiology research[END_REF]; [START_REF] Ditlov | Track theory and radiation effects[END_REF].

Unfortunately two main aspects of tumor growth are often missing in target and hit models: (i) the target reparation between two consecutive dose fractions and (ii) the heterogeneity of damages induced by radiations in the cancer cell population after each dose fraction. In [START_REF] Keinj | A multinomial model of tumour growth treated by radiotherapy[END_REF], an extension of hit models is proposed through a multinomial model that takes into account the different states of a surviving cell corresponding to the degree of DNA damage. Moreover, the proposed model is able to describe both target reparation between two consecutive dose fractions of the radiation schedule and the heterogeneity of damages induced by radiations in the cancer cell population after each dose fraction. In this multinomial model, a cell is supposed to contain m targets which must be all deactivated to produce cell death. The cell population is then split up into m + 1 categories. Like the majority of models used to measure the cell mortality, the multinomial model examines the number of surviving cells in a population. In this paper, we are interested in developing this model for the following objectives:

• to introduce the lifespan of a single cell, a tumor and a normal tissue and to study the randomness of these lifespans by determining their mean values, variances and their cumulative distribution functions; • to use the previous probabilistic features to provide original formulations of TCP and NTCP. A concrete application of these results is to estimate the suitable dose to both control the tumor and limit the damage in healthy tissues.

In Section 2 we give a reminder on the tumor growth multinomial model. The modeling of the cancer cell lifespan T is developed in Section 3. We then introduce the tumor and normal tissue lifespans L n 0 and L n, n0 in Section 4 and we give new expressions of TCP and NTCP. Before concluding, we also propose a new diagram that summarizes all the information needed to choose the appropriate treatment schedule.

MULTINOMIAL MODEL

In [START_REF] Keinj | A multinomial model of tumour growth treated by radiotherapy[END_REF], we have proposed a multinomial model of tumor growth relying on the target and hit modeling paradigm. This model is based on the following assumptions :

• a cell has m targets;

• each target may be made inactive by a single hit caused by a radiation particule; • the cell death happens when the m targets are deactivated;

• between two consecutive dose fractions, if the cell is still alive then a target repair mechanism occurs; • fractionated radiation doses are delivered at times k = 1, 2, 3, • • • , k is the number of treatment days.

Heterogeneity of cell states after radiations

We have m + 1 possible states for a cell:

• state i, the cell has i inactive targets, i ∈ {0, 1, . . . , m -1}, they are the m states of a surviving cell; • state m, the cell having m inactive targets is a dead cell.

Figure 1 shows the case of a 3-target cell and the corresponding cell states. Let Z k the random variable denoting the state of

Z k 0 1 2 3
Fig.

Π = PR,
(1) where P models the effects of dose fractions, and R describes repair mechanisms. In this paragraph, the first row and first column of a matrix will be noted by the index value 0.

Treatment effect modelling

Let P(i, j) the probability to deactivate j targets at time k + 1 when i targets are disabled at time k. We assume that

• an active target becomes inactive with probability q;

• a deactivated target remains inactive with probability 1;

• the disabling of targets in the cell are independent events. Therefore, if i targets are inactive then after applying a new fraction dose we obtain j inactive targets with probability

P(i, j) = ( m-i j-i )q j-i (1 -q) m-j i ≤ j 0 j < i.
(2)

Cellular reparation modelling

We introduce now repair mechanisms of deactivated targets which occur between two consecutive dose fractions. We suppose that

• a deactivated target in an alive cell can be repaired with probability r; • an active target remains active with probability 1;

• the behaviors of the targets are independent.

Therefore, a cell in state i < m at time k may switch at time k +1 to the state j ≤ i if ij targets among the i inactive targets are repaired and with probability

R(i, j) = ( i j )r i-j (1 -r) j j ≤ i < m 0 i < j. (3)
When i = m, i.e. the cell is dead, the only possible case is j = m, then R(m, m) = 1.

Probability distribution of Z k

Set

ν k = (ν 0 k , . . . , ν m k ) the probability distribution vector of Z k with ν i k = Pr(Z k = i).
From the transition matrix Π and the initial probability distribution ν 0 of Z 0 , we can determine ν k , for all k ∈ N * , using the Markov chain property,

ν k = ν 0 Π k .
(4) If we assume that all cells are initially in state i 0 (i 0 inactive targets), then Pr(Z 0 = i 0 ) = 1 and ν 0 = (ν 0 0 , . . . , ν i 0 0 , . . . , ν m 0 ) = (0, . . . , 1, . . . , 0). Therefore, we obtain

ν i k = Π k (i 0 , i) i ∈ {0, . . . , m}.
(5)

Multinomial model of the tumor

Consider a group of n 0 cells that behave independently and with the same dynamic. We suppose that each cell has m targets and admits alternative damage and repair phases as explained at the begining of this section. For i ∈ {0, . . . , m}, let the random variable X k (i) denote the number of cells in state i, among the n 0 initial cells, at time k:

X k (i) = n 0 ∑ l=1 1 {Z (l) k =i} , i ∈ {0, 1, • • • , m}, (6)
where

Z (l)
k is the number of inactive target(s) in the l th cell. From our assumptions, Z

(1) k , . . . , Z (n 0 ) k are i.i.d. and distributed as Z k . Consequently, the random vector

(X k (0), • • • , X k (m))
follows a multinomial distribution with parameters n 0 and ν k = (ν 0 k , . . . , ν m k ) given in (5). Therefore, Red circles show the 99.9% confidence interval on the population size.

P X k (0) = a 0 , • • • , X k (m) = a m = n 0 ! a 0 ! • • • a m ! ν 0 k a 1 . • • • . ν m k a m (7) where a 0 , • • • , a m ∈ {1, • • • , n 0 } and a 0 + • • • + a m = n 0 . The number X k (i) of cells in state i at time k follows the binomial distribution with parameters n 0 and ν i k = Π k (i 0 , i), X k (i) ∼ B n 0 , Π k (i 0 , i) . ( 8
! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! !

LIFESPAN T OF A CANCER CELL

Cell lifespan: definition and properties

In this section, we develop an original approach based on the lifespan of a cancer cell which has received a sequence of successive dose fractions of radiation. The cell lifespan is the first (random) instant of treatment time when the cell is dead. We would like to emphasize this new viewpoint of cell mortality modeling. Indeed, in the usual hit models and in the previous multinomial model, k is fixed and we only consider random variables at time k, e.g. the number Z k of deactivated targets. In the multinomial setting it has been proved that the distribution of Z k is given by Π k . Now, we consider a random variable T denoting the lifespan of a single cell composed of m targets. In other terms, T is the first instant k when the cell is in state m (death state) :

T = in f {i ≥ 1; Z i = m}. (9)
The knowledge of T allows to recover information at time k for instance T ≤ k is equivalent to Z k = m. Moreover T reveals the dynamical aspect of the cell lifespan, so

{T = k} = {Z 1 = m, • • • , Z k-1 = m, Z k = m} and therefore involves the sequence Z 1 , • • • , Z k instead of Z k only.
T gives naturally rise to Monte Carlo simulations.

The counterpart is that we have to deal with a random element. In Proposition 3.1 below, we calculate its mean time and the associated variance. The second parameter permits to control the variability of T around its mean value. This rough information can be completed by the cumulative distribution which highlights how T fluctuates. Proposition 3.1. If the cell is initially in state Z 0 = i 0 ∈ {0, . . . , m -1}, the random variable T of the cancer cell lifespan verifies the following properties :

(1) The cumulative distribution function F of T is related to the transition matrix Π:

F(k) = P(T ≤ k) = Π k (i 0 , m).
(2) For all n ∈ N * , the expected value E(T n ) is finite and

E(T ) = ∑ k≥1 k Π k (i 0 , m) -Π k-1 (i 0 , m) .
(3) The variance V (T ) is finite and given by

V (T ) = ∑ k≥1 k 2 Π k (i 0 , m) -Π k-1 (i 0 , m) -E(T ) 2 .
Proof 1. According to the definition of T (cf (9)) and the crucial fact that m is an absorbant state we have

P(T ≤ k) = P(Z k = m | Z 0 = i 0 ) ∀ k ∈ N * . ( 10
)
According to (5), we obtain

F(k) = P(T ≤ k) = Π k (i 0 , m). (11) 2.
Set

λ 1 = max{1 -Π(i, m), 0 ≤ i ≤ m -1} (12) = max{ ∑ j≤m-1 Π(i, j), 0 ≤ i ≤ m -1} (13)
λ 2 = max{Π(i, m), 0 ≤ i ≤ m -1} (14) λ = max(λ 1 , λ 2 ). ( 15
)
Note that, λ 1 < 1 and λ 2 < 1 then λ ∈]0, 1[. We prove that

P(T = k) ≤ λ 2 λ k-1 1 ≤ λ k
only when k = 3. The general proof being similar:

P(T = 3) = P(Z 1 = m, Z 2 = m, Z 3 = m) = ∑ i 1 ,i 2 ≤m-1 Π(i 0 , i 1 )Π(i 1 , i 2 )Π(i 2 , m) ≤ λ 2 ∑ i 1 ≤m-1 ∑ i 2 ≤m-1 Π(i 1 , i 2 ) Π(i 0 , i 1 ) ≤ λ 2 λ 1 ∑ i 1 ≤m-1 Π(i 0 , i 1 ) ≤ λ 2 λ 2 1 ≤ λ 3 .
As λ ∈]0, 1[, we deduce that the mean value E(T n ) < ∞, n ≥ 1:

E(T n ) = ∑ k≥1 k n P(T = k) (16) ≤ ∑ k≥1 k n λ k < ∞.
For n = 1, E(T n ) = E(T ) < ∞. Now, we express E(T ) in terms of the cumulative distribution function F. Thus,

E(T ) = ∑ k≥1 k P(T < k + 1) -P(T < k) (17) = ∑ k≥1 k F(k) -F(k -1) . = ∑ k≥1 k Π k (i 0 , m) -Π k-1 (i 0 , m) .
3. For n = 2, E(T 2 ) < ∞ and the variance of T is finite as it given by:

V (T ) = E(T 2 ) -E(T ) 2 (18) = ∑ k≥1 k 2 P(T < k + 1) -P(T < k) -E(T ) 2 = ∑ k≥1 k 2 Π k (i 0 , m) -Π k-1 (i 0 , m) -E(T ) 2 .

Confidence interval of T

We are interested in a most probable interval I i.e. I such that

P(T ∈ I) ≥ 1 -α, (19)
where 1α is the confidence level. We restrict ourselves to interval I of the type

I =]E(T ) -ε 1 , E(T ) + ε 2 ].
The positive real numbers ε 1 and ε 2 are determined by the conditions :

P T ≤ E(T ) -ε 1 ≤ β 1 and P T > E(T ) + ε 2 ≤ β 2 , (20)
where β 1 + β 2 = α. For instance, α = 0.05, β 1 = β 2 = α 2 = 0.025. The confidence interval I represents values of T for which the difference between T and E(T ) is not statistically significant at the α level. F being the cumulative distribution function of T , it is clear that ( 20) is equivalent to

F E(T ) -ε 1 ≤ β 1 and F E(T ) + ε 2 ≥ 1 -β 2 . (21)

Numerical analysis

We present here a set of numerical experiments to show the influence of the parameters q, r and the initial state Z 0 = i 0 , on the interval I in the case β 1 = β 2 = α 2 = 0.025. We fix two parameters and we vary the third one. We implemented the cumulative distribution function F on Matlab with m = 3. Table 1, figures 3 and 4 present all the numerical results.

We can observe and quantify numerically the monotonicity of the confidence interval I, specially when the efficiency of the treatment or the repair mechanisms are changed. It also appears that the interval I is smaller when the initial state of the cell is fragile and also when the treatment is more efficient. Conversely, the interval I becomes larger when the cellular repair mechanisms are stronger.

ORIGINAL FORMULATIONS OF TCP AND NTCP

The TCP is the probability that all cancer cells are dead in the irradiated region and the NTCP measures the sensitivity

Fixed parameters

Varied parameter E(T ) I q = 0.4

Z 0 = 0 5.76 [1, 17] r = 0.3 Z 0 = 2 3.63 [1, 14] r = 0.3 q = 0.6 2.90 [1, 7] Z 0 = 0 q = 0.8 1.70 [1, 3] q = 0.6 r = 0.2 2.75 [1, 6] Z 0 = 0 r = 0.6 3.47 [1, 10]
Table 1. Influence of the multinomial parameters q, r and i 0 on the confidence interval I of T with m = 3, n 0 = 10 4 and β 1 = β 2 = 0.025.

of the adjacent normal tissue for a given radiation treatment schedule. These probabilities are simultaneously used to optimize the quality of a treatment planning. The treatment goal is to maximizes the TCP and to keep the NTCP lower than a given tolerance threshold. In this section, we give an original approach based on the cell lifespan to choose the optimal dose fraction number for a treatment plan and we give new expressions of TCP and NTCP.

Tumor Control Probability

The tumor control probability at time k is the probability to destroy the tumor at this time. If L n 0 is the lifespan of the tumor, then TCP k is the probability that L n 0 is lower than k, Table 2. Influence of the initial cancer cell number n 0 on the expected value and confidence interval J of L n 0 .

TCP k = P(L n 0 ≤ k). ( 22
) n 0 E(L n 0 ) J n 0 =
In order to find an explicit expression of TCP k , we express L n 0 in terms of the single cell lifespans.

Lifespan L n 0 of the tumor. We consider a tumor initially composed of n 0 living cells. For each single cell j we associate the corresponding lifespan:

T ( j) = in f {i ≥ 1; Z ( j) i = m}.
(23) L n 0 corresponds to the time when the whole tumor is eliminated, i.e. all the cells are killed. L n 0 can then be expressed with respect to the variables {T ( j) } j∈{1,...,n 0 } as L n 0 = max{T (1) , . . . , T (n 0 ) }.

(24) Consequently,

P(L n 0 ≤ k) = P(T (1) ≤ k, . . . , T (n 0 ) ≤ k).
(25) Under the assumption that all cells behave independently and with the same dynamic, T (1) , . . ., T (n 0 ) are then i.i.d. random variables and have the same distribution as T (cf ( 11)). Therefore, we obtain

P(L n 0 ≤ k) = P(T (1) ≤ k) × . . . × P(T (n 0 ) ≤ k) (26) = P(T ≤ k) n 0 . (27)
According to Proposition 3.1, equation ( 17) and the assumption that all cells are initially in state i 0 , we have

TCP k = Π k (i 0 , m) n 0 (28)
Confidence interval of L n 0 . Let J the confidence interval of L n 0 of the type J =]E(L n 0 )κ 1 , E(L n 0 ) + κ 2 ] and we determine the positive real numbers κ 1 and κ 2 by the conditions :

P L n 0 ≤ E(L n 0 ) -κ 1 ≤ 0.025 (29) and P L n 0 > E(L n 0 ) + κ 2 ≤ 0.025. (30) These conditions are equivalent to G E(L n 0 ) -κ 1 ≤ 0.025 and G E(L n 0 ) + κ 2 ≥ 0.975, ( 31
) where G(•) denotes the cumulative distribution function of L n 0 , i.e. G(k) := P(L n 0 ≤ k).

Numerical analysis

We implemented the cumulative distribution function G and E(L n 0 ) in Matlab. Numerical results are given in table 4.1.2. The mean lifespan of the tumor increases with n 0 while the width of the confidence interval remains almost constant.

Normal Tissue Complication Probability

Normal Tissue Complication Probability (NTCP) is the probability that a complication appears in the adjacent normal tissue. A complication occurs when the number of dead cells in the normal tissue is larger than a given threshold number n.

We suppose that the normal cells behave similarly than the cancer cells, the difference coming from different values of the model parameters. A normal cell contains m targets, each of them is deactivated with a probability q and it is repaired with a probability r. For simplicity, we use the same notation, as for cancer cells, to denote the probability Πk ( ī0 , m) that a normal cell is in state m at time k (dead cell).

Lifespan T of a normal cell. We consider a single normal cell that contains m targets. Let Zk be the random variable denoting the state of the normal cell at time k and T its lifespan. Similary to the case of a cancer cell, T is defined as T = in f {i ≥ 1; Zi = m}, (32) and the corresponding cumulative distribution function is

F(k) = P( T ≤ k) = Πk ( ī0 , m). ( 33
)
Lifespan L n, n0 of the normal tissue. Consider now a normal tissue containing initially n0 cells and we define L n, n0 the random variable denoting the lifespan of the entire normal tissue. For each cell j we associate the corresponding lifespan, (35) where T (l: n0 ) denotes the lifespan of l normal cells among n0 . This sequence is called the order statistics of T (1) , T (2) , . . . , T ( n0 ) . T (l: n0 ) is called the l th order statistic.

T ( j) = in f {i ≥ 1; Z( j) i = m}. ( 34
Proposition 4.1. The normal tissue lifespan L n, n0 is the nth order statistic and the NTCP at time k is given by

NTCP k = P(L n, n0 ≤ k) = n0 ∑ i= n n0 i Fi (k) 1 -F(k) n0 -i . (36)
Proof Firstly, we determine the distribution of T ( n: n0 ) . As a matter of fact, if we put

A j (k) = { T ( j) ≤ k}, (37) then { T ( n: n0 ) ≤ k} = { n0 ∑ j=1 1 A j (k) ≥ n}. ( 38
)
As we have P A j (k) = F(k) and the events

A 1 (k), • • • , A n 0 (k) are independent then ∑ n0 j=1 1 A j (k) ∼ B( n0 , F(k)): P( T ( n: n0 ) ≤ k) = n0 ∑ i= n n0 i F(k) i 1 -F(k) n0 -i . ( 39
)
Finally, it is clear that NTCP at time k is given by

NTCP k = P(L n, n0 ≤ k). (40)

Determination of a suited number of dose fractions

The main task in radiotherapy is to choose the treatment schedule that permits to eliminate the tumor without causing significant damage on the adjacent normal tissue. This question may be reformulated as the choice of the appropriate treatment duration. In this diagram, the ideal therapeutic goal is the origin point O (0, 0) since it corresponds to a 100% efficient treatment and a null risk of complication. Note that k → TCP k and k → NTCP k are increasing.

Figure 5 shows the ECT curve for m = m = 3, q = 0.7, q = 0.6, r = 0.3, r = 0.2, n 0 = 10 4 , and n = n0 = 10 3 . We propose three strategies to select a suitable treatment duration:

• A first choice consists in limiting the probability of complication on normal tissues with a threshold fixed in this example at 0.05 (dotted vertical line). Beyond this level, the risk of complication is not acceptable. Subsequently, the most appropriate number of dose fractions to be applied is given by k = 13 dose fractions with a treatment non-efficiency probability estimated at 1 -TCP ≈ 0.3. • A second choice is to give priority to the therapeutic efficiency by using an non-efficiency threshold fixed in this example to 0.05 (dotted horizontal line). This choice leads to choose at least k = 16 days of treatment with a risk of complication estimated at NTCP ≈ 0.3. • A third possibility is a trade-off choice which consists in selecting the nearest point of the curve from the ideal point O (0, 0). In this case, the choice leads to select k = 14 days of treatment with a non-efficiency probability lower than 0.2 and a complication probability NTCP ≈ 0.05.

Limits

The current multinomial model used to calculate the tumor lifespan, TCP and NTCP has two main limitations:

• it does not take by-standing effects into account in the cell mortality process; • and numerical errors occur during Matlab computations when n 0 > 10 14 . In other terms, this numerical problem currently limits the application of our model to avascular tumors (< 1mm 3 ).

Nevertheless, studies are currently in progress to overcome these two computational limits.

CONCLUSION

We have proposed an original approach that expresses the probability distribution of the cancer and normal cells lifespans in terms of the treatment duration (or number of dose fractions) in radiotherapy. Our approach relies on a multinomial model, which is itself based on a discrete-time Markov chain. This model takes into account (i) the variety of cell responses according to their biological states and (ii) the repair mechanisms that occur between dose fractions. The main contribution is the introduction of two new response variables : the cancer cell and tumor lifespans, T and L n 0 respectively. A new expression of the TCP is proposed based on the cumulative distribution function of L n 0 . The (1α) confidence intervals of T and L n 0 were also determined. We have also introduced the normal tissue lifespan L n, n0 and provided an original formulation of the NTCP in terms of its cumulative distribution function G of T . Finally, we have proposed a new representation, called ECT (Efficiency-Complication Trade-off) which allows to optimize both efficiency and complication of the treatment. This characteristic could help the radiotherapeutist to choose the number of dose fractions suited to control the tumor while preserving the adjacent normal tissue.

  k [day]) Radiation schedule and living cell population

  Fig. 2. Fractionated radiation schedules (black bullets) with 5 dose fractions per week, and the mean size of the living cell population: n 0 -E(X k (m)) where k denotes the number of dose fractions in the treatment schedule and d is the discrete time based on a daily sampling rate. Red circles show the 99.9% confidence interval on the population size.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Sensitivity of the cell lifespan T to the target deactivation probability q

  ) Let us rearrange these random variables into a nondecreasing sequence T (1: n0 ) ≤ T (2: n0 ) ≤ . . . ≤ T ( n0 : n0 ) .

Fig. 5 .

 5 Fig. 5. ECT diagram (Efficiency Complication Trade-off)