Babatunde A Ogunnaike
email: ogunnaike@udel.edu

Claudio A Gelmi
email: cgelmi@ing.puc.cl

Jeremy S Edwards
email: jsedwards@salud.unm.edu

A probabilistic framework for microarray data analysis: Fundamental probability models and statistical inference

Keywords: Mixture models, Poisson distributions, Gamma distributions, Beta distributions, Gene Expression

Gene expression studies generate large quantities of data with the defining characteristic that the number of genes (whose expression profiles are to be determined) exceed the number of available replicates by several orders of magnitude. Standard spot-by-spot analysis still seeks to extract useful information for each gene on the basis of the number of available replicates, and thus plays to the weakness of microarrays. On the other hand, because of the data volume, treating the entire data set as an ensemble, and developing theoretical distributions for these ensembles provides a framework that plays instead to the strength of microarrays. We present theoretical results that under reasonable assumptions, the distribution of microarray intensities follows the Gamma model, with the biological interpretations of the model parameters emerging naturally. We subsequently establish that for each microarray data set, the fractional intensities can be represented as a mixture of Beta densities, and develop a procedure for using these results to draw statistical inference regarding differential gene expression. We illustrate the results with experimental data from gene expression studies on Deinococcus radiodurans following DNA damage using cDNA microarrays.

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

Standard statistical analysis involves estimating (and drawing statistical inference about)

p parameters using a data sample of size n, typically with n >> p. The precision of the estimates, as well as the power of the accompanying statistical tests of significance, are well-known to improve with increasing n. The converse is also true: reducing sample size n reduces estimation precision and the power of statistical hypotheses (other things being equal). The problem with microarray data analysis is well-documented: p (in this context the expression levels of thousands of genes to be estimated) is very large, while n, the number of independent replicates, from which these estimates are to be determined, is quite small. Many techniques have been proposed for extracting information about gene expression levels from microarray data, but this intrinsic characteristic of the data set (large p, small n) remains a major impediment (see, for example, [START_REF] Nadon | Statistical issues with microarrays: processing and analysis[END_REF].

To complicate matters further, genes are known to operate, not independently as individuals, but in networks, suggesting that microarray data sets represent the collective behavior of a population best studied jointly. Still, the inference problem of primary concern is the quantitative determination of the differential expression level of individual genes represented in the data set:

which ones are over-expressed, which are under-expressed-and by how much-and which are unchanged. None of the preceding factors change this primary focus, but they motivate us to consider analytical techniques that take into consideration the intrinsic characteristics of microarray data we have highlighted here. Our perspective in this paper is that by examining the entire data set as an ensemble and characterizing it as such, we are able to develop an analysis technique that actually plays to the strength of microarray technology. The development of techniques for spot-by-spot analysis of microarray data has been the subject of much research effort; and many of the key results and the still unresolved issues have been discussed in several reviews [START_REF] Allison | Microarray data analysis: from disarray to consolidation and consensus[END_REF][START_REF] Dharmadi | DNA microarrays: Experimental issues, data analysis, and application to bacterial systems[END_REF][START_REF] Nadon | Statistical issues with microarrays: processing and analysis[END_REF][START_REF] Quackenbush | Computational analysis of microarray data[END_REF][START_REF] Sebastiani | Statistical challenges in functional genomics[END_REF].

A multitude of techniques have been proposed for testing for differential gene expression, each incorporating varying degrees of statistical rigor, ranging from the very simple (for example, based on ratios, or the so-called "fold-change" criterion), to the more sophisticated and computationally demanding approaches. One of the simplest, and most intuitive techniques for analyzing microarray data uses the "fold change" criterion to decide if the changes observed in a gene's mRNA-expression level between two distinct experimental conditions are truly significant or not. It is based on the ratio of measured expression levels, with significance indicated when this ratio (the so-called "fold change") exceeds a pre-specified threshold, typically 2.0. "Fold-change" analysis is an example of a gene-by-gene approach where expression data for genes on a microarray are treated and analyzed as unrelated individual measurements. Analyses based on ANOVA [START_REF] Cui | Statistical tests for differential expression in cDNA microarray experiments[END_REF][START_REF] Kerr | Experimental design for gene expression microarrays[END_REF] and t-test [START_REF] Cui | Statistical tests for differential expression in cDNA microarray experiments[END_REF]Troyanskaya et al., 2002) are other examples of gene-bygene approaches.

An attractive alternative to the gene-by-gene approaches is provided by the class of techniques based on mixture models where the genes represented on the microarray are considered as consisting of two or more populations. In doing so, these techniques take advantage of the vast number of genes on any given array and explicitly recognize the fact that these genes are not A c c e p t e d m a n u s c r i p t its biological neighbor. Examples of mixture model approaches include mixture models for gene effect [START_REF] Lee | Importance of replication in microarray gene expression studies: statistical methods and evidence from repetitive cDNA hybridizations[END_REF]; empirical Bayes analysis [START_REF] Efron | Empirical Bayes analysis of a microarray experiment[END_REF][START_REF] Kendziorski | On parametric empirical Bayes methods for comparing multiple groups using replicated gene expression profiles[END_REF]Newton et al., 2001); fits of various probability densities to intensity ratios (Ghosh and Chinnaiyan, 2002), or to raw intensities [START_REF] Hoyle | Making sense of microarray data distributions[END_REF][START_REF] Steinhoff | Gaussian mixture density estimation applied to microarray data[END_REF]; and hierarchical mixture methods [START_REF] Ghosh | Mixture models for assessing differential expression in complex tissues using microarray data[END_REF][START_REF] Newton | differential variability of expression ratios: improving statistical inference about gene expression changes from microarray data[END_REF]. Also several authors have used Beta distributions to analyze microarray data [START_REF] Allison | A mixture model approach for the analysis of microarray gene expression data[END_REF][START_REF] Baggerly | Identifying differentially expressed genes in cDNA microarray experiments[END_REF][START_REF] Ji | Applications of beta-mixture models in bioinformatics[END_REF], mainly because of its mathematical convenience, not for reasons of biological plausibility. These methods are all based on empirical approximations fitted to observed data, with little or no possibility for mechanistic interpretation of the resulting distributions.

The objective of this paper is to present first-principles theoretical (as opposed to merely empirical) results for representing and characterizing microarray data as an ensemble. These results are then assembled into a theoretical framework for: (i) appropriately analyzing microarray data sets on the basis of this ensemble characterization; and (ii) drawing realistic inferences from the analysis. The proposed theoretical framework is them illustrated with an example experimental data set and in simulation. The experimental validation of the complete inference procedure will be presented in an upcoming publication.

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

Intensity measurements and gene expression levels

For each gene represented by a spot on a microarray, the expression level O, corresponding to the number of cDNA molecules, is considered fixed, but unknown. It is estimated by measurements of fluorescent dye intensity.

For a fixed expression level, repeated fluorescent signal intensity measurements taken n separate and independent times, I j ; j = 1,2,...,n, will differ because of inherent variability. For fluorescent intensity measurements acquired from cDNA microarrays, it is typical to represent this fact as follows [START_REF] Sebastiani | Statistical challenges in functional genomics[END_REF]:

O j j a I (1)
where the uncertainty in the measurement is considered as having arisen primarily from a j being a random proportionality constant, characteristic of the measurement device, relating the number of cDNA molecules to the observed intensity. In reality, it is more accurate to represent the raw intensity signal as:

j j j b a I O ' (2)
where b j is another random component corresponding to background noise; i.e. the intensity measurement values one would obtain even when there are no cDNA molecules present on the spot [START_REF] Rocke | A model for measurement error for gene expression arrays[END_REF]. Because of the non-negativity of intensity measurements, b j has a non-zero mean value that is typically subtracted off to obtain the background-corrected In what follows, we assume that the intensity signal has been appropriately background corrected so that Eq. ( 1) applies. Nevertheless, we will still pay due respect to the intensity magnitude.

Model development

Let I represent the observed fluorescent intensity corresponding to a total of z cDNA molecules; furthermore, we consider this intensity to be cumulative in the following sense:

z I I I I ! 2 1 (3)
where I (m) is the intensity contribution from the m th molecule. We now refer to z(I) as the number of molecules corresponding to the observed (background-corrected) total intensity measurement I. From this perspective, for a given observed I, z is an unknown random variable, an integer count.

The complementary element in this model development is the fluorescent intensity measurement obtained from k molecules, which, for purposes of this derivation, we will refer to as y(k). From this intensity measurement perspective, given a fixed number k, y is an unknown (continuous) random variable.

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

= 1/E, represents the mean number of molecules per unit observed intensity. These quantities, characteristic of the measurement device (akin to a calibration factor), can be used to characterize the random behavior of both the continuous y and the discrete z and relate one to the other as we now show.

The problem at hand is to develop a probability model to characterize y as a function of k, the number of molecules responsible for the observed intensity signal measurement. We adopt a 2step strategy: (i) begin with the complementary problem of characterizing z(I); the discrete count of the number of molecules corresponding to the observed intensity, and (ii) convert the integer count to the continuous variable, y.

Characterizing z(I).

First-principles analysis and experimental confirmation [START_REF] Ozbudak | Regulation of noise in the expression of a single gene[END_REF]Thattai and van Oudenaarden, 2001) have established that the steady state statistical properties of mRNA expression levels results in a Poisson distribution for z, i.e.

! z e z f z T T (4)
where T is the mean number of molecules. Within the context of our current formulation, in relation to the observed intensity I, and given the characteristic parameter K as described above, Eq. (4) becomes: which is now the probability model for characterizing z, the number of mRNA molecules in terms of I, the corresponding observed total intensity. Equation ( 5) provides an expression for computing the probability that z, the number of molecules responsible for generating the observed intensity I, takes on the integer values 0, 1, 2,…, f given the population parameter K; It is important to note that Newton et al. (2001), invoked this model, strictly for convenience, but with a constant shape parameter. Also, [START_REF] Sebastiani | Statistical challenges in functional genomics[END_REF], alluded to the fact that the Gamma distribution may be a more appropriate model for intensity data since a logarithmic transformation still leaves the data with unsatisfactory asymmetry.

A theoretical model of differential expression

At each spot on the microarray, the expression level of each represented gene is evaluated under two different conditions: 1 = "Test" conditions, 0 = "Reference" conditions (one labeled with Cy3 dye, the other with Cy5 dye). Let I i1 and I i0 be the fluorescent signal intensities measured from each spot i on a microarray, with I i1 obtained from the gene in question under test conditions, and I i0 from the gene under control conditions. Let O i1 and O i0 respectively represent the corresponding actual, but unknown, expression levels of the gene in question.

Then from Eq. (1), we have:

1 1 1 i i i a I O (8) 0 0 0 i i i a I O (9)
For various well-documented reasons [START_REF] Dudley | Measuring absolute expression with microarrays with a calibrated reference sample and an extended signal intensity range[END_REF][START_REF] Sebastiani | Statistical challenges in functional genomics[END_REF] [START_REF] Brenner | Gene expression analysis by massively parallel signature sequencing (MPSS) on microbead arrays[END_REF][START_REF] Mikkilineni | Digital quantitative measurements of gene expression[END_REF]) differential expression determination via standard microarray technology remains an important aspect of gene expression studies.

The customary measure of differential expression is the so-called fold change statistic defined as the ratio:

0 1 i i i O O U (10)
and the common practice is to base inference about U i on the measured intensity signal ratio r i defined as

0 1 / i i i I I r ( 11
)
It is important to observe from here however, that such a strategy is based on the (tacit) assumption that for both the Cy3 label signal intensity as the with the Cy5, the multiplicative factors are approximately the same. This may not be the case in general, especially since we have just shown that these multiplicative constants are related directly to the "mean intensity observed per molecule" a number that is clearly different for Cy3 and Cy5 fluorescent dye labels (This is potentially one source of the well-documented dye bias). Dye-swaps and data normalization are some of the tactics typically employed to minimize these problems.

Also, observe that if 2 a V is the variance associated with the random multipliers, then the observed variance in the signal intensity will be given by:

2 2 2 a I V O V (12)
A c c e p t e d m a n u s c r i p t

(and hence with intensity); or equivalently, the standard deviation varies linearly with intensity.

This is in agreement with actual microarray data [START_REF] Rocke | A model for measurement error for gene expression arrays[END_REF]. For example, Figure 1 shows a plot of experimental data from gene expression studies in Deinococcus radiodurans following DNA damage [START_REF] Tanaka | Analysis of Deinococcus radiodurans's transcriptional response to[END_REF] where the data variance is observed to increase with increasing intensity values. The implication for data analysis is that the statistical character of microarray measurements (especially measurement variance) is dependent on signal magnitude [START_REF] Chen | Ratio-based decisions and the quantitative analysis of cDNA microarray images[END_REF]Newton et al., 2001).

The most commonly employed strategy for addressing the preceding issue is a logarithmic transformation of the data to obtain:

0 0 0 log log log i i i a I O (13)
turning the expression into the more familiar additive error form:

i i i X Y H (14)
In most cases, it is now assumed that log(a i0 ) is normally distributed, and the analysis proceeds along these lines. However, it has been shown by various investigators [START_REF] Rocke | A model for measurement error for gene expression arrays[END_REF][START_REF] Sebastiani | Statistical challenges in functional genomics[END_REF] that this assumption, even though convenient, is not necessarily valid. It would be valid if a lognormal density is an appropriate model for the intensity distribution.

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

13 However, we have just established that the Gamma density is more appropriate, not the lognormal.

An alternative transformation

The fan-shaped characteristics of raw cDNA microarray intensity data (illustrated by the example in Figure 1) suggests a strategy where the original data in the form of the ordered pair of measured signal intensities (y i1 , y i0 ) for each spot i is converted from this "Cartesian" representation to the corresponding "polar" form ( R i , x i ) where R i , the vector magnitude defined by:

2 0 2 1 i i i y y R (15)
is the "intensity magnitude", and x i , the "fractional intensity" defined by: 0 1

1 i i i i y y y x ( 16
)
The following are important attributes of this coordinate transformation: (i) it naturally compensates for the inherent heteroskedasticity of microarray data since pie-shaped regions in Cartesian coordinates naturally map into rectangular regions in polar coordinates. For example,

Figure 2 shows the effect of this transformation on the data in Figure 1; (ii) it normalizes the data efficiently without losing any of the original information. Observe that while x i is dimensionless and naturally scaled between 0 and 1, R i retains the information about intensity magnitude, thus Even though the traditional intensity ratio r is related to x i according to the expression.

i i i i i x x y y r 1 0 1 (17)
we recommend x i as a more convenient variable on which to base inference about differential expression for the following reasons: a probabilistic model can be derived for x i from the model (as we will show shortly) while R i provides a rational metric for partitioning the data (allowing us to separate low magnitude data from higher magnitude data).

Probability model for fractional intensity

If we represent as y i1 , the measured intensities, corresponding to gene i (on spot i) under test conditions, and y i0 , the value obtained for the same gene i under reference conditions, then the implications of the results given above are that:

1 1 1 , ~E D i i Gamma y (18) 0 0 0 , ~E D i i Gamma y ( 19
)
allowing the equipment calibration parameter (mean intensity per molecule), E, to be dye dependent, but not gene dependent, so that we need only retain the subscripts 0 and 1 but not i.

A well-known result from probability theory (e.g., [START_REF] Evans | Statistical distributions[END_REF] is that if Eqs. ( 18) and ( 19) are true then:

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

1 , / 1 1 1 i i Gamma y D E (20) 1 , / 0 0 0 i i Gamma y D E (21)
and the fractional variable, i x ~ defined by

0 0 1 1 1 1 / / / ~E E E i i i i y y y x 0 1 1 i i i wy y y ( 22
)
where

0 1 E E w ( 23
)
possesses a Beta distribution, Beta(D j1 , D j0 ), with a probability density function given by:

1 1 0 1 1 0 0 1 1 * * * i i x x x f i i i i D D D D D D (24)
Some of the consequences of these results for practical data analysis are now discussed.

Nominal data transformation

A subtle but important point about Eqs. ( 22) and ( 23) is that raw data consists of y i1 , and y i0 pairs only; the parameters E 1 and E 0 are unknown calibration parameters. In fact, this latter point is the reason why absolute estimates of expression levels are not possible from A c c e p t e d m a n u s c r i p t 16 microarray data. As noted earlier, the rationale for using the intensity ratios as the basis for drawing inference about differential expression is predicated on the tacit assumption that:

1 0 E E (25)
We will designate these customary assumptions as representing what we refer to as nominal conditions. Note that in reality these conditions may in fact not hold. The implications of departure from these nominal conditions, how to detect and correct for such departures will be discussed shortly.

For now, note that under these conditions, w = 1 and Eq. ( 22) indicates that the nominal fractional intensity defined by: 0 1

1 i i i i y y y x ( 26
)
obtained directly from raw data, possesses a Beta distribution, so that for each gene, the probability model is as shown in Eq. ( 24).

When nominal conditions do not hold, then in reality the computed fractional intensity x i ought to be corrected to obtain i x ~. The relationship between the corrected variable and the nominal fractional intensity x i is easily established as the following dual expression: When nominal conditions do not hold, (i.e. w z 1) it is necessary to estimate w, and this is done as follows. Under conditions when there is no change in gene expression, the expected value of the nominal fractional intensity data, E(X) = P 0 should be theoretically equal 0.5.

i i i i x w x x x 1 ~ (27)
However, when w z 1 in ( 22), even though there is no change in gene expression, E(X) z 0.5, and hence any significant deviation of estimates of P 0 from 0.5 provides an indication that nominal conditions do not hold. Given a value of E(X) = P 0 z 0.5, we are able to use Eq. ( 28)

to estimate the value of w required so that the expectation of the corrected intensity i x ~ is 0.5. It is easy to show that under these conditions:

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

The theoretical implications of the results presented in Eqs. ( 24) and ( 26) are as follows: for a particular microarray dataset involving a total of N genes, each nominal fractional intensity x i = 1,2,…,N is a random sample from its own (possibly unique) Beta distribution. Consequently, f(x), the underlying theoretical distribution of (nominal) spot fractional intensities for the entire collection of genes, will be a mixture of overlapping Beta density functions consisting, in principle, of as many Beta densities as there are genes.

Because of the finite support [0,1] over which f(x) is defined and the large number of genes involved, but most importantly as a result of the necessary "coalescence" of the distributions of genes with theoretically identical-or statistically similar-expression levels, in practice, the actual number of distinct and truly discernible Beta densities required to describe any specific collection of microarray data will be relatively few. Practically therefore, this suggests a characterization strategy in which the histogram of the (nominal) fractional intensities is modeled as a mixture of just a few Beta densities.

In particular, the distributions of all genes showing no differential expression will coalesce into one clearly discernible distribution, f 0 (x). If, at the barest minimum, the distributions of genes showing lower differential expression also coalesce into a single discernible distribution and the same happens for the genes showing higher differential expression, we may thus characterize the histogram of fractional intensity data (at the barest minimum) as a mixture of (at least) three Beta distribution functions, f 1 (x) for the data from the collection of genes showing lower differential expression, f 0 (x) for those showing no differential expression, and f 2 (x) for those showing higher differential expression (Eq. ( 31)):

x f x f x f x f 2 2 0 0 1 1 I I I (31)
with 1 2 1 0

I I I .
Alternatively, and more generally, this exercise may be considered as akin to an expansion of f(x) in a set of basis functions, f i (x), that are beta densities. The total number of basis functions used to represent the true function depends on the level of desired "accuracy". The specific choice of 3 basis functions in Eq. ( 31) is based on the arguments advanced above: each of these functions can be related directly to the status of the gene in question.

Thus Finally, we note that [START_REF] Diaconis | Quantifying prior opinion[END_REF], have showed that any data distribution on the interval [0, 1] can be modeled as a finite mixture of Beta distributions, a theoretical result whose validity has been illustrated in a wide range of practical applications [START_REF] Allison | A mixture model approach for the analysis of microarray gene expression data[END_REF][START_REF] Parker | Identifying Important Results From Multiple Statistical Tests[END_REF][START_REF] Skliar | Interaction of 2,2,6,6tetramethyl-3,5-heptanedione with the Si(100)-2 x 1 surface: Scanning tunneling microscopy and density functional theory study[END_REF]. In the current context, the primary implication of this result is that the probabilistic framework we have presented thus far does in fact extend beyond cDNA microarray data: data from Affymetrix GeneChips, and any other two-channel or onechannel array technology (old or new), once expressed in the form of fractional intensities, can also be analyzed using the results presented here.

Statistical Inference: Probability of expression status

The mixture model in Eq. ( 31) can now be used as the reference distribution for drawing rigorous statistical inference about the expression status of each gene as follows. The probabilities that each gene g i with associated fractional intensity x i is up-regulated, not differentially expressed, or down-regulated may be computed from Eq. (1), respectively as:

P up,i = P(g i 2 | X=x i ) = ) ( ) ( 2 2 i i x f x f I (32) P non,i = P(g i 0 | X=x i ) = ) ( ) ( 0 0 i i x f x f I (33)
These probabilities provide quantitative measures of the evidence contained in the microarray data regarding the true expression status of each gene. However, since each probability is itself computed from data x i , the uncertainty inherent in the data is transmitted to these computed probabilities and must be quantified. As part of the inference framework we propose the following technique for quantifying these uncertainties and for converting these to a "confidence index" -a measure of the confidence associated with each computed probability.

Confidence index associated with the probability of expression status

Confidence limits around statistical estimates are generated from measures of uncertainty in the data employed for computing the estimates in the first instance; and naturally, such uncertainty measures are determined from replicates or other independent measures of pure variability. Microarray data are atypical in the sense that pure replicates, when available at all are often limited in number. Therefore, in recognition of the fact that there are very few situations where microarray data are generated from unreplicated experiments [START_REF] Reid | Determining gene expression on a single pair of microarrays[END_REF], we treat this case first before dealing with the most common and important case, where the data are available in (limited) replicates.

From unreplicated data: Under these circumstances, a measure of data uncertainty can be derived from considerations of the data generation technology itself. Under the assumption that

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

22 if the uncertainty in the intensity signal is proportional to the magnitude of the signal, the true but unknown values can be written as:

1 1 1 1 ~i i i y y y r J (35) 0 0 0 0 ~i i i y y y r J (36)
with y i1 as the actual measured intensity obtained from the gene in question under test conditions, and y i0 from the gene under control conditions. J 0 and J 1 are constants representing the extent of multiplicative uncertainty associated with the respective signal intensities. Such parameters are often supplied as part of instrument calibration characteristics.

From Eqs. ( 35) and ( 36), the true but unknown fractional intensity i

x is given by:

0 0 1 1 0 1 1 1 1 0 1 1 ~i i i i i i i i i i y y y y y y y y y x r r J J J (37)
From standard propagation-of-error analysis (Taylor, 1997), and under the reasonable assumption of equal multiplicative uncertainties in each signal channel (i.e., J 1 = J 0 = J ), the uncertainty propagation simplifies to:

i i i x x x r | J 2 (38)
Each computed nominal fractional intensity x i is thus located in the range

i i i i i x x x x x ' '
with 'x i = 2J•x i , from which corresponding probability ranges can be computed for each gene:

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

will be scaled between 0 and 1 (since the minimum value of the numerator term is 0 and the maximum value is 2). If we now define a confidence index as:

i i S c 1 (41)
then observe that 0 d c i d 1, and a low value of c i corresponds to low certainty regarding the true state of the expression status of the gene i, as indicated by the computed probabilities;

conversely, a high value of c i corresponds to a higher degree of confidence in what the computed probabilities imply about the expression status of the gene i.

From replicated data: For each gene g i , the fractional intensities computed from n replicates:

n i i i x x x , 2 , 1 , ..., , , (42)
are characterized by the usual average and the range (since typical microarray data seldom have enough replicates for determining a reliable estimate of standard deviation):

n i i i i x x x x , 2 , 1 , min , ,..., , min (43)
n i i i i x x x x , 2 , 1 , max , ,..., , max (44)
The nominal probabilities of expression status in Eqs. ( 32 thus scaled between 0 and 1, the confidence index, as before, is computed using Eq. ( 41).

We now illustrate the application of the theoretical results we have presented thus far with an example involving experimental data.

Illustrative examples

Example 1: Experimental Data

The data set shown first in Figure 1 is from studies of approximately 3000 genes from D. radiodurans 0.5 hours after ionizing radiation (following DNA damage), using cDNA microarrays.

Data representation and partitioning:

The polar coordinate transformed data is shown in Figure 2 where the intensity magnitude R is plotted against nominal fractional intensity x. The entire range of intensity magnitude is from 5.310 4 to 1.4210 7 . The data is sorted by intensity magnitude (from the lowest to the highest R value) and then partitioned into three groups containing an equal number of data points, (5.310 4 < R l ow < 6.6310 5 ); (6.6310 5 < R m e d < 1.2710 6 ) and (1.2710 6 < R high < 1.4210 7 ).

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

Nominal data characterization and model fits: A histogram of the data and the theoretical fit of mixture distributions are shown below. Figures 3 and4 show the data and model fit for "Low R" data; Figures 5 and6 for "Medium R", and Figures 7 and8 for "High R" data.

The probability model for the Low R group is: show that the probability that genes with x < 0.3 are down-regulated (dashed line) is almost 1.0, while the probability that genes for which x > 0.7 are up-regulated (thin solid line) is also close to 1.0. And when the fractional intensity for a gene is 0.50 (y i1 = y i0 ), not surprisingly the figures indicate that the probability that the gene in question is differentially expressed is extremely low, almost zero (for both lines).

x f x f x f x f
A typical output of the analysis procedure is displayed in Table 1: the fractional intensity x and the corresponding probabilities of expression status, are shown for the minimum, maximum and average values of x. The last column of the table is the confidence index (c i ).

There are several ways in which the researcher may use the computed probabilities and the confidence index. For example, Figure 12 shows a histogram of the confidence index values associated with the computed probabilities of expression status. This provides a means of quantifying the overall "quality" of the experimental data. For instance note that the figure

indicates that a substantial number (76.3%) of the analyzed genes have a confidence index greater than 0.80 associated with their computed probabilities of expression status -implying relatively low variability associated with the data in general.

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

29 However, the most important use of the computed results is in identifying and selecting differentially expressed (down-regulated or up-regulated) genes. By using as the cut-off criteria, threshold values for P down or P up appropriate to the study at hand, it is possible to select candidate genes satisfying the stipulated criteria. For example, to generate a list of down-regulated candidate genes, choosing a threshold value of 0.75 for P down (indicating a desire to select only those genes that have a probability greater than or equal to 0.75 of being down-regulated, regardless of the precision associated with the computed probability) results in a list containing 218 out of 3,094 genes.

By imposing a threshold value on the precision of the computed probabilities, we can further narrow down the search for potentially down-regulated genes. Thus, for example, imposing a cut-off value of c i t 0.90 reduces the number of candidate down-regulated genes from 218 to 129 genes. The researcher may now use a more precise, higher resolution experimental technique to determine the true expression status of this smaller subset of selected genes.

Example 2: A Simulation study

The purpose of this example is to complement Example 1 by assessing the performance of our approach in detecting differentially expressed genes when the "true" circumstances are known.

The study involves five microarray datasets simulated with statistical characteristics similar to the ones observed in real microarray data. Each simulated microarray dataset consisted of 3,000 genes with 750 genes differentially expressed: 300 genes were down-regulated (10% of the 2.

Based on the probabilities of expression status, the top 300 down-regulated and the top 450 upregulated genes were chosen as differentially expressed. Table 3 summarizes the total number of genes correctly identified as a function of the number of simulated microarrays (replicates).

Next, the simulated microarray data was also analyzed using the nonparametric method Significance Analysis of Microarrays (SAM) [START_REF] Efron | Empirical Bayes analysis of a microarray experiment[END_REF], a popular technique that does not require any assumption about the distribution of the microarray data. SAM does not produce probabilities; rather, it generates a sorted list of candidate genes given a specified False Discovery Rate (FDR). (The FDR associated with a decision rule is the proportion of false positives among all the genes initially identified as being differentially expressed (Wit and McClure, 2004)). The total number of genes correctly identified by SAM (with FDR = 3%) is also shown in Table 3 as a function of the number of available replicates.

We observe first, that the performance of the probabilistic approach is uniformly better than that obtained using SAM. This should come as no surprise: statistical inference based on distributionfree techniques general tend to be less powerful when the underlying data come from known distributions. Notwithstanding, the performance obtained with SAM is still quite good although it is unable to handle the case where no replicates are available. Second, and, again, not In particular, it is interesting to observe that with five replicates, the probabilistic approach identified 744 of the original 750 differentially expressed genes, a virtually-perfect precision of 99.6%; with SAM, 710 of the original 750 differentially expressed genes were identified, a stillimpressive precision of almost 95%. For example, for the case shown in Figures 3 and4, with the model fit shown in Eq. ( 31), the implications are the following: the contributions from the class of potentially down-regulated genes, is about 18%; the contribution from the class of genes whose expression levels are potentially unchanged is 79%; the remaining 3% is the complement, the contribution from potentially up-regulated genes.

In conclusion, we observe that just as classical statistical inference is based on theoretical reference distributions (such as the Gaussian, t-, Chi-square, or the F-distributions) we have developed a methodology for drawing statistical inference using the mixture distribution f(x) as its theoretical basis: it consists of (i) a probability statement that a gene belongs to a category (the down-regulated, the up-regulated, or the no-change) and (ii) a degree of confidence associated with such probability statements, determined from the variability estimated from replicates, or else by propagation-of-error techniques when there are no replicates. The final outcome is an ordered triplet of results for each gene: a raw computed fractional (or relative) change in expression level, an associated probability that this number indicates lower, higher, or no differential expression (a category-membership probability) and a measure of confidence associated with the stated result. This analysis technique, which has been illustrated with a real experimental data set and also via simulation, provides an alternative approach specifically tailored to the inherent characteristics of microarray data.

The key result at this point is the following identity for the RHS of Eq. ( 50) above:

³ ¦ f I k z z k z I dz z e k z I e K K K 1 1 0 ! 1 1 ! (51)
To establish this result, first let

³ f a k z dz z e k 1 W (52)
from where it is easy to show via integration by parts, that

1 1 1 k k e a k a k W W (53)
a recursion equation that leads immediately, upon successive substitution and rearrangement, to establishing the result in Eq. ( 51).

We may now use this result in Eq. ( 48) to obtain:

> @ > @ ³ f d ! I k z dz z e k I k y P I k y P K 1 ! 1 1 1 (55)
or, from the definition of the cumulative distribution F y (I),

  Different technologies employ different techniques for such background correction. We note here that for high and medium intensity measurements, the contributions from b j are relatively insignificant; appropriately accounting for the presence of this component is critical however at very low intensity measurements.

2.

  Convert count variable z(I) to continuous variable y. Let y(k) be the total observed intensity for k molecules, i.e. in Appendix A, the probability distribution function for y given the probability distribution function for z(I) is: this result are that the total fluorescent intensity measurement for each spot on the microarray possesses a Gamma distribution with the shape parameter D related directly to the number of cDNA molecules responsible for the observed intensity, y; the scale parameter E is the mean intensity per molecule.

  preserving the two dimensional character of the original data; (iii) the characteristics of x i are less than 0.5 indicate down regulation; x i = 0.5 indicates no change and x i > 0.5 indicates up-regulation.

  an estimate of the expected value of the nominal x for those genes showing no change in gene expression, we are able (i) to determine if nominal conditions hold and then (ii) correct the nominal fractional intensity for the entire data set using the estimate of w given above when nominal conditions do not hold.

  by fitting a mixture of Beta distributions to histogram data, we obtain a theoretical probabilistic characterization of the data where each contributing Beta distribution constitutes the population description of a category of genes with common differential expression attributes. Observe that the coefficients I 0 , I 1 and I 2 indicate the fraction of the entire population that belongs in each respective indicated category. In particular, for the fraction I 0 of genes for which there is truly no differential expression, x possesses a Beta distribution that is perfectly symmetric, with mean = mode = median = x* = 0.5. In some cases, it may be necessary to represent the data histogram with more than three Beta densities. Under such circumstances, the appropriate number of Beta densities required to represent the histogram structure adequately may be determined by employing such standard model selection criterion 's Information Criterion (AIC) or the Bayesian Information Criterion (BIC) in selecting among various possible models involving n = 4, 5, 6,… individual densities in the mixture model.

  i = P(g i 1 | X=x i ) = 1 -P up,i -P non,i

  )-(34) are now computed for the average fractional intensities, and the associated uncertainty determined from the probabilities computed at the extremes, m in , i x and m ax , i x. The magnitude of the uncertainty vector in this case is defined as: factor d depends on the number of replicates. It is easy to show that specifically, d = 2 if n = 2 (duplicates); otherwise d = 3 if n t 3 (triplicates or more). With S i

  Figures 3, 5 and 7 visually different; quantitatively, the parameter values of the contributing Beta

  450 genes were up-regulated (15% of the genome) and 2,250 genes were not differentially expressed (75% genome). The parameters of the Beta distributions were chosen to mimic very closely actual distributions of real microarray data (after the polar transformation) as observed in our laboratory. The main statistics of the three Beta distributions used to simulate the microarray data are summarized in Table

  techniques show improved performance as the number of replicates increase.

  resulting coefficients of the component Beta densities therefore indicate what fraction of the entire data set belong to the category represented by the Beta density in question.

  to I using Leibnitz's rule for differentiating under the integral,

  desiccation reveals novel proteins that contribute to extreme radioresistance. Genetics 168, 21-33.

FiguresFigure 1 :

 1 Figures

Figure 2 :

 2 Figure 2: Raw data of Figure 1 after the polar coordinates transformation.

Figure 3 :

 3 Figure 3: Histogram of Low R data and mixture model.

Figure 4 :

 4 Figure 4: Contributing Beta densities for mixture model fit to data in Figure 3.

Figure 5 :

 5 Figure 5: Histogram of Med R data and mixture model.

Figure 6 :

 6 Figure 6: Contributing Beta densities for mixture model fit to data in Figure 5.

Figure 7 :

 7 Figure 7: Histogram of High R data and mixture model.

Figure 8 :

 8 Figure 8: Contributing Beta densities for mixture model fit to data in Figure 7.

Figure 9 :

 9 Figure 9: Histogram of bias-corrected High R data and mixture model.

Figure 10 :

 10 Figure 10: Contributing Beta densities for mixture model fit to corrected data in

Figure 11 :

 11 Figure 11: Fitted mixture model (thick solid line) for the different intensity signal

Figure 12 :

 12 Figure 12: Histogram of the confidence index associated with the probabilities of

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

	it is a near-
	impossible task to determine the absolute value of the expression levels O i1 , O i0 strictly from the
	observed microarray intensity measurements. Even so, determining by how much O i1 differs
	from O i0 without necessarily knowing the absolute value of each number is still considered
	extremely valuable information about the gene in question. While methods for determining
	absolute gene expression levels are being perfected for lower cost and higher-throughput (for

  Taylor, J.R., 1997. An introduction to error analysis: the study of uncertainties in physical measurements. University Science Books, Sausalito, Calif.

	Thattai, M., and van Oudenaarden, A., 2001. Intrinsic noise in gene regulatory networks.
	Proc Natl Acad Sci USA 98, 8614-8619.
	Troyanskaya, O.G., Garber, M.E., Brown, P.O., Botstein, D., and Altman, R.B., 2002.
	Nonparametric methods for identifying differentially expressed genes in microarray
	data. Bioinformatics 18, 1454-1461.

Tusher, V.G.,

Tibshirani, R., and Chu, G., 2001

. Significance analysis of microarrays applied to the ionizing radiation response. Proceedings of The National Academy of Sciences of The United States of America 98, 5116-5121. Wit, E., and McClure, J.D., 2004. Statistics for microarrays: design, analysis, and inference. John Wiley & Sons, Chichester, England; Hoboken, NJ, USA.

Table 1 .

 1 Example of a typical output of the probabilistic framework.

	Fractional intensity	(x)
		Gene

P down P non P up Confidence index (c i ) Min Aver. Max Min Aver. Max Min Aver. Max Min Aver. Max

		0.97	1.00	1.00	0.91	0.96	1.00	0.96	1.00	0.61	0.59	0.95	0.88	1.00	1.00	0.95	1.00	0.97		0.91	0.79	0.51	1.00	1.00	0.95	0.55	0.76	0.82
		0.52 0.01 0.02 0.04 0.96 0.97 0.98 0.00 0.01 0.01	0.90 0.00 0.00 0.00 0.00 0.00 0.00 1.00 1.00 1.00	0.89 0.00 0.00 0.00 0.00 0.00 0.00 1.00 1.00 1.00	0.35 0.73 0.77 0.84 0.16 0.23 0.27 0.00 0.00 0.00	0.69 0.00 0.00 0.00 0.07 0.09 0.12 0.88 0.91 0.93	0.95 0.00 0.00 0.00 0.00 0.00 0.00 1.00 1.00 1.00	0.75 0.00 0.00 0.00 0.00 0.02 0.05 0.95 0.98 1.00	0.70 0.00 0.00 0.00 0.00 0.00 0.00 1.00 1.00 1.00	0.44 0.48 0.80 0.96 0.04 0.20 0.52 0.00 0.00 0.00	0.41 0.49 0.81 0.99 0.01 0.19 0.51 0.00 0.00 0.00	0.72 0.00 0.00 0.00 0.01 0.04 0.07 0.93 0.96 0.99	0.51 0.04 0.09 0.19 0.81 0.91 0.96 0.00 0.00 0.00	0.93 0.00 0.00 0.00 0.00 0.00 0.00 1.00 1.00 1.00	0.69 0.00 0.00 0.00 0.00 0.00 0.00 1.00 1.00 1.00	0.70 0.00 0.00 0.00 0.00 0.02 0.06 0.94 0.98 1.00	0.68 0.00 0.00 0.00 0.00 0.00 0.00 1.00 1.00 1.00	0.51 0.03 0.04 0.07 0.93 0.96 0.97 0.00 0.00 0.00		0.71 0.00 0.00 0.00 0.04 0.09 0.14 0.86 0.91 0.96	0.69 0.00 0.00 0.00 0.09 0.25 0.35 0.65 0.75 0.91	0.46 0.37 0.77 0.97 0.03 0.23 0.63 0.00 0.00 0.00	0.78 0.00 0.00 0.00 0.00 0.00 0.00 1.00 1.00 1.00	0.74 0.00 0.00 0.00 0.00 0.00 0.00 1.00 1.00 1.00	0.39 0.94 0.97 0.99 0.01 0.03 0.06 0.00 0.00 0.00	0.39 0.44 0.81 1.00 0.00 0.19 0.56 0.00 0.00 0.00	0.73 0.00 0.00 0.00 0.01 0.13 0.30 0.70 0.87 0.99	0.62 0.00 0.00 0.00 0.04 0.11 0.25 0.74 0.89 0.96
		0.51	0.90	0.88	0.35	0.69	0.94	0.72	0.69	0.41	0.40	0.71	0.49	0.92	0.68	0.66	0.66	0.49		0.69	0.67	0.42	0.78	0.73	0.38	0.38	0.70	0.61
	5. Table 1	arsenate reductase 0.49	cinA protein 0.89	citrate lyase, beta subunit 0.86	cytochrome P450-related protein 0.33	dihydroxy-acid dehydratase 0.68	DNA gyrase, subunit A 0.94	DNA helicase II 0.70	frnE protein 0.68	glutaryl-CoA dehydrogenase 0.39	glycerol uptake facilitator protein 0.38	glycerol-3-phosphate dehydrogenase 0.69	GMC oxidoreductase 0.47	Holliday junction DNA helicase 0.92	mocR protein 0.67	nitric oxide synthase-related protein 0.64	nitrogen regulatory protein P-II 0.64	nodulation protein-related protein 0.47	peptide methionine sulfoxide	reductase 0.68	PhoH-related protein 0.66	phosphoribosylamine--glycine ligase 0.40	tellurium resistance protein TerA 0.77	thiophene and furan oxidation protein 0.72	0.36 threonine synthase	0.38 tryptophan synthase, beta subunit	0.66 uroporphyrinogen decarboxylase	0.60 v-type ATP synthase, D subunit

Table 2 .

 2 Main statistics of the distributions used to simulate the microarray data.

	5. Table 2
	Status	Beta parameters (α, E E)	Mean fractional intensity (x)	Standard deviation	Mean fold-change	Log 2 (fold-change)
	Down-regulated	(38, 62)	0.38	0.048	0.62	-0.69
	Not
	differentially	(75, 75)	0.50	0.041	1.00	0.0
	expressed
	Up-regulated	(62, 38)	0.62	0.048	1.63	0.70

Table 3 .

 3 Total number of correctly identified genes by the probabilistic framework and SAM (the percentages of correctly identified genes are shown in parenthesis).

	5. Table 3
		1 array	2 arrays 3 arrays 4 arrays 5 arrays
	Probabilistic	603	688	729	743	747
	framework	(80.4%)	(91.5%)	(97.2%)	(98.9%)	(99.6%)
		Requires	668	699	701	710
	SAM
		replicates	(88.7%)	(93.2%)	(93.3%)	(94.7%)

Table 4 .

 4 Table of complete model results.

		Low R	Medium R	High R
	f 1 (D 11 , D 10 )	(64.0, 98.8)	(76.0, 126.1)	(15.2, 30.7)
	corrected	(69.3, 88.1)	(74.9, 100.3)	(16.9, 28.6)
	f 0 (D 01 , D 00 )	(98.8, 120.2)	(121.5, 148.8)	(31.1, 41.3)
	corrected	(106.9, 106.7)	(136.8, 136.1)	(30.1, 30.2)
	f 2 (D 21 , D 20 )	(137.8, 104.4)	(95.6, 72.5)	(29.4, 20.3)
	corrected	(150.3, 92.8)	(66.9, 42.4)	(50.3, 22.3)
	f 3 (D 31 , D 30 )			(61.6, 9.5)
	corrected			(96.3, 11.1)
		0.18	0.22	0.22
	I 1 corrected
		0.18	0.24	0.16
		0.79	0.73	0.68
	I 0 corrected
		0.79	0.71	0.77
		0.03	0.05	0.07
	I 2 corrected
		0.03	0.05	0.05
				0.03
	I 3 corrected
				0.02

5. Table

4

/

Acknowledgements

We would like to thank Dr. John R. Battista of the Louisiana State University for providing us the data on D. radiodurans. BAO and CAG acknowledge funding from the Delaware Biotechnology Institute; JSE acknowledges funding from the US Department of Energy Grant # DE-FG02-O1ER63200 and a subcontract grant from Genomatica.

. The superscripts '-' and '+' respectively denote the probabilities associated with the fractional intensity at the low end of the range, x x ' , and at the high end x x ' . The uncertainties in the fractional intensity data have thus been translated into uncertainties in the probabilities with ranges defined by , respectively.

Thus, associated with the nominal probabilities computed from Eqs. ( 32)-( 34) are the uncertainties represented by the ranges of the probabilities of expression status, for any given J.

Clearly, the wider these ranges are, the greater the degree of uncertainty associated with each computed probability, reducing our confidence in these probabilities. Conversely, narrow ranges imply less uncertainty and hence increased confidence.

Because the probabilities of certain event must add up to 1, it is easy to establish that the elements of ΔP, the vector of uncertainties associated with the probabilities, must satisfy the following fundamental constraint:

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

Summary and Conclusions

The three main results of this paper are now summarized:

1. The distribution of spot intensity follows a Gamma model as indicated in Eq. ( 7): the shape parameter, D, is related directly to the number of cDNA molecules responsible for the observed intensity, and the scale parameter, E, is the mean intensity per molecule an equipment calibration parameter.

2. For modeling differential expression, a polar transformation of the raw spot intensity data (as in Eqs. ( 15) and ( 16)) yields the intensity magnitude R i , (to be used for rationally segregating the data) and the fractional intensity x i , which possesses a Beta distribution as in Eq. ( 24); the distribution parameters, D i1 and D i0 are related directly to the number of cDNA molecules under test and reference conditions, respectively.

3. The distribution of entire microarray fractional intensities is characterized, theoretically, by a mixture of overlapping Beta densities; in practice (and as demonstrated with an example) as a result of the coalescence of distributions of genes with similar expression levels, the discernible features of the histogram of real data may be captured by as few as three overlapping Beta densities, as shown in Eq. ( 31).

This last result has the most direct impact on the quest for reliable determination of differential expression for thousands of genes from noisy, potentially biased, microarray data, with few or no replicates. It indicates that, f(x), the underlying probability density function of the fractional intensity, is a linear combination of Beta density basis functions providing the following important perspective of the data set: each contributing Beta density may be viewed as a

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

Appendix A: Derivation of the pdf in Eq. ( 7)

To derive a probability distribution function for y as defined in Eq. ( 6), given the probability distribution function for z(I), we begin by observing that z(I) and y(k) are related according to the following expression:

which, from the definition of the variables z and y, translates in words to: "if the z-count is not yet up to k then the observed intensity I must be less than y". In the language of probability, the event that z is less than k is equivalent to the event that the observed intensity I corresponding to z is less than y, and the probability of one event is identical to that of the other. From Eq. ( 48) we obtain

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

For each data group, segregated according to the intensity magnitude R, the parameters (D j1 , D j0 ) for the contributing probability densities, f j (x), and the associated mixture coefficient weights, I j (j = 0, 1, 2, 3) are shown in Table 4 for both nominal and corrected conditions. For each group, note the inequality of the nominal parameters (D 01 , D 00 ) associated with f 0 (x); the indicated bias is eliminated upon correction.

Note also that bias correction has little or no effect on the mixture coefficient weights for Low and Medium R data, the influence of bias correction on these weights is more noticeable for the High R group.


Dimitri Perrin
email: dperrin@computing.dcu.ie

Heather J Ruskin

Tohru Niwa

induced, aberrant DNA methylation in gastric cancer

Keywords: Cancer, Computational model, DNA methylation, Epigenetics

Epigenetic changes correspond to heritable modifications of the chromatin structure, which do not involve any alteration of the DNA sequence but nonetheless affect gene expression. These mechanisms play an important role in cell differentiation, but aberrant occurrences are also associated with a number of diseases, including cancer and neural development disorders. In particular, aberrant DNA methylation induced by H. Pylori has been found to be a significant risk factor in gastric cancer. To investigate the sensitivity of different genes and cell types to this infection, a computational model of methylation in gastric crypts is developed. In this article, we review existing results from physical experiments and outline their limitations, before presenting the computational model and investigating the influence of its parameters.

Introduction

Early advances in Genetics led to the all-genetic paradigm: the phenotype is a deterministic consequence of the genotype. Obvious counter-examples were outlined and this was later amended: visible characteristics of a living organism (phenotype, P ) combine hereditary genetic (genotype, G) and environmental factors (E). This is sometimes summarised as P = G + E.

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

However, differences between monozygotic twins can be seen even for high heritability diseases such as schizophrenia, while identification of environmental factors, e.g. smoking and air quality for lung cancer, [START_REF] Spurny | Chemical mixtures in atmospheric aerosols and their correlation to lung diseases and lung cancer occurence in the general population[END_REF][START_REF] Mucha | Meta-analysis of disease risk associated with smoking, by gender and intensity of smoking[END_REF], is relatively rare.

Early work by [START_REF] Waddington | The genetic control of development[END_REF], and more recently over the last decade, e.g. [START_REF] Bird | DNA methylation patterns and epigenetic memory[END_REF]; Wilkins (2005), has emphasised that genotype expression can be altered without changing DNA sequence itself, a phenomenon called Epigenetics:

P = G + E + EpiG.
Epigenetic mechanisms involve heritable alterations in chromatin structure, (e.g. DNA methylation and histone acetylation), amongst other epigenetic "signatures". In turn these regulate gene expression, but do not involve changes in DNA sequence. These "stable alterations" arise during proliferation and persist through cell division [START_REF] Feinberg | The history of cancer epigenetics[END_REF]. While information within the genetic material is not changed, instructions for its assembly and interpretation may be.

Such alterations have, in particular, been linked with tumor-suppressor inactivation and cancer initiation [START_REF] Esteller | Cancer epigenomics: DNA methylomes and histonemodification maps[END_REF]. In this study, we focus on aberrant DNA methylation induced by H. Pylori infection. This modification has recently been associated with higher cancer risk. We summarise the results obtained from physical experiments, and propose a computational model to complement these.

Infection-induced epigenetic perturbations and cancer initiation

DNA methylation

DNA methylation corresponds to the addition of a methyl group to a cytosine. In humans, only 1% of DNA bases undergo DNA methylation. In differentiated cells, DNA methylation is typically limited to CpG dinucleotides, while non-CpG methylation can be found in embryonic stem cells, [START_REF] Dodge | De novo methylation of MMLV provirus in embryonic stem cells: CpG versus non-CpG methylation[END_REF].

Of particular interest are CpG islands. These correspond to areas with higher proportion of CpG, and are formally defined as follows [START_REF] Gardiner-Garden | CpG islands in vertebrate genomes[END_REF]):

• Length of the considered region is at least 200 base pairs.

• GC percentage is greater than 50%, (i.e. more than half of amino-acids are cytosine or guanine).

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

• Observed/expected CpG ratio is greater than 60%.

In humans, these islands are found in or near to 70% of gene promoters, [START_REF] Saxonov | A genome-wide analysis of CpG dinucleotides in the human genome distinguishes two distinct classes of promoters[END_REF]. While most CpG dinucleotides are methylated over the genome, methylation of a CpG island is specifically associated to silencing of the associated gene. Aberrant changes in CpG island methylation are, therefore, linked with abnormal gene expression.

Infection-induced methylation in gastric cancers

Epigenetic alterations in non-disease tissues can be used as markers for disease risk and past exposure to some disease-inducing factors. In particular, current focus is on detection of aberrant DNA methylation in non-cancerous gastric mucosae, as the presence of such patterns can be used as a marker for both the risk of gastric cancers and for past exposure to Helicobacter pylori. A detailed presentation in the context of this study is proposed by [START_REF] Nakajima | DNA methylation: a marker for carcinogen exposure and cancer risk[END_REF]. For convenience, main points are summarised here.

In cancer development, aberrant DNA methylation is involved at two levels:

• Overall hypomethylation, which affects repetitive DNA sequences and causes both chromosomal instability, and, hence, tumours, [START_REF] Gaudet | Induction of tumors in mice by genomic hypomethylation[END_REF], as well as aberrant expression of normally methylated genes, [START_REF] Smet | DNA methylation is the primary silencing mechanism for a set of germ line-and tumor-specific genes with a CpG-rich promoter[END_REF]).

• Regional hypermethylation, most of which affects CpG islands and causes, (if these islands are located in gene promotion region), transcriptional silencing of the downstream gene. In the context of cancer, methylation affecting tumour suppressor genes is well documented, see e.g. [START_REF] Jones | DNA methylation and cancer[END_REF][START_REF] Baylin | Epigenetic gene silencing in cancer -a mechanism for early oncogenic pathway addiction?[END_REF]. This is sometimes referred to as driver methylation, because of causal involvement in carcinogenesis, (as opposed to passenger methylation, which refers to genes whose methylation is a consequence of cancer development, and therefore requires careful analysis of newly detected genes [START_REF] Ushijima | Detection and interpretation of altered methylation patterns in cancer cells[END_REF]).

In the context of gastric cancers, gene inactivation, (e.g. for tumour suppressor gene p16), is more frequently a consequence of aberrant promoter methylation than of defaults at the genetic level, [START_REF] Ushijima | Focus on gastric cancer[END_REF], but [START_REF] Kaneda | Identication of silencing of nine genes in human gastric cancers[END_REF] also observed that low levels of aberrant methylation occur in non-cancerous mucosae of cancer patients. This was

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

tested against tissues from healthy individuals, and it was found that methylation levels for these patients were one to two orders of magnitude higher in H. pylori -positive individuals than H. pylori-negative individuals [START_REF] Maekita | High levels of aberrant DNA methylation in Helicobacter pylori-infected gastric mucosae and its possible association with gastric cancer risk[END_REF]. This is a very significant finding, since H. pylori is known to be a major risk factor for gastric cancers.

It was also highlighted that part of this hypermethylation is temporary and will decrease after eradication of infection. This is not due to active demethylation, but simply to cell turnover. The structure of the gastric crypt, (one stem cell, multiple progenitor cells and many differentiated cells, as shown in Figure 1), is the cause for this two-part methylation, where:

• The permanent component is due to methylation of stem cells. Since methylation is conserved during cell division, progenitor and differentiated cells, obtained from an aberrantly methylated stem cell, will exhibit identical abnormal patterns.

• The temporary component is due to methylation in progenitor and differentiated cells. This, if the stem cell of the crypt is not methylated, will disappear because of cell turnover and creation of new, unmethylated, cells.

Experimental limitations and computational alternative

3.1. Issues: experimental cost and granularity of the results Investigation of infection-induced methylation relies on a time-consuming experimental protocol. Based on fixed infection dynamics, several lines of rodents are used: one line is infected for an extended period, (for instance from week 5 to week 55), a second is used as a control, (no infection), while a third is subjected to the eradication procedure only, (at week 55 in the current instance), in case this itself induces aberrant methylation. As the experiment progresses, samples are regularly taken from each population, and the average methylation level is obtained for the genes under investigation. Each sample contains an average of thirty cells, and requires termination of the sampled rodent.

This protocol provides a suitable proxy to investigate aberrant methylation in humans, as the behaviour for cell types of interest is similar, (see e.g. [START_REF] Nakajima | Persistence of a component of DNA methylation in gastric mucosae after helicobacter pylori eradication[END_REF]; [START_REF] Niwa | Causal role of helicobacter pylori infection in induction of aberrant DNA methylation in gastric epithelial cells, and critical involvement of inflammation[END_REF] for details). Yet, while this experimental protocol provides essential results, (see [START_REF] Nakajima | DNA methylation: a marker for carcinogen exposure and cancer risk[END_REF]

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

for a more comprehensive presentation), there are a number of limitations. These include reproducibility: each experiment lasts over a year, and requires large rodent populations. This is an obstacle to more comprehensive studies, for instance involving multiple infection dynamics. Investigating these would be particularly interesting, as the clearance rate then becomes crucial, but this is not immediately feasible, at least not in vivo.

The experimental contraints also limit the number of sampling timepoints. To satisfy statistical assumptions, each time-point requires multiple samples, taken from each sub-population, which implies a rapid increase, over the extended period, of the overall population required for the experiment. Given this, the sampling intervals is subject to practical limitation: typically every five to ten weeks under "stable" conditions, and every week at most during transitions, (initiation and eradication of the infection).

This, in addition to being able to obtain only the average methylation level, results in low granularity of the results. The "real-time" dynamics of epigenetic changes at the cell level, or even crypt level, are not accessible.

Gaining a better insight into these methylation dynamics, both during infection, (summarised in Figure 2a), and in the long term, (Figure 2b), is the objective of the prototype model we have developed. This is detailed in the following.

Computational model 3.2.1. Structure and typical iteration

To investigate the dynamics of infection-induced aberrant methylation in the crypt, we implement an object-oriented model of the entities involved, (samples, crypts, cells). Each entity corresponds to a computing object, and may inherit other objects, (e.g. a crypt contains cells). Each object has a number of attributes which characterise it further, (e.g. each cell has a type: stem cell, progenitor cell, or differentiated cell).

A sample, which mimics in silico the samples obtained in vivo during the physical experiments described above, is implemented as an array of crypts. Each crypt is initialised with one stem cell, six progenitor cells, and one hundred differentiated cells on each "wall" of the crypt. Initially, no cell is aberrantly methylated.

At each time step, (fixed at one minute for future inclusion of histone modifications), each crypt is updated, as follows. First, cells at the top of both "walls" are analysed and, as these have a finite life span, are removed if necessary. This is followed by update of the bottom of the crypt. A new cell may be created on each side, and progenitor cells may be replaced, (as these are limited in terms of number of differentiated cells produced prior

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

to replacement by a new progenitor cell). Newly-created cells inherit the epigenetic status of the cell from which they derive through the differentiation process.

Finally, the methylation status of the progenitor and stem cells is updated, depending on cell type and infection status.

The parameters used during this iterative process are described next.

Parameters

Key attributes in control of the cell renewal dynamics in the crypt are the life span of each cell type, and the probability that a new cell is created at the bottom of the crypt. In order to realistically reproduce observed dynamics, the life span, under normal conditions, is 3 days for differentiated cells (i.e. 4320 time steps). When the infection is established, this is reduced to 2.5 days, (3600 time steps). For progenitor cells, the life span is taken to be equivalent to the total number of differentiated cells they can produce, at unit intervals; set at 25 each. Stem cells in the model are considered to have an infinite life span.

To balance the limited life span of differentiated cells, new cells are produced at the bottom of the crypt. At each time step, this occurs with a probability p of 2.3 × 10 -3 under normal conditions. This maintains the crypt size to approximately 200 cells. During an infection, this size is increased by a factor ranging from 2 to 4, and this is reflected in the model by increasing p to a value between 7.0 × 10 -3 and 1.2 × 10 -2 , (randomly chosen for each crypt).

In order to investigate aberrant methylation, as described above, probabilities that this change occurs are required, both under normal conditions and during the infection, for both stem cells and progenitor cells. The influence of these parameters is analysed in Section 4.

Transition function

Most model parameters vary depending on the infection status. To account for a progressive transition between the two parameter values, (from normal conditions to established infection), a sigmoid function is used. Equation 1 represents a generic transition between values a and b. In the context of cell production, for instance, a would correspond to p = 2.3 × 10 -3 , and b to p = 7.0 × 10 -3 . A typical sigmoid function is shown in Figure 3.

Using this sigmoid function introduces two extra parameters, λ and T , (with λ < 0 and T ≥ 0). After eradication of the infection, initial values are restored, (Equation 2). In both equations, units for t are weeks.

f (t) = a + (b -a) × 1 1 + e λ(t-T )
(1)

f (t) = a + (b -a) × (1 - 1 1 + e λ(t-T ) )
( 2 )

Results and discussion

Reproduction of two-part methylation

The hypothesis, based on physical experimentation, is that, under normal conditions, aberrant stem cell methylation does not occur and that, during H. Pylori infection, the probability for methylation of these remains significantly lower than that for progenitor cells.

To validate this, a series of simulations was run, with a range of values for these probabilities, and infection dynamics were set to reproduce experimental conditions: infection, if included, starts at week 5, and is eradicated at week 55. Three typical runs are shown in Figure 4, and correspond to the three populations used during the physical experiments. These results clearly indicate that significant methylation levels can only be obtained through infection and that, after eradication, only part of this aberrant methylation remains. This permanent level is equal to the proportion of methylated stem cells in the sample.

The results from the physical experiments have been successfully reproduced, and the model hypothesis validated. The remainder of this section analyses the sensitivity of this model to its various parameters, and the implications of this sensitivity with respect to the observed phenomena. Because aberrant methylation is modelled as a stochastic process, this analysis is based on multiple simulations for each model configuration. In order to obtain statistically significant values, averages are calculated over twenty sim-

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

ulations. A similar protocol is used during the in vivo experiments, in that several rodents are used to obtain values for each time and each configuration.

Impact of cell-type dependency

The model is based on specific cell types having different sensitivity to H. Pylori infection. In particular, the main conclusion of the physical experiments is that infection-induced aberrant methylation is significantly more likely to occur on progenitor cells than on stem cells, and has a relatively short lifetime.

To account for this, four parameters are necessary. α and β represent the probability of aberrant methylation once the infection is established, for stem cells and progenitor cells respectively. Under normal conditions, the respective probabilities are denoted α 0 and β 0 . A first assumption, as mentioned, is that α β. Based on the experimental information, we can also assume that α 0 = 0.

Sensitivity to α is shown in Figure 5. This probability of stem cell methylation at each time step has two consequence for model behaviour: higher values of α increase both (i) the methylation level, at any given point during the infection, and (ii) the level of permanent methylation after eradication. In physical experiments, the level of permanent methylation for various genes was observed to be between 0.04 and 1.02%. This corresponds, in our model, to α between 0 and 0.9 × 10 -8 .

Sensitivity to β is shown in Figure 6. Increasing the probability of methylation for progenitor cells only affects the methylation level at a given timepoint during the infection. This is increased, but both the rate of increase during infection and the methylation level after eradication are unchanged.

These results are significant, as they imply that accurate values for α and β for any gene can be extracted from the model. Using data from physical experiments, we initially obtain α, by focusing on the level of permanent methylation after eradication of the infection. Then, once α is fixed, the methylation level at several time-points during the infection can be used to deduce β. Finally, variations about the average values both before and after the infection are sufficient to estimate β 0 .

Impact of infection dynamics 4.3.1. Sensitivity to transition periods

Neither the initial size of the crypt nor the size it reaches once the infection is established have any impact on the observed methylation level. This is A c c e p t a consequence of the methylation level being a proportion of methylation cells rather the number of such cells, (and has been verified experimentally). However, the transition between these two values does have a significant impact, which we investigate in detail.

:J= K L M N O P Q M R S T P L U L P V W X Y Z 8J< [ Y Z 8J> [ Y Z 8J@ [ Y Z 9J8 [ Y Z 9J: [ Y Z 9J< [ \]^_ \]^_ \]^_ \]^_ \]^_ \]^_
Two parameters are used in the sigmoid functions which describe this transition: λ and T . Two sets of values are used during each simulation: λ s and T s , (for changes with respect to cell renewal dynamics in the crypt); λ g and T g , (for modifications of the methylation probabilities). Figures 4 to 6 were shown for λ c = λ g = -4, and T c = T g = 2. Here, we analyse the impact of these values on the model behaviour.

Figure 7 summarises the effect of reducing λ g compared to λ s , (which must coincide with an increase of T g , so that the transition function remains realistic). This corresponds to a situation where the impact of the infection on aberrant methylation of a specific gene is established more slowly than the effect on crypt size, (through perturbed cell renewal dynamics). In such a case, the average methylation level over the sample is slower to rise, but the most significant difference occurs when the infection is eradicated. Cells are still likely to be methylated, since this transition is slower, but the crypt size is already decreasing. This results in a brief but sharp increase of the methylation level, as most destroyed cells are unmethylated, (due to low methylation level overall). Such an increase was observed for several genes during the physical experiments, and a difference between λ s and λ g is one plausible explanation of this phenomenon. As the observed increase is generally slow, this would correspond to a small difference between the two values.

Figure 8 summarises the effect of reducing λ s , (which, similarly to before, must coincide with an increase in T s , so as to maintain realism of the transition function). This occurs when aberrant methylation starts before the crypt size has reached its maximum infected level. The effect is, therefore, opposite to that of reducing λ s , and a temporary, sharp increase of the overall methylation level is observed when the H. Pylori infection begins.

Due to experimental limitations, it was possible to test in vivo only one infection configuration, (starting at week 5, and eradicated at week 55). However, shorter infections are not only possible, but can also occur repeatedly. These conditions give rise to another series of simulations. We first considered two infections, from weeks 5 to 25 and 35 to 55. Figure 9 summarises the results, for λ g < λ s , λ g = λ s , and λ g > λ s . We then simulated a succession of shorter infections, from weeks 5 to 13, 21 to 29, 37 to 45, and 53 to 61, (summarised in Figure 10 for similar configurations of the transition period). As expected from the earlier results, configurations where there is a significant difference between λ g and λ g create distortions of the methylation level during transition periods. However, if each occurrence of the infection is long enough, and if there is a sufficient time between these, the parameters related to the transition period have no impact on the longterm methylation level. This is of the same order as that for a single infection for which the length would be equal to the combined length of these shorter occurrences.

A significant observation, for short infections, is that the distortions created by differences between λ g and λ g may persist for most of the infection period. During physical experiments, (which have to yet to be performed for infection dynamics such as these), this would lead to a significant overestimation of the size of actual changes in methylation.

For clarity, it must be noted that, in Figures 9 and 10, the difference between the long-term methylation level of the three configurations is not a consequence of the values for λ g and λ s . In Figure 9, the lower methylation level for the first configuration is a consequence of the methylation of fewer stem cells during the first infection phase for some simulations. In Figure 10, a similar slight underestimation exists, and is a result of fewer methylations during the third infection phase. These differences occurred after the infection was fully established, and are an artefact of the finite number of simulations: the shorter the infection length, the greater the standard variation in the number of methylated stem cells between simulations, (because of very low probability, at each time step, to methylate a stem cell).

A similar increase in the variation can be expected for in vivo experiments. Given the limited number of samples at each time-point, this would further reduce the reliability of the measurements under these infection dynamics.

Increasing the number of simulations permits the reduction of this effect in the computational model, but is not a realistic for the in vivo counterpart, because of cost and space constraints associated with increasing the rodent population. From a single series of physical experiments, it is possible, with the computational model proposed, to extract parameters controlling the aberrant methylation of the genes under consideration. Using these parameters, we can estimate the level of aberrant methylation resulting from any infection pattern, without requiring a biological cell or a sample from the individual, for all genes investigated during the in vivo experiments. Further development of the model will focus on other genes of potential interest.

These results offer crucial insights towards a better understanding of the initiation of gastric cancers: methylation levels have been observed to be significantly higher in cases of gastric cancer compared to healthy volunteers [START_REF] Maekita | High levels of aberrant DNA methylation in Helicobacter pylori-infected gastric mucosae and its possible association with gastric cancer risk[END_REF], and significantly higher in patients with multiple gastric cancers compared to those with a single gastric cancer [START_REF] Nakajima | Higher methylation levels in gastric mucosae significantly correlate with higher risk of gastric cancers[END_REF].

The ability to model such features represents an important milestone in establishing DNA methylation as a marker for gastric cancer risk. This is likely to be generalised to other types of cancer and tissues: aberrant DNA methylation in non-cancerous tissues was also identified e.g. in the colon, [START_REF] Issa | Methylation of the oestrogen receptor CpG island links ageing and neoplasia in human colon[END_REF], the liver, [START_REF] Kondo | Genetic instability and aberrant DNA methylation in chronic hepatitis and cirrhosis -A comprehensive study of loss of heterozygosity and microsatellite instability at 39 loci and DNA hypermethylation on 8 CpG islands in microdissected specimens from patients with hepatocellular carcinoma[END_REF], and the stomach, (Waki et al., 2002).

A limitation of current in vivo and in vitro studies is that, while successful in investigating specific phenomena, such experiments so far fail to explain system-wide complex interactions. These may reflect overall system complexity, but technical constraints lead most research groups to focus on a single epigenetic change in a given context. The need for integration of these partial results is crucial to understanding the overall biological system, and computer-based modelling can provide a useful complementary framework in this, as in other fields, [START_REF] Dove | Virtual models best live cells at predicting biology[END_REF].

Conclusion

The model, developed here, is able to reproduce and complement in vivo experiments on cell type-dependent aberrant DNA methylation during H. Pylori infection. As such, it provides useful insights on the dynamics of infection-induced aberrant epigenetic changes, which have been linked to gastric cancer initiation.

Costly and time-consuming lab-testing is intrinsic to in vivo epigenetic research. This computer-based study demonstrates that in silico modelling can facilitate quantitative methylation analysis and the investigation of key factors influencing disease initiation. The model can thus play a significant role in developing better understanding of gastric cancers. Notably, as an important milestone in establishing DNA methylation as a marker for gastric cancer risk, it contributes to a better identification of at-risk patients. Furthermore, as the first to successfully reproduce epigenetic events in silico, the current version can be seen as a proof of concept for more advanced models, incorporating several epigenetic changes and the complex interactions between them. In particular, future model versions will take partial methylation levels into account, and will include histone modifications.

Figure 1 :

 1 Figure 1: Structure of a gastric crypt, with one stem cell (black), a few progenitor cells (grey), and, approximately, one hundred differentiated cells on each side

Figure 2 :

 2 Figure 2: Methylation dynamics in gastric crypts. Colour represents cell type, (black for stem cells, grey for progenitor cells, white for differentiated cells). Aberrant methylation is shown with a cross inside the cell.

Figure 3 :

 3 Figure 3: Typical sigmoid function, with a = 1, b = 5, λ = -2, and T = 4.

Figure 4 :

 4 Figure4: Simulated methylation level. Without infection, methylation is almost absent. H. Pylori infection leads to progenitor cell methylation. Unless the stem cell is methylated, this results in rapid but limited fluctuations of the methylation level, since these cells have a finite life span and are routinely replaced with unmethylated cells. The infection also leads, with lower probability, to stem cell methylation, which results in permanent and complete methylation of the corresponding crypt, (and a correspondingly sharp increase of the methylation level of the sample).

Figure 5 :

 5 Figure 5: Sensitivity to α. Each curve corresponds to the average methylation level over 20 simulations, for a specific value of α, (probability of stem cell methylation at each time step).

Figure 6 :

 6 Figure 6: Sensitivity to β. Each curve corresponds to the average methylation level over 20 simulations, for a specific value of β, (probability of progenitor cell methylation at each time step).

Figure 7 :

 7 Figure 7: Sensitivity to λ g and T g . Each curve corresponds to the average methylation level over 20 simulations, for specific values of λ g and T g .

Figure 9 :

 9 Figure 9: Sensitivity to λ g and λ s during repeated 20-week infections.

Figure 10 :

 10 Figure 10: Sensitivity to λ g and λ s during repeated 8-week infections.

Acknowledgements

The authors would like to thank Prof. Toshikazu Ushijima and his team at the National Cancer Center (Tokyo, Japan), for providing access to their results and participating in a number of thought-provoking discussions on Epigenetics.

Financial support from Science Foundation Ireland (Research Frontiers Programme 07/RFP/CMSR724) is warmly acknowledged.

The authors also wish to acknowledge the SFI/HEA Irish Centre for High-End Computing (ICHEC) for the provision of computational facilities and support.