text:
# Long range correlations in DNA sequences

## I INTRODUCTION

Recently, there has been considerable interest in the finding of long range correlations in genomic DNA sequences . A DNA sequence is a chain of sites, each occupied by either a purine (Adenine and Guanine) or a pyrimidine (Cytocyine and Thymine) group. For mathematical modeling, the DNA sequence might be considered as a string of symbols (G, A, T and C) whose correlation structure can be characterized completely by all possible base-base correlation functions or their corresponding power spectra. Different techniques including mutual information functions and power spectra analyses , auto correlation , DNA walk representation , wavelet analysis  and Zipf analysis  were used for statistical analyses of DNA sequences. But despite the effort spent, it is still an open question whether the long range correlation properties are different for protein coding (exonic) and non coding (intronic, intergenemic) sequences . One more fundamental ground, there is still continuing debate as to whether the reported long range correlations really mean a lack of independence at long distances or simply reflect the patchiness (bias in nucleotide composition) of DNA sequences. There have been attempts to eliminate local patchiness using methods such as min-max , detrended fluctuation analysis (DFA)  and wavelet analysis . In spite of its success in modeling the long range correlations observed in DNA sequences, as indicated by the power law increase in the variance and the inverse power law spectrum , the problem of the correct statistical interpretation of DNA walk is still unresolved and is attracting the attention of an increasing number of investigators. Since approaches based on different models predict different correlation structure, there is no unique measure of the degree of correlation in DNA sequences. Therefore, it is very important to investigate the correlations and extract the power law exponent $`\alpha `$ rather in a model independent way so that the interpretation of the data including the theoretical analysis becomes more meaningful. There is another confusion related to this study is the absence of a clear definition of the term ”long range”. Clearly, what is considered to be long is relative to what is considered to be short. To over come some of these problems, recently we have suggested a new method  to measure the degree of correlations using the variance analysis of the density distribution of a single or a group of nucleotides. We have also suggested a way to find out an approximate length scale above which all DNA sequences show strong long range correlations irrespective of their intron contents while below this, the correlation is relatively weak. Further, the density distribution which is nearly Gaussian at short distances shows significant deviations from the Gaussian statistics at large distances. In this paper, we present the details of the analyses and also extract the correlation parameter $`\alpha `$ for several intron containing and intronless sequences.

## II Density distribution and Factorial moments:

In the present method, we build the frequency spectrum of a single or a group of nucleotides by dividing the DNA sequence into many equal intervals of length $`l`$. For example, to build a purine spectrum, we compute

$$n=\underset{i=l_0}{\overset{l_0+l}{}}u_i$$
(1)
where $`u_i`$=1 if the site is occupied by a G or A and $`u_i`$=0 otherwise. Ideally, one can divide the entire DNA sequence of length $`L`$ into $`m`$ equal intervals of size $`l`$ $`(l=L/m)`$. The purine or GA spectrum can be built by computing $`n`$ from all the intervals. Alternatively, $`n`$ can be computed in any segment between $`l_0`$ and $`l_0+l`$ and the spectrum ($`n`$ distribution or $`P_n`$) is built by varying the starting position $`l_0`$ from 1, 2, 3 etc upto $`Ll`$ so as to cover the whole sequence ***At short distances, $`n`$ can be zero due to the non occurence of a given nucleotide. In such cases, the density spectrum can be built either including or excluding zeroth channel. In this analysis, we include zeroth channel also so that the complementarity is satisfied which is unlike the case when the zeroth channel is excluded. See appendix B for details. We adopt this second procedure for better statistics. Finally, the standard deviation (SD) of this $`P_n`$ distribution can be obtained from $`\sigma ^2=<n^2n_{0}^{}{}_{}{}^{2}>`$ which in general will depend on the interval or the window size $`l`$.

In addition to the standard deviation $`\sigma ^2`$, we also compute the factorial moments $`F_q`$’s of $`P_n`$. The normalized factorial moments of order q are written as

$$F_q=\frac{f_q}{f_1^q}$$
(2)
where

$$f_q=\underset{n=q}{\overset{\mathrm{}}{}}P_nn(n1)\mathrm{}..(nq+1)=\underset{n=q}{\overset{\mathrm{}}{}}\frac{n!}{(nq)!}P_n$$
(3)
As will be shown later, the factorial moment has the distinct advantage over the normal moments in identifying the genomic sequence from the random one. It may be mentioned here that for random Poisson distribution, the factorial moments for all q’s become unity i.e. for

$$P_n=\frac{a^ne^a}{n!}$$
(4)
the above factor for $`f_q`$ becomes

$$f_q=\underset{n=q}{\overset{\mathrm{}}{}}\frac{n!}{(nq)!}\frac{a^ne^a}{n!}=\underset{n=q}{\overset{\mathrm{}}{}}\frac{a^ne^a}{(nq)!}=\underset{m=0}{\overset{\mathrm{}}{}}\frac{a^{m+q}e^a}{m!}=a^q\underset{m=0}{\overset{\mathrm{}}{}}\frac{a^me^a}{m!}=a^q$$
(5)
which gives $`F_q`$=1.

In this work, we have applied the above factorial moment analysis (generally used to study the fluctuations during a phase transition ) to study the dynamical fluctuations present in the DNA sequences.

## III Principle of complimentarity

A general property noticed for all the genomic sequences (of statistically significant length) with a few exceptions is that the distributions of any single or group of nucleotides which has a probability of occurrence $`p`$ has the same variance $`\sigma `$ as that of its complimentary group that has the probability of occurrence $`(1p)`$, although both have different distribution functions. This would imply that even a single nucleotide distribution say $`G`$ distribution will have same variance as that of $`ATC`$ distribution or a $`GA`$ distribution will have identical variance as that of $`TC`$ distribution. Figure 1 shows $`\sigma `$ versus $`l`$ plots for $`G`$ and $`GA`$ distributions (solid curves) for two typical sequences of $`DROMHC`$ (Drosphilia Melanogaster, MHC, 22663 bps, $`20.5\%`$ $`G`$, $`30.3\%`$ $`A`$, $`25.4\%`$ $`T`$, $`23.8\%`$ $`C`$) and $`SC\mathrm{\_}MIT`$ (yeast mitochondrial DNA, $`9.1\%`$ $`G`$, $`42.2\%`$ $`A`$, $`40.7\%`$ $`T`$, $`8.0\%`$ $`C`$). As can be seen from the figure, the $`G`$ and $`GA`$ distributions have same $`\sigma `$ at all scale as that of $`ATC`$ and $`TC`$ distributions (filled circles) although the distribution functions of the two complimentary groups need not be identical. The above agreement is exact for most of the DNA sequences (with a few exceptions) as well as for the random sequences. For example, the $`\sigma `$ for $`G`$ and $`ATC`$ distributions of $`SC\mathrm{\_}MIT`$ and $`E.Coli:TN10`$ ($`E.Coli`$ with a $`TN10`$ mobile transposion (9147 bps) at location 22000 bps) show $`2\%`$ to $`3\%`$ deviations at all scale depending on the total length of the sequences where as for other DNA as well as random sequences, this agreement is exact. (This difference is not visible from figure 1 in case of $`SC\mathrm{\_}MIT`$ as the deviation is insignificant over a large distance).

Within the present formalism, we can also reproduce the result of random walk $`(RW)`$ model (See appendix for more detail) by assigning $`u_i=1`$ for purine group ($`G`$ and $`A`$) and $`u_i=1`$ for pyrimidine group ($`T`$ and $`C`$). However, unlike the random walk model of interpreting $`+1`$ and $`1`$ as the probability of step up and step down, $`P_n`$ can be considered as the frequency distribution of $`n`$ which gives the excess or deficit of purines over pyrimidines. The $`\sigma `$ versus $`l`$ as obtained from this assignment has also been shown in figure 1 (see the dotted curves labeled $`RW`$) for comparison. It is interesting to note that the $`RW`$ curves shows a parallel shift with respect to the $`GA`$ or $`TC`$ curves indicating that $`GA`$ or $`TC`$ distributions and $`RW`$ model have similar fluctuations at all scale. This is an interesting observations, as we can now use $`GA`$ or $`TC`$ distributions as alternatives to the DNA walk representation to study the correlation. The advantage is, since $`n`$ represents a sum, unlike the DNA walk model, the entire spectrum lies to the positive side of the coordinates which is essential to compute various higher moments like $`F_q`$ of the distributions.

It is also important to note that although the complimentary distributions have same $`\sigma `$ at all scale, the distribution functions need not be exactly identical. Figure 2 shows a typical normalized density distribution functions $`P_n`$ of two complimentary distributions $`G`$ and $`ATC`$ for the above two sequences ($`SC\mathrm{\_}MIT`$ and $`DROMHC`$) as a function of $`nn_0`$ (where $`n_0`$ is the average count ) at a typical length scale of $`l=150`$ (figures in left). The figures to the right shows $`P_n`$ distributions ($`x`$-axis is shifted by 100 for clarity) corresponding to the two purely random sequences having same length and nucleotide contents as that of $`DROMHC`$ and $`SC\mathrm{\_}MIT`$ sequences. It is interesting to note that although $`\sigma `$ versus $`l`$ plots are (nearly) identical $`i.e.`$, both distributions have same fluctuations at all scales, the distribution functions are not identical. This is an important characteristic of a DNA sequence which is not found in case of a random one.

## IV Extraction of slope parameter

The long range correlations are generally studied from the relation $`\sigma l^\alpha `$ where the parameter $`\alpha `$ is extracted from the $`\sigma `$ versus $`l`$ plot in the log-log scale. For the case of a completely random sequence, $`\alpha 0.5`$. The deviation of $`\alpha `$ from $`0.5`$ indicates presence of long range correlations. We have estimated $`\sigma `$ of $`G`$, $`A`$, $`T`$, $`C`$ and $`GA`$ distributions for several DNA sequences and found that $`\sigma `$ versus $`l`$ plot in the log-log scale is not linear over the entire length We consider only the $`G`$, $`A`$, $`T`$ and $`C`$ distributions to extract the correlation parameters for the individual nucelotides and $`GA`$ distributions to simulate the results of random walk model. Figure 3 shows $`\sigma `$ versus $`l`$ plot (bottom panel) for a typical $`E.Coli`$ sequence of length $`L=1.2`$ Mbps (solid curves) and $`L=30`$ Kbps (dotted curves) respectively. The top panel shows the factorial distributions of $`q`$=2, 3, 4 and 6 for a typical $`A`$ distributions, although similar plots can be obtained for other nucleotide distributions as well. A general feature of the factorial moments of the DNA sequence with a few exception is that at short distances, $`F_q<1.0`$ for all $`q^{}s`$ and exceeds unity at some point say at $`l_q`$. This behavior is not found in case of a purely random sequence where $`F_q`$ is always $`1.0`$. Further, all $`q`$’s do not cross unity exactly at the same point, $`l_q`$ being more for higher $`q`$ values. However, this variation is insignificant over a very large scale if we restrict to some of the lower moments say up to $`q=6`$.

From these plots and also from the several other studies, we make following few observations; (i) The $`\sigma `$ versus $`l`$ plot is not linear through out, rather starts bending around some region (say $`l_c`$, which could be different for different distributions) indicating a change of slope from $`\alpha _1`$ to $`\alpha _2`$, (ii) For most of the cases, while $`\alpha _1`$ shows weak deviation from $`0.5`$, $`\alpha _2`$ deviates significantly from $`0.5`$ and also depends on the sequence length $`L`$, (iii) The individual nucleotide distributions may have stronger correlations than any sum like $`GA`$ and $`TC`$ distributions or any other combinations.

Since $`\sigma `$ versus $`l`$ in the log-log plot starts bending around $`l_c`$, we can extract the slope by dividing the entire length into two segments; one for $`l<l_c`$ and the other one for $`l>l_c`$. This can be done by examining each case individually. However, we have noticed an approximate correlation between this bending region in $`\sigma `$ versus $`l`$ plot and the cross over points $`l_q`$ of the corresponding factorial moments i.e. the slope changes around the same region where the factorial moments become unity. This is a pure phenomenological observation which is found for several DNA sequences as listed in tables with a few exceptions which we will discuss below. It may be mentioned here that although, the two complimentary distributions have same fluctuations, both need not have identical factorial moments. Figure 4 shows the plots of $`F_q`$ versus $`l`$ for $`A`$ and $`GTC`$ distribution for a $`LAMCG`$ sequence. Since both are complimentary, they have identical fluctuations at all scales (hence same bending region), but the cross over regions in $`F_q`$ plots are different, being higher for $`ATC`$ distributions (due to large average values $`n_0`$ at all scales). While the $`l_q`$ value of the $`A`$ distribution shows an approximate correlation with the bending region of $`\sigma `$ versus $`l`$ plot where a possible slope change occurs, the $`l_q`$ values of $`GTC`$ distribution has no such correlations. This is true for any complementary distributions of $`G`$, $`A`$, $`T`$ and $`C`$ except for $`GA`$ and $`TC`$ distributions since both have nearly same overlapping cross over regions.

Therefore, only the $`l_q`$ values of the $`G`$, $`A`$, $`T`$, $`C`$ and $`GA`$ distributions are used as an approximate length scales $`(l_c)`$. The entire length of the sequence is divided into two parts one for $`0<l<l_{c1}`$ and other for $`l_{c2}<l<L_{max}`$ where $`l_{c1}`$ and $`l_{c2}`$ are the minimum and maximum of all the $`l_c`$ corresponding to $`G`$, $`A`$, $`T`$, $`C`$ and $`GA`$ distributions. The $`L_{max}=L/30`$, i.e. we have at least $`30`$ independent data sets so that the statistical analysis becomes meaningful. Therefore, excluding the region $`l_{c1}<l<l_{c2}`$, we have extracted $`\alpha _1`$ and $`\alpha _2`$ since the linearity in these two segments are found to be extremely good for most of the cases. The results are summarized in three tables which covers both intronless and intron containing sequences. The table shows the length of the sequence $`L`$ used in the analyses, the cross over values $`l_q`$ ( same as $`l_c`$), the slope parameters $`\alpha _1`$ and $`\alpha _2`$ and also the corresponding percentage of the nucleotide contents $`P`$. A general observation is that the sequence is weakly correlated at short distance with $`\alpha _1`$ which is quite close to $`0.5`$ where as for $`l>l_c`$, the correlation is relatively stronger with a larger value of $`\alpha _2`$. Now we discuss a few exceptions like in the case of $`SC\mathrm{\_}MIT`$ and $`PODOT7`$ ($`T7`$ bacteriophage, $`39936`$ bps). Figure 5 shows the factorial moments of a typical $`G`$ distributions. In both the cases, the factorial moments do not have any cross over point. In case of $`SC\mathrm{\_}MIT`$, the factorial moments are much higher than unity even at small distance and starts decreasing afterwards. The similar behavior is found for $`C`$ distribution also. However, the $`A`$, $`T`$ and $`GA`$ distributions do have $`l_c`$ points. Therefore, using $`l_{c1}`$ as $`36`$ and $`l_{c2}184`$, we estimated $`\alpha _1`$ and $`\alpha _2`$ for $`G`$, $`A`$, $`T`$, $`C`$ and $`GA`$ distributions which are listed in table III. The symbol $`{}_{}{}^{}_{}^{}`$ indicates absence of any critical value. It is interesting to note that $`\alpha _1`$ is quite large and in some cases $`\alpha _1>\alpha _2`$. On the other hand , the factorial moments of the sequence like $`PODOT7`$ do not reach unity at any scale. The absence of such type of scale has been indicated by the symbol $`{}_{}{}^{}_{}^{}`$ in table III. This type of sequences behave like a pure random one having $`\alpha `$ values quite close to $`0.5`$. We have listed a few such sequences with exceptions in table III.

Further, we would like to mention here that we have noticed that the factorial moments for many sequences starts decreasing at large distances. Also for a few cases, the factorial moments start decreasing even at a very short distances. Consequently, the slope also changes accordingly. However, we would not like to assign any reasons due to lack of enough statistics.

The slope with $`\alpha =0.5`$ corresponds to the case of a normal diffusion process of a random Brownian trajectory. The basic idea of a Brownian motion is that of a random walk having a Gaussian distribution probability for the position of the random walker after a time $`t`$ with the variance ($`\sigma ^2`$) proportional to $`t`$ ($`\sigma t^\alpha `$ where $`\alpha =0.5`$). This corresponds to the case of normal diffusion. However, nature shows enough examples of anomalous diffusion characterized by a variance which does not follow a linear growth in time . In such cases either the diffusion is accelerated if $`\alpha >0.5`$ or the growth is dispersive if $`\alpha <0.5`$. As found in the analyses (see tables I and II), $`\alpha _2>0.5`$ at large distances for most of the sequences irrespective of their intron contents. However, a few sequences as shown in table III, not only peculiar, may also have $`\alpha `$ which decreases at large distances. In such cases, $`\alpha <0.5`$ which may indicate the influence of dispersive dynamics. This aspect needs further investigations. Finally, we would like to add here that $`\alpha _1`$ is close to $`0.5`$ for most of the sequences at short distance (see tables I and II). Although, $`\alpha =0.5`$ would imply about a random behavior, it can not be told conclusively from the present analyses unless the short distance effects are taken into consideration .

## V Patchy sequences

In the following, we investigate whether the mosaic character of DNA consisting of patches of different composition can account for apparent long range correlations in DNA sequences. The Chargaff’s second parity rule states that in a single strand $`GC`$ and $`TA`$. However, asymmetries in base composition have been observed in many sequences. A quantitative estimate of the $`GC`$ and $`AT`$ skews can be obtained from the relation $`(GC)/(G+C)`$ (Excess of $`G`$ nucleotides over $`C`$ nucleotides) and $`(AT)/(A+T)`$ (Excess of $`A`$ nucleotides over $`T`$ nucleotides). This is, operationally equivalent to estimating $`n`$ as defined in Eq.(1) except $`n`$ now represents the count $`(GC)/(G+C)`$ for $`GC`$ skew and $`(AT)/(A+T)`$ for $`AT`$ skew in a fixed window size of $`(L/20)`$. We consider $`LAMCG`$ as an example and plot $`n`$ (defined appropriately) versus $`l_0`$ where the starting position of the sliding window $`l_0`$ varies from $`1`$, $`2`$, $`3`$ etc upto $`Ll`$. Figure 6 shows the plots of $`GC`$ and $`AT`$ skews as a function of the length for a typical $`LAMCG`$ sequence. The plots show a change in the direction of the slope with a change in sign of the skew. The quantity and quality of the skew can be assessed from the $`V`$ or from the inverted-$`V`$ shape of the curves.

From the above plots, we can identify three well known compositional domains of $`LAMCG`$ of size 22000 bps ($`GA`$ contents 0.54), 17000 bps ($`GA`$ contents 0.47) and 9000 bps ($`GA`$ contents 0.54). We also consider an artificially generated sequence by joining three random patches of size 22000 bps, 17000 bps and 9000 bps respectively with appropriate $`G`$, $`A`$, $`T`$ and $`C`$ contents. We also consider another heterogeneous sequence generated from $`E.Coli`$ DNA by a mobile insertion of TN10 at location 22000 bps. The corresponding random patches are of size 22000 bps, 9147 bps and 22000 bps respectively  Please note the distinction between the random sequence which is generated by joining three random patches of total length $`L`$ and a pure random one of length $`L`$. Although, both the sequence has same percentage of nucleotide contents in the length $`L`$, the former is random only patch wise.

Figure 7 shows the $`F_q`$ versus $`l`$ plot of a typical $`C`$ distribution for $`LAMCG`$ and for an artificially generated sequence (random only patch wise). Interestingly, the factorial moments for both the cases behave similarly. Figure 8 shows a similar $`\sigma (l)`$ versus $`l`$ plot both for real and artificially generated (from random patches) sequences. Although, in some cases both agree, in general they are not identical at the individual nucleotide levels particularly at large distances (Note that the scale is highly compressed). This deviation would mean that at large distances, the density distribution functions will have significant discrepancy due to different widths. So at a first look from the $`\sigma `$ versus $`l`$ plot, we can say that the actual DNA sequences and the RANDOM patches need not have identical slopes $`\alpha `$ (hence the width $`\sigma `$) at large distances for all the nucelotides although they agree in some cases. Even at short distances, although the DNA and the RANDOM sequences have nearly identical width $`\sigma `$, the full shape of the distributions need not be identical. To demonstrate this, we invoke the principle of complimentary which was mentioned before.

Figure 9(a) shows a $`G`$ and $`ATC`$ distribution (left most) for a $`LAMCG`$ sequence at $`l=300`$. Notice that although $`\sigma `$ versus $`l`$ plots are identical, i.e. both distributions have same fluctuations at all scales, the distribution functions are not same. Such differences are not found for a real random sequence (right most). The middle figure corresponds to the case of artificially generated random sequence. Although, the artificially generated sequence mimics the real sequence to some extent, it is not fully capable of reproducing the characteristic of a real sequence. Figure 9(b) shows another comparison for a $`E.Coli::TN10`$ sequence for $`A`$ and $`GTC`$ distributions. This discrepancy will be more prominent at higher $`l`$ values which the artificially generated sequence can not reproduce.

## VI Density distributions

In , we had demonstrated that the density distribution $`P_n`$ is Gaussian at short distances and starts deviating from it as the distance increases. Figure 10 shows another example where $`P_n`$ has been plotted for two complimentary distributions at $`l=25`$, $`100`$ and $`200`$ respectively. The complimentary distributions are nearly identical at short distance and coincide with the random distributions where as $`P_n`$ distributions for $`G`$, $`ATC`$ and pure random one are all different at larger distances.

Thus, irrespective of intron contents, most of the sequences follow Gaussian statistics at short distances. However, at large distances, the statistics deviates significantly from the Gaussian nature.

## VII Conclusions

In conclusion, we have extended our previous work to extract the slope parameter $`\alpha `$ for several intron containing and intron less DNA sequences. The advantage of the present method is that the variance analysis can be applied to any individual or group of nucleotides. We believe that the individual nucleotides provide a more fundamental measure of the correlation than any combination or group (like the DNA walk representation) where the effects may get reduced or washed out. Another interesting aspect is the (lower) factorial moments of most of the DNA sequences cross unity in a very narrow region in $`l`$ where the $`\sigma `$ versus $`l`$ plot in the log-log scale also shows a bending. Although, a formal justification to this correlation has not been provided, we have used this scale as an approximate measure to exclude the bending regions from the slope analyses. Based on this scale, we divide the DNA sequence into two segments to extract the slope parameters. It is found that below this scale, the correlation is weak and the DNA statistics is essentially Gaussian while above this all DNA sequences show strong long range correlations irrespective of their intron contents with a significant deviation from the Gaussian behavior. It may be mentioned here that the controversies that exist in this field of research are primarily due to different approaches that are adopted in various models. In this context, our analyses is model independent as it only involves the counting of an individual or a group of nucleotides in a given length to build the density distribution. In this work, we do not advocate for any specific model, although the extracted slope parameters indicate the presence of anomalous diffusion of both enhanced and dispersive nature. Instead, we provide an elegant tool to measure the degree of correlations unambiguously so that the interpretation of the data including theoretical analyses will become more meaningful. This work will also provide further impetus to develop models for the understanding of the DNA dynamics.

## Random walk model

The method of DNA walks, first suggested by Peng et al  is based on the rule that the walker either moves up $`(u_i=1)`$ or down $`u_i=1)`$ for each step $`i`$ of the walk. This is the case of a correlated random walk and differs from an uncorrelated walk where the direction of each step is independent of the previous steps. Further they assign $`u_i=1`$ if a pyrimidine occurs at the site $`i`$ whereas $`u_i=1`$ if the site contains a purine. The net displacement $`(y)`$ of the walker after $`l`$ steps is defined as

$$y(l)=\underset{i=1}{\overset{l}{}}u(i)$$
(6)
The standard deviation of the above quantity can be estimated from

$$\sigma ^2(l,L)=\frac{1}{Ll}\underset{l_0=1}{\overset{Ll}{}}(\mathrm{\Delta }y(l_0,l)\overline{\mathrm{\Delta }(l)})^2$$
(7)
where $`L`$ is the number of nucleotides in the entire sequence and

$$\overline{\mathrm{\Delta }y(l)}=\frac{1}{Ll}\underset{l_0=1}{\overset{Ll}{}}\mathrm{\Delta }y(l_0,l)$$
(8)
where $`\mathrm{\Delta }y(l_0,l)=y(l_0+l)y(l_0)`$. It was found  that the fluctuations can be approximated by

$$\sigma (l,L)l^\alpha $$
(9)
where $`\alpha `$ is the correlation exponents. For $`\alpha `$ close to $`0.5`$, there is no correlation or only short range correlation in the sequence. If $`\alpha `$ is significantly different from $`0.5`$, it indicates long range correlations.

## B

In the previos analyses, we account for the non-occurence of a particular nucleotide. This is operationally equivalent to building the density spectrum $`P_n`$ including $`n=0`$. If the nucleotide compositional asymmetry is quite large like $`SC\mathrm{\_}MIT`$, the occurence $`n`$ can be zero for some nucleotides particularly at short distances. Therefore, we can build $`P_n`$ distribution either including or excluding zeroth channel. The figure B1(a) shows the comparison of $`\sigma `$ versus $`l`$ plot for two complimentary distributions corresponding to a $`LAMCG`$ sequence both with (top panel where $`G`$ and $`ATC`$ distributions have identical slopes at all scales) and without (bottom panel) inclusion of $`n=0`$ channel in the $`P_n`$ spectra. Interestingly, absence of $`n=0`$ channel does not satisfy the complimentarity relation particularly at short distances. However, the difference does not exist at larger distances where always $`n>1`$. Figure B1(b) shows another example of $`F_q`$ versus $`l`$ plot for a typical $`SC\mathrm{\_}MIT`$ sequence. The spectrum with exclusion of $`n=0`$ channel behaves differently when zeroth channel is included (compare it with figure 5 where $`F_q`$ versus $`l`$ has no cross over).

Since the spectrum behaves differently when zeroth channel is not included, we have analysed the spectrum of three typical sequences listed in the table below. Notice now that while $`\alpha _2`$ values are essentially same as before, the $`\alpha _1`$ values are quite different. In fact, we have noticed a general trend where $`\alpha _1`$ is higer than the previous values although the corresponding density distributions do not deviate significantly from the Gaussian behavior at short distances. However, in the previous analysis, we alwyas include the zeroth channel so that the complimentarity properties is satisfied at all scales. Moreover, we also found a correlation between $`\alpha `$ and Gaussian statistics, namely the deviation of $`\alpha `$ from $`0.5`$ also shows a corresponding deviation of $`P_n`$ distribution from Gaussian behavior. For example, in case of $`SC\mathrm{\_}MIT`$, the $`\alpha `$ is quite large at a short distance. Accordingly, the $`P_n`$ distribution also shows strong deviation from the Gaussian statistics. However, this is not necessarilly true when complimentarity is not preserved while building the spectrum. At short distances, the deviation of $`\alpha `$ from $`0.5`$ does not always mean a strong deviation from the Gaussian statistics.


format: text

text:
# 1 Introduction

## 1 Introduction

There has been recently a renewed interest in the conditions under which $`\alpha `$-helices, a common structure in proteins, are formed or dissolved. It is long known that $`\alpha `$-helices undergo a sharp transition towards a random coil state when the temperature is increased. The characteristics of this so-called helix-coil transition have been studied extensively , most recently in Refs. . In Refs.  evidence was presented that the helix-coil transition in polyalanine exhibits a true thermodynamic phase transition when interactions between all atoms in the molecule are taken into account .

The later results were obtained from gas-phase simulations of poly-alanine. While there is some experimental evidence  supporting the numerical results of these gas-phase simulations, the question remains how these results relate to the biologically more relevant case of solvated molecules. First investigations of this question were described in Refs. where it was claimed that the transition temperature is lower in water than in vacuum. However, that investigation relies on a single representation of the protein-water interaction and the dependence of their results on the details of the solvation term is not clear.

In this paper, we have investigated how the characteristics of helix-coil transition change with the details of the solvation term. For this purpose, we have performed multicanonical simulations of polyalanine molecules of length 10. The protein-water interaction was included in two ways: either by a distance-dependent dielectric permittivity or by a term that is proportional to the solvent-accessible surface area of the peptide. For the later case we have considered four different parameter sets: OONS , JRF , W92  and SCH . Quantities such as the energy, helicity and susceptibility were calculated as function of temperature. Our result were compared with that of gas phase simulations. A strong dependence of the characteristics of the helix-coil transition from the details of the solvation term was found.

## 2 Methods

Our investigation of the helix-coil transition for polyalanine is based on a detailed, all-atom representation of that homopolymer. The interaction between the atoms was described by a standard force field, ECEPP/2, (as implemented in the program package SMMP ) and is given by:

$`E_{tot}`$ $`=`$ $`E_C+E_{LJ}+E_{HB}+E_{tor},`$ (1)
$`E_C`$ $`=`$ $`{\displaystyle \underset{(i,j)}{}}{\displaystyle \frac{332q_iq_j}{ϵr_{ij}}},`$ (2)
$`E_{LJ}`$ $`=`$ $`{\displaystyle \underset{(i,j)}{}}\left({\displaystyle \frac{A_{ij}}{r_{ij}^{12}}}{\displaystyle \frac{B_{ij}}{r_{ij}^6}}\right),`$ (3)
$`E_{HB}`$ $`=`$ $`{\displaystyle \underset{(i,j)}{}}\left({\displaystyle \frac{C_{ij}}{r_{ij}^{12}}}{\displaystyle \frac{D_{ij}}{r_{ij}^{10}}}\right),`$ (4)
$`E_{tor}`$ $`=`$ $`{\displaystyle \underset{l}{}}U_l\left(1\pm \mathrm{cos}(n_l\chi _l)\right).`$ (5)
Here, $`r_{ij}`$ (in Å) is the distance between the atoms $`i`$ and $`j`$, and $`\chi _l`$ is the $`l`$-th torsion angle. We have chosen ECEPP/2 instead of the newer ECEPP/3 differs because this choice allows a more easy comparison with our previous work. Both force fields differ from each other only in the way in which prolines and end groups are treated. In preliminary polyalanine simulations we found no qualitative differences in our results when ECEPP/3 was used instead of ECEPP/2 (data not shown).

The interactions between our homo-oligomer and water are approximated by means of two implicit water models. In the first model (DDE) the electrostatic interactions in the presence of water rely on a distance dependent electrostatic permittivity :

$$\epsilon (r)=D\frac{D2}{2}[(sr)^2+2sr+2]e^{sr}.$$
(6)
For the parameters $`D`$ and $`s`$ empirical values are chosen such that for large distances the permittivity takes the value of bulk water ($`\epsilon 80`$), and the value $`\epsilon =2`$ for short distances (protein interior space). Equation (6) is the result of interpolation of two types of interactions. For short distances it models the interaction of two charges placed in continuum medium, while over long distances it represents a Debye curve. This is clearly a gross over-simplification of protein-solvent interactions. However, approximating solvation effects by a distance-dependent dielectric permittivity was used by many authors to study the proteins and nucleic acids (e.g. ) since it does not significantly slow down protein simulations below that of simple in vacuo simulations.

In another common approximation of the protein-solvent interaction one assumes that the free energy contributions from atomic groups immersed in the protein interior differ from contributions of groups exposed to the water. It is commonly accepted  that this free-energy difference is proportional to the surface area of the atomic group which is exposed to the solvent. Within this approximation, the total solvation energy of a protein is given by the sum of contributions from each solvated atomic groups:

$$E_{sol}=\underset{i}{}\sigma _iA_i,$$
(7)
where $`E_{sol}`$ is the solvation energy, $`A_i`$ is the conformational dependent solvent accessible area of the surface of the $`ith`$ atom and $`\sigma _i`$ is the atomic solvation parameter for the atom $`i`$. The summation is extended over all atomic groups. The solvation parameters are evaluated experimentally by measuring the free energy needed to bring the group from a nonpolar environment (usually octanol or ethanol are used as convenient compounds) into water. Many sets of solvation parameters were evaluated by several authors with different methods, but unfortunately it is not always obvious which one is the most appropriate one. The sets we study here are named by us OONS , JRF , W92  and SCH , and are described in the respective references.

Simulations of detailed models of biological macromolecules are notoriously difficult. This is because the various competing interactions within the polymer lead to an energy landscape characterized by a multitude of local minima. Hence, in the low-temperature region, canonical Monte Carlo or molecular dynamics simulations will tend to get trapped in one of these minima and the simulation will not thermalize within the available CPU time. Only recently, with the introduction of new and sophisticated algorithms such as multicanonical sampling  and other generalized-ensemble techniques  was it possible to alleviate this problem in protein simulations . For polyalanine, both the failure of standard Monte Carlo techniques and the superior performance of the multicanonical algorithm are extensively documented in earlier work . For this reason, we use again this sophisticated simulation technique for our project.

In the multicanonical algorithm  conformations with energy $`E`$ are assigned a weight $`w_{mu}(E)1/n(E)`$. Here, $`n(E)`$ is the density of states. A simulation with this weight will lead to a uniform distribution of energy:

$$P_{mu}(E)n(E)w_{mu}(E)=\mathrm{const}.$$
(8)
This is because the simulation generates a 1D random walk in the energy space, allowing itself to escape from any local minimum. Since a large range of energies are sampled, one can use the reweighting techniques  to calculate thermodynamic quantities over a wide range of temperatures $`T`$ by

$$<𝒜>_T=\frac{𝑑x𝒜(x)w^1(E(x))e^{\beta E(x)}}{𝑑xw^1(E(x))e^{\beta E(x)}},$$
(9)
where $`x`$ stands for configurations.

Unlike in the case of canonical simulations the weights

$$w(E)=n^1(E)=e^{S(E)}$$
(10)
are not a priori known. Instead estimators for these weights have to be determined. This is often done by an iterative procedure in which for reasons of numerical stability Eq. 10 is replaced by

$$w(E)=e^{\beta (E)E\alpha (E)}.$$
(11)
The multicanonical parameters $`\beta (E)`$ and $`\alpha (E)`$ are defined through

$$\beta (E)=\frac{S(E^{})S(E)}{E^{}E}\mathrm{and}\alpha (E)=\{\begin{array}{cc}0\hfill & ,EE_{max}\hfill \\ \alpha (E^{})+(\beta (E^{})\beta (E))E^{}\hfill & ,E<E_{max}\hfill \end{array}$$
(12)
with $`E`$ and $`E^{}`$ adjacent bins in the array $`S(E)`$. The $`\beta (E)`$ are then iteratively updated  by the relation

$$\beta ^{i+1}(E)=\beta ^i(E)+g_0(E)\left(\mathrm{ln}H^i(E^{})\mathrm{ln}H^i(E)\right)/(E^{}E),$$
(13)
in which $`H^i(E)`$ is the histogram of the $`i`$-th run (and $`H(E)1`$). In Ref.  the factor $`g_0(E)`$ in Eq. 13 was defined through

$$g_0(E)=\frac{\widehat{g}^i(E)}{\widehat{g}^i(E)+_j^{i1}\widehat{g}^j(E)}\mathrm{with}\widehat{g}^i(E)=\frac{H^i(E^{})H^i(E)}{H^i(E^{})+H^i(E)}.$$
(14)
The above relation assumes that the histogram $`H(E)`$ counts independent events which is in general not true. Hence, it is more appropriate and leads to a faster convergence of $`\beta (E)`$ if the array $`\widehat{g}^i(E)`$ in Eq. 14 is instead defined by

$$\widehat{g}^i(E)=\frac{K^i(E^{})K^i(E)}{K^i(E^{})+K^i(E)}$$
(15)
where the auxiliary array $`K(E)`$ now counts only the number of independent visits at energy $`E`$.

With the above described iterative procedure we needed 200,000 sweeps for the weight factor calculations. All thermodynamic quantities were then estimated from one production run of $`1,000,000`$ Monte Carlo sweeps starting from a random initial conformation, i.e. without introducing any bias.

## 3 Results and Discussion

In previous gas-phase simulations of poly-alanine  we observed at $`T=430K`$ a pronounced transition between a high-temperature phase dominated by disordered coil structures and an ordered phase with single, extended helices. A natural order parameter for this helix-coil transition is the average number $`<n_H(T)>`$ of residues in the oligomer which are part of an $`\alpha `$helix. Following earlier work  we define a residue as helical if the pair of backbone dihedral angles $`\varphi ,\psi `$ takes a value in the range $`(70\pm 20,37\pm 20)`$. In Fig. 1a this order parameter is displayed as function of temperature for a gas-phase simulation (GP) of Ala10 and simulations with the various solvation terms. Fig. 1b shows the corresponding plots for the susceptibility $`\chi (T)`$ defined by

$$\chi (T)=<n_H^2(T)><n_H(T)>^2.$$
(16)
In Fig. 1a and 1b the curves, representing the various simulations, fall into three groups. For the case where the protein-solvent interaction was approximated by a distance-dependent permittivity (DDE), both $`<n_H>`$ and $`\chi `$ have a similar temperature dependence than is observed for poly-alanine in gas-phase simulations (GP). However, the transition temperature $`T_c`$ is shifted from $`T=435\pm 20`$ K (gas-phase) to a higher value $`T=495\pm 20`$. This temperature was determined from the maximum of the susceptibility $`\chi (T)`$ in Fig. 1b and is listed in table 1. To the same group belong the simulations in which the solvation energy was approximated by a solvent accessible surface term with either the OONS  or SCH  parameter set. In both cases susceptibility $`\chi `$ and order parameter $`<n_H(T)>`$ show also a temperature dependence similar to the one of gas-phase simulations. Only now, the transition temperature $`T_c`$ is shifted to lower temperatures. The corresponding transition temperatures can be again determined from the positions of the maximum in $`\chi (T)`$ and are also listed in table 1. The shift towards lower temperatures was one of the main results reported in Refs.  for simulations with the OONS solvation energy, and our $`T_c=345\pm 20`$ K agrees well with their value $`T_c=340`$ K (no errors quoted) in Refs. .

A somehow different behavior is observed in the simulation where the protein-water interaction was approximated by a solvent accessible surface term relying on the W92  parameter set. Here, the form of $`<n_H>`$ indicates only partial helix formation and occurs only at much lower temperatures. The susceptibility $`\chi (T)`$ in Fig. 1b gives no indication for a helix-coil transition. For this reason no value of $`T_c`$ is listed for the W92 parameter set in table 1. Instead, we observe in Fig. 2 for this case at low temperatures even the appearance of residues whose backbone dihedral angles $`\varphi ,\psi `$ take values typical for a $`\beta `$-sheet $`(150\pm 30,150\pm 30)`$.

Yet another behavior is observed in simulations where the solvation energy of Eq. 7 is evaluated by means of the JRF parameter set. No formation of helices or sheets is observed in Figs. 1 and 2. Since no transition temperature can be determined, we do not list a value of $`T_c`$ for the JRF parameter set in table 1.

The same grouping can be found in Fig. 3a-f where we display various energy terms as a function of temperature. In these figures we have shifted the solvation energies and the partial ECEPP/2 energies $`E_C,E_{LJ},E_{HB}`$ and $`E_{tor}`$ of Eq. 5 by a constant term such that we have for all solvation models at $`T=1000`$ K $`E_{sol}=0`$ and $`E_C=E_{LJ}=E_{HB}=E_{tor}=0`$. Such a shift by an irrelevant constant allows a better comparison of the different simulations. The average total energy $`<E_{tot}>`$ which is the sum of intramolecular potential energy $`E_{ECEEP/2}`$ and the solvation energy $`E_{sol}`$, is displayed in Fig. 3a. We observe again that simulations with the parameters sets OONS and SCH, and such with distant dependent permittivity (DDE), have a similar temperature dependence as gas phase simulations (GP). On the other hand, in simulations relying on the W92 parameter set, the energy varies less with temperature and is at low temperatures considerably higher than in the simulations with other solvation energy terms. Finally, the energy in simulations with the JRF parameter set is an almost linear function of temperature and is especially at high temperatures much lower than the energies found in gas phase simulations.

The dissimilar behavior of energy for simulations with different solvation terms is even more obvious in Fig. 3b where the average intramolecular energy $`E_{ECEPP/2}`$ is drawn. While this energy term decreases between $`1000`$ K and $`150`$ K by $`50`$ Kcal/mol (with most of that change, $`30`$ Kcal/mol, happing around the respective transition temperature $`T_c`$) in gas-phase simulations (GP) and in simulations with OONS, SCH and DDE solvation terms, it changes in the same temperature interval only by $`20`$ Kcal/mol in simulations utilizing the JRF or W92 parameter sets. Since for these two parameter sets also no or only little helix formation was observed it seems likely that the formation of helices is related to the large gain in potential energy observed for GP, OONS,SCH and DDE simulations. This gain in potential energy is in part due to the formation of hydrogen bonds between a residue and the fourth following one in the polypeptide chain which stabilize an $`\alpha `$-helix. Fig. 3c displays the average hydrogen-bonding energy $`<E_{HB}>`$ of Eq. 5 as a function of temperature and one can clearly see the gain in energy for the GP, DDE, OONS and SCH simulations at the respective helix-coil transition temperatures of table 1. No such gain is observed in W92 and JRF simulations where also no helix formation was found. A similar gain in energy with helix formations in gas-phase and simulations with DDE,OONS and SCH solvent representations is also observed for the average Lennard-Jones energy $`<E_{LJ}>`$ and the electrostatic energy $`<E_C>`$ displayed in Fig. 3d and 3e, respectively. Note also in Fig. 3e the large gain in $`E_C`$ for DDE at the helix-coil transition temperature which additionally stabilizes the $`\alpha `$-helix in this model.

A complementary picture is found in Fig. 3f where the solvation energy $`E_{sol}`$ is shown as a function of temperature. The observed helix formation in gas phase simulations and such with OONS,SCH and DDE solvent representations is correlated with an increases of the solvation energies by $`5`$ Kcal/mol. On the other hand, in simulations with the W92 and JRF parameter sets, for which no helix-formation was observed in Fig. 1, $`E_{sol}`$ decreases with temperature. This decrease is only $`5`$ kcal/mol for W92, but it is much larger (of order $`30`$ kcal/mol) in the case of JRF where the solvation energy is the dominant term.

The effects of the dominant solvation term in simulations with the JRF parameter set can also be seen in Fig. 4. In this figure the average radius of gyration, a measure for the compactness of configurations, is shown as a function temperature. One can see that this quantity changes little with temperature for the JRF data. However, its value is over the whole temperature range considerably smaller than observed in the other simulations. This indicates that the JRF solvation term favors already at high temperatures compact configurations, and that the pressure towards compact structure is such that the more elongated helices cannot be formed. Note however, that the tendency towards compact configurations does not lead to a lower Lennard-Jones energy $`E_{LJ}`$ as one would expect. Fig. 3d indicates that $`<E_{LJ}>`$ is at low temperatures in JRF simulations even larger than in GP, DDE, OONS and SCG simulations where helix-formation was observed. The tendency towards compact structures in JRF simulations may be due to the fact that JRF parameter set was developed from minimum energy (i.e. compact) conformations of peptides (the low-energy conformations of 13 tetrapeptides derived by NMR studies ), and therefore this parameter set may have an intrinsic bias towards compact structures.

On the other hand, the W92 parameter set was developed from measurements of free energies of amino acid side chain analogs from vapor to water . The parameters for this set are negative for all atoms except carbon meaning that the nitrogen, oxygen and sulfur atoms are considered hydrophilic, i.e. favoring solvent exposure. This explains not only the small solvation energies observed for this parameter set in Fig. 3f, but also why in Fig. 4 the radius of gyration is consistently larger for this parameter set than for the others indicating that extended configurations are favored with this parameter set. This bias towards extended structures limits again the formation of $`\alpha `$-helices.

While the OONS parameter set was derived from experimental free energies of gas-to-water transfer of small aliphatic and aromatic molecules, the SCH is not directly based on experimental free energy data. Instead, it was developed as an optimized parameter set to complement the CHARMM force field . In both parameter sets the hydrophobic character of the carbon atoms is increased and the hydrophilic character of uncharge oxygen and nitrogen atoms decreased resulting into the large solvation energies of these two parameter sets (when compared with the one of the W92 parameter set) that one observes in Fig. 3f. The OONS and SCH solvation energies again favor extended structures (the radius of gyration has larger values than found in gas-phase simulations), however, the interplay of solvation energies and intramolecular ECEPP/2 energy is such that the radius of gyrations (and consequently the compactness) of polyalanine configurations as a function of temperature shows a similar behavior as the gas-phase simulation. The same is true for the DDE simulation where the protein-solvent interaction was approximated by a distance-dependent permittivity.

Our results demonstrate that the helix formation is due to the gain in potential (intramolecular) energy while (with the exception of the JRF parameter set) the solvent-accessible surface terms favor extended peptide configurations. Table 2 summarizes the differences in total energy $`\mathrm{\Delta }E_{tot}`$ , solvation energy $`\mathrm{\Delta }E_{sol}`$, potential energy $`\mathrm{\Delta }E_{ECEEP/2}`$, and the partial energies $`\mathrm{\Delta }E_C`$, $`\mathrm{\Delta }E_{LJ}`$, $`\mathrm{\Delta }E_{HB}`$ and $`\mathrm{\Delta }E_{tor}`$ between complete helical configurations (all residues with exception of the terminal ones are part of an $`\alpha `$-helix) and coil configurations at temperature $`T=300`$ K for gas-phase, DDE OONS and SCH simulations. Note, that the intramolecular energy differences $`\mathrm{\Delta }E_{ECEPP/2}`$ of gas-phase, OONS and SCH simulations have within their error bars the same values. For simulations with the W92 parameter set the longest found helix consists of 6 consecutive residues. Hence, we measured for this case only the energy difference between configurations with at least three consecutive helical residues (i.e. one turn of an $`\alpha `$-helix) and coil configurations. This modified definition of the energy differences is also the reason for the smaller value of $`\mathrm{\Delta }E_{ECEPP/2}`$ listed for W92 in table 2. We do not list energy differences for the JRF parameter set since no helices were found in simulations utilizing this parameter set.

Note that in simulations with distant dependent permittivity (DDE) helices are energetically more favored than in the gas-phase simulations. This is due to the increased contribution from the Coulomb term $`E_C`$ as one can also see in Fig. 3e. The larger energy gap between helical and coil conformations (when compared with gas-phase simulations) explains why the transition temperature is higher in DDE simulations than in gas-phase simulations.

For the OONS and the SCH parameter set the solvation energy difference $`\mathrm{\Delta }E_{sol}`$ is positive (indicating that coil structures are energetically favored), but its magnitude is only approximately half that of the potential energy difference $`\mathrm{\Delta }E_{ECEPP/2}`$. Hence, there is still an overall energetic gain connected with helix formation. However, in both cases the total energy difference between helical and coil configurations is reduced by the solvation energy when compared with the gas-phase simulation. This reduction of the energy gap leads to the lower transition temperatures observed in OONS and SCH simulations.

On the other hand, for the W92 parameter set we find that $`\mathrm{\Delta }E_{ECCP/2}`$ and $`\mathrm{\Delta }E_{sol}`$ are of same magnitude so that helical configurations are not or only weakly energetically favored. This is consistent with our results in Fig. 1a and 1b where we find at $`T=280`$ K a high average helicity in OONS and SCH simulations but only a small value of $`<n_H>`$ and no indications for a helix-coil transition in W92 simulations. An evaluation of energy differences was not possible for simulations with the JRF parameter set since no helices were found.

The above results indicate that the existence and characteristics of the helix-coil transition in polyalanine depend strong on the details of the solvent representation. In order to evaluate the validity of the different solvent models one has to compare the numerical results with experimental data. For this purpose we have calculated the helix propagation parameter $`s`$ which was also determined by experiments . According to the Zimm-Bragg model  the average number of helical residues $`<n>`$ and the average length $`<\mathrm{}>`$ of a helical segment are given for large number of residues $`N`$ by

$`{\displaystyle \frac{<n>}{N}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{1s}{2\sqrt{(1s)^2+4s\sigma }}},`$ (17)
$`<\mathrm{}>`$ $`=`$ $`1+{\displaystyle \frac{2s}{1s+\sqrt{(1s)^2+4s\sigma }}},`$ (18)
where $`s`$ is the helix propagation parameter and $`\sigma `$ the nucleation parameter of the Zimm-Bragg model. From these equations with the values of $`<n>/N`$ and $`<\mathrm{}>`$ calculated from the multicanonical production runs, we have calculated $`s`$ at temperature $`T=280`$ K for gas-phase and the different solvation models. Our values are summarized in table 3 which also lists our $`\sigma `$ values. Our results for gas-phase, DDE and OONS simulations are in agreement with the experimental results of Ref.  where they list values of $`s`$(Ala) between $`1.5`$ and $`2.19`$. On the other hand, the $`s`$ value obtained in the SCH simulation agrees well with the one obtained by the host-guest technique of Ref. . However, the $`s`$ values which were obtained in W92 or JRF simulations do not agree with either of the experimental data. Hence, we conclude that the W92 and JRF parameter sets are not appropriate solvation models in simulations of polyalanine. Otherwise, the variation in the experimental data is too large to give indications whether one of the remaining solvent representations (DDE, OONS, SCH or even no solvent at all (GP)) is preferable over the others.

## 4 Conclusions

We have performed multicanonical simulations of polyalanine. The intramolecular forces were modeled by the ECEPP/2 force field and various approximations for the solvation energy were studied. We observed that whether a helix-coil transition is observed for poly-alanine, and at what temperature, depends strongly on the chosen approximation for the protein-solvent interaction. Our results demonstrate both the importance (and need) of including solvation terms into protein simulations and the difficulties in chosing an adequate representation of the protein-water interactions. Especially when using the solvent-accessible surface approach, it seems necessary to carefully choose a parameter set that is suitable for the problem under consideration. Use of a specific parameter set without further justification could otherwise generate miss-leading results.

## Acknowledgement

U. Hansmann gratefully acknowledges support by a research grant from the National Science Foundation (CHE-9981874). This article was written in part while U.H. was visitor at the Department of Physics at University of Central Florida. He thanks Alfons Schulte, Weili Luo, Aniket Bhattacharya and Brian Tonner for their kind hospitality during his stay in Orlando.

Tables:

Figure Captions:

Temperature dependence of (a) the average number $`<n_H>`$ of helical residues and (b) the susceptibility $`\chi (T)`$ for ALA10 as calculated from a gas-phase simulation and from simulations with various solvation energy terms. All results rely on multicanonical simulations of 1,000,000 Monte Carlo sweeps each.

Temperature dependence of the average number $`<n_B>`$ of residues whose backbone dihedral angles $`\varphi ,\psi `$ take values as typically found in $`\beta `$-sheets. Results from a gas-phase simulation and such with various solvation terms are presented for ALA10. All data rely on multicanonical simulations of 1,000,000 Monte Carlo sweeps.

Temperature dependence of (a) the total energy $`<E_{total}=E_{ECEEP/2}+E_{sol}>`$, (b) the intramolecular energy $`<E_{ECEPP/2}>`$, (c) the hydrogen-bonding energy $`<E_{HB}>`$, (d) Lennard-Jones energy $`<E_{LJ}>`$, (e) Coulomb energy $`<E_C>`$, and (f) the solvation energy $`<E_{sol}>`$ as calculated from a gas-phase simulation and from simulations with various solvation energy terms. All results rely on multicanonical simulations of ALA10 with 1,000,000 Monte Carlo sweeps for each case.

Temperature dependence of the average radius-of-gyration $`<R_{gy}>`$ as measured in gas-phase simulations and simulations with various solvent representations. All data rely on multicanonical simulations of 1,000,000 Monte Carlo sweeps.


text:
# Detectability of the Supernova Relic Neutrinos and Neutrino Oscillation

## 1 Introduction

A core-collapse supernova explosion produces a number of neutrinos and 99% of the gravitational energy is transformed to neutrinos. It is generally believed that the core-collapse supernova explosions have traced the star formation history in the universe and have emitted a great number of neutrinos, which should make a diffuse background. This supernova relic neutrino (SRN) background is one of the targets of the currently working large neutrino detectors, SuperKamiokande (SK) and Sudbury Neutrino Observatory (SNO). Comparing the predicted SRN spectrum with the observations by these detectors provides us potentially valuable information on the nature of neutrinos as well as the star formation history in the universe. This SRN background has been discussed in a number of previous papers -. The work after Totani et al.  takes into account the realistic star formation history inferred from various observations and theoretical modeling of galaxy formation, to calculate the SRN flux. Totani et al. calculated the energy flux of SRN at SK detector and compared it with neutrinos emitted by other sources (solar, atmospheric, and reactor neutrinos and so on). Then they concluded that the visible event rate of SRN at SK is 1.2 yr-1 in the energy range from 15 to 40 MeV. On the other hand, Kaplinghat et al.  calculated the upper limit of SRN and discussed the possibility of the detection at SK detector comparing to the other backgrounds. Their result is that at the energy range from 15 to 40 MeV, where no significant background had been considered to exist in previous studies, there is a huge background of the invisible muon decay and the detection of SRN is difficult. They also discussed the effects of neutrino oscillation and briefly mentioned that if the SRN flux is in the vicinity of their upper bound and all three flavors are maximally mixed, it may be detectable as a distortion of the expected muon background.

In this paper we calculate the SRN flux and the event rate at SK, and discuss about the detectability of SRN, with the following new aspects compared with previous studies: 1) realistic neutrino oscillation parameters are incorporated based on the recent solar and atmospheric neutrino experiments, 2) a realistic neutrino spectrum from one supernova explosion is used, which is obtained from a numerical simulation by the Lawrence Livermore group, and 3) we have examined other contaminating background events against the detection of SRN, in more detail than previous studies.

Recent experiments of SK and SNO support neutrino oscillation by solar - and atmospheric neutrino  data. Although some previous studies (e.g. ) mentioned the possibility of neutrino oscillation, no quantitative calculation of SRN has been made incorporating the realistic oscillation parameters. Here we consider four neutrino oscillation models, which satisfy the solar and atmospheric neutrino data, and investigate the dependence of the positron spectra at SK on the oscillation models. If neutrino oscillation occurs, $`\overline{\nu }_{\mu ,\tau }`$’s are converted into $`\overline{\nu }_e`$’s which are mainly detected at SK detector. Because $`\overline{\nu }_{\mu ,\tau }`$’s interact with matter only through the neutral-current reactions in supernovae, they are weakly coupled with matter compared to $`\overline{\nu }_e`$. Thus the neutrino sphere of $`\overline{\nu }_{\mu ,\tau }`$’s is deeper in the core than $`\overline{\nu }_e`$’s and their temperatures are higher than $`\overline{\nu }_e`$’s. Therefore neutrino oscillation enhances the mean $`\overline{\nu }_e`$’s energy and enhances event rate at SK detector.

In all of the past studies, neutrinos emitted from a supernova are assumed to obey the Fermi-Dirac distribution with zero chemical potential. However, since neutrinos are not in the thermal equilibrium states in supernovae, the real spectrum should be different from the pure Fermi-Dirac distribution. Thus we use in this paper a realistic supernova model established by the Lawrence Livermore group  which in fact shows clear difference from the Fermi-Dirac distribution .

There are several background events which hinder the detection of SRN. These includes atmospheric and solar neutrinos, anti-neutrinos from nuclear reactors, and decay electrons from invisible muons. We should find the energy region which is not contaminated by these background events and then calculate the detectable event rate of SRN. By careful examination of these events, we found that there is a narrow energy window at which the detection of the SRN event might be possible.

In addition we use three models of the cosmic star formation history, i.e., the evolution of star formation rate (SFR) density of high redshift galaxies (e.g., ), which is inferred from several star formation indicators such as rest-frame UV luminosity, H$`\alpha `$ lines, or dust emission in far-infrared or submillimeter wavebands. It should be noted that SFR models contain various uncertainties, especially at high redshift regions. (e.g., the faint ends of luminosity functions have not been well established at high redshifts and the uncertainties in dust extinction are large at all redshifts .) This is the reason we applied several SFR evolution models described below.

Throughout this paper, we consider only electron antineutrinos ($`\overline{\nu }_e`$’s) from collapse-driven supernovae, because $`\overline{\nu }_e`$’s are most easily detected in a water Čherenkov detector like SK.

This paper is organized as follows: In Section 2, we illustrate the neutrino oscillation models, the supernova model, and the supernova rate models considered in this paper and discuss about these models. Formulations for calculation of flux and event rate at SK are given in Section 3. In Section 4, we show the calculated flux and event rate. Detailed discussions especially on other background events are presented in Section 5.

## 2 Models of Neutrinos and Supernova Rate

### 2.1 Supernova Model

We use a realistic model of a collapse-driven supernova calculated by the Lawrence Livermore group . We show in Fig. 1 the time-integrated spectrum for $`\overline{\nu }_e`$’s (e.g., see Ref.  for detail). We also show the Fermi-Dirac (FD) distribution with zero chemical potential for comparison. Comparing with the FD model used in the other studies, the deficit of both low- and high-energy neutrinos can be seen.

In this paper, we assume that this numerical model represents all of the past supernovae, although in this model the progenitor of the supernova has been assumed to have the mass of $`20M_{}`$. It is clearly an oversimplification, but we note that the mean mass of progenitor stars of type II supernovae (above $`8M_{}`$) weighted by number is about $`15M_{}`$ when a typical initial mass function is applied .

Another assumption made in the simulation is the isotropic radiation of neutrinos. All supernova progenitors are rotating, and it may have significant effect on the degree of isotropy of neutrino emission , although the effect is very difficult to estimate quantitatively. However, since our interest is in the SRN background, which is the sum of all past supernova neutrinos, the rotation effect is expected to be small.

### 2.2 Neutrino Oscillation Models

We adopt the four parameter sets for the neutrino mixing (see Table 1) and the normal mass hierarchy. These parameter sets are introduced to explain the observations of the solar and the atmospheric neutrinos -. The neutrino spectra emitted by each supernova are calculated numerically by Takahashi et al.  for the same parameter sets assuming the normal mass hierarchy. (Supernova model they have used is the same one discussed in the previous subsection.) We use their calculated $`\overline{\nu }_e`$ spectrum as that from any supernova.

In Table 1, “LMA” and “SMA” indicate MSW solution of the solar neutrino problem (see Ref. , for the review of MSW effect). Recent SK and SNO observations  show that LMA solution are more favorable but SMA solution cannot be rejected at more than $`3\sigma `$ level (see Ref.  and references therein). Therefore we try the SMA solution as well. The suffixes “-L” and “-S” indicate whether $`\theta _{13}`$ is large or small. A large (small) $`\theta _{13}`$ means “higher resonance” is adiabatic (nonadiabatic). (There are two resonance points where neutrino mass eigenstates may flip in the supernova matter, and we call them lower and higher resonance points. The lower resonance point is at lower density, and the higher resonance point at higher density.) The adiabatic higher resonance enhances the energy of the electron neutrinos and enhances the event rate of $`\nu _e`$ scattering, but does not affect the electron anti-neutrinos . In this paper, since we deal with only the reaction of $`\overline{\nu }_e`$’s, this “-L” and “-S” hardly influence the result.

Now, we semi-qualitatively present a simple illustration of the expected spectral shape of $`\overline{\nu }_e`$’s for LMA and SMA models. We can naively deal with anti-neutrino oscillation effect as vacuum oscillation, since $`\overline{\nu }_e`$’s are not affected by the resonance. Further we take two generation formalism, for simplicity (this approximation is justified when $`\theta _{13}`$ is sufficiently small). The conversion probability from $`\overline{\nu }_\mu `$ to $`\overline{\nu }_e`$ (and its inverse) averaged over distance is

$$P(\overline{\nu }_\mu \overline{\nu }_e)=P(\overline{\nu }_e\overline{\nu }_\mu )=\frac{1}{2}\mathrm{sin}^22\theta _{12}.$$
(1)
This probability depends only on $`\theta _{12}`$. Then, $`P(\overline{\nu }_\mu \overline{\nu }_e)=P(\overline{\nu }_e\overline{\nu }_\mu )=0.44(\mathrm{LMA}),2.5\times 10^3(\mathrm{SMA})`$. Using these conversion probabilities we can present $`\overline{\nu }_e`$ flux as

$`F_{\overline{\nu }_e}`$ $`=`$ $`[1P(\overline{\nu }_e\overline{\nu }_\mu )]F_{\overline{\nu }_e}^0+P(\overline{\nu }_\mu \overline{\nu }_e)F_{\overline{\nu }_\mu }^0`$ (2)
$``$ $`\{\begin{array}{cc}0.57F_{\overline{\nu }_e}^0+0.44F_{\overline{\nu }_\mu }^0\hfill & (\mathrm{LMA})\hfill \\ F_{\overline{\nu }_e}^0\hfill & (\mathrm{SMA})\hfill \end{array},`$ (5)
where $`F_{\overline{\nu }}^0`$ means neutrino flux for no oscillation model. This equation shows that in LMA models the $`\overline{\nu }_e`$ flux at higher energy region is enhanced, on the other hand that at lower energy region is suppressed, because the original mean $`\overline{\nu }_\mu `$’s energy is higher than that of $`\overline{\nu }_e`$’s as mentioned in Section 1.

However, in the case of the inverted hierarchy, this situation changes dramatically although we do not give calculation for this case. In this case, the higher resonance point is in the anti-neutrino sector , and hence resonant conversion of $`\overline{\nu }_e`$’s and $`\overline{\nu }_{\mu ,\tau }`$’s is possible. The conversion probability would be dependent on the neutrino energy. Another possibility of such resonant conversion effect is the spin-flip oscillation of $`\overline{\nu }_e`$’s by a flavor-changing neutrino magnetic moment .

### 2.3 Supernova Rate

A number of studies have modeled the expected evolution of the cosmic SFR with redshift. The cosmic SFR can now be traced to $`z4`$ observationally, although some details remain controversial. We use three SFRs per unit comoving volume which are also used in Ref. , where all of them are derived following the same method employed by Madau et al.  first. They use the intergalactic absorption to identify high redshift galaxies in broadband multicolor survey, and convert the UV luminosity to SFR and the metal ejection rate assuming a typical initial mass function (IMF) of stars. The supernova rate is expected to be proportional to SFR since the lifetime of progenitors of core-collapse supernovae is much shorter than the cosmological time scale. We show in Fig. 2 the supernova rates corresponding to these three SFR models that are explained below.

In the Einstein-de Sitter (EdS) universe ($`\mathrm{\Omega }_m=1.0,\mathrm{\Omega }_\lambda =0.0`$), the first (hereafter SF1) is taken from Ref. :

$$R_{SF1}(z)=0.3h_{65}\frac{\mathrm{exp}(3.4z)}{\mathrm{exp}(3.8z)+45}M_{}\mathrm{yr}^1\mathrm{Mpc}^3,$$
(6)
where $`h_{65}=H_0/65`$km s-1 Mpc-1. This SFR increases rapidly between $`z=0`$ and $`z=1`$, peaks between $`z=1`$ and $`z=2`$, and gently declines at higher redshifts. This includes an upward correction for dust reddening of $`A_{1500}=1.2`$ mag. The original SFR (before correction) should be interpreted as lower limits to the real values, since a significant fraction of UV light from young stars could be absorbed and re-emitted in far-infrared band. In the second model, the SFR remains instead roughly constant at $`z2`$ (SF2) , as:

$$R_{SF2}(z)=0.15h_{65}\frac{\mathrm{exp}(3.4z)}{\mathrm{exp}(3.4z)+22}M_{}\mathrm{yr}^1\mathrm{Mpc}^3.$$
(7)
We should consider SF2 in addition to SF1, because of the uncertainties associated with the incompleteness of the data sets and the amount of dust extinction at early epochs. The third SFR (SF3) ,

$$R_{SF3}(z)=0.2h_{65}\frac{\mathrm{exp}(3.05z0.4)}{\mathrm{exp}(2.93z)+15}M_{}\mathrm{yr}^1\mathrm{Mpc}^3$$
(8)
represents even more star formation at early epochs. This SFR is based on the studies which suggest that the evolution of the SFR up to $`z1`$ may have been overestimated , while the rates at high-$`z`$ may have been severely underestimated due to large amounts of dust extinction .

These SFR evolutions should change when a different cosmological model is assumed, since they are based on the observed data of high-$`z`$ galaxies. The correction to other cosmological models can be written as (e.g., ):

$`R_{SF}(z;\mathrm{\Omega }_m,\mathrm{\Omega }_\lambda ,h_{65})`$ $`=`$ $`h_{65}{\displaystyle \frac{\sqrt{(1+\mathrm{\Omega }_mz)(1+z)^2\mathrm{\Omega }_\lambda (2z+z^2)}}{(1+z)^{3/2}}}`$ (9)
$`\times R_{SF}(z;1,0,1).`$
To obtain the supernova rate ($`R_{SN}`$), we multiply the SFRs by the coefficient

$$\frac{_8^{125}𝑑m\varphi (m)}{_{0.01}^{125}𝑑mm\varphi (m)}=0.0122M_{}^1,$$
(10)
where $`\varphi (m)`$ is the Salpeter IMF ($`\varphi (m)m^{2.35}`$) and $`m`$ is the stellar mass in solar units . (Here, we assume that all stars whose mass is greater than $`8M_{}`$ explode as core-collapse supernovae.) This resulting rates agree with the local observed value. We label these models “SN1”, “SN2”, and “SN3”, respectively.

## 3 Formulation of Flux Calculations

In section 2, we calculated the number of emitted $`\overline{\nu }_e`$’s from a supernova per unit energy $`q`$, $`dN_\nu /dq`$, and supernova rate, $`R_{SN}`$. We remark that $`R_{SN}(z)`$ is supernova rate per comoving volume, and hence we should multiply the factor $`(1+z)^3`$ to obtain the rate per physical volume at that time. The present number density of $`\overline{\nu }_e`$’s, whose energy is in the interval of the $`qq+dq`$, emitted in the interval of the redshift $`zz+dz`$ is given by

$`dn_\nu (q)`$ $`=`$ $`R_{SN}(z)(1+z)^3{\displaystyle \frac{dt}{dz}}dz{\displaystyle \frac{dN_\nu ((1+z)q)}{dq}}(1+z)dq(1+z)^3`$ (11)
$`=`$ $`R_{SN}(z){\displaystyle \frac{dt}{dz}}dz{\displaystyle \frac{dN_\nu ((1+z)q)}{dq}}(1+z)dq,`$ (12)
where the factor $`(1+z)^3`$ comes from the expansion of the universe. The Friedmann equation gives the relation between $`t`$ and $`z`$ as follows:

$$\frac{dz}{dt}=H_0(1+z)\sqrt{(1+\mathrm{\Omega }_mz)(1+z)^2\mathrm{\Omega }_\lambda (2z+z^2)}.$$
(13)
We now obtain the differential number flux of SRN, $`dF_\nu (q)/dq`$, using the relation $`dF_\nu (q)/dq=c[dn_\nu (q)/dq]`$:

$`{\displaystyle \frac{dF_\nu }{dq}}`$ $`=`$ $`{\displaystyle \frac{c}{H_0}}{\displaystyle _0^{z_{\mathrm{max}}}}R_{SN}(z){\displaystyle \frac{dN_\nu ((1+z)q)}{dq}}`$ (14)
$`\times {\displaystyle \frac{dz}{\sqrt{(1+\mathrm{\Omega }_mz)(1+z)^2\mathrm{\Omega }_\lambda (2z+z^2)}}},`$
where we assume that gravitational collapses began at the redshift parameter $`z_{\mathrm{max}}=5`$.

It can be seen that, from equations (9) and (14), the flux does not depend on the cosmological parameters such as $`\mathrm{\Omega }_m,\mathrm{\Omega }_\lambda ,`$ and $`H_0`$, while in Ref. , the flux depend on these parameters. The illustration of this difference is as follows: Totani et al. used supernova rate evolution derived from their theoretical model of galaxy evolution which reproduces various properties of present-day galaxies. On the other hand, our supernova rate is based on the observational estimate of luminosity densities of high-$`z`$ galaxies, and the dependence of the cosmological volume element on the parameters such as $`\mathrm{\Omega }_m`$ and $`\mathrm{\Omega }_\lambda `$ is cancelled out. Therefore, our supernova rate depends on the cosmological parameters while the neutrino flux does not depend on them.

## 4 Results

### 4.1 Calculation of SRN Flux

We calculate the flux of the SRN using the formula eq. (14) for various models. In Fig. 3, the flux for the three supernova rate models are shown assuming no oscillation model. The fluxes of these models are almost the same above $`8`$ MeV and the models with more SFR at early epochs have a higher peak at lower neutrino energy. These properties come from the following effects: The energy of neutrinos which were emitted at redshift $`z`$ is reduced by a factor $`(1+z)^1`$ when we observe. Then at observation, the high energy tail ($`8\mathrm{MeV}`$) is mainly contributed by low redshift supernovae. Since at low redshift three supernova rate models are almost the same, the fluxes at $`E_{\overline{\nu }_e}8\mathrm{MeV}`$ are not much different. Similarly, because high redshift supernova neutrinos contribute more at low energy region when we observe, we can see model dependence at low energy clearly.

In Fig. 4, we show the flux of various neutrino oscillation models assuming SN1 model. In SMA and no oscillation models we can see higher peak around $`5`$ MeV than in LMA models, because in LMA models the low energy $`\overline{\nu }_e`$’s are deficit due to the conversion into $`\overline{\nu }_{\mu ,\tau }`$’s. But above $`10`$ MeV, in LMA models we can see more flux, because $`\overline{\nu }_{\mu ,\tau }`$’s which have higher mean energy at production have more changed into $`\overline{\nu }_e`$’s than in the case of SMA or no oscillation. (See also qualitative explanation in Section 2.2.)

Integrated flux over the entire neutrino energy range is shown in Table 2.

### 4.2 Event Rate at SuperKamiokande Detector

The SK detector is a water Čerenkov detector whose fiducial mass is 22,500 ton. The detector efficiency is 100$`\%`$ for electrons (positrons) whose energy is above 5 MeV, and 50$`\%`$ at 4.2 MeV. We only consider the reaction $`\overline{\nu }_e+pe^++n`$, because the cross section \[$`9.52\times 10^{44}(E_e/\mathrm{MeV})(p_e/\mathrm{MeV})\mathrm{cm}^2`$ \] is much larger than the other reactions.

In Fig. 5, the event rate at SK of the three supernova rate models are shown assuming no oscillation model. Because of the detector threshold, we see only the positrons whose energies are above $`5`$ MeV, so that the differences between models are small as shown also in Fig. 3 where the model dependence is also small for neutrinos above $`7`$ MeV. Then, although our three supernova rate models have very different properties at high redshift regions, our results are not influenced by the behavior of SFRs at those high $`z`$ regions. (These results are hardly changed when we assume the other LMA or SMA models, instead of the no oscillation model.)

In Fig. 6, we show the event rate of various neutrino oscillation models assuming SN1 model. We can see clear difference of the LMA model from SMA or no oscillation models, especially at the high energy tail. This property results from the flux dependence on oscillation models (see also Fig. 4).

As a result, if we can detect SRN events above $`10`$ MeV, we can discriminate the LMA model from the SMA or no oscillation models in any supernova rate models and any sets of cosmological parameters.

## 5 Discussion

### 5.1 Background Events against the Detection

We discuss in this subsection about neutrinos from other sources which may become an obstacle to the SRN detection. They are atmospheric and solar neutrinos, anti-neutrinos from nuclear reactors, and decay electrons from invisible muons. We show in Fig. 7 the number flux of SRN and these background events, and we discuss the each background event below.

The flux of the atmospheric neutrinos is usually calculated using Monte Carlo method including various relevant effects (flux of primary cosmic rays, solar modulation, geomagnetic field, interaction of cosmic rays in the air, and so on), and in that simulation one-dimensional approximation is used, i.e., after the interaction of primary cosmic ray particles with air nuclei, all the particles are assumed to be moving along the line of the momentum vector of the primary cosmic ray particles . ( Recently preliminary results of three-dimensional flux calculations have been reported .) There are many authors who calculated the atmospheric neutrino flux (see, e.g.,  for a recent result). We use in this paper the flux calculated by Gaisser et al.  and Barr et al. . Although their calculations are rather old, the fluxes of low energy neutrinos, in which we are interested, are also calculated, while most of other papers show only higher energy region ($`>1\mathrm{GeV}`$).

Solar neutrino flux is dominant at energy range 10–20 MeV. We use the flux predicted by the standard solar model (SSM) in Fig. 7 . Since the solar neutrinos are not $`\overline{\nu }_e`$’s but $`\nu _e`$’s, the cross section for solar neutrinos is about two order smaller than that for $`\overline{\nu }_e`$’s. Furthermore recoil electrons scattered by solar neutrinos strongly concentrate to the opposite direction of the Sun, in contrast to the isotropic distribution of $`\overline{\nu }_e`$ events. Therefore the solar neutrino is an avoidable background. In fact, we can demonstrate how we can avoid the solar neutrino background, as follows. The solar neutrino event number is about $`5\times 10^3`$ $`\mathrm{yr}^1(22.5\mathrm{kton})^1`$ , about 3 or 4 orders of magnitude larger than our expected event rate $`1`$ $`\mathrm{yr}^1(22.5\mathrm{kton})^1`$. The recoil electrons scatter obeying the Gaussian of one-sigma error $`\sigma 25^{}`$ for 10 MeV electrons . Then, 99.98% of solar neutrino events within $`3.7\sigma `$ (corresponding to $`90^{}`$) can be avoided. Therefore, restricting our discussion on the nearby side of the hemisphere from the Sun, we can ignore the solar neutrino events.

The third background which we must consider is anti-neutrinos from nuclear reactors. In each nuclear reactor, almost all the power comes from the fissions of the four isotopes, 235U ($`75\%`$), 238U ($`7\%`$), 239Pu ($`15\%`$), and 241Pu ($`3\%`$) . Each isotope produces a unique electron anti-neutrino spectrum through the decay of its fission fragments and their daughters. The $`\overline{\nu }_e`$ spectrum from 235U, 239Pu, and 241Pu can be derived using the semi-empirical formula with which we fit data of detected $`\beta `$-spectrum from fission by thermal neutrons . (238U undergoes only fast neutron fission and hence electron spectrum from 238U cannot by measured by this kind of experiment.) Above 7 MeV, the number of $`\beta `$ counts drops dramatically and fitting error becomes large. In addition, with this method, as we determine the maximum $`\beta `$ energy and derives the energy distribution below that energy, it is difficult to estimate the errors on the high energy range . While the $`\overline{\nu }_e`$ spectra in Ref.  are given as tables, we use for simplicity somewhat less accurate analytical approximation given in Ref. . As a normalization factor we use energy-integrated $`\overline{\nu }_e`$ flux at Kamioka, $`1.34\times 10^6\mathrm{cm}^2\mathrm{s}^1`$, which are the summation of the flux from various nuclear reactors in Japan and Korea . This fit is not valid above $`8`$ MeV. However, there is an estimation that we would get few (10 or less) events per year above 10 MeV from reactors.

With three backgrounds we discussed above, we expect the energy window of SRN events opening from 10 MeV to 27 MeV. However, according to Kaplighat et al.  electrons or positrons from invisible muons are the largest background in the energy window from 19 to 35 MeV. This invisible muon event is illustrated as follows. The atmospheric neutrinos produce muons by interaction with the nucleons (both free and bound) in the fiducial volume. If these muons are produced with energies below Čherenkov radiation threshold (kinetic energy less than 53 MeV), then they will not be detected (“invisible muons”), but their decay-produced electrons and positrons will be. Since the muon decay signal will mimic the $`\overline{\nu }_epne^+`$ process in SK, it is difficult to distinguish SRN from these events. The energy spectrum of this invisible muon events is obtained by the stopped muon decay spectrum

$$\frac{dN}{dE_e}=\frac{G_F^2}{12\pi ^3}m_\mu ^2E_e^2\left(3\frac{4E_e}{m_\mu }\right),$$
(15)
(Michel spectrum ), where $`G_F`$ is the Fermi constant, $`m_\mu `$ the muon mass, and $`E_e`$ the electron (positron) energy. This equation is valid for $`E_e<m_\mu /2`$. From the observation of Kamiokande II detector the estimated event rate from these muon decays is around unity for 0.58 kton yr exposure and this is forming the principal source of background after the various cuts had been implemented . Kaplinghat et al.  used this value corrected for SK by multiplying the volume ratio and concluded that it is impossible to detect SRN unless their upper limit realizes and maximum neutrino oscillation occurs. This event rate is now inferred directly from the new data of the SK, which is about 100 per 1,258 days per 22.5 kton fiducial volume above 18 MeV . This value corresponds to about 20% of the Kamiokande II data. In Fig. 8 we show SRN event rate compared to the invisible muon events. From this figure we conclude that SRN events can be seen only below about 12 MeV.

In practice, there is another serious background, i.e., spallation products induced by cosmic ray muons. Ultra high energy cosmic ray muons spall oxygens in the detector, and radioactive decay processes of these spalled nuclei occur. The event rate of the spallation background is several hundred per day per 22.5 kton. Although most of them can be rejected by the information of preceding muons, even a small fraction can not be. (Roughly, this spallation products produces about 200,000 events per a year. Because expected SRN event rate is less than 1 per a year, we should reject all these 200,000 events. For future detectors this problem is also quite difficult to solve.) Then, this makes a serious background at the energy range below the maximum energy of beta spectrum of spallation products, 16 MeV . From these discussions presented in this subsection we conclude that there is no energy window of SRN.

### 5.2 Calculating the detectable event rate at SK

In the previous subsection we have found that there is no energy region where SRN is dominant. However, we can detect the SRN events by subtracting the other background events from total detected events. We consider the energy range $`17<(T_e/\mathrm{MeV})<25`$, where $`T_e`$ is positron kinetic energy. This range corresponds to $`19<(E_{\overline{\nu }_e}/\mathrm{MeV})<27`$ by the simple relation, $`E_{\overline{\nu }_e}=T_e+1.8\mathrm{MeV}`$. We find two advantages in using this energy region. First, SRN event rate is rather large, and second, the background (invisible muon) event rate is fairly well known by SK observation. We show in Table 3 the SRN event rate at SK in this energy range, i.e., $`0.40.8\mathrm{yr}^1`$. In contrast, the event rate of the invisible muon over the same energy range is $`3.4\mathrm{yr}^1`$. When SRN event rate is larger than the statistical error of background event rate, we can conclude that the SRN is detectable as a distortion of the expected invisible muon background event. Unfortunately, only one year observation does not provide any useful information about SRN. However, we can expect that ten-year observation provides several statistically meaningful results. The statistical error of invisible muon events in ten years is $`\sqrt{34}=5.8`$, which is smaller than the event rate of LMA models and is larger than that of SMA and no oscillation models. Then we conclude that these neutrino oscillation models can be distinguished by the observation of the event rate of invisible muon events. (If there is a discrepancy from expected event rate, this is due to SRN events and LMA models are favored.)

In future, it is expected that next generation of water Čherenkov detectors have much larger volume than that of SK. For example HyperKamiokande project is now under consideration. HyperKamiokande detector is planed to be a water Čherenkov detector whose mass is about 1,000,000 tons (about 20 times larger than SK), and its location is near SK detector. We expect that the SRN event becomes about 10 per one year for this detector, and statistically sufficient discussion of SRN is possible even using only one year data.

### 5.3 SNO and KamLAND detectors

In this subsection we discuss SRN detectability with SNO  and KamLAND . An advantage of these detectors is that we are able to identify $`\overline{\nu }_e`$ events using delayed coincidence signals. For this reason, we can remove other backgrounds from non-$`\overline{\nu }_e`$ origin (solar neutrinos, invisible muon decay products, and spallation products).

At SNO detector, used reaction for $`\overline{\nu }_e`$ detection is

$$\overline{\nu }_e+de^++n+n.$$
(16)
Then these neutrons react with surrounding nuclei through

$`n+d{}_{}{}^{3}\mathrm{H}+\gamma (E_\gamma =6.3\mathrm{MeV},\mathrm{efficiency}:24\%),`$ (17)
$`n+{}_{}{}^{35}\mathrm{Cl}{}_{}{}^{36}\mathrm{Cl}+\gamma (E_\gamma =8.6\mathrm{MeV},\mathrm{efficiency}:83\%),`$ (18)
where NaCl is added to efficiently capture neutrons. This delayed signal with the preceding Čherenkov radiation from $`e^+`$ in eq. (16) shows that the detected neutrino is $`\overline{\nu }_e`$. Using this criterion we can reject backgrounds from non-$`\overline{\nu }_e`$ origin. These are solar neutrinos, invisible muon decay products, and spallation products, which are great obstacle for SRN detection at SK. Thus, we can find the energy window below $`T_e23\mathrm{M}\mathrm{e}\mathrm{V}`$. Actual calculation shows that the expected event rate at SNO is $`0.03\mathrm{yr}^1`$. (We used cross section calculated in Ref. .)

We can use a similar criterion at KamLAND. $`\overline{\nu }_e`$’s are detected through below reactions:

$$\overline{\nu }_e+pe^++n,$$
(19)
$$n+pd+\gamma (E_\gamma =2.2\mathrm{MeV}).$$
(20)
The energy window for SRN is from 10 MeV to 25 MeV. (Below 10 MeV reactor $`\overline{\nu }_e`$’s are large.) In that range the calculated event rate is $`0.1\mathrm{yr}^1`$.

Unfortunately these values (0.03 yr-1 for SNO; 0.1 yr-1 for KamLAND) are quite small, since the fiducial volume of these detectors (1 kton) is much smaller than that of SK. However, the future same kind of detector of larger volume, if ever built, might detect SRN.

## 6 Acknowledgments

We would like to thank A. Suzuki, M. Nakahata, and Y. Fukuda for useful discussions and also would like to thank K. Takahashi for preparing the numerical data of the neutrino oscillation and for useful discussions. S. Ando also would like to thank S. Nagataki for useful discussions. This work was supported in part by Grants-in-Aid for Scientific Research provided by the Ministry of Education, Science and Culture of Japan through Research Grant No.07CE2002.