Antti Häkkinen
email: antti.hakkinen@tut.fi

Shannon Healy
email: shannon.healy@tut.fi

Howard T Jacobs
email: howard.t.jacobs@uta.fi

Andre S Ribeiro
email: andre.ribeiro@tut.fi

Genome wide study of NF-Y type CCAAT boxes in unidirectional and bidirectional promoters in human and mouse

A subset of CCAAT boxes are known binding sites for the transcription factor NF-Y. We characterize their number, mismatches to the consensus sequence, and locations in bidirectional and unidirectional promoter sequences in human and mouse. We confront the findings with an analytical null model of DNA sequences and find that NF-Y type CCAAT boxes play key, but distinct roles in the two types of promoters. They are found above chance in both, but in unidirectional only when having few mismatches. In bidirectional, the relative positions of multiple boxes differ from what is expected by chance, suggesting the need for contiguity. In agreement, when there are four boxes (four-box configurations), these have much lower number of mismatches than expected in bidirectional promoters alone. Positioning of the first box differs in the two types of promoters and the null model, and mismatches and positioning are found to be correlated. Finally, four-box configurations are conserved between human and mouse, supporting the relevance of the findings. We conclude that bidirectional and unidirectional promoters, while sharing some similarities, appear to possess distinct regulatory mechanisms at the sequence level.

DNA sequences, Transcription Factor, Transcription Factor Binding Site

Background

In mammalian genomes many pairs of genes are co-regulated by bidirectional promoters [START_REF] Engstrom | Complex loci in human and mouse genomes[END_REF]. These pairs are organized in a head-to-head configuration with transcription start sites (TSS) less than 1 kb apart [START_REF] Trinklein | An abundance of bidirectional promoters in the human genome[END_REF]. The transcriptional activity of the genes in such pairs is, in general, co-regulated [START_REF] Zanotto | The bidirectional promoter of two genes for the mitochondrial translational apparatus in mouse is regulated by an array of ccaat boxes interacting with the transcription factor nf-y[END_REF].

A recent study [START_REF] Zanotto | Modulation of mrps12/sarsm promoter activity in response to mitochondrial stress[END_REF] characterized the bidirectional promoter of two mouse nuclear genes encoding components of the mitochondrial translational apparatus, mitoribosomal protein S12 (Mrps12) and mitochondrial seryl-tRNA ligase (Sarsm). Their TSSs are less than 200 bp apart, the intervening region containing four CCAAT boxes, which are transcription factor binding sites (TFBS) of the transcription factor (TF) NF-Y that facilitate efficient bidirectional transcription. These boxes play an identical role in the homologous human promoter. Altering the sequences of any of the boxes, e.g. by deletion, it is possible to change the relative rates of transcription of the two genes while not necessarily changing their combined transcriptional activity, or vice versa [START_REF] Zanotto | Nf-y influences directionality of transcription from the bidirectional mrps12/sarsm promoter in both mouse and human cells[END_REF].

A similar array of CCAAT boxes interacting with NF-Y regulates the human SOX3 promoter [START_REF] Krstic | Regulation of sox3 gene expression is driven by multiple nf-y binding elements[END_REF], which is unidirectional, suggesting that these boxes may also allow the transcriptional regulation of unidirectional promoters. The existence of multiple NF-Y type CCAAT boxes in promoters is common in cell cycle-regulated genes [START_REF] Colter | Regulation of the human sox9 promoter by the ccaat-binding factor[END_REF][START_REF] Li | Xenopus nf-y pre-sets chromatin to potentiate p300 and acetylation-responsive transcription from the xenopus hsp70 promoter in vivo[END_REF][START_REF] Muller | Chimpanzee, orangutan, mouse, and human cell cycle promoters exempt ccaat boxes and chr elements from interspecies differences[END_REF]. The boxes appear to recruit histone acetylase complexes such as p300 [START_REF] Caretti | Dynamic recruitment of nf-y and histone acetyltransferases on cell-cycle promoters[END_REF][START_REF] Salsi | Interactions between p300 and multiple nf-y trimers govern cyclin b2 promoter function[END_REF] and other TFs [START_REF] Roder | Fire3 in the promoter of the rat fatty acid synthase (fas) gene binds the ubiquitous transcription factors cbf and usf but does not mediate an insulin response in a rat hepatoma cell line[END_REF][START_REF] Imbriano | Direct p53 transcriptional repression: in vivo analysis of ccaat-containing g2/m promoters[END_REF]. Defects or malfunctioning of these CCAAT boxes have been linked to several diseases, including some forms of cancer [START_REF] Farina | Down-regulation of cyclin b1 gene transcription in terminally differentiated skeletal muscle cells is associated with loss of functional ccaat-binding nf-y complex[END_REF][START_REF] Pabst | Dominantnegative mutations of cebpa, encoding ccaat/enhancer binding proteinalpha (c/ebpalpha), in acutemyeloid leukemia[END_REF].

Previous studies suggest that NF-Y type CCAAT boxes are common regulators of bidirectional mammalian promoters [START_REF] Lin | Transcription factor binding and modified histones in human bidirectional promoters[END_REF][START_REF] Zanotto | Nf-y influences directionality of transcription from the bidirectional mrps12/sarsm promoter in both mouse and human cells[END_REF]. However, several questions remain regarding their role as regulators of gene expression. Are there any common patterns in the number, sequence, orientation or spacing of these boxes between different classes of promoters? Are these box arrays also present in unidirectional promoters? If so, is their function the same as in bidirectional promoters? What are the mismatch distributions of these TFBSs in bidirectional and unidirectional promoters and are these mismatches informative of their functioning? Finally, are the sequences equally conserved in both type of promoters, and if not, why?

These questions can only be addressed using a proper "null model", since only by comparing the findings on numbers, positioning, and variability with a null model it is possible to determine whether these features are under selection, or are what is expected by chance.

Here we present the first genome-wide study of the occurrence of NF-Y type CCAAT boxes in bidirectional and unidirectional promoters combined with an exact null model to rigorously evaluate the statistical significance of the findings. Given the importance of NF-Y type CCAAT boxes in the regulation of many genes, it is of interest to analyze their presence at the genomewide level. We perform this analysis in the promoter sequences of human and mouse, by searching and characterizing the number, the mismatches, and locations of these boxes in bidirectional and unidirectional promoters. We then confront the results with a null model here proposed. Our procedure can be used to search for any sequence with a predefined number of possible mismatches, and quantify precisely whether their number and location differ significantly from what is expected by chance.

Methods

Search of NF-Y type CCAAT boxes

The DNA sequences of the complete genome of human (Homo sapiens) and mouse (Mus musculus), were extracted from UCSC Genome Browser (http://genome.ucsc.edu/) database (December 1st, 2008). Promoter sequences were extracted according to the annotations of TSSs.

To extract bidirectional promoters we define a consecutive sequence of nucleotides (substring) of less than 1 kb (kilobase) in length as a bidirectional promoter if and only if there is a TSS downstream of the sequence in both strands and the substring is not a substring of another promoter determined by this criteria [START_REF] Trinklein | An abundance of bidirectional promoters in the human genome[END_REF]. Conversely, a substring of 1 kb in length is defined as a unidirectional promoter if and only if a TSS is downstream of the substring and one (or both) of the following conditions are met: i) the same strand has another TSS annotation within 1 kb, or, ii) the reverse complementary strand has no annotated TSS within 10 kb.

The consensus sequence for NF-Y binding CCAAT box was defined to be 5'-YRRCCAATCA-3' [START_REF] Bi | Dna binding specificity of the ccaat-binding factor cbf/nf-y[END_REF][START_REF] Koessler | Human replication-dependent histone h3 genes are activated by a tandemly arranged pair of two ccaat boxes[END_REF] and also its re-verse complement (i.e. 5'-TGATTGGYYR-3'), according to [START_REF] Mantovani | A survey of 178 nf-y binding ccaat boxes[END_REF], as previously used in [START_REF] Zanotto | Nf-y influences directionality of transcription from the bidirectional mrps12/sarsm promoter in both mouse and human cells[END_REF]. In previous studies, the reverse orientation of the NF-Y binding CCAAT boxes was found to be slightly favored over the forward orientation [START_REF] Mantovani | A survey of 178 nf-y binding ccaat boxes[END_REF]. We note that this is only applies to unidirectional promoters, since in bidirectional promoters the orientation of the box is ambiguous. As certain mismatches in the sequences do not affect the binding significantly [START_REF] Mantovani | A survey of 178 nf-y binding ccaat boxes[END_REF][START_REF] Zanotto | Nf-y influences directionality of transcription from the bidirectional mrps12/sarsm promoter in both mouse and human cells[END_REF], we allow mismatches to some extent in all positions but the CCAAT box, as evidence suggest that this box is strictly required for NF-Y to bind [START_REF] Dorn | A multiplicity of ccaat box-binding proteins[END_REF][START_REF] Van Huijsduijnen | Properties of a ccaat box-binding protein[END_REF]. In case of overlapping matches, only the first is recorded. Overlapping is possible and is more likely as the number of allowed mismatches increases. Such possibility is accounted for in the null model.

Generalized null model

The significance of the findings is accessed by statistical hypothesis testing. Let Σ be a finite alphabet and let there be γ ∈ N random strings, denoted by R i i . Let the length of the ith string be iid (independent and identically distributed) i ∼ D , and let each character be iid c ik ∼ D c . The D and D c are arbitrary univariate discrete distributions with support of ∈ N and c ∈ Σ, respectively. To determine the probability that a random string R i i contains one or more sequences S 1 , . . . , S k as a substring, say P[R i i ∈ L], we find the regular language L and build a deterministic finite state automaton (DFA) that recognizes that language.

Let M = (Q, q 0 , A, Σ, δ) be the DFA, Q being the finite set of states, q 0 ∈ Q the initial state, A ⊆ Q the accepting states, Σ the alphabet, and δ(q, c) : Q × Σ → Q the transition function. Let P be the stochastic matrix of size |Q| × |Q| corresponding to the DFA M , i.e. if I(q) : Q → {1, . . . , |Q|} is an arbitrary bijective mapping P is defined as follows:

∀ q, q ∈ Q : [P] I(q),I(q ) = q =δ(q,c)

P [c ik = c] ,
(1)

where P[c ik = c] is the probability of a random character c ik ∈ Σ being c.

The probability for the DFA M to transit from state i to state j through exactly k transitions is given by [P k ] i,j , thus:

P[R i i ∈ L] = e I(q 0 ) * P i q∈A e I(q) , ( 2
)
where e i is the unit vector having the ith element equal to 1. We construct the DFA M using the Aho-Corasick algorithm [START_REF] Aho | Efficient string matching: An aid to bibliographic search[END_REF], for which |Q| is proportional to the sum of the length of the strings. To obtain the DFAs for multiple matches, we concatenate several of the DFAs. Next, e.g. given that k out of γ sequences met a certain condition, we use the two-tailed hypothesis testing to determine if the findings differ from the by-chance findings suggested by the null model. The p-value is:

p =1 + F (E[K] -|k -E[K]|) -F (E[K] + |k -E[K]|) , ( 3
)
where K ∼ D K is a random variable representing the number of promoters meeting the criterion under the null hypothesis H 0 , E[K] is its expected value, and

F (k) = P[K ≤ k] is the cumulative distribution function (CDF) of the distribution D K . It follows that the distribution is a binomial Bin(γ, P[R D ∈ L]
), where:

P[R D ∈ L] = ∞ j=0 P[ = j] P[R j i ∈ L] , (4)
is the probability mass of the mixture density of the probability mass functions (PMF) of different promoter lengths.

To determine the spatial probability that the nth match is at position k, we construct languages L i such that L i recognizes i (or more) matches to the consensus. Denoting a random string of length by R , the probability is:

P[•] = P[k ≤ ] P[R k ∈ L n ] -P[R k-1 ∈ L n ] ,
(5) for a given length distribution D . The position k denotes the ending position of a match. To obtain the distance between (n -1)th and nth match, we modify the length distribution s.t. instead of D , we use the distribution of the remaining lengths after exactly (n-1) matches, D . Eq. 5 can be readily evaluated using the PMF obtained using the convolution:

f (k ) = ∞ k=k +1 f (k) P[R k-k ∈ L n-1 ] -P[R k-k -1 ∈ L n-1 ] , (6)
Several kinds of null hypotheses about DNA sequences have been used [START_REF] Fitch | Random sequences[END_REF]. Typically, one or more sequences of non-related DNA are used, either as-is or randomly shuffled. Random shuffling guarantees the independence between nucleotides but this is not strictly required, as our method can be extended for Markov chains of order up to min i |S i |, and, with some modifications, even higher.

It is known that in many promoters in human and mouse, some dinucleotides are more common than others [START_REF] Bajic | Mice and men: Their promoter properties[END_REF]. The different multinucleotide distributions can be imposed to the null model (e.g. using first-or second-order Markov chains). Instead, we use the alphabet Σ = {A, C, G, T}, and the experimentally observed PMFs for the nucleotide distribution D c and the promoter length distribution D . This is because it is unknown whether and to what extent multinucleotide distributions in mouse and human DNA sequences are due to physical-chemical constrains or due to evolutionary pressure. Since we aim to use a null model whose sequence is absent of any evolutionary pressure, we opted for not imposing any dependence in the distribution of nucleotides.

In that regard, for example, while many promoters have CpG islands and TATA boxes, many others do not [START_REF] Bajic | Mice and men: Their promoter properties[END_REF]. This implies that they do not have to be necessarily present due to any physical-chemical constraints, and that selection can remove or add these motifs and alter the ratio of dinucleotides. Therefore, our null model should not assume a priori these distributions, if it is to be used to determine if a sequence has been selected in or out.

Results and Discussion

NF-Y type CCAAT box occurrence

First, we searched for NF-Y type CCAAT boxes in human and mouse promoters, as NF-Y (also called CBF) is currently considered the principal CCAAT-box-binding activator [START_REF] Mantovani | A survey of 178 nf-y binding ccaat boxes[END_REF]. The search for the consensus-binding sequence of NF-Y was made in both strands, allowing varying number of mismatches in specific nucleotides, as described. The findings, their corresponding expected probabilities given by the null model, and the hypothesis testing are shown in Table 1.

The occurrence of NF-Y type CCAAT boxes is similar in human and mouse although substantially different from their null models. In bidirectional promoters, the boxes were found above chance. This finding is in agreement with a study using chromatin-immunoprecipitation (ChIP) on 118 regions in the human genome that showed that CCAAT boxes, among a few other small motifs, are over-represented in bidirectional promoters, in comparison to unidirectional ones [START_REF] Lin | Transcription factor binding and modified histones in human bidirectional promoters[END_REF]. Interestingly, our analysis of a much larger set of promoters further shows that in unidirectional promoters, the boxes are also found above chance, provided that they have a low number of mismatches from consensus. However, unlike in bidirectional, they were found below chance when a large number of mismatches from consensus were allowed. The striking differences between our findings in bidirectional and unidirectional promoters suggest that the boxes perform different functions in the two types of promoters.

We hypothesize that, in general but not exclusively, each box functions as an independent TFBS in unidirectional promoters, while in bidirectional there are usually arrays of multiple boxes, as in Mrps12/Sarsm in both mouse and human [START_REF] Zanotto | The bidirectional promoter of two genes for the mitochondrial translational apparatus in mouse is regulated by an array of ccaat boxes interacting with the transcription factor nf-y[END_REF][START_REF] Zanotto | Nf-y influences directionality of transcription from the bidirectional mrps12/sarsm promoter in both mouse and human cells[END_REF] which possesses four such boxes. NF-Y binding is likely to determine the ratio of expression of the two genes in such pairs, depending on how it interacts with the boxes, which depends on their number, mismatches, and positioning in the template. This was shown to be the case in Mrps12/Sarsm promoter, as deleting any of the boxes alters the ratio between the transcription rates of the two genes [START_REF] Zanotto | Nf-y influences directionality of transcription from the bidirectional mrps12/sarsm promoter in both mouse and human cells[END_REF]. Further, the deletion of boxes usually also affects the ratio by which each of the two genes is expressed, rather than their combined expression rate. Importantly, no single deletion blocked transcription completely [START_REF] Zanotto | The bidirectional promoter of two genes for the mitochondrial translational apparatus in mouse is regulated by an array of ccaat boxes interacting with the transcription factor nf-y[END_REF][START_REF] Zanotto | Modulation of mrps12/sarsm promoter activity in response to mitochondrial stress[END_REF]. This suggests that one can expect several configurations of these boxes in promoters, instead of a single fixed configuration.

These observations are supported by our finding that the boxes in bidirectional promoters can sustain more mismatches from consensus. The occurrence, above chance, of boxes with many mismatches indicates that they maintain some functionality (otherwise they would not exist beyond what is expected by chance), and act co-operatively. On the other hand, since in unidirectional promoters boxes with high number of mismatches are "selected out" (i.e. are found below chance), the same is likely not to occur in unidirectional promoters. Interestingly, if boxes with many mismatches were totally non-functional in unidirectional promoters, their amount should be similar to what is predicted by the null model. Instead, they were found below chance. We suggest that such boxes with mismatches allow a TF, other than NF-Y, to bind, creating a co-operative array which would be pernicious in some way to the function of unidirectional promoters.

First box position

The boxes form co-operative arrays in some bidirectional promoters. Thus, we studied the spatial locations of the boxes in the two types of promoters, starting with the location of the first box found upstream of each TSS in human and mouse promoters.

Results for human and mouse are shown in Fig. 1, and are almost identical. The probability distributions for the locations of NF-Y type CCAAT boxes in bidirectional promoters are shown in the top plot of Fig. 1. There is a region adjacent to the 5' end of the TSS (from -50 to 0) that is less prone than in the null model to contain the first box, suggesting that, during evolution, the boxes have been selectively eliminated from this region; alternatively, some other class of sequence element is constrained to be located in this region. Finally, another possible explanation is that placing a CCAAT box in this region causes a change in the selection of the TSS as reported in [START_REF] Kabe | Nf-y is essential for the recruitment of rna polymerase ii and inducible transcription of several ccaat box-containing genes[END_REF], since NF-Y plays a role in the recruitment of RNA polymerase II. NF-Y type CCAAT boxes were most common in the subsequent region, between -120 and -50 (particularly around position -85). Beyond that, several small peaks were found, e.g. around -200, -350 (in human only), and -580. Beyond -600, NF-Y type CCAAT box occurrence was the same as expected by chance, implying that such boxes located at this distance or greater from the TSS might not be functional.

As in the bidirectional case, in unidirectional promoters (bottom plot in Fig. 1) there was also a region where the first box was preferentially located, again around -85, with a comparatively barren region closer to the TSS. Beyond that and until approximately -800, the findings appear to match the predictions of the null model. Next, there is a larger peak centered on -930. Specific reasons for the existence of a peak in this location are unknown. It is noted that the y-axis scales on the upper and lower figures are not the same, and that the findings in unidirectional thus differ less from the null model than in the case of bidirectional promoters.

Subsequent box positions

We determined the relative positions of each subsequent NF-Y type CCAAT box (when present) up to the fourth box relative to the position of the first such element. The positions of these boxes for bidirectional promoters are presented in Fig. 2. The results for human and mouse are similar and show significant "non-random" positioning. In general, each subsequent NF-Y type CCAAT box appeared to be preferentially located at approximately -40 from the previous box. The region immediately adjacent to the previous box was strongly disfavored. These features were found to be slightly more robust in mouse than in human. The -40 peak was also found in unidirectional promoters, but was far less pronounced peak. Thus, we conclude that, at least in bidirectional promoters, NF-Y type CCAAT boxes function in combination, in sets of two to four boxes, with very specific distances between them. Relative positioning is thus inferred to be a condition for functionality.

It may be noted that the chance of finding two, three or four boxes arrayed with this precise spacing is increasingly small. Thus, the relative positioning of multiple boxes, when present, appears to be under stronger selection than the position of the first box. In unidirectional this is less frequent, as depicted by Fig. 3. Interestingly, in bidirectional, the strongest correlation in position is between the third and the fourth boxes, whereas in unidirectional the probability of precise spacing is similar between all boxes.

We found that in mouse, on average, the distance between the boxes is slightly smaller than in human, especially the relative positioning of boxes 2, 3 and 4. In general, the box positions approximately match the configuration in the Mrps12/Sarsm bidirectional promoter [START_REF] Zanotto | Nf-y influences directionality of transcription from the bidirectional mrps12/sarsm promoter in both mouse and human cells[END_REF], with four consecutive, equidistant boxes (≈ 40 nucleotides apart), suggesting that NF-Y binding to a precisely spaced array of boxes might commonly regulate the expression of genes driven by bidirectional promoters.

Four-box configurations

We analyzed promoters with exactly four NF-Y type CCAAT boxes, allowing up to 5 mismatches from consensus in each such element. Previous studies of the Mrps12/Sarsm bidirectional promoter [START_REF] Zanotto | Nf-y influences directionality of transcription from the bidirectional mrps12/sarsm promoter in both mouse and human cells[END_REF] showed that NF-Y binds to a specific spatial configuration of four consecutive boxes, approximately 40 nucleotides apart from each other. However, measurement of the transcriptional activity in promoters with boxes deleted showed that most single and double deletions do not completely hamper expression [START_REF] Zanotto | The bidirectional promoter of two genes for the mitochondrial translational apparatus in mouse is regulated by an array of ccaat boxes interacting with the transcription factor nf-y[END_REF][START_REF] Zanotto | Modulation of mrps12/sarsm promoter activity in response to mitochondrial stress[END_REF][START_REF] Zanotto | Nf-y influences directionality of transcription from the bidirectional mrps12/sarsm promoter in both mouse and human cells[END_REF], although they can change the ratio of expression levels of the two genes of the pair. Thus, it is likely that multiple variations of the box configuration (in both number and position of boxes) in the Mrps12/Sarsm promoter might still be under regulation by NF-Y.

The locations of boxes in bidirectional promoters with four boxes are shown in Fig. 4. The relative positions, in general, differ significantly from what would be expected by chance. In general, they are contained within small regions as expected given the preferential distance between consecutive boxes (Fig. 2). As an example, 50% of the sets of four boxes in human and 38% in mouse are within a range of 200 nucleotides or less, significantly above from what would be expected by chance (Table 2).

In unidirectional promoters (figure not shown), the higher variability of distances between the boxes suggests that their relative position does not follow a specific pattern. This implies that, in general, they are independent TFBSs, which are nevertheless likely to be functional as they were found to occur above chance.

We measured the mean μ and variance σ 2 of distances between consecutive boxes and the correlation ρ between their positions (Table 2). In bidirectional promoters the locations of the boxes was highly correlated, above 0.9. The only exception was between the first and second boxes in human, which are weakly correlated. One possible explanation is that, in human, three-box configurations might be more common than four-box ones, compared with mouse, with any additional (fourth) box functionally superfluous. However, it is noted that such proportions would change slightly if one had considered independently cases with five or more boxes. Such promoters might have four "active" boxes and "superfluous" ones.

Interestingly, the correlation between the box positions in unidirectional promoters is slightly higher than in the null model (but much lower than in bidirectional promoters). One possibility is that some apparently unidirectional promoters are in fact bidirectional, with the "other" direction used to transcribe e.g. a short non-coding RNA that has not yet been mapped.

Generally, the findings suggest that, in unidirectional promoters, NF-Y type CCAAT boxes do not function in the same manner as in bidirectional promoters since, in the unidirectional case, there is much less conservation across promoters and correlation between the positions of different boxes in each promoter.

We studied the degree of conservation of the four-box configurations between human and mouse. Although there is some evidence suggesting that certain dinucleotides might mutate at faster rates than others, it is unknown if such higher rates result in faster changes on these sites in long time scales, as that depends on if the changes are conserved or not. Thus, rates of mutation are not considered when computing significance of findings. The significance of findings is made by confronting the conservation rate of the boxes with the mean conservation rate of all nucleotides of the promoters.

We first located each human (mouse) bidirectional promoter and their neighboring genes. Then, using Ensembl (http://www.ensembl.org/), we obtained the protein sequences produced by the sense and antisense genes, and mapped them to the best matching mouse (human) protein sequences using BLAST [START_REF] Altschul | A basic local alignment search tool[END_REF], thus obtaining the corresponding genes. The intervening region between the two genes was selected for further analysis. We found that 18 out of 26 cases in human and 14 out of 21 cases in Next, for the cases that could be mapped, we obtained the optimal alignment between nucleotides of the pairs of promoters using the Smith-Waterman algorithm [START_REF] Smith | Identification of common molecular subsequences[END_REF], and from there, the differences between the sequences (fraction of deletions, insertions and conserved nucleotides). Results are shown in Table 3.

We performed an exact multinomial test to obtain a p-value to test if the number of deletions, insertions and conservations in the boxes is similar to that of all nucleotides in the promoter or differs significantly, i.e. to test whether the boxes were more strongly conserved than average in cases i) and ii). The small p-values indicate that the differences between fractions of deleted, inserted and conserved nucleotides between nucleotides that form the boxes and all other nucleotides are significant (significance level α = 0.01). From the comparison, in both human and mouse, NF-Y type CCAAT boxes are strongly conserved.

Table 3: Fraction of nucleotides and CCAAT boxes deleted, inserted, and conserved (%) in the aligned subregions when mapping i) the 18 human bidirectional promoters containing exactly four NF-Y type CCAAT boxes to mouse bidirectional promoters and ii) the 14 mouse promoters to human promoters, and the p-values for the distributions to match.

Number of mismatches

We investigated the distribution of the number of mismatches. First, we studied these distributions as a function of the position relative to the 5' end of the TSS. In Fig. 5 are the results for bidirectional and in Fig. 6 for unidirectional promoters.

In both types of promoters, the mismatch distributions differ significantly from what is expected by chance. However, the spatial location also seems relevant. Both in unidirectional and bidirectional promoters there were fewer mismatches in the regions where the boxes are preferred, especially in unidirectional (Fig. 6). In bidirectional this was also true but the result is less robust as the sample size is smaller, especially in the regions far from the TSS. Nevertheless, it appears that the regions where the boxes are found above expectation is also where the number of mismatches from consensus is lower, consistent with selection for functionality.

In bidirectional promoters, the distributions differ up to 5 mismatches. However, in unidirectional promoters the distributions for more than 3 mismatches are very similar to the null prediction. This supports our hypothesis, that due to the boxes functioning in combination in bidirectional promoters, more mismatches are allowed. In unidirectional, it appears that for more than 3 mismatches, the boxes are non-functional, as they do not function in co-operation with other boxes.

We analyzed the number of mismatches within four-box configurations. Assuming that, in many cases, the boxes work in combination in bidirectional promoters, but not in unidirectional promoters, we hypothesized that it would be unlikely to find promoters with four boxes with the same mis- match distributions in the two types of promoters. Since in unidirectional promoters the boxes are proposed to function as independent TFBSs, we expected results similar to the null model. Conversely, we hypothesized that the mismatch distribution in NF-Y type CCAAT boxes in bidirectional promoters with four boxes would be far from what is expected by chance at the single box level, as we expect the mismatches to be evenly distributed by the four boxes, which is not what the null model predicts. The results of this test are shown in Fig. 7 and confirm this prediction.

The results support the hypothesis. In bidirectional promoters the number of sets of four boxes where all boxes have 3 or less mismatches was significantly above chance or, equivalently, the cases where most boxes have a large number of mismatches was significantly below chance, whereas in unidirectional no clear difference from the null model was found.

Specificity of the findings to NF-Y versus other TFBSs

We also performed additional analysis to study to what extent the findings presented are specific to NF-Y, compared to other similar TFBSs. While allowing for mismatches is necessary to determine the possible binding sites of NF-Y [START_REF] Zanotto | Nf-y influences directionality of transcription from the bidirectional mrps12/sarsm promoter in both mouse and human cells[END_REF], it may also allow the nucleotide sequences identified as binding sites for NF-Y to also be binding sites for another TF, in the case where the two TFBS are similar. Such an eventuality becomes more likely, the higher the number of mismatches allowed in the search. Due to this, it is necessary to estimate the probability of this eventuality.

We verified that only a small fraction of our findings are potential binding sites of other TFs. This was done by comparing the findings in human and mouse genomes using the NF-Y binding site consensus with overlapping findings of two other TFBSs, namely the binding sites of transcription factors C/EBP and CTF/NF-I. Their consensus sequences were defined to be 5'-RTTGCGYAAY-3' for C/EBP [START_REF] Osada | Dna binding specificity of the ccaat/enhancer-binding protein transcription factor family[END_REF], and 5'-TTGGCNNNNNGCCAA-3' for CTF/NF-I [START_REF] Roulet | Experimental analysis and computer prediction of ctf/nfi transcription factor dna binding sites[END_REF].

The results in Table 4 indicate that regardless of the number of mismatches allowed to the NF-Y consensus sequence, there is only a small probability that an overlapping match with low number of mismatches to a consensus sequence of other type is present. Namely, for any class of the NF-Y findings, there is less than 2% chance for an overlapping finding of the other TFBSs with no more than 1 mismatch, and less than 10% chance for no more than 2 mismatches. Thus, we expect that very few of the findings reported by a match to the NF-Y consensus defined, are functional TFBSs of TFs other than NF-Y.

Comparison with alternative TFBS mapping techniques

Instead of the consensus-with-mismatches method used here, it is also possible to use other methods to locate putative TFBSs. Position weight matrix (PWM) based techniques (e.g. [START_REF] Staden | Computer methods to locate signals in nucleic acid sequences[END_REF] are one such alternative, allowing fuzzy or probabilistic matching, at the expense of higher computational cost. To find if such a matrix is well approximated by our choice of consensus sequence, we searched the promoters also with a PWM constructed from the findings reported in [START_REF] Bi | Dna binding specificity of the ccaat-binding factor cbf/nf-y[END_REF]. The PWM was constructed such that the PWM score corresponds to the log-odds of the respective sequence over the background nucleotide distribution present in each set of promoters. For each possible value of threshold we then obtained the 90% confidence interval for the match count, and tested if one of our findings with consensuswith-mismatch technique lies within that interval. This was done by using the Beta-Binomial distribution. The set of ranges of thresholds that the consensus-with-mismatch findings correspond to are presented in Table 5. We found that for each number of mismatches, there is a range of PWM thresholds for which the results of the two methods do not differ significanly.

Conclusions

From a genome wide survey of promoter sequences in human and mouse we found that the amount of NF-Y type CCAAT boxes in both unidirectional and bidirectional promoters differ significantly from what is expected by chance, indicating that they play a key role. Further, since the boxes occupy specific locations in bidirectional, more than in unidirectional promoters, their performance is likely to depend on distance to the TSS.

Strikingly, in bidirectional promoters there is a strong correlation between the locations of boxes when multiple boxes exist. This correlation is far less pronounced in unidirectional promoters although still slightly higher than if they were randomly located. Further, in bidirectional promoters, many boxes have an arbitrary number of mismatches, while in unidirectional, only boxes with a low number of mismatches are found significantly above chance.

The findings support the hypothesis that in bidirectional promoters the boxes function as a co-operative array, while in unidirectional they are independent TFBSs (whose position is under selection but not relative to one another). This suggests that these boxes function inherently differently in unidirectional versus bidirectional promoters, and that in the latter case the binding requires multiple boxes. This would explain why the number of mismatches from consensus, per box, is evenly distributed and higher in bidirectional promoters, since due to co-operation, mismatches on one box affect less the total binding strength between TF and TFBS, than if only one box was required for the binding.

The above conclusions are similar in human and mouse, suggesting similar functions in both, in agreement with the strong conservation observed in sets of four boxes in bidirectional promoters in human and mouse. The results on the location, number and mismatch tolerance in bidirectional promoters, its similarity in mouse and human, and the strikingly different results in unidirectional promoters, suggest that these two types of promoters are likely to be subject to distinct regulatory mechanisms.

The variability in positioning, mismatches and number of boxes in bidirectional promoters with multiple boxes, in both human and mouse, suggests that this regulatory mechanism has plasticity and multi-functionality, i.e. it allows different pairs of genes to differ in total and relative expression of the two genes and might facilitate interactions with other TFs, as suggested in [START_REF] Dorn | A multiplicity of ccaat box-binding proteins[END_REF].

In the methods, we mentioned that our null model should be free from any constrain in the ordering of the existing nucleotides, as there is no known physical constrain that would impede a certain sequence. The only imposed condition to our null model was the ratio of each nucleotide, which was set to be identical to those in human and mouse. It is of importance to note that we did not necessarily had to impose this as well. In fact, most likely, these ratios were also imposed by selection, as it selected for certain sequences more than others, instead of being due to physical-chemical constraints. For that reason, we also confronted the sequences of human and mouse with a null model where all nucleotides exist in equal ratios. The results do not differ, qualitatively, from those here presented.

Finally, the conservation of sets of four boxes between mouse and human and the overall number of boxes in both genomes confirms their widespread use. This strengthens the conclusions of previous studies [START_REF] Kabe | Nf-y is essential for the recruitment of rna polymerase ii and inducible transcription of several ccaat box-containing genes[END_REF] about the critical role of NF-Y as a regulator of transcription of many bidirectional promoters [START_REF] Mantovani | A survey of 178 nf-y binding ccaat boxes[END_REF], and of CCAAT boxes as binding sites for several other TFs [START_REF] Dorn | A multiplicity of ccaat box-binding proteins[END_REF].

Figure 1 :

 1 Figure 1: First box position relative to the TSS, up to 3 mismatches. Upper figure bidirectional promoters, lower figure unidirectional promoters. Experimental results were smoothed using moving average window of 10.

Figure 2 :

 2 Figure 2: Second (top), third (middle), and fourth (bottom) box positions relative the previous box, up to 3 mismatches, in bidirectional promoters. Results from the searches in mouse and human sequences were smoothed using a moving average window of 10.

Figure 3 :

 3 Figure 3: Second (top), third (middle), and fourth (bottom) box positions relative the previous box, up to 3 mismatches, in unidirectional promoters. Results from the searches in mouse and human sequences were smoothed using a moving average window of 10.

Figure 4 :

 4 Figure 4: All four-box configurations found in bidirectional promoters, sorted by the position of the first box. Human (upper plot) and mouse (lower plot). First box in black and the next boxes are in increasingly lighter grey.

Figure 5 :

 5 Figure 5: Distribution of mismatches in bidirectional promoters in each position. Upper figure for human, and lower for mouse. Results were smoothed using moving average window of 10. Bar on the right is the null prediction. Black: no mismatches, lightest gray: 5 mismatches.

Figure 6 :

 6 Figure 6: Distribution of mismatches in unidirectional promoters in each position. Upper figure for human, and lower for mouse. Results were smoothed using moving average window of 10. Bar on the right is the null prediction. Black: no mismatches, lightest gray: 5 mismatches.

Figure 7 :

 7 Figure 7: Distribution of number of boxes with up to 3 mismatches in bidirectional (top) and unidirectional (bottom) promoters with four boxes.

Table 1 :

 1 Percentage of promoters meeting each criterion in human, mouse, and null model promoter sequences, and the results of hypothesis testing at significance level α = 0.01. Percentages of findings above chance in dark gray background, below chance in light gray background and, if the null hypothesis H 0 is not rejected, the background is white. Criteria of type k-N-ε denote k or more NF-Y type CCAAT boxes, up to ε mismatches each.

		Human	Mouse	Human (null) Mouse (null)
	Criterion	Bi	Uni	Bi	Uni	Bi	Uni	Bi	Uni
	γ	1718 37201 1447 31407 N/A N/A N/A N/A
	1-N-0	4.77	2.91	7.81	4.22	0.29	1.48	0.32 1.54
	1-N-1	16.82 13.95 21.15 16.04 3.06 13.84 3.28 14.41
	1-N-2	27.65 36.45 32.41 39.14 11.53 44.94 12.26 46.40
	1-N-3	34.81 57.42 39.74 61.75 22.30 72.33 23.56 73.91
	1-N-4	38.77 68.23 42.78 72.99 28.62 83.05 30.16 84.38
	1-N-5	39.46 70.44 43.33 75.22 30.17 85.16 31.77 86.41
	2-N-0	0.41	0.21	0.62	0.25	0.00	0.01	0.00 0.01
	2-N-1	4.54	2.20	6.57	2.50	0.07	0.99	0.08 1.07
	2-N-2	10.77 10.37 13.06 10.95 1.08 11.87 1.21 12.76
	2-N-3	16.18 24.84 19.49 28.14 4.15 36.30 4.60 38.35
	2-N-4	17.99 35.81 22.11 40.66 6.97 52.40 7.65 54.81
	2-N-5	18.86 38.47 22.60 43.63 7.77 56.27 8.52 58.71
	3-N-0	0.06	0.04	0	0.04	0.00	0.00	0.00 0.00
	3-N-1	1.40	0.58	1.31	0.41	0.00	0.05	0.00 0.05
	3-N-2	4.89	2.88	5.74	2.98	0.08	2.16	0.10 2.42
	3-N-3	7.45	9.10	9.26 10.40 0.67 13.32 0.80 14.64
	3-N-4	8.85 15.98 10.44 18.67 1.49 25.27 1.75 27.40
	3-N-5	8.96 17.93 10.78 20.90 1.77 28.75 2.07 31.08
	4-N-0	0	0.01	0	0.01	0.00	0.00	0.00 0.00
	4-N-1	0.70	0.17	0.55	0.11	0.00	0.00	0.00 0.00
	4-N-2	1.98	0.99	2.28	0.85	0.01	0.29	0.01 0.34
	4-N-3	3.73	3.17	4.22	3.36	0.09	3.77	0.12 4.32
	4-N-4	4.42	6.38	5.04	7.34	0.27	9.60	0.35 10.83
	4-N-5	4.42	7.45	5.11	8.59	0.35 11.64 0.44 13.09

Table 2 :

 2 Spatial features of four-box configurations in unidirectional and bidirectional promoters along with the null predictions. Null prediction was obtained using Monte-Carlo simulations rather than analytically. X i denote random variables of the position of ith box and P[•] = P[X 4 -X 1 ≤ 200] the fraction of four-box promoters where boxes separated by 200 bp or less.

		Human	Mouse	Human (null) Mouse (null)
	Measure	Bi	Uni	Bi	Uni	Bi	Uni	Bi	Uni
	P[•]	0.50 0.10 0.38 0.08 0.06	0.03	0.03 0.02
	ρ X 1 ,X 2	0.67 0.70 0.95 0.72 0.41	0.61	0.66 0.60
	ρ X 2 ,X 3	0.99 0.72 0.95 0.78 0.77	0.66	0.73 0.66
	ρ X 3 ,X 4	0.95 0.71 0.92 0.67 0.74	0.60	0.68 0.63
	μ X 2 -X 1 140.5 180.3 74.0 181.8 191.8 199.1 173.7 202.5
	μ X 3 -X 2	65.1 187.1 93.9 178.8 150.3 204.1 150.1 201.5
	μ X 4 -X 3	72.2 175.5 81.5 190.1 151.2 200.8 184.0 193.1
	σ X 2 -X 1 181.0 170.1 43.3 158.5 150.6 160.8 121.7 164.1
	σ X 3 -X 2	30.4 174.1 75.4 149.8 116.2 161.9 118.2 160.8
	σ X 4 -X 3	78.3 165.5 68.2 175.1 130.9 160.3 135.2 152.9
	mouse (i.e. 69.2% and 66.7%, respectively) could be mapped to the other
	genome.

Table 4 :

 4 Distributions of the minimum number of mismatches to any overlapping C/EBP or CTF/NF-I consensus sequences, in percentage, for findings with various number of mismatches in NF-Y consensus. N-ε denotes that up to ε mismatches are allowed to the NF-Y consensus sequence.

	Criterion			Mismatches to C/EBP consensus
		0	1	2	3	4	5	6
	N-0	0	0.21 7.52 41.11 47.03 4.13	0
	N-1	0	0.42 7.56 39.59 47.44 4.98 0.01
	N-2	0.01 0.51 7.68 37.88 48.34 5.57 0.00
	N-3	0.00 0.59 7.44 35.84 49.48 6.64 0.01
	N-4	0.00 0.48 6.76 34.99 50.42 7.33 0.01
	N-5	0.00 0.46 6.50 34.37 51.03 7.63 0.01
	Criterion		Mismatches to CTF/NF-I consensus
		0	1	2	3	4	5	6	7
	N-0	0.32 0.46 9.19 24.19 46.49 18.56 0.78	0
	N-1	0.55 0.60 6.98 24.85 44.95 20.98 1.08 0.01
	N-2	0.47 0.68 5.96 23.15 44.43 23.80 1.52 0.00
	N-3	0.44 0.57 4.93 21.33 44.55 26.47 1.71 0.00
	N-4	0.42 0.48 4.40 19.95 44.68 28.16 1.90 0.01
	N-5	0.42 0.46 4.26 19.66 44.60 28.69 1.92 0.01

Table 5 :

 5 Intervals of PWM thresholds for which the results of the consensus-with-mismatch findings lie within the 90% confidence interval. N-ε denotes that up to ε mismatches are allowed to the NF-Y consensus sequence.

	Criterion PWM threshold range
	N-0	[9.9381, 10.0121]
	N-1	[7.7843, 7.8172]
	N-2	[5.3332, 5.3931]
	N-3	[2.5715, 2.6928]
	N-4	[-0.5868, -0.3383]
	N-5	(-∞, -2.6060]

Acknowledgements

We thank the Academy of Finland, TEKES, the Sigrid Juselius Foundation, and Tampere Univ. Hospital Medical Research Fund.


Vicenç Méndez
email: vicenc.mendez@uab.es

Isaac Llopis

Daniel Campos

Werner Horsthemke

Effect of environmental fluctuations on invasion fronts

Keywords: Invasion front, Environmental fluctuations, Allee effect

We determine the density profile and velocity of invasion fronts in one-dimensional infinite habitats in the presence of environmental fluctuations. The population dynamics is reformulated in terms of a stochastic reaction-diffusion equation and is reduced to a deterministic equation that incorporates the systematic contributions of the noise. We obtain analytical expressions for the front profile and velocity by constructing a variational principle. The effect of the noise differs, depending on whether it affects the density-independent growth rate, the intraspecific competition term or the Allee threshold. Fluctuations in the density-independent growth rate increase the invasion velocity and the population density of the invaded area.

Fluctuations in the competition term also change the population density of the invaded area, but modify the invasion velocity only for certain initial conditions.

Fluctuations in the Allee threshold can induce pulled or pushed invasion fronts as well as invasion failure. We compare our analytical results with numerical solutions of the stochastic partial differential equations and show that our procedure proves useful in dealing with reaction-diffusion equations with multiplicative noise.

INTRODUCTION

Biological invasions correspond to the introduction and spread of exotic organisms into regions outside of their native range. They have emerged as a major environmental, economic and public health problem tied to the rapid, ongoing expansion of international trade and travel and have been studied for many decades [START_REF] Fisher | The Wave of Advance of Advantageous Genes[END_REF][START_REF] Kolmogorov | Study of the equation of diffusion equation combined with an increase of matter and its application to a biological problem[END_REF][START_REF] Skellam | Random Dispersal in Theoretical Populations[END_REF][START_REF] Hengeveld | Dynamics of Biological Invasions[END_REF]Williamson, 1996;[START_REF] Shigesada | Biological Invasions: Theory and Practice[END_REF][START_REF] Murray | Mathematical Biology[END_REF]. The invasion process can be regarded as a traveling wave front that propagates through a homogeneous medium, and it can be described mathematically by reactiondiffusion equations [START_REF] Shigesada | Biological Invasions: Theory and Practice[END_REF][START_REF] Murray | Mathematical Biology[END_REF][START_REF] Ortega-Cejas | Role of the delay time in the modelling of biological range expansions[END_REF][START_REF] Méndez | Reaction-Transport Systems[END_REF]. These equations can be generalized to take into account memory effects in the animal motion for logistic dynamics [START_REF] Ortega-Cejas | Role of the delay time in the modelling of biological range expansions[END_REF]) and dynamics with Allee effect [START_REF] Méndez | Wavefronts in bistable hyperbolic reactiondiffusion systems[END_REF].

Another important topic in ecology is the role of fluctuations in population dynamics [START_REF] Horsthemke | Noise-Induced Transitions[END_REF]. Populations with a small number of individuals experience internal fluctuations and require a stochastic description. Even if these fluctuations can be neglected because the population contains a large number of individuals, other stochastic factors are usually present and can affect significantly the dynamic behavior of the population. Examples are demographic stochasticity, environmental stochasticity, demographic heterogeneity and stochastic sex determination. Various authors have investigated the effect of random fluctuations on population dynamics.

Marcus studied a stochastic partial differential equation (SPDE), namely a reaction-diffusion equation with logistic growth and multiplicative spatiotem-poral Gaussian white noise [START_REF] Marcus | A stochastic logistic diffusion equation[END_REF]. He carried out a small noise expansion and established the existence and uniqueness of a weak solution of the SPDE. More recently, Malchow and coworkers have studied numerically the effect of external noise on pattern formation in reaction-diffusion models of plankton dynamics (Malchow et al, 2004a;Malchow et al, 2004b). Scarsoglio and coworkers have used analytical and numerical tools to investigate the effects of external fluctuations on vegetation patterns [START_REF] Scarsoglio | Spatial pattern formation induced by Gaussian white noise[END_REF].

Already [START_REF] Skellam | Random Dispersal in Theoretical Populations[END_REF] analyzed the effects of random dispersal on population spread and applied his results to the spread of oak in Britain and the muskrat in Central Europe. Mollison (1991) used linear stochastic models to estimate invasion velocities. The author stressed the limitations of linear models and the need for a better understanding of nonlinear stochastic models.

Lewis similarly stressed the importance of stochastic factors for biological invasions [START_REF] Lewis | Spread rate for a nonlinear stochastic invasion[END_REF] and used a discrete-time continuous-space model, an integrodifference equation, to study the effect of intrinsic stochastic factors.

Deriving equations for the first two spatial moments of the population density via moment closure approximations, he obtained an upper bound for the invasion speed via comparison methods. Neubert and Caswell (2000) combined integrodifference equations for dispersal with matrix models for population growth to analyze invasion velocity in populations with stage structure.

Though they did not incorporate stochastic effects in their description, they acknowledged the importance of including environmental stochasticity and demographic stochasticity in their discussion of directions for future research.

They also acknowledged the importance of including Allee effects [START_REF] Allee | The Social Life of Animals[END_REF].

The issue of environmental stochasticity has been addressed by Schreiber and Ryan, who studied the effect of a fluctuating environment in integrodifference matrix models (Schreiber and Ryan, 2011). Melbourne and Hastings have reported experimental studies to assess the effects of various stochastic factors on extinction risk [START_REF] Melbourne | Extinction risk depends strongly on factors contributing to stochasticity[END_REF] and the effects of demographic stochasticity on spread rates in populations for the red flour beetle (Tribolium castaneum) [START_REF] Melbourne | Highly variable spread rates in replicated biological invasions: fundamental limits to predictability[END_REF]).

Here we investigate the effect of environmental stochasticity on the invasion velocity of two types of population dynamics, namely logistic growth and growth with Allee effect. For logistic growth, we consider an external noise that gives rise to fluctuations in either the density-independent growth rate or the competition term. For growth with Allee effect, we consider an external noise that gives rise to fluctuations in the Allee threshold. This requires us to address the open problem of how to include stochasticity in an analytically tractable model as stated by different authors (see, e.g., Mollison, 1991;[START_REF] Melbourne | Extinction risk depends strongly on factors contributing to stochasticity[END_REF].

We adopt a continuous-time continuous-space description of the population, namely reaction-diffusion equations with multiplicative, spatially correlated, temporally white Gaussian noise. The effect of external noise on propagating fronts has been studied in physical and chemical reaction-diffusion systems, mostly for the case of bistable and excitable systems (see, e.g., [START_REF] Schimansky-Geier | Kink propagation induced by multiplicative noise[END_REF][START_REF] Armero | External Fluctuations in Front Propagation[END_REF][START_REF] García-Ojalvo | Noise in Spatially Extended Systems[END_REF][START_REF] Hizanidis | Noise-Induced Front Motion: Signature of a Global Bifurcation[END_REF]. Our approach consists in decomposing the population density into a systematic (or deterministic) part and a fluctuating part. We derive analytically the reaction-diffusion equation for the deterministic part and construct a variational principle that provides the invasion velocity and the front profile by maximizing a functional. This allows us to study analyti-cally the effect of environmental stochasticity on the invasion velocity and on the front profile.

We find for logistic growth that if fluctuations are present in the densityindependent growth rate, then the final state after invasion experiences an increase in the population density and the invasion velocity increases with the intensity of fluctuations. If fluctuations affect the competition term, then two types of invasion fronts can occur. The first type corresponds to the front observed in the deterministic system. It occurs for certain initial conditions.

The population density of the final stable state after invasion increases in the presence of fluctuations. The invasion front, however, travels with the Fisher velocity, that is, the invasion velocity does not depend on the intensity of fluctuations. For other initial conditions, and if the intensity of the fluctuations is low enough, a noise-induced front occurs that corresponds to the collapse of a noise-induced high-density state to the effective carrying capacity of the system. This front travels with a velocity that decreases with the intensity of the fluctuations. In the cases mentioned above, the front is pulled-its dynamics is governed by the tail of the front-in the absence of fluctuations and remains pulled in the presence of fluctuations. For growth with strong Allee effect, fluctuations in the Allee threshold give rise to a richer behavior.

In this case, the invasion front is pushed-its dynamics is governed by the interior part of the front-in the absence of fluctuations. In their presence, the invasion front can be pulled, pushed or even disappear all together, depending on the intensity of the fluctuations.

EQUATION FOR THE SYSTEMATIC POPULATION DEN-

SITY

Consider a population with density φ(x, t). For simplicity, we assume that the population lives in a one dimensional system and evolves dynamically according to a reaction-diffusion equation,

∂φ(x, t) ∂t = D ∂ 2 φ(x, t) ∂x 2 + F [φ(x, t), t] . (1)
We consider two types of population dynamics, i.e., two forms of F (φ), namely logistic population growth,

F (φ, t) = r 1 φ(x, t) -r 2 φ(x, t) 2 , (2)
and growth with Allee effect,

F (φ, t) = rφ (1 -φ) (φ -a) . (3)
For logistic growth, r 1 is the density-independent growth rate, or the coefficient of increase to use Skellam's terminology [START_REF] Skellam | Random Dispersal in Theoretical Populations[END_REF]. The coefficient r 2 is the limiting coefficient in Skellam's terminology; it describes intraspecific competition or the struggle for resources. For growth with Allee effect, r is the growth rate and a the Allee threshold. We consider only the case a > 0, i.e., a strong Allee effect. For logistic growth, one usually writes r 1 = r and r 2 = r/K, where K is the carrying capacity of the environment. In general, we can consider these parameters to be functions of space and time, e.g., r 1 (x, t) and r 2 (x, t), though they are usually treated as constants in the literature [START_REF] Murray | Mathematical Biology[END_REF]. Equation (1) does not include fluctuations explicitly. Environmental fluctuations, i.e., external noise, turns the parameters r, K and a into stochastic quantities and (1) into a SPDE, whose solution will be a functional of the noise [START_REF] San Miguel | Instabilities and Nonequilibrium Structures VI[END_REF],

∂φ ∂t = D ∂ 2 φ ∂x 2 + F (φ) + ε 1/2 g(φ)η(x, t), (4)
where η(x, t) represents the noise and ε measures the noise strength. Note that the external noise appears in a multiplicative way in (4) if g(φ) is not constant. Motivated by the Central Limit Theorem, we consider the case that the noise is Gaussian with zero mean, η(x, t) = 0. We further assume that the fluctuations occur on a faster time scale than that of the systematic population growth. This implies that we can consider the noise to be white in time, i.e., it has a correlation function given by

η(x, t)η(x , t ) = 2C(|x -x | / )δ(t -t ). (5)
Here C(x) is the spatial correlation function of the noise and its characteristic length. In the numerical simulations C(x) has been chosen to be Gaussian with variance 1 and C(0) = 10. Since the noise is white in time, the population density is a Markovian stochastic variable.

The mean value of the noise term g(φ)η(x, t) does not vanish if g(φ) is not constant, which gives rise to a systematic contribution to the dynamics of the density [START_REF] García-Ojalvo | Noise in Spatially Extended Systems[END_REF]. Adopting the Stratonovich interpretation of the SPDE (4), we evaluate g(φ)η(x, t) using Novikov's theorem [START_REF] Novikov | Functionals and the random-force method in turbulence theory (Euler velocity field described by random forces method, using Lagrangian representation of turbulence)[END_REF], see Appendix A for details,

g(φ)η(x, t) = ε 1/2 C(0) dg(φ) dφ g(φ) . ( 6
)
This result allows us to rewrite (4) in a form where the noise term has zero mean,

∂φ ∂t = D ∂ 2 φ ∂x 2 + F eff (φ) + ε 1/2 ζ(φ, x, t). ( 7
)
The effective reaction term and the new noise term are given by

F eff (φ) = F (φ) + εC(0) dg(φ) dφ g(φ), (8)
and

ζ(φ, x, t) = g(φ)η(x, t) -ε 1/2 C(0) dg(φ) dφ g(φ), (9)
respectively. We decompose the population density into a systematic, or deterministic, part and a fluctuating part by writing the density as

φ(x, t) = ρ(x, t) + ∞ n=1 ε n/2 ρ n (x, t). ( 10
)
The first term ρ(x, t) is deterministic and includes the systematic contributions. The sum represents the fluctuating part of the density. This expansion represents an extension of the standard small-noise expansion procedure (see, e.g., [START_REF] Marcus | A stochastic logistic diffusion equation[END_REF]; the form (10) keeps the relevant, systematic contributions of the multiplicative noise. It allows us to separate the systematic effects from the purely stochastic effects of the external noise on the population density.

We assume that the deterministic contribution dominates, while the stochastic part plays only a small role. We will show by comparing the analytical results with numerical simulations that focusing solely on the deterministic density represents an excellent approximation.

Substitution of ( 10) into (7) yields to the lowest order the equation for the systematic (deterministic) contribution,

∂ρ ∂t = D ∂ 2 ρ ∂x 2 + F eff (ρ). ( 11
)
The equation for the next order, the first stochastic contribution, is given by

∂ρ 1 ∂t = D ∂ 2 ρ 1 ∂x 2 + dF eff (ρ) dρ ρ 1 + ζ(ρ, x, t). ( 12
)
By construction, see ( 9), the mean value of the noise ζ(ρ, x, t) is zero, and its correlation function is given by

ζ(ρ, x, t)ζ(ρ, x , t ) = g(ρ(x, t))g(ρ(x , t )) η(x, t)η(x , t ) . ( 13
)
The equation for the systematic density, the dominant contribution, can be written as

∂ρ ∂t = D ∂ 2 ρ ∂x 2 + F (ρ) + ε 0 dg(ρ) dρ g(ρ), (14)
where ε 0 ≡ εC(0).

We investigate this equation explicitly for three cases: (i) fluctuations in the density-independent growth rate, (ii) fluctuations in the carrying capacity, and (iii) fluctuations in the Allee threshold. Stochastic variations in environmental variables, such as temperature and precipitation for example, lead to random changes in the vital rates of populations (Schreiber and Ryan, 2011).

Fluctuations in the density-independent growth rate [START_REF] Tuckwell | A study of some diffusion models of population growth[END_REF] occur because the difference between births and deaths changes with time and fluctuates around a mean value, due to external stochastic factors. Fluctuations in the competition term represent temporal variations in the availability of resources [START_REF] Feldman | A population's stationary distribution and chance of extinction in a stochastic environment with remarks on the theory of species packing[END_REF]. Finally, fluctuations in the Allee threshold occur because the minimum number of individuals necessary for the population to grow fluctuates in space and time due to environmental factors.

Fluctuating density-independent growth rate

If the density-independent growth rate fluctuates around a mean value r m , we write r 1 (x, t) = r m 1 + ε 1/2 η(x, t) and r 2 = r m /K, with r m and K constants.

Introducing the expression for the fluctuations of the density-independent growth rate into (1), we obtain the following SPDE:

∂φ ∂t = D ∂ 2 φ ∂x 2 + r m φ 1 - φ K + ε 1/2 r m φη(x, t). (15)
In this case g(φ) = r m φ, as can be seen by comparing ( 15) and (4). According to ( 14), the systematic density obeys the reaction-diffusion equation

∂ρ ∂t = D ∂ 2 ρ ∂x 2 + r m (1 + r m ε 0 )ρ - r m K ρ 2 . ( 16
)
The systematic effects of the fluctuating net growth rate are represented by the term r 2 m ε 0 ρ. In the following we assume K = 1 for simplicity. Note that fluctuations in the density-independent growth rate lead to a larger effective coefficient of increase.

Fluctuating competition term

To model fluctuations in the carrying capacity, we consider (2) with r 2 (x, t) = r Km 1 + ε 1/2 η(x, t) and r 1 (x, t) = r. Equation (1) turns into the SPDE

∂φ ∂t = D ∂ 2 φ ∂x 2 + rφ 1 - φ K m - r K m ε 1/2 φ 2 η(x, t). (17)
In this case g(φ) = -rφ 2 /K m . According to (14), we obtain the following equation for the systematic population density,

∂ρ ∂t = D ∂ 2 ρ ∂x 2 + rρ 1 - ρ K m + 2r 2 K 2 m ε 0 ρ 3 . ( 18
)
Note that noise in the competition term introduces a singular perturbation in the logistic growth model. The reaction function of the logistic model is quadratic, while the systematic contributions of the noise give rise to a cubic term in the effective reaction function. In other words, noise in the limiting term has the effect of introducing a mutualistic term due to its Gaussian nature. For densities above a critical value, see below, this cubic term will overcome the quadratic competition term and lead to unbounded growth. In actual populations, higher-order limiting term, e.g., terms of order φ 4 , will come into play at these high densities and prevent a population explosion.

Consequently, the range of validity of ( 18) is restricted to densities below the critical value and to sufficiently small noise strengths.

This does not affect our studies, since the expansion ( 10) is essentially a smallnoise expansion and self-consistency requires sufficiently small noise strengths.

In the following we assume K m = 1 for simplicity.

Fluctuating Allee threshold

If the Allee threshold fluctuates around the mean value a m , we write a(x, t) = a m 1 + ε 1/2 η(x, t) , r 1 (x, t) = r 2 (x, t) = r and K = 1 and obtain the SPDE

∂φ ∂t = D ∂ 2 φ ∂x 2 + rφ (1 -φ) (φ -a m ) -ra m ε 1/2 φ (1 -φ) η(x, t). (19)
Here g(φ) = -ra m φ (1φ), and according to ( 14) we find

∂ρ ∂t = D ∂ 2 ρ ∂x 2 + rρ (1 -ρ) 1 -2rε 0 a 2 m ρ -a m (1 -rε 0 a m ) . ( 20
)
Note that the third factor of the effective reaction term is very similar to the case without noise (ε 0 = 0) but with different coefficients, which can be either positive or negative giving rise a wide range of different front solutions. In other words, external noise can decrease the effective Allee threshold and even render it negative, i.e., turn a strong Allee effect into a weak Allee effect.

INVASION FRONTS

Consider the situation that the reaction term F (φ) in (1) has two steady states, F (φ) = 0, given by φ u = 0 (unstable) and φ s = 1 (stable). The velocity of the invasion front joining states 0 and 1 can be obtained from the variational formula,

v 2 = 2D max 0<λ≤2 λΓ(4) Γ(2 -λ)Γ(2 + λ) 1 0 F (ρ)ρ -λ (1 -ρ) λ dρ , (21)
where λ parameterizes the set of trial functions {ρ -λ (1ρ) λ }. See Appendix B for details of the derivation of the variational formula (21).

Invasion velocities

In the following, we determine the invasion velocity for the three cases of external noise described in Sect. 2. To do so, we apply the variational formula (21)

to the evolution equation for the systematic part of the density, namely ( 16), ( 18) and ( 20), replacing F (φ) by F eff (ρ) and renormalizing the density where necessary, so that the steady states are given by 0 and 1. We also compare our analytical results with numerical calculations. Details of the numerical algorithm to calculate the invasion velocities and the front profile are provided in Appendix C. In the following sections we find fronts invading both unstable and metastable states. The main difference is that in the first case the front velocity depends on the tail of the initial condition and in the second case the front velocity is independent. For more information about pulled and pushed fronts and fronts propagating into unstable and metastable states see the works by van Saarloos, 2003 and by [START_REF] Méndez | Reaction-Transport Systems[END_REF].

Fluctuating density-independent growth rate

If the density-independent growth rate fluctuates, the effective reaction term is given by F eff (ρ) = r m (1 + r m ε 0 )ρr m ρ 2 according to ( 16). The effective reaction term has two steady states, namely ρ = 0 and ρ = 1 + r m ε 0 . The first steady state is unstable to small density perturbations, while the second one is stable. To be able to apply (21) we introduce the renormalized density u = ρ/(1 + r m ε 0 ). This change of variables turns ( 16) into the Fisher equation with a rate constant given by r m (1 + r m ε 0 ). Consequently, the front velocity is given by

v = 2 Dr m (r m ε 0 + 1). ( 22
)
We confirm this result by numerical integration of ( 1) and (2), see Figure 1.

The velocity grows sublinearly with the stochastic noise strength as ε 1/2recall that ε 0 = εC(0). Moreover, the velocity becomes the Fisher front velocity, v F , for ε → 0, as expected.

Fluctuating competition term

If the carrying capacity fluctuates, the effective reaction term is given by F eff (ρ) = rρrρ 2 + 2r 2 ε 0 ρ 3 according to (18). This reaction term possesses three different steady states: ρ = 0 and

ρ ± = 1 ± √ 1 -8ε 0 r 4ε 0 r . ( 23
)
The nontrivial solutions, ρ ± , exist only if ε 0 < 1/8r. This is a reflection of the fact that the validity of ( 18) is restricted to the regime of small noise strength, as discussed in Sect. 2.2. The state ρ -is stable, and the states 0 and ρ + are unstable. The latter is the critical value of the population density, see Sect. 2.2, above which the population will grow without bound. Writing the square root as a power series, we find that

ρ -= 1 + ε 0 r/4 + O(ε 2 0 ), (24)
ρ + = 1 2ε 0 r -1 + O(ε 0 ). ( 25
)
The state ρ -correspond to the deterministic case of the population being at the carrying capacity of the system, φ = K m = 1. Noise in the competition term increases the effective carrying capacity of the system. The other state, ρ + , is a pure noise-induced state. Its origin lies in the fact that noise in the competition term, i.e., in the limiting coefficient, gives rise to a systematic effect of the opposite nature, see Sect. 2.2.

We exclude here situations where the initial densities are larger than ρ + , i.e., we exclude all cases of unbounded population growth. Then there are two possible pulled fronts propagating into unstable states, one invading the state 0 and connecting it to the state ρ -and another invading the state ρ + and connecting it to ρ -. To obtain both fronts, one has to prepare the system with appropriate initial conditions. Below we indicate which initial condition we employed. Equation ( 21) allows us to calculate the front velocities of the two different fronts, and we compare the results with numerical results.

3.1.2.1

Pulled front connecting 0 to ρ -. This case corresponds to the one that occurs in the deterministic system, (1) with (2). Since the noise increases the effective carrying capacity, we must again renormalize the density, u = ρ/ρ -, to be able to apply the variational formula (21). The front velocity is calculated from

v 2 = 4Dr 2 ε 0 ρ + ρ -max 0<λ≤2 λ 1 - 1 4 ρ - ρ + (2 -λ) . ( 26
)
Since the function to be maximized is monotonically increasing with λ, the maximum is attained at λ = 2, which results in the Fisher velocity, v F = 2 √ Dr. Note that in this case the front velocity does not depend on the noise strength ε 0 . This result is depicted in Figure 2a, where we compare our analytical result with results from numerical simulations. For the numerical calculations we chose the initial conditions to be ρ(x,

0) = ρ -[1 -θ(x -x 0 )],
where θ(•) is the Heaviside function and x 0 > 0.

3.1.2.2 Pulled front connecting ρ + to ρ -. Such a front can occur only in the presence of external noise; it has no deterministic counterpart. To obtain the front velocity using ( 21), we apply the change of variables u = (ρ +ρ)/(ρ +ρ -) and find

v 2 = 4Dr 2 ε 0 (ρ + -ρ -) ρ + max 0<λ≤2 λ 1 - 1 4 ρ + -ρ - ρ + (2 -λ) . (27)
Again the function of λ is monotonically increasing, the maximum is at λ = 2, and the front velocity expression reads,

v = D ε 0 1 -8ε 0 r + √ 1 -8ε 0 r. (28)
This front velocity decreases with the noise amplitude as ε -1/2 and is valid only if ε 0 < 1/8r. This result is depicted in Figure 2b, where we compare our analytical result with results from numerical simulations. For the numerical calculations we chose the initial conditions to be ρ(x, 0) = ρ -+(ρ + -ρ -)θ(xx 0 ) with x 0 = 10. The front corresponds to a collapse of the hight-density state ρ + to the "deterministic" state ρ -. The fact that the velocity decreases with the noise amplitude reflects the nature of the state ρ + . It is a noise-induced or noise-sustained state, and its survival time increases with increasing noise strength.

Fluctuating Allee threshold

Applying ( 21) to ( 20) we obtain

v 2 = 2Dr m max 0<λ≤2 1 -2rε 0 a 2 m 4 λ(2 -λ) + a m (rε 0 a m -1)λ . ( 29
)
The maximum must be determined carefully because it depends strongly on the values for the parameters a m and ε 0 . We can distinguish up to five different and connecting ρ = 0 to ρ = 1. In this case, the front can be pulled (the front dynamics is governed by the leading edge, i.e., the nonlinear term is not important) or pushed (i.e. the front dynamics is governed by its interior part and nonlinear term becomes more important than linear term) and can occur only for a m < 1/2. Using the variational formula (29) we can determine analytically where this transition takes place. Equation ( 29) implies that

v I = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ rD 2 1-2am √ 1-2rε 0 a 2 m , 1 am < ε 0 r ≤ 1+2am 4a 2 m (pushed), 2 Dra m (a m rε 0 -1), 1+2am 4a 2 m ≤ ε 0 r < 1 2a 2 m (pulled). (30)
In Figure 4 we plot the invasion velocity versus rε 0 (30) (solid line) and compare it with results from numerical solutions (circles). The transition between the pulled and pushed regime is located at rε 0 = (1 + 2a m )/(4a 2 m ).

Region II. Front propagating into the unstable state ρ = 0

and connecting ρ = 0 to ρ * = a m (rε 0 a m -1)/(2rε 0 a 2 m -1). Note the ρ * is the effective Allee threshold in the presence of external noise. We again need to renormalize the density, u = ρ/ρ * . In this case, the front is always pulled and the region is defined by

ε 0 r > 1/2a 2 m if a m < 1/2 or by ε 0 r > 1/a m if a m > 1/2.
According to (29), the velocity of the front is given by v II = 2 Dra m (a m rε 0 -1). ( 31)

Region III. Front propagating into the metastable state

ρ = 0 and connecting ρ = 0 to ρ = 1. This region is defined by ε 0 r < 1/a m for a m < 1/2. In this case, the front is always pushed and it travels with the velocity

v III = rD 2 1 -2a m 1 -2rε 0 a 2 m . ( 32
)
If ε 0 r > 1/a m and a m < 1/2 there is no front. In Figure 3 we depict these regions and compare our theoretical results with numerical simulations.

To construct numerically the region boundaries we fix a m and then vary rε 0 .

For each value of rε 0 we compute numerically the front profile and check in which region we are. For example, to find the boundary between regions I and II in Figure 3 we begin with a high value of rε 0 (belonging to region I) and determine if the front velocity depends on the shape of initial condition. To this end we take ρ(x, 0) = 1θ(x)e μx and observe that the front velocity depends explicitly on μ. This indicates that the front invades the unstable state ρ = 0 and we are in region I. When we decrease progressively rε 0 and it crosses the boundary to enter in region III, the front velocity ceases to depend on μ which indicates that the front invades the metastable state ρ = 0 and we are in region III. To compute the boundary between regions I and II we must observe which is the final stable state, since in both regions the front invades the unstable state ρ = 0. To this end, we take the initial condition

ρ(x, 0) = 0.1[1 -θ(x)].
If the stable state (state behind the front) is ρ = 1 then we are in region I. Otherwise, we are in region II. Finally, the boundary between region II is trivial to find numerically because below the boundary no front is formed.

Front profiles

Besides the front velocity, front profiles are also of interest. The variational principle allows us also to predict the profile shape, i.e., the analytical expression for ρ(x, t), when the front is fully developed and the velocity has reached its stationary value. If the trial function is chosen to be g

(ρ) = ρ -λ (1 -ρ) λ ,
and using the fact that (see Appendix B)

- d dρ ln g = v Dp , (33)
we find that

p = v λD ρ(1 -ρ) ≡ - dρ dz . (34)
Integrating over z, we obtain the expression for the density

ρ(z) = 1 1 + b -1 exp -vz λD , (35)
which is a sigmoidal function. The integration constant is b

-1 = (1 -ρ * )/ρ * ,
where ρ * is the density at z = 0. The asymptotic front profile given by ( 35) is the same for each of the three cases we have studied. To specify the solution for each case, we must substitute the value of λ for which the maximum in ( 21)

is attained and the corresponding velocity v. For example, when the densityindependent growth rate fluctuates, the velocity is given by ( 22). The density profile reads

ρ(x, t) = 1 1 + b -1 exp -rm(1+rmε 0 ) D x -2 Dr m (r m ε 0 + 1)t , (36)
which is plotted in Figure 5 and compared with numerical integration of the dynamic equations. The other cases are straightforward.

CONCLUSIONS

Variations in temperature and other climatic factors can produce environmental fluctuations that can modify the invasion patterns of populations. The question of how such external stochastic factors affect population dynamics will become more crucial with the expected increase in climatic variability (Schreiber and Ryan, 2011). We have obtained the density profile and the invasion velocity in the presence of such external fluctuations for reaction-diffusion models of population growth. We have considered the situations where they give rise to noise in the density-independent growth rate, the competition term or the Allee threshold. We have been able to address this problem in an analytically tractable model by reducing the SPDE to an effective deterministic equation for the mean density that incorporates the systematic effect of the noise. We have constructed a variational principle that is able to provide the mean density profile and the invasion velocity for both pulled and pushed fronts.

If fluctuations affect the density-independent growth rate, the mean population density of the invaded area and the invasion velocity experience an increase with the noise intensity. Both increases have exactly the same cause, namely an increase in the effective density-independent growth rate due to the systematic contribution of the multiplicative external noise. This result differs from previous studies that either considered intrinsic stochastic factors or adopted discrete-time integrodifference descriptions of population growth.

Lewis analyzed the effect of internal fluctuations in a discrete-time model [START_REF] Lewis | Spread rate for a nonlinear stochastic invasion[END_REF]. He found, in agreement with a long-time conjecture in the literature, that "intrinsic stochastic factors associated with interacting individuals can slow the spread of a population or disease, even in a uniform

(φ)η(x, t) = 2ε 1/2 dg(φ) dφ g(φ) dx C (|x -x | /λ) δ(x -x) = ε 1/2 C(0) dg(φ) dφ g(φ) , (40)
where in the last equality we have made use of the the Stratonovich interpretation of stochastic calculus ( dx δ(xx) = 1/2). Equation ( 40) is Novikov's theorem for spatially extended systems [START_REF] Novikov | Functionals and the random-force method in turbulence theory (Euler velocity field described by random forces method, using Lagrangian representation of turbulence)[END_REF].

Appendix B. Variational principle

In order to formulate a variational principle we start with the reaction-diffusion

equation ∂ρ(x, t) ∂t = D ∂ 2 ρ(x, t) ∂x 2 + F (ρ). (41)
The time evolution of an initial condition ρ(x, 0) has been studied both for F (ρ) monostable or bistable. It was proven [START_REF] Aronson | Multidimensional nonlinear diffusion arising in population genetics[END_REF] that suitable initial conditions evolve into a monotonic traveling front ρ(x-vt),

joining the stable state ρ = 1 to ρ = 0. If F is monostable, e.g., logistic growth, there exists a continuum of values of v for which a monotonic (pulled) front occurs. If the initial condition has compact support, the front with minimal velocity is selected. If F is bistable, e.g., growth with Allee effect, there is a single isolated value of the velocity for which the (pushed) front exists. In this appendix we show how to derive the velocity for pushed fronts and the minimal velocity in the case of pulled fronts from a variational principle following the method developed by [START_REF] Benguria | Speed of fronts of the reaction-diffusion equation[END_REF] To this end, consider (41) transformed to the frame commoving with the front

D d 2 ρ dz 2 + v dρ dz + F (ρ) = 0, ( 42
)
where v is the asymptotic front velocity and z = xvt is the commoving coordinate. For simplicity, we assume that the front connects the state ρ = 0 (unstable) to the state ρ = 1 (stable). One can define the front slope as p(ρ) = -dρ/dz, which is a positive function for every value of the density.

Therefore, (42) can be written as,

Dp dp dρ -vp + F (ρ) = 0, (43)
with the conditions p(0) = p(1) = 0. Defining g(ρ) as a positive trial function on (0, 1), such that h = -dg/dρ > 0, multiplying (43) by g(ρ), and integrating by parts we obtain

1 0 F gdρ = 1 0 vpg - 1 2 Dhp 2 dρ. (44)
For fixed ρ, the functional

Φ[p] ≡ vpg - 1 2 Dhp 2 (45)
has a maximum at p max = vg/(Dh), since v, p, g and h are positive. Therefore Φ[p] ≤ Φ[p max ] = v 2 g 2 /(2Dh) for any value of ρ. Equations ( 44) and (45) lead to the following condition, which is valid for any reaction function F (ρ):

v 2 ≥ 2D 1 0 F gdρ 1 0 (g 2 /h)dρ . (46)
The equality in (46) holds if

- d dρ (ln g * ) = v Dp , (47)
which can be integrated to yield

g * (ρ) = exp ρ 0 ρ v Dp dρ ( 48
)
where 0 < ρ 0 < 1 comes from the integration constant. Obviously, g * (ρ) is a continuous, positive and decreasing function on (0, 1) with g * (1) = 0. Near ρ = 0, g * diverges. We linearize (43) to find from (48) that g * diverges as

g * (ρ) ∼ ρ -v/(Dμ + ) (49) near ρ = 0, with μ + = v + v 2 -4DF (0) /2D. Consequently, F g ∼ ρg ∼ ρ 1-v/(Dμ + ) near ρ = 0. Analogously g 2 /h ∼ ρ 1-v/(Dμ + ) near ρ = 0, so the integrals in (46) exist if v/(Dμ + ) < 2.
Since this condition is always satisfied, we can ensure that the equality in ( 46) is well established. Thus

v 2 = max g 2D 1 0 F gdρ 1 0 (g 2 /h)dρ , (50)
where the maximum is taken over all positive, decreasing trial functions on (0,1). Since (49) implies that g must behave like ρ -λ , with λ > 0, near ρ = 0

and must be such that g(1) = 0, we can choose the trial function g(ρ) = ρ -λ (1ρ) λ . With this choice, (50) has to be maximized over the parameter λ, but only for values such that the integrals in the variational principle exist,

v 2 = 2D max 0<λ≤2 λΓ(4) Γ(2 -λ)Γ(2 + λ) 1 0 F (ρ)ρ -λ (1 -ρ) λ dρ . ( 51
)
To confirm that the variational formula provides the correct front velocities, we apply it to a simple case. If the reaction term is logistic, F (ρ) = rρ(1ρ),

  situations, which correspond to different regions of the parameter space (rε 0 versus a m ), see Figure 3. (The first two situations are shown as one region, Region I, in the figure.) 3.1.3.1 Region I. Front propagating into the unstable state ρ = 0

  environment." Schreiber and Ryan employed an integrodifference matrix population model to investigate the effect of external stochastic factors and also found that increased variability in vital rates leads to slower invasion speeds(Schreiber and Ryan, 2011). Clearly, the nature of the fluctuations, intrinsic versus external, plays an important role for the invasion speed. So does the nature of the dynamics, discrete-time integrodifference equations versus continuous-time reaction-diffusion equations.If fluctuations affect the competition term, two different fronts travelling with different velocities appear depending on the initial conditions. One corresponds to the deterministic invasion front. (The initial conditions for this type of front correspond to the usual conditions of a non-native species being introduced locally at small to moderate densities.) Here the noise raises the carrying capacity of the system, while having no effect on the invasion speed. The other front represents a noise-induced phenomenon, namely the "collapse" of a noise-sustained high-density state. (The initial conditions for this type of front correspond to a high-density state that arose due to some catastrophic event triggering a sudden explosive growth in the population.) The front speed is strongly dependent on the noise strength and decreases with increasing noise amplitude. In other words, stronger noise extends the lifetime of the noise-induced high density state.In the case of logistic growth, the systematic effect of the external noise in either the coefficient of increase or the limiting coefficient does not change the pulled character of the front. If the population experiences the Allee effect, the invasion front is pushed in the absence of noise. In this case, the invasion velocity cannot be calculated by the linear conjecture, because the front dynamics is governed by the nonlinear terms. However, in the presence of fluctuations in the Allee threshold, both pulled and pushed fronts are possible, and we have been able to predict analytically the transition diagram. The dynamics of the system is very rich and the velocity and success or failure of the invasion depends on the balance between the average Allee threshold and the noise strength. Our results agree with numerical simulations and can be generalized to two dimensional invasions and the case of colored noises.Our results for the logistic growth model show that fluctuations in the densityindependent growth rate are more "dangerous" than fluctuations in the competition term. The former increase the velocity of invasion, while the latter have no effect. This suggests that if the effects of external variability on the vital rates of a non-native species can be mitigated, then the effort should be directed towards minimizing the noise strength of the density-independent growth rate. It should be possible to test our predictions in experiments of the type conducted by Melbourne and Hastings. The authors were interested in the effects of endogenous variability on invasion velocity in those experiments and tightly controlled the ambient environment[START_REF] Melbourne | Highly variable spread rates in replicated biological invasions: fundamental limits to predictability[END_REF]. It should be possible to repeat the experiments and introduce deliberate and controlled variability of the environment. A comparison of the two types of experiments should reveal the effects of external fluctuations on the invasion g

  without noise, the above equation leads to the Fisher front velocity,v ≡ v F = 2 √ Dr. (52)van Saarloos, W., 2003. Front propagation into unstable states. Physics Reports 386, 29-222. Williamson, M., 1996. Biological Invasions. Chapman and Hall, London.

Figure captions Figure 1 .

 captions1 Figure captionsFigure 1. Invasion velocity as a function of the noise amplitude if r 1 fluctuates; C(0) = 10, r m = 0.1, D = 1. Solid lines correspond to the analytical expression (22) and circles to numerical solutions. Inset: Logarithmic plot to show that the dependence of the velocity on the noise amplitude is square root like (slope 1/2).

Figure 2 .

 2 Figure 2. Invasion velocity as a function of the noise amplitude if K fluctuates, with K m = 1. Solid lines correspond to the analytical expression (28) and circles to numerical solutions. a) Invasion velocity for a front that connects ρ -with ρ = 0. b) Invasion velocity for a front that connects ρ -with ρ + . C(0) = 10, r = 0.1, D = 1. In both cases x 0 = 10.

Figure 3 .

 3 Figure 3. Diagram of the different invasion front regimes if the Allee threshold fluctuates. Solid lines correspond to the analytical expressions and circles to numerical solutions.

Figure 4 .

 4 Figure 4. Invasion velocity for a front connecting ρ = 0 to ρ = 1 (region I). The boundary between the regions of pulled and pushed fronts occurs at rε 0 = 6; a m = 1/4.The initial condition taken is ρ(x, 0) = 1θ(x).

Figure 5 .

 5 Figure 5. Front profiles at different times for a fluctuating intrinsic net growth rate; C(0) = 10, r m = 0.1, D = 1 and L = 5000. The final stable state is ρ 0 = 1 + r m ε 0 = 2. Analytical and numerical results are plotted. The analytical profile calculated from (36) is drawn for t = 500 and ρ * = 4/3.

Acknowledgments

The work has been funded by Direcció General de Recerca-Generalitat de Catalunya under grant 2009 SGR-164 (VM and DC) and by MICINN Grant No. FIS 2006-12296-C02-01 (VM).

velocity.

Appendix A. Mean value of the noise term

The mean value can be evaluated as follows: g(φ)η(x, t) = t 0 dt dx η(x, t)η(x , t ) δg(φ(x, t)) δη(x , t

Equation ( 4) can be integrated formally to obtain

and (38) implies that

Substituting ( 39) into (37), we obtain

Appendix C. Numerical algorithm

The numerical integration of ( 15), ( 17) and ( 19) was performed using the general approach for SPDE [START_REF] San Miguel | Instabilities and Nonequilibrium Structures VI[END_REF]. The habitat is simulated as a grid consisting of N lattice sites of size Δx, such that L = NΔx is the habitat size. Every lattice site i is identified by its position x i , which is a discrete variable. The continuum density, φ(x, t), is replaced by its analogue in discrete space, φ(x i , t), and the external noise is simulated via the generation of Gaussian numbers, where we replace the white noise η(x, t) by 1 √ Δx η i (t), as a result of the relation between the Dirac and the Kronecker delta functions.

The correlation function is given by η i (t)η j (t ) = δ ij δ(tt ). Moreover, the noise amplitude is εC(0) = ε/Δx [START_REF] Santos | Noise-induced fronts[END_REF].

Equations ( 15), ( 17) and ( 7) can be written as,

where G(φ) = ε 1/2 g(φ). Assuming uniform time steps of size Δt and integrating ( 15) and ( 17) using the Euler algorithm, we obtain a dynamic equation that is discrete in both space and time,

where Δt 1/2 u t,i = η i (t). Simulations have been performed with Δx = 0.1, Δt = 0.0025 and D = 1. We have varied the parameters r, K, ε and L in different simulations.