Peter Van Der Gulik

Serge Massar

Dimitri Gilis

Harry Buhrman

Marianne Rooman
email: mrooman@ulb.ac.be

The first peptides: the evolutionary transition between prebiotic amino acids and early proteins

Keywords: origin of life, prebiotic peptides, protein ancestor, peptide structure and function

The issues we attempt to tackle here are what the first peptides did look like when they emerged on the primitive earth, and what simple catalytic activities they fulfilled. We conjecture that the early functional peptides were short (3 to 8 amino acids long), were made of those amino acids, Gly, Ala, Val and Asp, that are abundantly produced in many prebiotic synthesis experiments and observed in meteorites, and that the neutralization of Asp's negative charge is achieved by metal ions. We further assume that some traces of these prebiotic peptides still exist, in the form of active sites in present-day proteins. Searching these proteins for prebiotic peptide candidates led us to identify three main classes of motifs, bound mainly to Mg 2+ ions: D(F/Y)DGD corresponding to the active site in RNA polymerases, DGD(G/A)D present in some kinds of mutases, and DAKVGDGD in dihydroxyacetone kinase. All three motifs contain a DGD submotif, which is suggested to be the common ancestor of all active peptides. Moreover, all three manipulate phosphate groups, which was probably a very important biological function in the very first stages of life. The statistical significance of our results is supported by the frequency of these motifs in today's proteins, which is three times higher than expected by chance, with a P-value of 3 10 -2 . The implications of our findings in the context of the appearance of life and the possibility of an experimental validation are discussed.

Introduction

The development of molecular biology over the past half century has transformed the question of how life appeared from the level of philosophical speculations to the level of rigorous scientific enquiry. A number of key ideas have emerged which govern our thinking about this question, among which the synthesis of some amino acids and sugars by prebiotic synthesis [START_REF] Miller | Which Organic Compounds Could Have Occurred on the Prebiotic Earth? Cold Spring Harbor Symp Quan[END_REF][START_REF] Ricardo | Borate Minerals Stabilize Ribose[END_REF]), the RNA-world in which RNA catalyses its own duplication [START_REF] Gilbert | The RNA World[END_REF][START_REF] Johnston | RNA-Catalyzed RNA Polymerization: Accurate and General RNA-Templated Primer Extension[END_REF], and the role of lipid vesicles to limit the spatial extent of the cell precursor [START_REF] Deamer | Boundary Structures Are Formed by Organic Components of the Murchison Carbonaceous Chondrite[END_REF][START_REF] Chen | The Emergence of Cells During the Origin of Life[END_REF][START_REF] Mansy | Template-directed synthesis of a genetic polymer in a model protocell[END_REF]. At some point the controlled synthesis of proteins emerged, and proteins then took over many functions, presumably because of their great specificity and efficiency.

In this letter we address an important question, namely what did the first functional peptides and proteins look like when they emerged during very early life? Answering these questions is a difficult and unsolved issue. Indeed, the smallest natural proteins that can take a stable structure by themselves, without interacting with other biomolecules, are composed of around 20 amino acids (aa's). Some synthetic constructs are even smaller: for example, chignolin, a synthetic protein of 10 aa's, has been shown to have a stable structure in water [START_REF] Honda | 10-residue folded peptide designed by segment statistics[END_REF]. But it is difficult to imagine that functional proteins 10-20 aa's long suddenly appeared out of the blue. Indeed the sequences of these proteins are very specific, and reliably making functional proteins this long requires an efficient code for the aa sequence, and an efficient translation mechanism from the code into the protein (see the Discussion section).

Here we propose a solution to this "chicken and egg" paradox by suggesting that specific short peptides 3 to 8 aa's long could have served as catalysts during very early life.

Longer proteins would then have gradually evolved from these early precursors. The idea that very short peptides could have had a useful role in very early life has already been put A c c e p t e d m a n u s c r i p t 4 forward by Shimizu [START_REF] Shimizu | Specific Aminoacylation of C4N Hairpin RNAs with the Cognate Aminoacyl-Adenylates in the Presence of a Dipeptide: Origin of the Genetic Code[END_REF][START_REF] Shimizu | Histidine and Its Anticodon GpUpG are Similar Metabolic Reaction Rate Enhancers: Molecular Origin of the Genetic Code[END_REF][START_REF] Shimizu | Amino Acid and Anticodon Enhance Metabolic Reaction Rates Weakly but Specifically: Genetic Code World[END_REF] who showed that single amino acids and dipeptides could slightly enhance the rates of certain chemical reactions. But obviously some intermediate steps are required between Shimizu's dipeptides and the smallest of today's functional enzymes.

An additional constraint on any theory of how the first functional peptides emerged is that these should be composed exclusively or almost exclusively by the aa's that are efficiently produced by prebiotic synthesis. In most prebiotic synthesis experiments (see next section) the most efficiently produced amino acids are Gly, Ala, Val and Asp (or G, A, V and D in one letter code) [START_REF] Miller | Which Organic Compounds Could Have Occurred on the Prebiotic Earth? Cold Spring Harbor Symp Quan[END_REF]. Of these, the first three are neutral, and Asp is negatively charged. At first sight this is problematic: the absence of positive charges compensating the negative Asp's are likely to limit their ability to form stable structures and to carry out a catalytic activity.

We propose to resolve these conundrums by suggesting that the first peptides were composed of short chains of prebiotic aa's bound to (one or more) positively charged metal ions. Supposing the first peptides to be bound to metal ions solves two problems at once.

Namely, the metal ion(s) provide(s) an anchor around which the peptide can organize itself, thereby (at least partially) stabilizing its structure, and secondly it provides a positive charge which would be very useful for catalytic activity.

Later, as the coding and translation mechanisms improved, these very first peptides gradually lengthened, thereby improving the efficiency and specificity of their biological activity. We further conjecture that some of these first peptides, composed of prebiotic aa's bound to metal ions, have been conserved across evolution. The fact that active sites are better conserved than all other protein regions speaks in favor of this conjecture. If this idea is correct, it should be possible to find today the memory of the very first functional peptides in the active sites of some present-day proteins.

A c c e p t e d m a n u s c r i p t 5

The idea of finding in present-day proteins traces of very early life is not entirely new.

For instance it has been argued that today's aa abundances reflect the order in which they were introduced in the genetic code (see Zuckerkandl et al., 1971;[START_REF] Jordan | A universal trend of amino acid gain and loss in protein evolution[END_REF], and the criticism of [START_REF] Hurst | Causes of trends in amino-acid gain and loss[END_REF].

To explore the validity of our conjectures, we carried out a search in the protein structure database to identify all active sites composed almost exclusively of the abundant prebiotic aa's bound to metal ions. In the following we review the abundances of amino acids produced in prebiotic synthesis experiments, then present our precise search methodology, and discuss the classes of prebiotic peptide candidates revealed by the search. These candidates all correspond to active sites in their host protein, and perform catalytic functions likely to be important during very early life. Furthermore there is strong evidence that some of these active sites have been conserved at least since the Last Universal Common Ancestor (LUCA) 3.5 billion years ago. We conclude this letter by discussing in more detail how our suggestions fit into the wider picture of the appearance of life, and how they could be tested experimentally.

State of the art: Prebiotic amino acids

In this section we review the relative abundances of aa's that are produced in prebiotic synthesis experiments such as those originally designed by [START_REF] Miller | Production of some organic compounds under possible primitive earth conditions[END_REF], or that are observed in meteorites, in order to determine the most plausible composition of the first peptides that appeared on the primitive earth. These aa abundances are summarized in Table 1. Note that we do not review the determinations of aa's in the interstellar medium, as this issue is under discussion [START_REF] Kuan | Interstellar glycine[END_REF]Snyder et al. 2005;[START_REF] Jones | A search for biomolecules in Sagittarius B2 (LMH) with the Australia telescope compact array[END_REF]Zaia et al. 2008).

Miller and coworkers performed synthesis experiments in several environments attempting to mimic possible prebiotic circumstances, as the real environment is not agreed [START_REF] Miller | Production of some organic compounds under possible primitive earth conditions[END_REF], and recently reanalyzed by [START_REF] Johnson | The Miller volcanic spark discharge experiment[END_REF], consist of a highly reducing mixture of CH 4 , NH 3 , H 2 O and H 2 , which is sparked in three different ways. They all produced Gly, Ala, Asp, Glu and Ser in significant amounts. Moreover, in the experimental set-up which is nowadays considered the most interesting one (see [START_REF] Johnson | The Miller volcanic spark discharge experiment[END_REF], Val was produced in about the same amount as Glu.

In later work of Miller and coworkers [START_REF] Ring | Prebiotic synthesis of hydrophobic and protein amino acids[END_REF], ammonia was only present in trace amounts and the major source of nitrogen was N 2 . In this case ten proteinaceous aa's were produced, in order of abundance: Ala, Gly, Asp, Val, Leu, Glu, Ser, Ile, Pro and Thr. [START_REF] Schlesinger | Prebiotic sunthesis in atmospheres containing CH 4 , CO, and CO 2 . I. Amino acids[END_REF] investigated the effect of changing the carbon source (CH 4 , CO, or CO 2 ). They found that with CO or CO 2 as carbon source, Gly was almost the only amino acid found, with small amounts of Ala and Ser, and trace amounts of Asp and Glu, and, in a single case, Val. Recently, however, Miller and coworkers showed, contrary to previous reports, significant production of Ala, Ser and Glu in neutral atmospheres [START_REF] Cleaves | A reassessment of prebiotic organic synthesis in neutral planetary atmospheres[END_REF]. They conjectured that nitrite and nitrate, which are also produced by spark discharge in neutral atmospheres, destroyed the amino acids in previous work. Addition of pH buffer and oxidation inhibitor prevents this destruction. They concluded that neutral atmospheres may have been productive in prebiotic synthesis of amino acids, provided the early oceans were buffered sufficiently with respect to pH and redox balance.

Another series of experiments, performed by [START_REF] Plankensteiner | Amino acids on the rampant primordial Earth: electric discharges and the hot salty ocean[END_REF], assumed the prebiotic atmosphere to be neutral and composed of CO 2 , N 2 , H 2 O. The aa's that were systematically produced were Gly, Ala and Val. In addition, Ser, Pro, and Lys were detected in significant amounts in several experiments. Surprisingly, Asp and Glu were not detected.

An interesting observation is the similarity between the proteinaceous aa's produced in prebiotic synthesis experiments and those present in carbonaceous meteorites, as first noticed by [START_REF] Kvenvolden | The origin of life and evolutionary biochemistry[END_REF] and [START_REF] Miller | The atmosphere of the primitive earth and the prebiotic synthesis of amino acids[END_REF]. A detailed analysis of the composition of several Clearly, some discrepancies occur among the amino acids produced, according to the experimental setups, the environments meant to represent the primitive earth, or the type of meteorites. However, general tendencies are noted: Gly and Ala always appear as the most abundant aa's, and Asp, Val, Ser, and Glu often appear in significant amounts.

The choice we made in this paper of which were the abundant prebiotic aa's was motivated by the prebiotic synthesis experiments performed by [START_REF] Ring | Prebiotic synthesis of hydrophobic and protein amino acids[END_REF]. These authors found ten proteinaceous aa's as well as several non-proteinaceous ones in detectable amounts, with a spectrum showing strong similarities with that of the Murchison meteorite. In this experiment, Gly, Ala, Asp and Val were the four most abundant aa's, and we chose them as the presumed components of the first active peptides.

Methods: Search for traces of prebiotic peptides in present-day proteins

We searched for traces of prebiotic peptides in the RCSB Protein structure DataBase (PDB) [START_REF] Berman | The Protein Data Bank[END_REF]. For the purpose of selecting all PDB entries that contain a protein structure interacting with one metal ion at least and of identifying all interactions between ions and aa's, we used the PDBsum summary and analysis tool [START_REF] Laskowski | PDBsum more: new summaries and analyses of the known 3D structures of proteins and nucleic acids[END_REF]. The metals considered are Li, Be, Na, Mg, K, Ca, Cr, Mn, Fe, Co, Ni, Cu, Zn, and Se.

The prebiotic peptide candidates we searched for in the PDB were assumed to be built almost exclusively from the aa's that were presumably the most abundant in early times, i.e.

Gly, Ala, Val and Asp (see previous section). Moreover, these peptide candidates were required to be bound to metal ions, present in the prebiotic environment, in order to neutralize the negative charge of Asp and to allow the formation of well defined tertiary structures and catalytic functions. The exact criteria that we used to characterize prebiotic peptide candidates are the following. Firstly, these peptides must correspond to sequence regions of maximum length 8, the longest prebiotic peptides envisaged. Note that doubling this maximum length did not change the results of the search. Secondly, they must be bound to a metal ion, and this ion may not interact with any other part of the protein. Thirdly, the ion must be bound to the abundant prebiotic aa's only, among which three or more Asp residues. Requiring three Asp's at least increases the strength and specificity of the peptide-ion binding and thus its chance to be part of the protein active site in addition to stabilizing the structure. Indeed, when relaxing the criterion and accepting two Asp's bound to an ion, many spurious fragments were selected. Fourthly, the sequence starting at the first and ending at the last residue bound to the ion, and including start and end residues, should be composed of the abundant prebiotic amino acids for 80% at least. Lastly, if this sequence is flanked in the host protein by abundant prebiotic aa's, these aa's are added to the peptide.

To estimate the statistical significance of finding these prebiotic peptide candidates in the PDB, we computed E, the number of such peptides expected to occur by chance, as follows. Let E ij be the expected number of sequences of length i with j residues binding an ion. We compute this as E ij = P ij N ij , where P ij is the probability that a motif matches our prebiotic criteria, and N ij is the number of motifs of length i with j residues attached to the ion present in the PDB. The probabilities P ij are estimated as: possibilities to choose the location of the remaining j-2 prebiotic aa's that also bind the ion. Among these we have to choose at least three Asp's, giving rise to the sum from k=3 to j. Last we have to choose, out of the i-j remaining positions that do not bind the ion, at least 0.8i j ª º « » prebiotic aa's in order to have at least 80% such aa's in the motif. The total expected number of prebiotic motifs is given by:

3 0.8 2 ( | ) ( / / | ) 2 ( / / / | ) ( / / / | ) , (1)
j k j k ij k i j m i j m m i j i j P PDi o n PAGVi o n j k i j P D A G V ion P D A G V ion m ª º « » § • § • ¨¸¨ © ¹ © ¹ § • ¨© ¹ ¦ ¦
8 8 3 3 3 3 . ( 2
) i i ij ij ij i j i j E E PN

¦¦ ¦¦

However, this formula contains several overcountings. First, the PDB contains many highly similar, homologous, proteins. This implies that when a motif is found in one of them, it is almost automatically found in the others. To correct for this bias, we replace N ij by cor ij N , which counts only the motifs that have at least one different aa, or that have at least 20% different aa's in a sequence stretch of 10 residues around the middle residue of the motif.

Another correction is due to the fact that two consecutive Asp's binding an ion are very rarely observed. Indeed, the joint probability of finding two Asp's bound to an ion at consecutive positions along the sequence is much lower than the product of the independent probabilities of finding an Asp bound to an ion; the ratio of these probabilities is only equal to 0.05. This effect is due to geometrical constraints imposed by the polypeptide chain, and the unfavorableness of having two negative charges close together. This effect is also, but to a lower extent, observed for other consecutive aa's bound to an ion. To take this effect into account, we have to distinguish between the different ion-binding motifs, which we label by the index v. For example, for a motif of length i=5 with j=3 residues bound to an ion, the

2 2 i j § • ¨ © ¹ =3 possible
2 ij i i cor cor cor cor ij ijv ijv ij i j i j v i E E P N j D D § • ¨ © ¹ ¦¦ ¦¦¦
where the corrected probability cor ijv P is given by:

3 1 0.8 ( | , ) ( | , ) ( | ) ( | , ) ( | , ) ( | ) ( / | ) ( / | ) , with . (4
E ª º « » § • § • ¨¸¨© ¹ © ¹ ¦¦ ¦
where a=A/G/V, and ( | , ) P y z ion is the probability that an aa bound to an ion is of type y,

given that the next aa along the sequence is bound to the same ion and is of type z. The index To objectively compare E (or E cor ), the expected number of prebiotic peptide candidates, and N (or N cor ), their actual number observed in the PDB, we estimate the P-value, defined as the probability of finding at least N (or N cor ) prebiotic peptide candidates if the null hypothesis was true. For this, we assume that E (or E cor ) follow a Poisson distribution, as it is a discrete probability distribution only defined for positive values. The P-value is then given by: Using the above definition of prebiotic peptide candidates, we screened systematically the whole PDB, and found two different ion binding motifs: DxDxD and DxxxxDxD, where x stands for any aa. These two motifs differ in the relative positions along the sequence of the three Asp residues binding the ion. The first binding motif further divides into four groups, according to the aa's between and flanking the Asp's, as summarized in Table 2.

w
The first two DxDxD binding motifs correspond to active sites of RNA polymerases, and are both bound to Mg 2+ , or sometimes Mn 2+ ions. The first motif, with strictly conserved sequence ADFDGD, or even NADFDGD if one includes the non abundant prebiotic residue N, is found in DNA-dependent RNA polymerase, the normal RNA polymerase engaged in protein-coding gene expression. More precisely, it is found in the second domain of RNA polymerase Rpb1, using the nomenclature of the protein families database Pfam [START_REF] Finn | The Pfam protein families database[END_REF].

It is worthwhile to linger a few moments on the age of this NADFDGD sequence. It corresponds to the active site in the mRNA-producing enzyme of many organisms, such as Homo sapiens, the nematode worm Caenorhabditis elegans, but also in the yeast Saccharomyces cerevisiae and the plant Arabidopsis thaliana [START_REF] Iyer | Evolutionary connection between the catalytic subunits of DNA-dependent RNA polymerases and eukaryotic RNA-dependent RNA polymerases and the origin of RNA polymerases[END_REF]. This sequence seems thus to be present in all eukaryotes, which populate the earth since two billion years at least. Moreover, it is present in other organisms, which share a common ancestor with eukaryotes in an even more remote past, in particular in the bacterial genera Thermus, Deinococcus, Escherichia, Agrobacterium, Mesorhizobium, Helicobacter, Ureaplasma, Bacillus, Thermotoga, Synechocystis, Aquifex, Chlamydia, Mycoplasma, Mycobacterium and The orthodox way in molecular biology to interpret the almost completely universal appearance of some characteristic in a class of organisms, is that their common ancestor already did have this characteristic. The sequence NADFDGD is found in the mRNAsynthesizing RNA polymerase of all living cellular organisms, and is thus assumed to be present in their common ancestor, which lived approximately 3.5 billion years ago. This is thus the age of the sequence NADFDGD, or at least its minimum age given that there is still a long history from the origin of life up to the time of LUCA.

The second motif, with sequence GGDYDGD, corresponds to the active site of a cellencoded RNA-dependent RNA polymerase of Neurospora crassa, which is part of the RdRP family according to the Pfam database. Its role is to produce double-stranded RNA from aberrant single-stranded RNA as part of an RNA silencing response [START_REF] Salgado | The Structure of an RNAi Polymerase Links RNA Silencing and Transcription[END_REF].

Although this kind of protein is not present in insects, vertebrates and the yeast Saccharomyces cerevisiae, it is a kind of protein that is a normal constituent of an eukaryotic cell. For example, related sequences are known from the nematode worm Caenorhabditis elegans, the plants Arabidopsis thaliana, Nicotiana tabacum and Oryza sativa, and the unicellular eukaryotes Dictyostelium discoideum and Giardia intestinalis [START_REF] Iyer | Evolutionary connection between the catalytic subunits of DNA-dependent RNA polymerases and eukaryotic RNA-dependent RNA polymerases and the origin of RNA polymerases[END_REF].

The GGDYDGD motif is not totally conserved among these sequences, but has the following (NTP) that is going to be built into the RNA polymer, as illustrated in Figs. 2a-b.

What is the evolutionary relationship between these two RNA polymerase active sites, which share the same geometry, the same mechanism and the same Asp residues and Mg 2+ or Mn 2+ ions? [START_REF] Iyer | Evolutionary connection between the catalytic subunits of DNA-dependent RNA polymerases and eukaryotic RNA-dependent RNA polymerases and the origin of RNA polymerases[END_REF] tend to view these common characteristics as a sign of very ancient common ancestry, dating from more than 3.5 billion years ago. In this view, the RNAdependent RNA polymerase would originally have been the mRNA-producing enzyme of an organism that lived before the last common ancestor of all contemporary cells, and would have found its way to the eukaryote by a pathway of horizontal gene transfer.

The Asp residues and the Gly residue in the ancestral RNA polymerase consensus sequence xxDxDGD can be expected to play such an important role in the polymerase function that they cannot be changed, Asp because of its specific interaction with the metal ions, and Gly because of its ability to make main chain turns impossible with other amino acids. This does not exclude other ways to make an RNA polymerase; however, when starting from this particular sequence, function is lost upon mutation. Note that the extreme conservation of the other aa's in the NADFDGD motif of DNA-dependent RNA polymerase cannot be attributed to the conservation of the function, but is probably related to their interaction with other parts of the protein and other parts of the transcription machinery.

A c c e p t e d m a n u s c r i p t 15

The DxxxxDxD motif is represented in our search by just one sequence, DAKVGDGD, which binds to two Mg 2+ ions and turned out to be the active site of dihydroxyacetone (DHA) kinase (Fig. 1c). Because this enzyme was producing the glycolytic intermediate dihydroxyacetone-phosphate (DHAP), we consider this enzymatic activity as interesting from the viewpoint of the transition of prebiotic to early biotic. Like the situation we came across with the RNA polymerases, the active site consists of two Mg 2+ ions, coordinated by three Asp residues, and these Asp residues and Mg 2+ ions are essential for catalysis (Siebold et al.

2003).

As with the RNA polymerase active sites, the phosphate groups of an NTP are complexed with the Mg 2+ ions (Fig. 2d). The J-phosphate is used to phosphorylate the substrate of this kinase enzyme. The substrate spectrum may be broader than the name DHA kinase suggests, because the activity of flavin mononucleotide (FMN) cyclase was found to be identical to that of human DHA kinase [START_REF] Erni | Small substrate, big surprise: fold, function and phylogeny of dihydroxyacetone kinases[END_REF].

The form of the active site of DHA kinase is different from that of the RNA polymerase and mutase classes (Fig. 1). The catalytic action presents similarities with that of RNA polymerases, but differs when looking in detail: in this case the bond between the second and the third phosphate is broken, while in the case of the RNA polymerases, it is the bond between the first and the second phosphate that is being broken.

Though the three binding motifs D(F/Y/M/L)DGD, DGD(G/A)D and DAKVGDGD have different sequences and structures, they share a common submotif, DGD, with the two Asp's binding the metal ion. The conformation adopted by this submotif differs according to their host motifs, which is not surprising as such a small peptide is necessarily flexible.

However, DGD could be assumed to be the ancestor of all three motifs, which could have evolved over time to fulfill more specific biological functions.

To assess the statistical significance of finding these sequence motifs in present-day proteins and proposing them as prebiotic peptide candidates, we compute the number E of

A c c e p t e d m a n u s c r i p t

16 such motifs that are expected to be found in the PDB by chance, as defined by eqs (1)-( 2). We found E=30, whereas the number of motifs indeed present in the PDB is equal to N=67 (Table 2); the associated P-value (see eq. ( 5)) is 4 10 -9 . When taking into account out all biases, in particular the presence of homologous proteins in the PDB, the geometrical constraints that tend to prevent successive residues to be bound to an ion, and the preference for an Asp to be bound to an ion rather than to be close but not bound to it, we find using eqs

(3)-( 4) that the expected number of occurrences is E cor =1.8, whereas its actual number is N cor =5. Thus, we found 2.8 times more prebiotic peptide candidates than expected by chance, leading to a P-value of 3 10 -2 . We would moreover like to stress that all these peptides are part of the active site of their host protein, although this was not among the selection criteria.

This strongly improves the statistical significance of our results, even though this improvement cannot be estimated quantitatively at the moment.

As described in the state-of-the-art section, it is not easy to determine exactly which were the most abundant aa's on the primitive earth, and therefore our choice of Ala, Val, Asp and Gly as prebiotic aa's may give rise to criticisms. Indeed, Asp and Val do not come out in all prebiotic synthesis experiments (Table 1). Asp cannot be removed because its negative charge is indispensable for binding tightly to the metal ions, but Val could be dropped.

However, this choice only removes a single motif from the prebiotic peptide candidates, i.e.

DxxxxDxD. The Ala-Gly-Asp choice of abundant prebiotic aa's thus does not change the statistical significance level of our results by much, as we have E cor =1.3 and N cor =4, with a Pvalue of 5 10 -2 .

As an additional test of our procedure, we investigate the effect of replacing Asp by His among the abundant prebiotic aa's, given that histidines also bind very efficiently to metal ions. This led to detect N cor =2 prebiotic peptide candidates, whereas the expected number is E cor =0.1. In these two hits, however, the structure of the host protein is incomplete: the

A c c e p t e d m a n u s c r i p t

17 coordinates of 5 and 11 residues, respectively, are lacking in the immediate vicinity of the metal ion. These matches have thus to be considered as spurious. Finally, when considering as prebiotic aa's all 10 aa's produced in the experiments of [START_REF] Ring | Prebiotic synthesis of hydrophobic and protein amino acids[END_REF], i.e. Val, Gly, Asp, Ala, Ile, Pro, Ser, Thr, Leu and Glu (Table 1), the number of expected and observed peptide candidates is of course higher: N cor =22 and E cor =20. However, given that this set of aa's contain the four abundant prebiotic aa's Val, Gly, Asp, and Ala, for which we found N cor =5

and E cor =1.8, we may conclude that the six additional aa's yield less peptide candidates than expected by chance.

In conclusion, the most statically significant results come out when considering Gly, Ala and Asp, and possibly Val, as abundant prebiotic aa's. This can be viewed as an a posteriori indication that these are really the true abundant prebiotic aa's.

Discussion: the origin of life and the first peptides

How do our results fit into a wider picture of the appearance of life? The present-day picture of how life appeared is punctuated by a number of established facts, in a background of unknown or dubious stages and links. We make here an attempt to review the results that are well verified or supported, in their tentative order of occurrence:

(1) Many essential molecules for life can be synthesized in a prebiotic environment: these include amino acids [START_REF] Miller | Which Organic Compounds Could Have Occurred on the Prebiotic Earth? Cold Spring Harbor Symp Quan[END_REF]), but also pentoses [START_REF] Ricardo | Borate Minerals Stabilize Ribose[END_REF]), nitrogenous bases [START_REF] Cleaves | The prebiotic synthesis of pyrimidines in frozen solution[END_REF], and vesicle-producing organic compounds [START_REF] Deamer | Boundary Structures Are Formed by Organic Components of the Murchison Carbonaceous Chondrite[END_REF].

(2) Still in a prebiotic environment, more complex molecules can be built from these basic building blocks: adenine and ribose can be condensed into a nucleoside by the presence of borate in aqueous solution [START_REF] Prieur | Étude de l'activité prébiotique potentielle de l'acide borique[END_REF]; polyphosphate, likely to be a source of energy and phosphate for producing nucleotides, can be produced by volcanic processes (Yamagata et al. 1991); montmorillonite clay has been discovered to catalyze

A c c e p t e d m a n u s c r i p t

18 oligonucleotide formation [START_REF] Ferris | Montmorillonite-catalysed formation of RNA-oligomers: the possible role of catalysis in the origins of life[END_REF]. Of particular importance to the present article, dipeptides can be produced in a dehydrating salt solution [START_REF] Schwendinger | Possible Role of Copper and Sodium Chloride in Prebiotic Evolution of Peptides[END_REF], and longer peptides can be synthesized in the presence of montmorrillonite [START_REF] Rode | The combination of salt induced peptide formation reaction and clay catalysis: a way to higher peptides under primitive earth conditions[END_REF].

(3) Vesicles, the ancestors of cells, have been shown to have a prebiotic origin [START_REF] Deamer | Boundary Structures Are Formed by Organic Components of the Murchison Carbonaceous Chondrite[END_REF] and to be catalyzed by montmorillonite [START_REF] Ferris | Montmorillonite-catalysed formation of RNA-oligomers: the possible role of catalysis in the origins of life[END_REF]. The fact that vesicles containing oligonucleotides have a higher fitness as compared with empty vesicles (Chen et al., 2004), and the fact that vesicles have a ribose-collecting power [START_REF] Sacerdote | Semipermeable lipid bilayers exhibit diastereoselectivity favoring ribose[END_REF] are indications of the early involvement of vesicles in the process of the origin of life. Also important is the fact that a membrane potential is not necessarily coupled to cells, but can be found in much simpler vesicles (Chen & Szostak 2004).

(4) Homochiral molecules were selected out of the racemic mixture, by a mechanism which is still not elucidated.

(5) At some stage RNA and peptides became part of the same biochemical system. However, which of the RNA or peptides came first is highly controversial. The arguments in favor of the RNA world include the existence of ribozymes, the essential role of RNA in transcription and translation, and the fact that aa's are activated by ATP before they are coupled to their cognate tRNAs. An argument against the RNA world is the huge size of the RNA-dependent RNA-polymerizing ribozyme [START_REF] Johnston | RNA-Catalyzed RNA Polymerization: Accurate and General RNA-Templated Primer Extension[END_REF]. The arguments in favor of the peptide world include the fact that peptides are more easily synthesized and more resistant to degradation (Plankensteiner et al., 2005). Moreover some peptides have been shown to catalyze their own synthesis [START_REF] Lee | A selfreplicating peptide[END_REF] and have catalytic activity (Kochavi et al., 1997). A hybrid picture between these two extremes would be an RNA-peptide world, in which both types of molecules evolved together. Exactly how these different more or less well established facts are related is the subject of much debate and speculation, and is outside the scope of the present paper. Here we addressed the issue of the form and functionality of the first peptides. We believe it should be possible to answer the latter question to a large extent independently of how points (1) to ( 6) fit together.

It can also to some extent be addressed independently of how the first peptides were synthesized and how their sequences were coded, and of whether a peptide, RNA, or RNApeptide world appeared first. Concerning the last point there are two conflicting results which cannot be resolved at present: on the one hand small peptides can be synthesized in prebiotic environments as mentioned in point (2) above, but on the other hand peptide synthesis as it occurs in the cell is deeply connected to the RNA world, and the traces of this connection can still be seen in the essential role in translation of mRNA, rRNA, and tRNA. However this may not have been the primitive situation.

Two main facts can be stated with confidence about the first functional peptides: first of all the first peptides were presumably made of the amino acids that are efficiently produced in prebiotic synthesis, as these are the chemically simplest aa's. A supporting argument is that the present genetic code, with its block structure, could easily have evolved from a more primitive code that used fewer amino acids: from time to time some codons were reassigned to other aa's without affecting the other codons.

Secondly, the first peptides probably appeared when a coding mechanism for the aa sequence was inexistent, or at least highly imperfect and error prone. Indeed the present (highly complex) coding, with its essentially error free mechanisms for duplication, transcription and translation, was almost certainly not present initially, but presumably

A c c e p t e d m a n u s c r i p t

20 evolved in response to the increasingly critical role of complex proteins in early life. This is supported by the fact that the genetic code seems to have been optimized against transcription and translation errors [START_REF] Freeland | The genetic code is one in a million[END_REF] A first question to address when trying to understand what were the first functional peptides is to decide whether they had lengths comparable to present-day protein chains (about 20 to 10 4 aa's long), or whether they were much shorter. We believe that the first possibility is highly unlikely. Indeed long random peptides will not have a well defined 3dimensional structure and certainly will not have any useful functionality. On the other hand random synthesis of very short peptides could have produced useful ones with reasonably high yield. For instance the common submotif of our prebiotic peptide candidates, DGD, could be one of those randomly selected peptides. In this case the success rate of the synthesis would have been 4 -3 =1/64 (assuming that only the four abundant prebiotic residues were used and were equally abundant). The fraction of active versus possible peptides drops very rapidly as the length increases (at least if all positions have to be occupied by specific residues): for 8 residues it is as small 4 -8 . However, if the synthesis occurred in the presence of metal ions, which tend to attract free Asp's, the amount of meaningful peptides could be much higher.

The prebiotic active sites of 5-8 residues long, such as those we identified in this paper, seem to form the limit between randomly synthesized and template-directed peptides.

To discover what the original functional peptides were, we have adopted here a theoretical approach. We have identified, using the method described above, active sites of present-day proteins that are made of prebiotic amino acids, carry out functionalities that are essential for early life, and can be traced back to LUCA.

We then take the bold step of supposing that these sequences in fact go back much earlier, and constitute some of the earliest useful peptides. This is of course highly conjectural. If we accept this conjecture, we obtain a remarkably precise picture of the Remarkable is the fact that ribozymes also need specifically Mg 2+ ions for both structure and function [START_REF] Glasner | Metal ions requirements for structure and catalysis of an RNA ligase ribozyme[END_REF]Takagi et al., 2004;[START_REF] Robertson | The Structural Basis of Ribozyme-Catalyzed RNA Assembly[END_REF].

Mg 2+ ions seem thus to have a particularly important role in handling phosphate groups and therefore in the development of life.

Moreover, in the peptides we identified, the biological functions performed by these active sites, i.e. RNA polymerization, glycolytic intermediate production and sugar phosphorylation, can be assumed as prebiotic, or at the very least, very early biotic. The existence of these early peptides would of course have required the existence of a stable environment, protected from toxic chemical species. This environment could have been provided by the vesicles that correspond to the ancestors of cells (point (3) above). At present we do not have enough information to decide whether these oligopeptide motifs appeared before, during or after the RNA world. In the latter case they could have constituted an interesting improvement to the phosphate chemistry available in the purely RNA world.

More generally, we can conjecture that, once small random peptides became useful to early life, mechanisms emerged to select for the active peptides. Thus gradually a coding and translation mechanism appeared, which in turn would allow for the synthesis of longer peptides, and for the exploration of the space of possible aa sequences, progressively selecting more active and specific oligopeptides and proteins. Therefore, the appearance of prebiotic oligopeptides is tightly connected with the appearance of the genetic code. Note however that we have no idea about how the transition occurred from the production of interesting compounds by noncoded synthesis of random peptides to the production of coded peptides by a ribosome, and about the gradual evolutionary construction of the complex system of tRNAs,

Conclusion: how to validate our findings?

We identified in nowaday protein structures five short segments bound mainly to Mg 2+ ions, and built for 80% at least from the four aa's Asp, Gly, Val, Ala, thought to be the most abundant aa's in prebiotic times. These motifs appear as statistically significant as they are 2.8 times more frequent than expected by chance, with a P-value of 3 10 -2 . Moreover, they all form the active sites of their host protein: DFDGD and DYDGD in RNA polymerases, DGDGD and DGDAD in mutases, and DAKVGDGD in dihydroxyacetone kinase, where they manipulate phosphate groups, thought to be an important biological function in the very first stages of life. We thus conjecture that these motifs could correspond to the first functional peptides, at the earliest stages of life, before or at the beginning of the setup of the genetic coding mechanism, and after the purely RNA world or during the RNA-peptide world.

These peptides could be viewed as transitional fossils, constituting a missing link between the prebiotic amino acids and the coded proteins.

Note that this conjecture must considered with caution, as there could be alternative explanations for the prominent presence of Ala, Val, Asp and Gly in small, metal-binding, universal active sites. For example, these aa's could be indispensable for achieving a specific biological function. This argument is however dubious, given that other RNA polymerases are known to exist with active sites consisting of different aa's. Another explanation is that an intense process of selection within the first organisms led to the small, metal-binding active sites notwithstanding the presence of many more proteinaceous aa's in the prebiotic environment and other possible ways of making the active site. In summary, the present work pursues a new approach to studying very early life.

Through the analysis of the sequence and structure of present-day proteins, more and more of which are available in public databases, we try to push the well established method of phylogeny to even remoter periods. In this way we formulate new hypotheses about the origin of life. Table 1. Relative frequencies of proteinaceous aa's in prebiotic synthesis experiments and in meteorites. 1 [START_REF] Miller | Production of some organic compounds under possible primitive earth conditions[END_REF]; [START_REF] Ring | Prebiotic synthesis of hydrophobic and protein amino acids[END_REF]; [START_REF] Schlesinger | Prebiotic sunthesis in atmospheres containing CH 4 , CO, and CO 2 . I. Amino acids[END_REF]; [START_REF] Johnson | The Miller volcanic spark discharge experiment[END_REF]. The residues in bold are bound to the metal ion. 2

The central structure of the group is that for which the rms deviation with respect of all other members of the group is minimum. 3

The rms deviation is computed between all heavy atoms for the main chain and side chain and the metal ion, after coordinate superimposition; the mean is taken between all pairs in the class. 4

The rms deviation is computed between the equivalent parts of the central structures of two classes, in particular between DFDGD and DYDGD (considering the common atoms in F and Y), and between DGDGD and DGDAD (considering the common atoms in G and A). 5

The consensus sequence is obtained from sequence alignments in homologous protein families.

  first three environments considered by

  performed by[START_REF] Botta | Relative amino acid concentrations as a signature for parent body processes of carbonaceous chondrites[END_REF], and showed the presence of Gly, Ala, Asp, Val, Ser and Glu.

  where '/' means 'or', and ( | ) P y ion [ ( | ) P y ion ] is the fraction of aa's of type y, given that they are bound [not bound] to one of the considered ions; these fractions are computed from . To obtain this formula, we first fixed the outer aa's binding the ion,

  motifs are iixxi, ixixi and ixxii, where i and x denote aa's bound to the ion, respectively; cor ijv N with i=5, j=3 and v=1 corresponds thus to the number of non-homologous matches of motifs of type iixxi in the PDB. This leads us to the corrected expected value E cor :

  labels the E jk motifs of type v; for example, for j=4 and k=3, the E jk =4 ion-binding motifs of type iiixi are DDDxa, DDaxD, DaDxD and aDDxD. The integers yz vw n denote the number of times, in a motif of type v (e.g. iixi) and subtype w(e.g. DDaxD), two successive aa's bound to an ion, the first being of type y and the second of type z. In the example DDaxD,

  the archaeal genera Sulfolobus, Aeropyrum, Methanosarcina and Halobacterium.

  variations: G(G/S/A)D(Y/F/M/L)DGD, where the '/' symbol means 'or'; it is thus very close to the NADFDGD motif.The catalytic domains of these two types of RNA polymerases, Rpb1 domain 2 and RdRP, are structurally similar, whereas their other domains are completely different. The geometries of their active sites bound to Mg 2+ or Mn 2+ ions are particularly alike: they present a root mean square (rms) deviation of 0.67Å (see Table2), and their resemblance is clearly .1. This leads us to merge these two motifs and to define a single RNA polymerase motif, of consensus sequence D(F/Y/M/L)DGD, with an invariant Gly and three invariant Asp residues. Note that the only non abundant prebiotic residue of the motif, F/Y/M/L, is not well conserved. Moreover, it points towards the protein core rather than to the ions. It has thus probably a structural role, but is unimportant for the catalytic function; this leads to the conjecture that this position was occupied by an abundant prebiotic residue A, G or V in early times. The three conserved Asp's of the motif are involved in the binding of one or two Mg 2+ or Mn 2+ ions, which coordinate the phosphates of the nucleotide triphosphate

  of thought, xxDxDGD can be viewed as one of the most profound molecular fossils in life, and possibly one of the first coded peptide sequences. In these very early times, more than 3.5 billion years ago, they were possibly built from Val, Ala, Asp and Gly residues only, and might have been originally composed of amino acids with a prebiotic rather than a biogenic origin.The last two DxDxD motifs presented in Table2correspond to the two sequences DGDAD and DGDGD, bound to Mg 2+ , Zn 2+ , or Ni 2+ ions. They have very similar structures, as shown in Fig.1b; their rms deviation is equal to 0.77Å. They are also part of the same homologous family according to Pfam[START_REF] Finn | The Pfam protein families database[END_REF], named PGM_PMM_II, with consensus sequence DGD(G/A/F)D. These motifs correspond to the active sites in enzymes with mutase activity (see Fig 2.c). These enzymes can usually bind several metal ions, but Mg 2+ yields systematically the maximum activity[START_REF] Regni | Crystal structure of PMM/PGM: an enzyme in the biosynthetic pathway of P. aeruginosa virulence factors[END_REF]. Some, like rabbit phosphoglucomutase, show a slight activity with Ni 2+[START_REF] Regni | Crystal structure of PMM/PGM: an enzyme in the biosynthetic pathway of P. aeruginosa virulence factors[END_REF], whereas Zn 2+ is usually observed to inhibit the enzyme(Regni et al. 2006a). As an example, phosphomannomutase/ phosphoglucomutase of the bacterium Pseudomonas aeruginosa is changing the position on which a phosphorylated hexose sugar molecule is carrying the phosphate group. The Mg 2+ ion, which is essential for catalysis(Shackleford et al. 2006), is coordinated by the three Asp's. The catalytic process, in which a biphosphorylated sugar intermediate has to engage in a dramatic, 180º reorientation, is described in detail inRegni et al. (2006b). Of course, a complex reaction mechanism like this is a long way from simple SUHELRWLF FLUFXPVWDQFHV 1HYHUWKHOHVV WKH SURWHLQV RI WKH Į-D-phosphohexomutase enzyme superfamily could have inherited the simple DGD(G/A)D Mg 2+ ion-binding motif from a more remote past, possibly from a very early biotic age, where this motif fulfilled simplified catalytic activities.

  The coding of genetic information in DNA was probably, given the difficulty of synthesizing DNA, one of the last steps before the appearance of life as we know it. But this issue is far from settled.

  the first, very short, oligopeptides could have played during very early life. In all cases we have identified, phosphate groups are manipulated, mainly by Mg 2+ ions, coordinated by oligopeptide motifs consisting of three Asp residues and at least one Gly residue.

  rRNAs and mRNAs; all this occurred moreover in the presence of the various toxic substances that must have existed on a primeval earth. These issues, however, are beyond the scope of the present paper.

  suggestions be tested experimentally? The easiest point to analyze is whether the short prebiotic peptide candidates are soluble in aqueous solution, and whether they form more or less unique and stable structures with Mg 2+ ions. The second point to test is whether these short, metal ion-complexed peptides do still have some kind of enzymatic activity, even very inefficient, when they are isolated from their normal protein context. This would possibly require introducing other environments, present in the prebiotic era, which would favor the assembly of the different partners required for enzymatic function. Finally, we could consider peptides built from the abundant prebiotic residues only, for example ADVDGD and DAAVGDGD, which could represent even earlier peptides, and test if their properties are comparable to those of ADFDGD and DAKVGDGD. A further area of research is the interaction of these kinds of oligopeptides with RNA molecules.

  ., Arnold, I., Garcia-Alles, L.F., Baumann, U. & Erni, B. 2003 Crystal structure of the Citrobacter freundii Dihydroxyacetone Kinase Reveals an Eight-VWUDQGHG Į-

Figure 1 .

 1 Figure legends

Figure 2 .

 2 Figure 2. Biological function of the three prebiotic peptide candidates. The figures are drawn

  Prebiotic peptide candidates found in the protein structure database. 1

  Figure 1.a

Figure

  Figure 2.d

  )

	cor ijv P	j	E	jk	DD vw n P D D ion P D a ion P D ion Da vw n	DD Da vw vw k n n	aD vw n P a D ion P a a ion P a ion aa vw n	vw j k n n vw aD aa
		k w
					m	i j i j	i j m	m P D a ion P D a ion	i j m	jk	j k

  HelicalBarrel ATP-binding Domain. J Biol Chem 278, 48236-48244.

	Snyder, L.E., Lovas, F.J., Hollis, J.M., Friedel, D.N., Jewell, P.R., Remijan, A., Ilyushin,
	V.V., Alekseev, E.A. & Dyubko, S.F. 2005 A rigorous attempt to verify interstellar
	glycine. Astrophys. J. 619, 914-930.
	Takagi, Y., Ikeda, Y. & Taira, K. 2004 Ribozyme mechanisms, Top. Curr. Chem. 232, 213-
	251.
	Yamagata, Y., Watanabe, H., Saitoh, M. & Namba, T. 1991 Volcanic production of
	polyphosphates and its relevance to prebiotic evolution. Nature 352, 516-519.

Zaia, D.A.M., Zaia, C.T.B.V. & De Santana, H. 2008 Which amino acids should be used in prebiotic chemistry studies? Orig. Life Evol. Biosp. 38, 469-488.

Zuckerkandl, E., Derancourt, J., and Vogel, H.1971 Mutational trends and random processes in the evolution of informational macromolecules. J.

Mol. Biol. 59,[473][474][475][476][477][478][479][480][481][482][483][484][485][486][487][488][489][490]

Acknowledgments

We thank Roman Laskowski for his kind help with PDBsum. MR and DG acknowledge support from the Belgian State Science Policy Office through an Interuniversity Attraction Poles Program (DYSCO); MR, DG and SM acknowledge support from the Fund for Scientific Research (FRS) through FRFC projects; PvdG and HB acknowledge support from the Bsik research grant BRICKS, and HB from an NWO VICI grant. SM and MR are Research Directors at the FRS.


Alessandro Bertuzzi
email: bertuzzi@iasi.cnr.itgandolfi@iasi.cnr.itsinisgalli@iasi.cnr.it

Antonio Fasano
email: fasano@math.unifi.it

Alberto Gandolfi

Carmela Sinisgalli

Necrotic core in EMT6/Ro tumour spheroids: Is it caused by an ATP deficit?

Keywords: tumour spheroid, necrosis, energy metabolism, ATP, diffusive transport, free boundary problems 1

Although commonly related to nutrient deprivation, the cause of the formation of the necrotic core in the multicellular tumour spheroids is still a controversial issue. We propose a simple model for the cell ATP production that assumes glucose and lactate as the only fuel substrates, and describes the main reactions occurring in the glycolytic and the oxidative pathways. Under the key assumption that cell death occurs when ATP production falls to a critical level, we formulate a multiscale model that integrates the energy metabolism at the cellular level with the diffusive transport of the metabolites in the spheroid mass. The model has been tested by predicting the measurements of the necrotic radius obtained by Freyer and Sutherland (1986a) in EMT6/Ro spheroids under different concentrations of glucose and oxygen in the culture medium. The results appear to be in agreement with the hypothesis that necrosis is caused by ATP deficit.

A c c e p t e d m a n u s c r i p t

Necrotic core in EMT6/Ro tumour spheroids: is it caused by an ATP deficit?

A c c e p t e d m a n u s c r i p t 1 Introduction

Poorly vascularized tumours show the occurrence of necrotic regions that have been mainly correlated to the deprivation of oxygen [START_REF] Thomlinson | The histologic structure of some human lung cancers and the possible implications for radiotherapy[END_REF][START_REF] Tannock | The relation between cell proliferation and the vascular system in a transplanted mouse mammary tumour[END_REF], Höckel and Vaupel, 2001). Investigations on in vitro growing tumour spheroids have been useful to characterize the dependence of necrosis formation on the external environment [START_REF] Sutherland | Cell and environment interactions in tumor microregions: the multicell spheroid model[END_REF]. During the growth of the spheroid the fraction of proliferating cells decreases, the cells in the inner region become deprived of oxygen, glucose and other nutrients whereas metabolic waste accumulates, and the formation of a necrotic core is observed. In an advanced stage of growth, indeed, spheroids show an outer viable rim (whose thickness ranges from about 100 to 250 μm) that surrounds the necrotic core. Because of the degradation of dead cells in the necrotic core [START_REF] Folkman | Cell-regulation of growth in three dimensions[END_REF], the spheroid eventually attains a limiting size with a final diameter of 1-3 mm. The spheroid growth has been described by a variety of mathematical models using either a continuum or a discrete approach (see [START_REF] Araujo | A history of the study of solid tumour growth: the contribution of mathematical modelling[END_REF], for a comprehensive review). Almost all these models assumed that cell proliferation and cell death depend on the concentration of a single critical chemical that diffuses from the external medium into the cell aggregate. The diffusion of glucose, oxygen, lactate and other substrates has been included in a spheroid model proposed by Casciari et al. (1992). In that model, the cell proliferation rate was taken as an empirical function of glucose, oxygen and hydrogen ion concentrations. Both glucose and oxygen transport and consumption were included in discrete models (Schaller andMeyer-Hermann, 2005, Jiang et al., 2005) and in a continuum model of the spheroid growth [START_REF] Schaller | Continuum versus discrete model: a comparison for multicellular tumour spheroids[END_REF].

The hypothesis that the formation of necrotic regions is caused by a deficit in ATP production was adopted in a model of cell population growth in a tridimensional (sandwich) culture [START_REF] Hlatky | Joint oxygen-glucose deprivation as the cause of necrosis in a tumor analog[END_REF]. More recently the cell energy metabolism, i.e. the intracellular ATP production, has been incorporated in continuous models of spheroid [START_REF] Venkatasubramanian | Incorporating energy metabolism into a growth model of multicellular tumor spheroids[END_REF] and tumour cord growth [START_REF] Astanin | Mathematical model of tumour cord growth along the source of nutrient[END_REF], as well as in cellular automata tumour growth models (Smallbone et al., 2007, Gerlee and[START_REF] Gerlee | A hybrid cellular automaton model of clonal evolution in cancer: the emergence of the glycolytic phenotype[END_REF]. [START_REF] Bertuzzi | ATP production and necrosis formation in a tumour spheroid model[END_REF] proposed a simple mathematical model of cell ATP production that accounts for the main reactions occurring in the glycolytic and the oxidative pathways. The formation of the necrotic core in spheroids was described in that paper by assuming that cell death occurs when the ATP production rate falls below a critical threshold. This assumption led us to consider the diffusion free boundary problem for the concentrations of glucose, lactate and oxygen inside the spheroid viable rim.

In the present paper, our previous model of cell energy production is modified by simplifying the description of lactate transport and by adopting a more realistic assumption on pyruvate flux toward the Krebs cycle. Moreover, the diffusion of the nutrients in an unstirred layer surrounding the spheroid surface has been taken into account. The model of necrosis formation, based on the threshold on ATP production rate, has been applied to predict the measurements of the necrotic radius obtained by Freyer and Sutherland (1986a) by culturing EMT6/Ro spheroids under different concentrations of glucose and oxygen in the culture medium.

A c c e p t e d m a n u s c r i p t 2 A simple mathematical model of ATP production

Although the in vivo cell energy metabolism involves glucose, fatty acids and amino acids, we only consider here glucose and lactate, so disregarding other energetic substrates possibly contained in the medium. The basic biochemical reactions of glucose metabolism are depicted in Fig. 1 (see e.g. Stryer, 1988) and are briefly described in the Appendix A, where the relationships (1) and ( 2) below are derived.

In particular, when the intracellular concentrations are constant in time, and under the hypothesis that glucose, lactate and oxygen are the only metabolites exchanged with the extracellular space, the stoichiometry of the metabolic reactions imposes the following relationship among the consumption rates (mol•cell -1 •s -1 ) of these compounds:

f O = 6f G + 3f L , ( 1
)
where f O > 0 and f G > 0 are the consumption rates of O 2 and glucose, respectively, and f L is the consumption rate of lactate. Note that f L can also be negative because, according to the reaction (A.2), there can be a net production of lactate. The production rate of ATP 1 , f AT P , can be expressed in terms of the above rates. Denoting by η 1 and η 2 the efficiency of the phosphorylation of ADP into ATP related to the oxidation of NADH and, respectively, FADH 2 , the relation

f AT P = 2f G + (5η 1 + η 2 + 1)f O /3 ( 2 )
follows (see the Appendix A).

Equations ( 1)-( 2) give the relations among the metabolic fluxes that are implied by the stoichiometry of reactions, but say nothing on the extent of these fluxes. To our aims, moreover, it is necessary to express the consumption/production rates in terms of the extracellular concentrations. Let us consider a cell surrounded by a medium in which glucose, oxygen and lactate have constant concentrations σ G , σ O and, respectively, σ L . The respective intracellular concentrations will be denoted as σG , σO , and σL . We assume: (i) glucose uptake follows a Michaelis-Menten kinetics and all glucose taken up by cells is consumed (Casciari et al., 1992[START_REF] Venkatasubramanian | Incorporating energy metabolism into a growth model of multicellular tumor spheroids[END_REF]; (ii) the exchange of pyruvate with the extracellular medium is disregarded; (iii) transmembrane transport of lactate is non-saturable in the concentration range of interest, and is driven by the gradient of lactate concentration; (iv) because of the high permeability of cell membrane to oxygen [START_REF] Alberts | Molecular Biology of the Cell[END_REF], oxygen concentration inside the cell is equal to the extracellular concentration; (v) the pyruvate flow toward the Krebs cycle is a saturable function of pyruvate and oxygen concentrations (see also [START_REF] Aubert | Interaction between astrocytes and neurons studied using a matematical model of compartmentalized energy metabolism[END_REF].

Experimental data [START_REF] Freyer | A reduction in the in situ rates of oxygen and glucose consumption of cells in EMT6/Ro spheroids during growth[END_REF], Casciari et al., 1992[START_REF] Schroeder | Spatial heterogeneity and oxygen dependence of glucose consumption in R3230Ac and fibrosarcomas of the Fischer 344 rat[END_REF] have shown that glucose consumption increases when oxygen concentration decreases (the so-called Pasteur effect) and that, at physiological pH, oxygen consumption increases when glucose concentration decreases (the so-called Crabtree effect). Expressions for the glucose and oxygen consumption including the above effects were proposed on an empirical basis (Casciari et al., 1992). These effects seem to depend on the pH of cell environment: for 1 Abbreviations: ATP, adenosine triphosphate; ADP, adenosine diphosphate; P i , inorganic orthophosphate; GTP, guanosine triphosphate; GDP, guanosine diphosphate; NAD + , nicotinamide adenine dinucleotide (oxidized form); NADH, nicotinamide adenine dinucleotide (reduced form); FAD, flavin adenine dinucleotide (oxidized form); FADH 2 , flavin adenine dinucleotide (reduced form); Ac-CoA, acetyl-CoA.

A c c e p t e d m a n u s c r i p t

instance, the Crabtree effect disappears in an acidic environment. Although the production of H + ions is related to lactic acid production, the description of pH changes is complicated by the presence of an intracellular buffer (see the mathematical models in Casciari et al., 1992, and[START_REF] Webb | Mathematical modelling of tumour acidity: regulation of intracellular pH[END_REF]. Since a mechanistic representation of these phenomena appears to be rather complex, they are neglected in our model.

For lactate transport, it is known that most of its influx and efflux occur through monocarboxylate transporters (MCTs) which exhibit saturable kinetics [START_REF] Poole | Transport of lactate and other monocarboxylates across mammalian plasma membranes[END_REF]. In Ehrlich tumour cells the Michaelis-Menten constant, K m , of transport kinetics has been found to be 20 mM at 37 o C and pH 7.2 [START_REF] Spencer | L-lactate transport in Ehrlich ascites-tumour cells[END_REF] and, with a different technique, 5.34 mM in the same conditions [START_REF] Carpenter | The kinetics, substrate and inhibitor specificity of the lactate transporter of Ehrlich-Lettre tumour cells studied with the intracellular pH indicator BCECF[END_REF].

In the mammary carcinoma MDA-MB231, high affinity (MCT2) and low affinity (MCT4) transporters have been evidenced [START_REF] Wang | The role of monocarboxylate transporter 2 and 4 in the transport of γ-hydroxybutyric acid in mammalian cells[END_REF], with a prevalent flux through MCT4 (K m = 21.3 mM). The K m values have been found to be highly pH dependent, with K m decreasing with acidity [START_REF] Dimmer | The low-affinity monocarboxylate transporter MCT4 ia adapted to the export of lactate in highly glycolytic cells[END_REF]. So, in the presence of a difference between extracellular and intracellular pH, an asymmetric flux could occur. In the present paper we neglect this possible asymmetry, at variance with our previous paper [START_REF] Bertuzzi | ATP production and necrosis formation in a tumour spheroid model[END_REF], and the high K m values reported above lead us to assume a simple non-saturable kinetics for lactate transport.

On the basis of the above assumptions, the rate of glucose consumption is given by

f G (σ G ) = F G σ G K G + σ G , (3)
where F G is the maximal rate of uptake and K G is the Michaelis-Menten constant. For the intracellular concentrations of pyruvate and lactate we write the following equations:

v dσ P dt = 2f G (σ G ) -ψ(σ P , σL ) -φ P (σ P , σ O ) , (4) v dσ L dt = ψ(σ P , σL ) + v h(σ L -σL ) , ( 5
)
where v is the cell volume, ψ is the net flow of pyruvate to lactate, φ P is the pyruvate flow toward the Krebs cycle, and h is the transfer coefficient of lactate transmembrane transport.

Although the reversible reaction between pyruvate and lactate is an enzymatic reaction, we simply represent it by a first-order reaction, writing

ψ(σ P , σL ) = v(k + σP -k -σL ) , (6)
where k + and k -are the forward and, respectively, the backward rate constants. The complex chain of reactions in the Krebs cycle and the electron transport system that lead to pyruvate oxidation, is summarized in the simplest way according to assumption (v):

φ P (σ P , σ O ) = F P σP K P + σP σ O K O + σ O . ( 7
)
Expression ( 7) modifies the previous expression for φ P adopted in [START_REF] Bertuzzi | ATP production and necrosis formation in a tumour spheroid model[END_REF], where the pyruvate flow toward the Krebs cycle was taken proportional to the product of pyruvate and oxygen concentrations. The expression (7) accounts for the saturative behaviour of enzymatic reactions in the pyruvate oxidation pathway.

A c c e p t e d m a n u s c r i p t

Using Eqs. ( 3)-( 7), the intracellular concentration of pyruvate at the steady state can be derived in terms of the extracellular concentrations. From ( 5)-( 6) we derive:

σL = k + k -+ h σP + h k -+ h σ L . ( 8
)
To write the following equations in a more compact form, we define the functions

M (σ O ) = σ O K O + σ O , N (σ G , σ L ) = 2f G (σ G ) + p 1 σ L ,
and the parameter combinations

p 1 = k -h k -+ h v , p 2 = k + h k -+ h vK P .
From ( 4), taking into account ( 5), ( 7) and ( 8), we then obtain

σP = - K P 2p 2 [p 2 + F P M (σ O ) -N (σ G , σ L )] + K P 2p 2 p 2 + F P M (σ O ) -N (σ G , σ L ) 2 + 4p 2 N (σ G , σ L ) 1/2 . ( 9
)
As seen in Appendix A, the stoichiometry of the reactions at the steady state gives for the oxygen consumption rate, f O , the equation

f O = 3φ P .
Therefore, dropping the arguments of the functions M (σ O ) and N (σ G , σ L ), we can write:

f O (σ G , σ L , σ O ) = 3 2 p 2 + F P M + N -[p 2 + F P M + N ] 2 -4F P MN 1/2 . ( 10
) Note that f O is zero for σ O = 0 and df O dσ O σ O =0 = 3 N p 2 + N F P K O .
For σ O increasing, f O attains the saturation value:

f Omax (σ G , σ L ) = 3 2 p 2 + F P + N -[p 2 + F P + N ] 2 -4F P N 1/2 . ( 11
)
The oxygen consumption rate increases with N , and then with the glucose consumption and lactate uptake, but it is always less than 3F P as expected from (7). When 4F P N (p 2 + F P + N ) 2 and p 2 + N > F P we have 4F P MN (p 2 + F P M + N ) 2 for any σ O > 0, and Eq. ( 10) can be approximated by the expression

f O (σ G , σ L , σ O ) 3 F P MN p 2 + F P M + N ,

A c c e p t e d m a n u s c r i p t

that gives an easier insight of ( 10) and [START_REF] Araujo | A history of the study of solid tumour growth: the contribution of mathematical modelling[END_REF]. Moreover, if p 2 F P + N , the previous equation can be further approximated as

f O (σ G , σ L , σ O ) γN M = γ 2f G (σ G ) + p 1 σ L σ O K O + σ O , ( 12
)
where γ = 3F P /p 2 . The equation for the rate of lactate consumption,

f L (σ G , σ L , σ O ), follows from (1) writ- ten as f L = 1 3 f O -2f G , ( 13
)
with f G and f O defined above. f L may be positive or negative according to the oxygen concentration: when σ O → 0 we will have f L < 0, that is, lactate production.

In view of the expressions of M (σ O ) and N (σ G , σ L ), when the equation ( 10) for f O is adopted, the consumption rates f G , f L and f O depend on the following parameters and parameter combinations: F G , K G , F P , K O , p 1 and p 2 . When the approximation ( 12) is used, the model parameters reduce to F G , K G , γ, K O and p 1 . Figure 2 gives an example of the profile of the oxygen consumption rate as a function of oxygen concentration, according to (10), for an assigned concentration of extracellular lactate and three concentrations of glucose. The behaviour of the maximal value of f O as a function of the extracellular glucose and lactate concentrations is depicted by Figure 3. The profile of lactate production/consumption rate as a function of oxygen concentration is shown in Fig. 4 for two values of extracellular glucose and three values of lactate concentration. In the absence of extracellular glucose (left panel), lactate is consumed with a rate that increases with its concentration. When glucose is present, lactate can be either produced or consumed depending on the extracellular lactate concentration. The right panel exemplifies this behaviour for σ G = 5 μM. However, when oxygen concentration is small enough, lactate is produced in the presence of glucose irrespective of the value of σ L . At

σ O = 0, in fact, it is f L = -2f G .
Experimental evidences have suggested that energy metabolism of tumour cells mainly relies on the glycolytic pathway for the ATP production, thus exhibiting enhanced lactate production and the formation of an acidic environment. These observations led O. Warburg to the notion of the glycolytic phenotype of cancer cells [START_REF] Warburg | On the origin of cancer cells[END_REF]. Recently, however, this view has been questioned by stressing that cells may rely on glycolysis because of the prevailing hypoxic conditions in tumours [START_REF] Guppy | Contribution by different fuels and metabolic pathways to the total ATP turnover of proliferating MCF-7 breast cancer cells[END_REF]. In our model, the balance between the glycolytic and the aerobic ATP production is governed by the parameters k + , F P and K P . From Eqs. ( 2) and ( 10), we see that if the product k + K P is large and F P is small, the aerobic production of ATP is small even in the presence of high oxygen concentrations, and thus the cells would express a glycolytic phenotype. The fraction of the total ATP produced by glycolysis at high oxygen concentration (M (σ O ) = 1) is plotted in Fig. 5 as a function of F P for three values of parameter p 2 . Values of this fraction close to the unity indicate a glycolytic phenotype. The curves show the increase of this fraction, and then the transition to the glycolytic phenotype, as F P decreases and the product k + K P increases.

A c c e p t e d m a n u s c r i p t 3 Development of central necrosis in multicellular spheroids

We consider a spheroid of radius R, in which the nutrient consumption rates of all viable cells depend on the concentrations of the chemicals according to the expressions given in the previous section. Since the growth dynamics of the spheroid is much slower than the process of diffusion of the nutrients, we will consider the diffusion process in spherical symmetry to be at the quasi-steady state at any time. Of course all the derivatives σ i , i = G, L, O, of the extracellular concentrations σ G (r), σ L (r), σ O (r) vanish for r = 0 by symmetry (whether or not a necrotic core is present).

We assume that all cells die when f AT P , as given by ( 2), falls to a threshold value f N , and we suppose that the impairment of the ATP production is the only cause of formation of the necrotic core. According to this hypothesis, the viable region of the spheroid is defined by the inequality

f AT P (r) > f N .
We obviously assume that the concentrations at the boundary R are such that f AT P (R) > f N . In the absence of a necrotic core, the viable region will extend from r = 0 to r = R.

In the viable region the concentrations σ G , σ L , σ O obey the equations

D i 1 r 2 ∂ ∂r r 2 ∂σ i ∂r = nf i , i = G, L, O , ( 14
)
where D i denote the effective diffusivities and n = ν * /v, with ν * the cell volume fraction assumed to be constant, is the number of cells per unit volume. The quantities f i , i = G, L, O, in ( 14) are the consumption/production rates per cell defined in the previous section by Eqs. ( 3), ( 13) and (10) (or ( 12)), respectively. For simplicity, the possible changes in the consumption rates between proliferating and quiescent cells are disregarded. We remark that the values of the diffusivities (see Table 1) guarantee that the contribution of the extracellular liquid flow to metabolite transport is negligible.

Experimental observations that evidenced the presence of concentration gradients in a layer surrounding the surface of tumour spheroids (Mueller-Klieser and [START_REF] Mueller-Klieser | Oxygen tensions in multicell spheroids of two cell lines[END_REF][START_REF] Mueller-Klieser | Method for the determination of oxygen consumption rates and diffusion coefficients in multicellular spheroids[END_REF][START_REF] Groebe | Distributions of oxygen, nutrient, and metabolic waste concentrations in multicellular spheroids and their dependence on spheroid parameters[END_REF], led us to prescribe the concentrations in the culture medium at a fixed distance L from the spheroid surface. These concentrations obey the following equations:

D w i 1 r 2 ∂ ∂r r 2 ∂σ i ∂r = 0 , i = G, L, O , R < r < R + L , ( 15
)
σ i (R + L) = σ * i > 0 , i = G, O , σ L (R + L) = σ * L ≥ 0 , ( 16
)
σ i (R -) = σ i (R + ) , D i σ i (R -) = D w i σ i (R + ) , i = G, L, O , ( 17
)
where D w i denote the diffusivities in the medium, and continuity of concentrations and fluxes is imposed at r = R.

In the presence of a necrotic core, we denote by ρ N the radius of the necrotic region. In this case the viable region reduces to a rim ρ N < r < R. At the necrotic boundary r = ρ N , we will have the conditions

f AT P (ρ N ) = f N , ( 18
)
σ i (ρ N ) = 0, i= G, L, O . ( 19
)

A c c e p t e d m a n u s c r i p t

In the necrotic core the nutrient concentrations are constant and f AT P = 0.

In [START_REF] Bertuzzi | ATP production and necrosis formation in a tumour spheroid model[END_REF], the concentrations σ * i in the culture medium were simply prescribed at the boundary R:

σ i (R) = σ * i > 0 , i = G, O , σ L (R) = σ * L ≥ 0 . ( 20
)
In that case we proved the following results:

• For R sufficiently large, the inequality f AT P (r) > f N cannot hold in the whole interval (0, R).

• For R sufficiently large, there exists one and only one solution σ G (r), σ L (r), σ O (r), ρ N of the free boundary problem defined by Eqs. ( 14), ( 18), ( 19) and ( 20).

Therefore, given the concentrations σ * i and the threshold f N , an increasing function

ρ N = Φ(R) is defined for R > R N > 0, with ρ N = 0 for R = R N ,
where R N is the spheroid radius at which the central necrosis starts to appear.

We note that the existence proofs in [START_REF] Bertuzzi | ATP production and necrosis formation in a tumour spheroid model[END_REF] remain valid even if the consumption rates are not precisely the ones there proposed, but just possess the same qualitative properties as those proposed here. Moreover, it can be easily seen that the above existence results still hold for the free boundary problem considered in the present paper, which is defined by Eqs. ( 14), ( 18), ( 19) and (15-17).

Necrotic core in EMT6/Ro spheroids

The growth characteristics and the metabolism of the EMT6/Ro spheroid (mouse mammary carcinoma cells) have been extensively investigated [START_REF] Mueller-Klieser | Oxygen tensions in multicell spheroids of two cell lines[END_REF][START_REF] Mueller-Klieser | Method for the determination of oxygen consumption rates and diffusion coefficients in multicellular spheroids[END_REF][START_REF] Freyer | A reduction in the in situ rates of oxygen and glucose consumption of cells in EMT6/Ro spheroids during growth[END_REF], Freyer and Sutherland, 1986a[START_REF] Casciari | Glucose diffusivity in multicellular tumor spheroids[END_REF][START_REF] Freyer | Role of necrosis in regulating the growth saturation in multicellular spheroids[END_REF][START_REF] Walenta | ATP concentrations in multicellular tumor spheroids assessed by single photon imaging and quantitative bioluminescence[END_REF][START_REF] Freyer | Cellular energetics measured by Phosphorous Nuclear Magnetic Resonance spectroscopy are not correlated with chronic nutrient deficiency in multicellular tumor spheroids[END_REF], Casciari et al., 1992, Casciari et al., 1992a). In Freyer and Sutherland (1986a) in particular, the radius of the necrotic core for spheroids of different size and under different oxygen and glucose concentrations in the medium were reported. These data have been used to evaluate the predictive capability of our model.

The values of a subset of model parameters were assigned according to the literature data (see Table 1). [START_REF] Freyer | A reduction in the in situ rates of oxygen and glucose consumption of cells in EMT6/Ro spheroids during growth[END_REF] reported, for an exponentially growing monolayer culture of EMT6/Ro cells with 5.5 mM glucose, the following values of the glucose consumption rate: 18 × 10 -17 mol/(cell•s) when oxygen concentration was 0.28 mM and 25 × 10 -17 mol/(cell•s) when oxygen concentration was 0.07 mM. This finding suggested to take for F G the mean of the above values. In the same paper, the oxygen consumption rate per cell is also reported. In particular, for an oxygen concentration of 0.28 mM and glucose concentration of 5.5 mM, the value of 8.3 × 10 -17 mol/(cell•s) has been found. We used this data as a constraint on the value of the unknown parameters. The thickness of the unstirred layer around the spheroid has been found to range from 60 -80 μm [START_REF] Groebe | Distributions of oxygen, nutrient, and metabolic waste concentrations in multicellular spheroids and their dependence on spheroid parameters[END_REF] to 100 -200 μm (Mueller-Klieser, 1984), hence we assumed L = 80 μm.

The other parameters and parameter combinations of the biochemical model and the ATP threshold, f N , were instead determined by minimizing the deviation from the experimental data. Let θ be the vector of the parameters to be estimated. We denote i , Y j i , i = 1, . . . , N j , the i-th measured values of the necrotic core radius, ρ N , and respectively of the spheroid radius, R, in the j-th experiment. By Φ -1 j (ρ N ) we denote the function that gives the spheroid radius, R, for which the radius of the necrotic core is ρ N under the experimental condition j. This function is well defined since the function Φ is increasing. The optimal fitting was achieved by minimizing the index

J(θ) = 8 j=1 N j i=1 (Φ -1 j (X j i ; θ) -Y j i ) 2 . ( 21
)
The choice of the index (21) deserves some comments. A more intuitive choice would be that of fitting the measurements of the necrotic radius by the function Φ defined as

Φ(θ) = 0 , R≤ R N (θ) Φ(R; θ) , R > R N (θ) .
In this way, however, the quadratic deviations for the measurements with R values smaller than R N (θ) have zero derivative with respect to θ, and then do not contribute to the selection of the descent direction in the minimization. This problem has been avoided by using the index ( 21). For the minimization, we used a global-search derivative-free algorithm, reported in [START_REF] Liuzzi | Multi-objective optimization techniques for the design of induction motors[END_REF], that searches the minimum in a rectangular box of the parameter space whose boundaries must be a priori assigned. The function Φ -1 j for a given ρ N was computed according to the steps briefly described in Appendix B.

Firstly, adopting Eq. ( 10) for the oxygen consumption rate, the estimation of the parameters F P , p 1 , p 2 , and of the threshold f N , was attempted. The minimization was performed under the constraint that the prediction, q O (θ), of the oxygen consumption rate for σ G = 5.5 mM, σ L = 0, σ O = 0.28 mM met the measured value Z O = 8.3 × 10 -17 mol/(cell•s) [START_REF] Freyer | A reduction in the in situ rates of oxygen and glucose consumption of cells in EMT6/Ro spheroids during growth[END_REF] within a prescribed tolerance. However, the parameters F P , p 1 , p 2 could not be reliably estimated, as different parameter sets were found to give similar fittings. In particular, the estimated values of p 2 were large enough to make approximation (12) feasible, so we used that approximation to fit the data. With the approximation (12) the estimation problem is considerably simplified since the quantity q O only depends on the unknown parameter γ. By setting q O = Z O , we obtained for γ the value 0.198 (see Table 2). Therefore, the only parameters to be estimated were p 1 and f N , and the fitting so obtained was found to be substantially equal to that obtained when using Eq. ( 10). The covariance matrix of the estimates was computed by the linearization method that yields

Cov( θ) = s 2 dg dθ T θ dg dθ θ -1
, where θ is the vector of the estimated parameters, g(θ) is the vector with components Φ -1 j (X j i ; θ) for all i and j, and s 2 is the estimate of the (common) variance of the measurements, given by

s 2 = J( θ) 8 j=1 N j -m ,

A c c e p t e d m a n u s c r i p t

with m = 2 being the number of unknown parameters.

The optimal prediction of the radius of the necrotic core as a function of the spheroid radius is shown in Fig. 6, together with the experimental data. It can be noted that the model reasonably predicts the clear dependence of the necrotic radius on the oxygen value in the medium at glucose 5.5 mM, and the small sensitivity to oxygen when glucose attains the extreme values of 0.8 and 16.5 mM. Less satisfactory is the fitting of the data at glucose 1.7 mM, particularly for the oxygen concentration of 0.07 mM when the data at glucose 0.8 and 1.7 mM are practically indistinguishable while the model exhibits a monotone dependence on σ * G . The optimal parameter values are reported in Table 2. As a comment to these values we note the following. The experimental data of lactate transport in Ehrlich tumour cells [START_REF] Carpenter | The kinetics, substrate and inhibitor specificity of the lactate transporter of Ehrlich-Lettre tumour cells studied with the intracellular pH indicator BCECF[END_REF] indicate, for the ratio between maximal transport rate and Michaelis constant, the value 0.13 s -1 . Assuming that the coefficient h in Eq. ( 5) is equal to this ratio, we obtain an upper limit of p 1 , given by hv, equal to 37.0×10 -14 liter/(cell•s), having taken v = 2.86×10 -12 liter/cell (v vas computed as ν * /n taking n as in Table 1 and ν * = 0.6). Thus, the obtained value of p 1 satisfies the above upper bound. The ATP production rate at the surface of spheroids cultured with 5.5 mM glucose and 0.28 mM oxygen was computed, giving a range from 85 to 95×10 -17 mol/(cell•s) at different R values. Thus the estimated value of the threshold f N is close to half of this value.

Figure 7 reports the predicted profiles of metabolite concentrations and of the ATP production rate across a spheroid of radius 400 μm cultured with 5.5 mM glucose and 0.28 mM oxygen. The metabolite concentrations in the unstirred layer are also shown. A marked increase in ATP production rate in the middle of the viable rim can be noted. This maximum is due to the increase in the oxygen consumption rate caused by the increase in lactate concentration, see the definition of N (σ G , σ L ). The presence of a maximum of f AT P internal to the viable rim matches the existence of an internal maximum of ATP concentration in EMT6/Ro spheroids as measured by bioluminescence [START_REF] Walenta | ATP concentrations in multicellular tumor spheroids assessed by single photon imaging and quantitative bioluminescence[END_REF]. The mean ATP production rate per cell over the viable rim was also computed as

f AT P = 3 R 3 -ρ N (R) 3 R ρ N (R) r 2 f AT P (r) dr .
For R increasing from 250 to 1000 μm, the mean production rate was found to slightly decrease with R, ranging from 120 to 114 × 10 -17 mol/(cell•s).

Concluding remarks

The model we propose appears to give a reasonable prediction of the radius of the necrotic core as measured in EMT6/Ro spheroids cultured in different conditions of oxygen and glucose (Freyer and Sutherland, 1986a). The main feature of the model is the assumption that the onset of necrosis is caused by an insufficient ATP production rate and this hypothesis, which is controversially discussed in the literature [START_REF] Walenta | Metabolic imaging in multicellular spheroids of oncogene-transfected fibroblasts[END_REF], appears thus in agreement with the data considered. Attempts to predict these data by using a threshold on a single nutrient, oxygen or glucose, gave far worse results (not shown). Nuclear magnetic resonance measurements in EMT6/Ro spheroids showed that high energy phosphates remained at a normal steady level irrespective of the stage of spheroid growth A c c e p t e d m a n u s c r i p t [START_REF] Freyer | Cellular energetics measured by Phosphorous Nuclear Magnetic Resonance spectroscopy are not correlated with chronic nutrient deficiency in multicellular tumor spheroids[END_REF]. Under the assumption that the ATP content per cell is proportional to the ATP production rate, our model meets this observation as it predicts that the mean ATP production rate per cell is only slightly decreasing with the spheroid radius.

A limitation of the model is that the presence in the spheroids of quiescent cells (Freyer and Sutherland, 1986b), possibly with different rates of nutrient consumption, is disregarded. In principle, the model could be extended to take into account a quiescent cell population. However, since the transition to quiescence should be related to the availability of nutrients or, more strictly, to the energy production, this extension would require additional assumptions and parameters thus making the model validation more uncertain.

As mentioned above, experimental investigations on the EMT6/Ro cell line have shown the presence of the Pasteur and Crabtree effects, that is, the dependence of glucose and, respectively, oxygen consumption rates on the concentration of the other metabolite. The complex mechanisms that underlie these phenomena and their dependence on pH, however, have prevented us from including these effects in a relatively simple model.

Modelling the pH profile in the spheroid would thus be a valuable improvement, even in view of the possible role of acidity in inducing cell death, that suggests different hypotheses on the formation of the necrotic core. The relevance of pH has been pointed out, in fact, not only in connection with the invasive capacity of the tumour [START_REF] Gatenby | A reaction-diffusion model for cancer invasion[END_REF][START_REF] Smallbone | The role of acidity in solid tumour growth and invasion[END_REF], Fasano et al., 2009), but also in connection with the formation of the tumour necrotic regions (Smallbone et al., 2005, Bianchini andFasano, 2009).

A c c e p t e d m a n u s c r i p t

In the above reaction, that refers to the oxidation of 1 glucose molecule, 2 NADH molecules derive from the pyruvate→acetyl-CoA step, 6 NADH from the Krebs cycle, and the remaining 2 molecules from the glycolytic step. We note that this reaction, together with (A.2), provides all the NAD + and FAD molecules required in reactions (A.1) and (A.3). The reactions in the ETS yield the energy necessary for the phosphorylation of ADP into ATP. The efficiency of this process appears to vary depending on cell type and cell condition.

The transformation of 1 NADH molecule in NAD + provides energy for the formation, in average, of about 2.5 ATP molecules, whereas 1 FADH 2 generates about 1.5 ATP molecules [START_REF] Hinkle | Mechanistic stoichiometry of mitochondrial oxidative phosphorylation[END_REF]. In addition, we remark that the GTP molecule produced in the Krebs cycle is equivalent to an ATP molecule.

Let us consider a cell at the steady state and assume that glucose, lactate and oxygen are the only metabolites exchanged with the extracellular space. Let φ G , φ P , φ + L , and φ - L denote the fluxes (mol

f O = 5 2 φ P + φ G - 1 2 φ + L + 1 2 φ - L .
As the flux balance at the steady state give

φ P = 2φ G -φ + L + φ - L , (A.5) it follows that f O = 3φ P . (A.6)
Moreover, taking into account that

f G = φ G , f L = φ - L -φ + L ,
from (A.5) and (A.6) we obtain equation (1):

f O = 6f G + 3f L .
Denoting by η 1 and η 2 the mean number of ATP molecules generated as a consequence of the oxidation of 1 NADH and 1 FADH 2 molecules, respectively, from reactions (A.1)-(A.3) we can write the rate of ATP production, f AT P , as

f AT P = 2φ G + φ P + η 1 (4φ P + 2φ G -φ + L + φ - L ) + η 2 φ P .
Note that the second term in the right-hand side represents the production of GTP, considered equivalent to ATP. Since φ P = f O /3, the above equation can be rewritten giving equation (2): The computation of the function Φ -1 j for a given ρ N was performed according to the steps briefly described in the following. Preliminarily, each equation in ( 14) is rewritten as a system of two first-order differential equations for σ i and σ i .

f AT P = 2f G + (5η 1 + η 2 + 1)f O /3 .
1. For any experimental value of ρ N , fix an initial guess for R = Φ -1 j (ρ N ) equal to (ρ N + R max )/2, with R max a suitable large value.

2. Find the values σi , i = G, L, O, such that the Cauchy problem defined by Eqs. ( 14), ( 15), ( 17), with the initial conditions σ i (ρ N ) = σi and ( 19), has a solution such that σ i (R + L) = σ * i . This "shooting" step is performed in the following way.

2a. Given a triple σi , i = G, L, O, integrate Eqs. ( 14) with conditions σ i (ρ N ) = σi and ( 19) on [ρ N , R], and continue the solutions up to the radius R + L according to ( 15) and ( 17).

2b. Minimize with respect to σi the index

I(σ G , σL , σO ) = i=G,L,O (σ i (R + L) -σ * i ) 2 .
The minimization is accomplished by the Levemberg-Marquardt algorithm.

3. Compute the ATP production rate at ρ N , f AT P (ρ N ; R), using the σi values obtained at step 2b.

4. Find the solution with respect to R of the equation f AT P (ρ N ; R) = f N by means of the bisection method, iterating steps 2 and 3 a suitable number of times (n=30). Figure 2. Oxygen consumption rate at σ L = 1 mM according to Eq. ( 10). Parameter values: F G , K G and K O as in Table 1, p 1 = 10 • 10 -14 liter/(cell•s), p 2 = 500 • 10 -17 mol/(cell•s), F P = 50 • 10 -17 mol/(cell•s). 2, other parameters reported in Table 1.

A c c e p t e d m a n u s c r i p t

  j = 1, . . . , 8, the different conditions in which the spheroids were grown inFreyer and Sutherland (1986a): j = 1, . . . , 4 refer to experiments with σ * O = 0.07 mM and σ * G equal to 0.8, 1.7, 5.5 and 16.5 mM, respectively; j = 5, . . . , 8 indicate experiments with σ * O = 0.28 mM and σ * G as above. For all the experiments we assumed σ * L = 0. Let us define by X j

Figure 1 .

 1 Figure 1. Simplified scheme of glucose metabolism. The inner box represents the mitochondria.

Figure 3 .

 3 Figure 3. Maximal value of the oxygen consumption rate as a function of σ G and σ L . Parameter values as in Fig. 2.

Figure 4 .

 4 Figure 4. Lactate consumption/production rate as a function of oxygen concentration for three values of lactate concentration. Left panel, σ G = 0; right panel, σ G = 5 μM. Other parameters as in Fig. 2.

Figure 5 .

 5 Figure 5. Fraction of ATP production (computed with σ G = 5 μM, σ O = 0.28 μM and σ L = 1 μM) due to the glycolytic step as a function of F P for three values of p 2 . Other parameters as in Fig. 2.

Figure 6 .

 6 Figure6. Optimal prediction of the radius of the necrotic core as a function of the spheroid radius, together with the experimental measurements reported inFreyer and Sutherland, 1986 (courtesy of Dr. J.P. Freyer). Left panel, oxygen concentration in medium 0.07 mM; right panel, oxygen concentration in medium 0.28 mM. Glucose concentrations are reported in the inset. The estimated values of the parameters p 1 , γ and f N are given in Table2, other parameters reported in Table1.

Figure 7 .

 7 Figure 7. Predicted profiles of glucose, lactate and oxygen concentrations and of ATP production rate as a function of the radial distance in a spheroid of radius R = 400 μm cultured with 5.5 mM glucose and 0.28 mM oxygen. The concentrations in the unstirred layer (L = 80 μm) are shown. The necrotic radius is ρ N = 209.38 μm. The values of p 1 , γ and f N are given in Table2, other parameters reported in Table1.

  Figure 1:

  Figure 2:

  Figure 4:

  Figure 6:

  •cell -1 •s -1 ) of glucose to pyruvate, pyruvate to acetyl-CoA, pyruvate to lactate and, respectively, lactate to pyruvate. Recalling (A.1)-(A.2), the net production rate of NADH in the reactions Glucose → Pyruvate Lactate is given by 2φ Gφ + L + φ - L , because when pyruvate is converted into lactate, 1 NADH molecule is consumed. Since the cumulative production rate of NADH and FADH 2 from pyruvate oxidation is 5φ P , and each NADH or FADH 2 molecule requires (1/2)O 2 to be oxidized (see reaction (A.4)), for the oxygen consumption rate, f O , we have:

Table 1 :

 1 Model parameter values

	Symbol Value	Source
	F G	21.5 × 10 -17 mol/(cell•s) assumed
	K G K O	0.04 mM 4.64 × 10 -3 mM	(Casciari et al., 1992) (Casciari et al., 1992)
	η 1	2.5	(Hinkle et al., 1991)
	η 2 D G D L D O D w G D w L D w O n	1.5 1.05 × 10 -6 cm 2 /s 1.78 × 10 -6 cm 2 /s 1.82 × 10 -5 cm 2 /s 9.25 × 10 -6 cm 2 /s 1.20 × 10 -5 cm 2 /s 3.0 × 10 -5 cm 2 /s 2.1 × 10 8 cells/cm 3	(Hinkle et al., 1991) (Casciari et al., 1988) (Casciari et al., 1992) (Mueller-Klieser, 1984) (Li, 1982) (Swabb et al., 1974) (Thews and Hutten, 1983) (Freyer and Sutherland, 1986a)
	L	80 μm	assumed

Table 2 :

 2 Parameter estimates and coefficients of variation (CV)

	Symbol	Estimate	CV (%)
	γ	0.198	-
	p 1 f N	18.78 × 10 -14 liter/(cell•s) 42.82 × 10 -17 mol/(cell•s)	19.91 0.90

  Table 2, other parameters reported in Table 1.

		Cell membrane
	Glucose
	ATP
	Pyruvate	Lactate
	A c c e p t e d Krebs cycle GTP NADH m FADH 2 CO 2 O 2 H 2 O ETS a n u s c r i p t ATP Ac-CoA NADH

Acknowledgements

The authors wish to thank Dr. G. Liuzzi for the precious help given in solving the optimization problem. The authors also acknowledge the work done by Drs. P. Nanci and P. Pasanisi in the software development. Finally, we thank the anonymous Referees for their constructive comments that allowed us to improve the paper.

Appendix A. Glucose energy metabolism: main biochemical reactions and metabolic fluxes

As shown by Fig. 1 glycolysis is the first step of glucose metabolism, in which a sequence of enzymatic reactions converts 1 molecule of glucose into 2 molecules of pyruvate with the simultaneous net production of 2 ATP molecules and the conversion of 2 NAD + into 2 NADH molecules:

Pyruvate may either be reversibly converted into lactate or undergo complete oxidation in the Krebs cycle and the associated electron transport system (ETS). The former reaction is given by Pyruvate

Pyruvate oxidation occurs through multiple steps, but the net reaction may be written as follows:

The NADH produced in the glycolytic step, the pyruvate→acetyl-CoA step and the Krebs cycle, together with the FADH 2 produced in the Krebs cycle, are oxidized in the so-called electron