# Genetic Correlations in Mutation Processes

## I Introduction

Biological evolution is influenced by a number of processes including population growth, mutation, extinction, and interaction with the environment, to name a few . Genetic sequences are strongly affected by such processes and thus provide an important clue to their nature. The ongoing effort of reconstructing evolution histories given the incomplete set of mapped sequences constitutes much of our current understanding of biological evolution.

However, this challenge is extraordinary as it involves an inverse problem with an enormous number of degrees of freedom. Statistical methods such as maximum likelihood techniques coupled with simplifying assumptions on the nature of the evolution process are typically used to infer the structure of the underlying evolutionary tree, i.e., the phylogeny .

Genetic sequences such as RNA/DNA or amino acid sequences can be seen as words with letters taken from an alphabet of 4 or 20 symbols, respectively. Generally, there are nontrivial intra-sequence correlations that influence the evolution of the entire sequence. Additionally, the structure of the evolutionary tree plays a role in this process as one generally expects that the closer sequences are on this tree, the more correlated they are . In this study, we are interested in describing the influence of the latter aspect, namely the phylogeny, on the evolution of sequences. Specifically, we examine correlations between sequences, thereby complementing related studies on changes in fluctuations and entropy due to the phylogeny . To this end, we consider particularly simple sequences and focus on a model that mimics the competition between the fundamental processes of mutation and duplication.

The rest of this paper is organized as follows. In Sec. II, the model is introduced, and the main result is demonstrated using the pair correlations. Correlations of arbitrary order are obtained and analyzed asymptotically in Sec. III. To examine the range of validity of the results, generalizations to stochastic tree morphologies and sequences with larger alphabets are briefly discussed in Sections IV and V. Sec. VI discusses implications for multiple site correlations in sequences with independently evolving sites. We conclude with a summary and a discussion in Sec. VII.

## II Pair Correlations

Let us formulate the model first. The sequences are taken to be of unit length and the corresponding alphabet consists of two letters. The numeric values $`\sigma =\pm 1`$ are conveniently assigned to these letters. We will focus on binary trees where the number of children equals two. This structure is deterministic in that both the number of children and the generation lifetime are fixed. Nevertheless, the results apply qualitatively to stochastic tree morphologies as well. Finally, the mutation process is implemented as follows: with probability $`1p`$ a child equals his predecessor while with probability $`p`$ a mutation occurs, as illustrated in Fig. 1. The mutation process is invariant under the transformation $`\sigma \sigma `$ and $`p1p`$, and we restrict our attention to the case $`0p1/2`$ without loss of generality.

A natural question is how correlated are the various leafs (or nodes) of a tree in a given generation (or equivalently, time)? Consider $`G_2(k)`$ the average correlation between two nodes at the $`k`$th generation

$$G_2(k)=\sigma _i\sigma _j.$$
(1)   
The first average should be taken over all realizations for a fixed pair of nodes $`ij`$, while the second average is taken over all different pairs belonging to the same generation. For example, consider this quantity at the second generation (see Fig. 1), $`G_2(2)=[\sigma _3\sigma _4+\sigma _3\sigma _5+\sigma _3\sigma _6]/3`$. One index ($`i=3`$) may be fixed since all nodes in a given generation are equivalent.

To evaluate averages, it is useful to assign a multiplicative random variable $`\tau _i=\pm 1`$ to every branch of the tree such that $`\sigma _i=\sigma _j\tau _i`$ with $`j`$ the predecessor of $`i`$. One has $`\tau _i=1`$ ($`1`$) with probability $`1p`$ ($`p`$), and consequently,

$$\tau \tau _i=12p.$$
(2)   
Pair correlations are readily calculated using the $`\tau `$ variables: writing $`\sigma _3=\sigma _0\tau _1\tau _3`$ and similarly for $`\sigma _4`$ gives $`\sigma _3\sigma _4=\sigma _0\tau _1\tau _3\sigma _0\tau _1\tau _4=\sigma _0^2\tau _1^2\tau _3\tau _4`$. Since $`\sigma _i^2=\tau _i^2=1`$, this correlation simplifies, $`\sigma _3\sigma _4=\tau _3\tau _4`$. Furthermore, mutation processes on different branches are independent and consequently $`\tau _i\tau _j=\tau _i\tau _j`$ when $`ij`$. Thus, $`\sigma _3\sigma _4=\tau ^2`$ and similarly $`\sigma _3\sigma _5=\sigma _3\sigma _6=\tau ^4`$. The overall picture becomes clear: when calculating two-point correlations, the path to the tree root is traced for each node. As $`\tau ^2=1`$, doubly counted branches cancel. Only branches that trace the path to the first common ancestor are relevant. In other words,

$$\sigma _i\sigma _j=\tau ^{d_{i,j}}$$
(3)   
with $`d_{i,j}`$ the “genetic distance” between two points, the minimal number of branches that connect two nodes. Indeed, at the second generation $`d_{3,4}=2`$, $`d_{3,5}=d_{3,6}=4`$ and consequently $`G_2(2)=(\alpha ^2+2\alpha ^4)/3`$ with the shorthand notation $`\alpha =\tau =12p`$. This generalizes into a geometric series $`G_2(k)=(\alpha ^2+2\alpha ^4+\mathrm{}+2^{k1}\alpha ^{2k})/(2^k1)`$. Evaluating this sum gives the pair correlation

$$G_2(k)=\frac{\alpha ^2}{2\alpha ^21}\frac{(2\alpha ^2)^k1}{2^k1}.$$
(4)   
Interestingly, pair correlations are not affected by the initial state, i.e., the value of the tree root.

For sufficiently large generation numbers, the leading order of the pair correlation decays exponentially with the generation number. However, different constants characterize this decay, depending on the mutation probability

$$G_2(k)\{\begin{array}{cc}\frac{\alpha ^2}{2\alpha ^21}\alpha ^{2k}\hfill & p<p_c\text{;}\hfill \\ \frac{\alpha ^2}{12\alpha ^2}\mathrm{\hspace{0.17em}2}^k\hfill & p>p_c\text{.}\hfill \end{array}$$
(5)   
As seen from Eq. (4), the transition between the two different behaviors occurs when $`2\alpha ^2=1`$ or alternatively at the following mutation probability

$$p_c=\frac{1}{2}\left(1\frac{1}{\sqrt{2}}\right).$$
(6)   
Although in general correlations decay exponentially $`G_2(k)\beta ^{2k}`$, the decay constant $`\beta `$ exhibits two distinct behaviors which depend on the mutation probability $`\alpha `$. When the mutation probability is smaller than the critical one $`p<p_c`$ then $`\beta =\alpha `$ while in the complementary case $`\beta =1/\sqrt{2}`$.

As a reference, it is useful to consider the decay of the average node value $`G_1(k)=\sigma `$. At the $`k`$th generation, the path to each node involves $`k`$ branches and thus, $`G_1(k)=G_1(0)\alpha ^k`$ with $`G_1(0)=\sigma _0`$. Writing $`G_1(k)\beta ^k`$ then $`\beta =\alpha `$ for all mutation probabilities, in contrast with the asymptotic behavior of $`G_2(k)`$. Below the critical mutation rate, $`G_2(k)[G_1(k)/G_1(0)]^2`$, indicating that knowledge of the one-point average suffices to characterize correlations.

In fact, the above behavior can be attributed to the tree morphology. To see that, it is useful to consider a structureless morphology where the only ancestor shared by two nodes is the tree root itself (see Fig. 2). Using the notation $`G^{}`$ to denote correlations on this “star” morphology, we see that the average remains unchanged $`G_1(k)=G_1^{}(k)=G_1^{}(0)\alpha ^k`$. The star morphology is trivial in that all genetic distances are equal: $`d_{i,j}=2k`$ when $`ij`$. Thus, pair correlations are immediately obtained from the average $`G_2^{}(k)=[G_1^{}(k)/G_1^{}(0)]^2=\alpha ^{2k}`$. As branches in the star morphology do not interact, no correlations develop.

In contrast, nontrivial phylogenies do induce correlations. Indeed, $`G_2(k)>G_2^{}(k)`$ when $`p>0`$. Interestingly, when $`p<p_c`$, merely the asymptotic prefactor $`\alpha ^2/(2\alpha ^21)>1`$ in Eq. (5) is enhanced and $`G_2(k)G_2^{}(k)`$. As the critical point is approached, this constant diverges thereby signaling the transition into a second regime. When $`p>p_c`$, the decay constant itself is enhanced and the ratio $`G_2(k)/G_2^{}(k)`$ grows exponentially. The mutation probability affects only the asymptotic prefactor, and the decay constant $`\beta =1/\sqrt{2}`$ is determined by the tree morphology. We conclude that the nontrivial phylogeny generates significant correlations for larger than critical mutation probabilities.

This behavior can be understood and partially rederived using a heuristic argument. Genetically close nodes are highly correlated, while distant pairs are weakly correlated, as indicated by Eq. (3). On the other hand, distant pairs are more numerous. Both effects are magnified exponentially for large generation numbers, and their competition results in a critical point. Different mechanisms dominate on different sides of this point. Specifically, the number of minimal genetic distance pairs ($`d=2`$) is $`2^{k1}`$, while the number of maximal distance pairs ($`d=2k`$) is $`2^{2(k1)}`$. The rule (3) gives the relative contributions of these two terms to the overall two-point correlation: $`2^{k1}\alpha ^2`$ versus $`2^{2(k1)}\alpha ^{2k}`$. These are simply the first and last terms in the geometric series that led to Eq. (4). Comparing these two terms in the limit $`k\mathrm{}`$ correctly reproduces the most relevant aspects, i.e., the location of the critical point (6) and the decay constants of Eq. (5). We conclude that competition between the multiplicity and the degree of correlation of close and distant nodes underlies the transition.

## III Higher Order Correlations

The above analysis gives useful intuition for the overall qualitative behavior. Yet, it can be generalized into a more complete treatment that addresses correlations of arbitrary order. This set of quantities is helpful in determining the extent to which this picture applies, and in particular, whether the transition is actually a phase transition.

Multiple point correlations obey a rule similar to Eq. (3). For example, consider the four-node average $`\sigma _3\sigma _4\sigma _5\sigma _6`$ in Fig. 1. Using the $`\tau `$ variables, we rewrite $`\sigma _3\sigma _4\sigma _5\sigma _6=\sigma _0^4\tau _1^2\tau _2^2\tau _3\tau _4\tau _5\tau _6`$, and since $`\sigma ^2=\tau ^2=1`$ we get $`\sigma _3\sigma _4\sigma _5\sigma _6=\tau _3\tau _4\tau _5\tau _6=\tau ^4`$ or $`\sigma _3\sigma _4\sigma _5\sigma _6=\sigma _3\sigma _4\sigma _5\sigma _6`$. The four point average equals a product of two-point averages with the indices chosen as to minimize the total number of branches. This can also be seen by tracing the path of each node to the tree root and canceling doubly counted branches. Thus, Eq. (3) generalizes as follows:

$$\sigma _i\sigma _j\sigma _k\sigma _l=\tau ^{d_{i,j,k,l}},$$
(7)   
with the four-point genetic distance

$$d_{i,j,k,l}=\mathrm{min}\{d_{i,j}+d_{k,l},d_{i,k}+d_{j,l},d_{i,l}+d_{j,k}\}.$$
(8)   
Similarly, the law for arbitrary order averages is $`\tau `$ raised to a power equal to the $`n`$-point genetic distance. Such distance is obtained by considering all possible decompositions into pairs of nodes. The genetic distance is the minimal sum of the corresponding pair distances. Averages over an odd number of nodes can be obtained by adding a “pseudo” node at the root of the tree and using the convention $`d_{i,\mathrm{root}}=k`$ when $`i`$ belongs to the $`k`$th generation. The average $`\sigma _0`$ is generated by the root and this factor multiplies all odd order correlation. Since even order correlations are independent of the root value, and odd correlations are simply proportional to $`\sigma _0`$, we set $`\sigma _0=1`$ in what follows without loss of generality.

The average $`n`$-point correlation is defined as follows

$$G_n(k)=\sigma _{i_1}\sigma _{i_2}\mathrm{}\sigma _{i_n},$$
(9)   
where the averages are taken over all realizations and over all possible choices of $`n`$ distinct nodes at the $`k`$th generation. For the trivial star phylogeny, the $`n`$-point genetic distance is constant and equals a product of the correlation order and the generation number, $`d=nk`$. Consequently, all averages are trivial as knowledge of the one-point average immediately gives all higher-order averages, $`G_n^{}(k)=[G_1^{}(k)]^n`$, or explicitly

$$G_n^{}(k)=\alpha ^{nk}.$$
(10)   
When the tree morphology is nontrivial, the minimal-sum rules (7)-(8) imply that such factorization no longer holds. For binary trees, it is possible to obtain these correlations recursively. Let us assign the indices $`1,2,\mathrm{},2^k`$ to the $`k`$th generation nodes and order them as follows $`1i_1<i_2<\mathrm{}<i_n2^k`$. As the average over the realizations is performed first, the average correlation requires a summation over all possible choices of nodes

$$F_n(k)=\underset{1i_1<i_2<\mathrm{}<i_n2^k}{}\sigma _{i_1}\sigma _{i_2}\mathrm{}\sigma _{i_n}.$$
(11)   
Proper normalization gives the $`n`$-node correlation

$$G_n(k)=F_n(k)/\left(\genfrac{}{}{0pt}{}{2^k}{n}\right).$$
(12)   
Consider a group of $`n`$ nodes taken from the $`k`$th generation. They all share the tree root as a common ancestor. The two first generation nodes naturally divide this group into two independently evolving subgroups. This partitioning procedure allows a recursive calculation of the correlations. Formally, a given choice of nodes $`1i_1<i_2<\mathrm{}<i_n2^k`$ is partitioned into two subgroups as follows $`1i_1<\mathrm{}<i_m2^{k1}`$ and $`2^{k1}+1i_{m+1}<\mathrm{}<i_n2^{k1}+2^{k1}`$. These subgroups involve different $`\tau `$ variables, so their correlations factorize

$$\sigma _{i_1}\mathrm{}\sigma _{i_n}\sigma _{i_1}\mathrm{}\sigma _{i_m}\sigma _{i_{m+1}}\mathrm{}\sigma _{i_n}.$$
(13)   
The proportionality constant depends upon the parity of $`m`$ and $`nm`$. Even correlations are independent of the tree root, while odd correlations are proportional to the average value of the tree root. This extends to sub-trees as well, and since $`\sigma _0=1`$, the average value of the root of both sub-trees is $`\tau `$. This factor accompanies all odd correlations. Substituting Eq. (13) into Eq. (11) shows that the summation factorizes as well. Using $`F_m(k1)=_{1i_1<\mathrm{}<i_m2^{k1}}\sigma _{i_1}\mathrm{}\sigma _{i_m}`$ reduces the problem to two sub-trees that are one generation shorter, and a recursion relation for $`F_n(k)`$ emerges

$$F_n(k)=\underset{m=0}{\overset{n}{}}F_m(k1)B_mF_{nm}(k1)B_{nm},$$
(14)   
with the boundary conditions $`F_n(0)=\delta _{n,0}+\delta _{n,1}`$. The summation corresponds to the $`n+1`$ possible partitions of a group of $`n`$ nodes into two subgroups. The weight of the odd correlations is accounted for by $`B_n`$

$$B_n=\{\begin{array}{cc}1\hfill & n=2r\text{;}\hfill \\ \tau \hfill & n=2r+1\text{.}\hfill \end{array}$$
(15)   
Using the definition (11), the sums $`F_n(k)`$ vanish whenever $`n>2^k`$. This behavior emerges from the recursion relations as well. Additionally, one can check that the sums are properly normalized in the no mutation case ($`\alpha =1`$), $`F_n(k)=\left(\genfrac{}{}{0pt}{}{2^k}{n}\right)`$ when $`n2^k`$.

For sufficiently small $`n`$, it is possible to evaluate the sums explicitly using Eqs. (14). The average correlations are then found using Eq. (12)

$`G_0(k)`$ $`=`$ $`1,`$ (16)   
$`G_1(k)`$ $`=`$ $`\alpha ^k,`$ (17)   
$`G_2(k)`$ $`=`$ $`{\displaystyle \frac{\alpha ^2}{2\alpha ^21}}{\displaystyle \frac{(2\alpha ^2)^k1}{2^k1}},`$ (18)   
$`G_3(k)`$ $`=`$ $`{\displaystyle \frac{3\alpha ^{k+2}}{2\alpha ^21}}{\displaystyle \frac{\frac{(4\alpha ^2)^k(4\alpha ^2)}{4\alpha ^21}(2^k2)}{(2^k1)(2^k2)}}.`$ (19)   
Indeed, these quantities agree with the previous results for $`n=1`$, $`2`$ and equal unity when $`p=0`$. We see that correlations involve a sum of exponentials. Furthermore, it appears that the condition $`2\alpha ^2=1`$ still separates two different regimes of behaviors. However, calculating higher correlations explicitly is not feasible as the expressions are involved for large $`n`$. Instead, we perform an asymptotic analysis that more clearly exposes the leading large generation number behavior.

Let us consider first the regime $`p<p_c`$ or equivalently $`2\alpha ^2>1`$. From Eq. (16), we see that the leading large $`k`$ behavior of the average correlation satisfies $`G_n(k)\alpha ^{nk}`$ for $`n=0,1,2`$, and $`3`$. We will show below that this behavior extends to higher order correlations, i.e.,

$$G_n(k)g_n\alpha ^{nk}.$$
(20)   
In other words, the following limit $`\alpha =lim_k\mathrm{}[G_n(k)]^{1/nk}`$ exists and is independent of $`n`$. As correlations are larger when the phylogeny is nontrivial, one expects that $`G_n(k)G_n^{}(k)`$ or in terms of the prefactors, $`g_ng_n^{}=1`$. Combining Eq. (12) with the leading behavior of the combinatorial normalization constant $`\left(\genfrac{}{}{0pt}{}{2^k}{n}\right)2^{nk}/n!`$ gives the the asymptotic behavior of the sums

$$F_n(k)=f_n(2\alpha )^{nk},\mathrm{with}f_n=\frac{g_n}{n!}.$$
(21)   
Substituting Eq. (21) into the recursion relation Eq. (14) eliminates the dependence on the generation number $`k`$, and a recursion relation for coefficients $`f_n`$ is found

$$f_n(2\alpha )^n=\underset{m=0}{\overset{r}{}}f_mB_mf_{nm}B_{nm},$$
(22)   
with $`B_n`$ of Eq. (15). These recursion relations are consistent with the conditions $`f_0=f_1=1`$. The case $`n=2`$ reproduces the coefficient $`f_2=\alpha ^2/[(2\alpha )^22]`$. The divergence at $`2\alpha ^2=1`$ indicates that the ansatz (20) breaks down at the critical point. To show that the ansatz holds in the entire range $`0p<p_c`$, one has to show that the coefficients $`f_n`$ are positive and finite for all $`n`$. Rewriting the recursion (22) explicitly $`f_n[(2\alpha )^n2B_n]=_{m=1}^{r1}f_mB_mf_{nm}B_{nm}`$ allows us to prove this. Since $`f_0=1>0`$, then to complete a proof by induction one needs to show that a positive $`f_{n1}`$ implies a positive $`f_n`$. The right hand side of the recursion is clearly positive and thus the positivity of $`f_n`$ hinges on the positivity of the term $`(2\alpha )^n2B_n`$. When $`2\alpha ^2>1`$, then $`\alpha >1/\sqrt{2}`$ and certainly $`2\alpha >1`$. Combining this with the inequality $`(2\alpha )^2>2>2B_n`$ shows that $`(2\alpha )^n2B_n>0`$ when $`n2`$. Hence $`f_n`$ is positive and finite for all $`n`$, which validates the ansatz (20) in the regime $`p<p_c`$.

In principle, the coefficients can be found by introducing the generating functions

$$f(z)=\underset{n}{}f_nz^n.$$
(23)   
Multiplying Eq. (22) by $`z^n`$ and summing over $`n`$ yields the following equation for the generating functions

$$f(2\alpha z)=\left[\frac{f(z)+f(z)}{2}+\alpha \frac{f(z)f(z)}{2}\right]^2.$$
(24)   
This equation reflects the structure of the recursion relations. A factor $`\alpha `$ is generated by each odd-index coefficient and as a results, the odd part of the generating functions $`[f(z)f(z)]/2=f_1z+f_3z^3+\mathrm{}`$ is multiplied by $`\alpha `$. Although a general solution of this equation appears rather difficult, it is still possible to obtain results in the limiting cases. It is useful to check that when $`\alpha =1`$, the above equation reads $`f(2z)=f^2(z)`$ which together with the boundary conditions $`f_0=f_1=1`$ gives $`f(z)=\mathrm{exp}(z)`$ or $`f_n=\frac{1}{n!}`$. As $`g_n1`$, the trivial correlations are recovered, $`G_nG_n^{}`$ indicating that role played by the tree morphology diminishes in the no mutation limit.

In the limit $`pp_c^{}`$ it is possible to extract the leading behavior of the asymptotic prefactors. Here, it is sufficient to keep only the highest powers of the diverging term $`1/(2\alpha ^21)`$. The calculation in this case is identical to the one detailed below for the case $`p>p_c`$ and we simply quote the results

$$G_n(k)\{\begin{array}{cc}\frac{2r!}{r!}\left[\frac{\alpha ^2}{2(2\alpha ^21)}\right]^r\alpha ^{nk}\hfill & n=2r\text{;}\hfill \\ \frac{(2r+1)!}{r!}\left[\frac{\alpha ^2}{2(2\alpha ^21)}\right]^r\alpha ^{nk}\hfill & n=2r+1\text{.}\hfill \end{array}$$
(25)   
In this limit, the odd order correlations simply follow from their even counterparts and for example $`f_{2r+1}=f_{2r}`$.

In the complementary case $`p>p_c`$, it proves useful to rewrite the recursion relations (22) for the even and odd correlations separately

$`F_{2r}(k)`$ $`=`$ $`{\displaystyle \underset{s=0}{\overset{r}{}}}F_{2s}(k1)F_{2r2s}(k1)`$ (26)   
$`+`$ $`\alpha ^2{\displaystyle \underset{s=0}{\overset{r1}{}}}F_{2s+1}(k1)F_{2r2s1}(k1)`$ (27)   
$`F_{2r+1}(k)`$ $`=`$ $`2\alpha {\displaystyle \underset{s=0}{\overset{r}{}}}F_{2s}(k1)F_{2r2s+1}(k1).`$ (28)   
The leading asymptotic behavior of Eq. (16) implies $`F_0(k)=f_0`$, $`F_1(k)f_0(2\alpha )^k`$, $`F_2(k)f_22^k`$, and $`F_3(k)f_22^k(2\alpha )^k`$ with $`f_0=1`$ and $`f_2=\alpha ^2/[2(2\alpha )^2]`$. Let us assume that this even-odd pattern is general

$`F_{2r}(k)`$ $`=`$ $`f_{2r}2^{rk},`$ (29)   
$`F_{2r+1}(k)`$ $`=`$ $`f_{2r}2^{rk}(2\alpha )^k.`$ (30)   
Substituting this ansatz into Eq. (26) shows that the second summation in the recursion for the even correlations is negligible asymptotically. Both equations reduce to

$$f_{2r}2^r=\underset{s=0}{\overset{r}{}}f_{2s}f_{2r2s},$$
(31)   
and therefore the pattern (29) holds when $`p>p_c`$. It is seen that odd correlators are enslaved to the even ones.

To obtain the coefficients, we introduce the generating functions $`f(z)=_rf_{2r}z^{2r}`$ which satisfies $`f(0)=1`$, $`f^{}(0)=0`$ and $`f^{\prime \prime }(0)=f_2=\alpha ^2/[2(12\alpha ^2)]`$. The recursion relation translates into the following equation for $`f(z)`$

$$f\left(\sqrt{2}z\right)=[f(z)]^2.$$
(32)   
Its solution is $`f(z)=\mathrm{exp}[(\alpha z)^2/2(12\alpha ^2)]`$. Thus, $`f_{2r}=\frac{1}{r!}[f_2]^r`$. From Eqs. (20)-(21), the leading asymptotic behavior in the regime $`p_c<p<1/2`$ is found

$$G_n(k)\{\begin{array}{cc}\frac{2r!}{r!}\left[\frac{\alpha ^2}{2(12\alpha ^2)}\right]^r2^{kr}\hfill & n=2r\text{;}\hfill \\ \frac{(2r+1)!}{r!}\left[\frac{\alpha ^2}{2(12\alpha ^2)}\right]^r\alpha ^k2^{kr}\hfill & n=2r+1\text{.}\hfill \end{array}$$
(33)   
Using the Stirling formula $`n!\sqrt{2\pi n}n^ne^n`$ it is seen that the coefficients $`g_{2r}`$ have nontrivial $`r`$ behavior as $`g_{2r}=g_{2r+1}/(2r+1)\sqrt{2}[2\alpha ^2/(12\alpha ^2)]^rr^r`$.

The even order correlations have identical asymptotic behavior to the two point correlation: $`lim_k\mathrm{}[G_{2r}(k)]^{1/2rk}=\frac{1}{\sqrt{2}}`$ for all $`r`$. The odd order correlations behave differently, however, as this limit depends on the correlation order: $`lim_k\mathrm{}[G_{2r+1}(k)]^{1/(2r+1)k}=\frac{1}{\sqrt{2}}(\sqrt{2}\alpha )^{1/2r+1}`$. Thus, only in the limit $`r\mathrm{}`$ do the even and odd order correlations agree. However, this conclusion is misleading since the decay rate of the (properly normalized) odd order correlations $`G_{2r+1}(k)/G_1(k)G_{2r}(k)`$ is identical to that of the even order correlations. We conclude that the decay rate of two-point correlations characterizes the decay of all higher order correlations.

From Eqs. (25) and (33), we see that the coefficients diverge according to

$$f_{2r}=f_{2r+1}|p_cp|^r$$
(34)   
as the critical point is approached, $`pp_c`$. Since the correlations must remain finite, this indicates that the purely exponential behavior must be modified when $`p=p_c`$. Indeed, evaluating Eq. (16) at $`p=p_c`$ yields $`F_2(k)f_22^k`$ and $`F_3(k)=f_22^{3k/2}`$ with $`f_2=k/4`$, i.e., the even-odd pattern of Eq. (29) is reproduced. Furthermore, the value of $`f_2`$ shows that the diverging quantity $`1/|12\alpha ^2|`$ is simply replaced by $`k`$. This implies that the coefficients become generation dependent, $`f_nf_n(k)`$. Assuming the pattern Eq. (29), substituting it into Eq. (25), and following the steps that led to Eq. (33) yields the critical behavior

$$G_n(k)\{\begin{array}{cc}\frac{2r!}{r!}\left[\frac{k}{4}\right]^r2^{kr}\hfill & n=2r\text{;}\hfill \\ \frac{(2r+1)!}{r!}\left[\frac{k}{4}\right]^r2^{k(r+1/2)}\hfill & n=2r+1\text{.}\hfill \end{array}$$
(35)   
Generally, the diverging quantity $`1/|12\alpha ^2|`$ is replaced with the finite (but ever growing) quantity $`k`$. The algebraic modification to the leading exponential behavior in Eq. (35) is reminiscent of the logarithmic corrections that typically characterize critical behavior in second order phase transitions .

## IV Stochastic Tree Morphologies

The question arises: how general is the behavior described above? The binary tree considered was particularly simple as it involved a fixed number of children and a fixed generation lifetime. Below we show that relaxing either of these conditions does not affect the nature of the results.

Let us first consider tree morphologies with a varying number of children, i.e., the trees are generated by a stochastic branching process where with probability $`P_r`$ there are $`r`$ children. This probability sums to unity $`_rP_r=1`$, and the average number of children is given by $`r=_rrP_r`$. As a result, the average number of nodes at the $`k`$th generation is $`r^k`$, indicating that the tree “survives” only if $`r>1`$, a classical result of branching processes theory . The rule (3) is independent of the tree morphology, and therefore, one can repeat the heuristic argument in Sec. II. The extreme contributions to the average pair correlations have the relative weights $`r^{k1}\alpha ^2`$ and $`r^{2(k1)}\alpha ^{2k}`$. Comparing these two terms asymptotically shows that the critical point is a simple generalization of Eq. (6)

$$p_c=\frac{1}{2}\left(1\sqrt{\frac{1}{r}}\right).$$
(36)   
The critical mutation rate varies from $`0`$ to $`1/2`$ as the average ancestry size varies between $`1`$ and $`\mathrm{}`$. This indicates that correlations are significant over a larger range of mutation rates for smaller trees. The heuristic argument also gives the decay constant $`\beta `$, and the leading asymptotic behavior of Eq. (5) is generalized by simply replacing $`2`$ with $`r`$. A more complete treatment of this problem is actually possible and closely follows Eq. (4). Again, the ancestry size $`r`$ replaces the deterministic value $`2`$. As both the results and the overall behavior closely follow the deterministic case, we do not detail them here.

A second possible generalization is to morphologies with a varying generation lifetime. Such tree morphologies can be realized by considering a continuous time variable. Branching is assumed to occur with a constant rate $`\nu `$. For such tree morphologies, the number of nodes $`n(t)`$ obeys $`\dot{n}(t)=\nu n(t)`$ which gives an exponential growth $`n(t)=e^{\nu t}`$. Similarly, the mutation process is assumed to occur with a constant rate $`\gamma `$. A useful characteristic of this process is the autocorrelation $`A(t)=\sigma (0)\sigma (t)`$. To evaluate it’s evolution, we note that $`A(t+dt)=(1\gamma dt)A(t)\gamma dtA(t)`$ when $`dt0`$. Therefore, $`\dot{A}(t)=2\gamma A(t)`$ and one finds $`A(t)=e^{2\gamma t}`$. The quantities $`n(t)`$ and $`A(t)`$ allow calculation of the average pair correlation.

Let us pick two nodes at time $`t`$ and denote their values by $`\sigma _i(t)`$ and $`\sigma _j(t)`$, and let the genetic distance between these two nodes be $`\tau `$. Using their first common ancestor $`\sigma _c(t\tau )=\sigma _i(t\tau )=\sigma _j(t\tau )`$ and the identity $`\sigma ^2=1`$, their correlation can be evaluated as follows $`\sigma _i(t)\sigma _j(t)=\sigma _i(t)\sigma _c(t\tau )\sigma _c(t\tau )\sigma _j(t)=\sigma _i(t)\sigma _i(t\tau )\sigma _j(t)\sigma _j(t\tau )=A^2(\tau )`$. Integrating over all possible genetic distances gives the average pair correlation

$$G_2(t)=\frac{_0^t𝑑\tau n(\tau )A^2(\tau )}{_0^t𝑑\tau n(\tau )}.$$
(37)   
The factor $`n(\tau )/_0^t𝑑\tau n(\tau )`$ accounts for the multiplicity of pairs with genetic distance $`\tau `$. Using $`A(t)=e^{4\gamma t}`$ and $`n(\tau )=e^{\nu t}`$, the average pair correlation is evaluated

$$G_2(t)=\frac{\nu }{\nu 4\gamma }\frac{e^{(\nu 4\gamma )t}1}{e^{\nu t}1}.$$
(38)   
For the star phylogeny the genetic distance is always $`t`$ and therefore $`G_2^{}(t)=e^{4\gamma t}`$. Here the relevant parameter is the normalized mutation rate $`\omega =\gamma /\nu `$. Again, there exists a critical point $`\omega _c=1/4`$. For smaller than critical mutation rates, $`\omega <\omega _c`$, correlations due to the tree morphology are not pronounced, $`G_2(t)G_2^{}(t)`$. On the other hand, when $`\omega >\omega _c`$, strong correlations are generated and $`G_2(k)e^{\nu t}`$ is exponentially larger than $`G_2^{}(t)`$. We conclude that the behavior found for the deterministic case is robust.

## V Multistate Sequences

We now consider larger alphabets. Previously, the two states satisfied $`\sigma ^2=1`$. A natural generalization is to $`\sigma ^n=1`$, i.e, the $`n`$th order roots of unity $`\sigma =e^{i2\pi l/n}`$ with $`l=0,1,\mathrm{},n1`$. Previously, with probability $`p`$ the mutation $`\sigma \tau \sigma `$ occurred with $`\tau =e^{i\theta }`$ and $`\theta =\pi `$. We thus impose the same transition but with $`\theta =2\pi /n`$. This can be viewed as a clockwise rotation in the complex plane by an angle $`\theta `$. Since the states are now complex, the definition of the pair correlation is now

$$G_2(k)=\overline{\sigma }_i\sigma _j,$$
(39)   
with $`\overline{\sigma }`$ the complex conjugate of $`\sigma `$. The real part of $`\overline{\sigma }_i\sigma _j`$ gives the inner product of the two-dimensional vectors corresponding to $`\sigma _i`$ and $`\sigma _j`$, respectively.

Consider the average $`\overline{\sigma }_3\sigma _4`$ in Fig. 1. Using the $`\tau `$ variables and $`\overline{\tau }\tau =\overline{\sigma }\sigma =1`$ one has $`\overline{\sigma }_3\sigma _4=\overline{\sigma }_0\overline{\tau }_1\overline{\tau }_3\sigma _0\tau _1\tau _4=\overline{\tau }_3\tau _4=\overline{\tau }_3\tau _3=\overline{\tau }\tau =|\tau |^2`$. All of our previous results hold if one replaces the average $`\tau `$ with its magnitude $`\alpha =|\tau |=|1p(1e^{i\theta })|=\sqrt{12p(1p)(1\mathrm{cos}\theta )}`$. Furthermore, it is sensible to consider arbitrary phase shifts $`0<\theta <2\pi `$ since the identity $`\overline{\sigma }\sigma =\overline{\tau }\tau =1`$ rather than $`\sigma ^n=\tau ^n=1`$ was used to evaluate correlations.

The critical point is determined from the condition $`r\alpha ^2=1`$. This equation has a physical solution only when $`2\phi <\theta <2(\pi \phi )`$ with the shorthand notation

$$\phi =\mathrm{cos}^1\sqrt{\frac{1}{r}}.$$
(40)   
In terms of the number of states, this translates to

$$\frac{\pi }{\pi \phi }<n<\frac{\pi }{\phi }.$$
(41)   
Hence, the transition may or may not exist depending on the details of the model, in this particular “clock” model case, the number of states. As we have seen before, correlations become less pronounced when the number of ancestors increases. Indeed, the transition always exists in the limit $`r1`$, while the transition is eliminated in the other extreme $`r\mathrm{}`$. When the transition do occur, the following critical mutation probability is found

$$p_c=\frac{1}{2}\left[1\sqrt{1\frac{\mathrm{sin}^2\phi }{\mathrm{sin}^2\frac{\theta }{2}}}\right].$$
(42)   
Indeed, Eq. (36) is reproduced in the two state case ($`\theta =\pi `$). This turns out to be the minimal critical point, $`p_c(1\sqrt{1/r})/2`$, reflecting the fact that transition $`\sigma \sigma `$ provides the most effective mutation mechanism.

Interestingly the transition is restored when both the mutation and the duplication processes occur continuously in time. In this continuous description duplication occurs with rate $`\nu `$ and the mutation $`\sigma e^{i\theta }\sigma `$ occurs with rate $`\gamma `$. The autocorrelation $`A(t)=\sigma (0)\overline{\sigma }(t)=\mathrm{exp}\left[\gamma (1e^{i\theta })t\right]`$ is found from its time evolution $`\dot{A}(t)=\gamma (1e^{i\theta })A(t)`$. It can be easily shown from the definition of the pair correlation (39) that $`A^2(\tau )`$ should be replace with $`|A(\tau )|^2=\mathrm{exp}\left[2\gamma (1\mathrm{cos}\theta )\tau \right]`$ in the integral (37). Comparing with the previous section results, we see that the effective mutation rate is now $`\gamma (1\mathrm{cos}\theta )/2`$. As a result, the location of the critical point is increased by a factor $`2/(1\mathrm{cos}\theta )`$. Using the normalized mutation rate $`\omega =\gamma /\nu `$ one finds

$$\omega _c=\frac{1}{2(1\mathrm{cos}\theta )}.$$
(43)   
This critical point increases with the number of states, and it diverges according to $`\omega _c(n/2\pi )^2`$ when $`n\mathrm{}`$. This behavior is intuitive as one expects that mutations between a large number of states diminishes correlations and consequently, phylogenetic effects.

## VI Two-site Correlations

When sequences are not of unit length, i.e., when there are two or more sites per sequence, the results can be used to characterize a correlation measure quantifying the interaction between sites. Assume there are two or more sites per sequence and that the sites evolve independently of each other. Denote the state of position $`a`$ in sequence $`i`$ as $`\sigma _i^a`$, and similarly denote the state of position $`b`$ in sequence $`i`$ as $`\sigma _i^b`$. If the sequences were not related by a phylogenetic tree, but instead were independent samples drawn from a given distribution, then the following quantity defined on a finite set of $`N=2^k`$ samples specifies a two-site correlation measure:

$$\rho =\frac{1}{N}\underset{i}{}\sigma _i^a\sigma _i^b\frac{1}{N^2}\underset{i}{}\sigma _i^a\underset{j}{}\sigma _j^b.$$
(44)   
Correlation between sites $`a`$ and $`b`$ is indicated by a non-zero value of $`\rho `$.

The quantity $`\rho `$ is well defined also when the sequences are related by a phylogenetic tree. Due to the assumption of independent positions, the mean of $`\rho `$ over all realizations vanishes $`\rho =0`$. This behavior is independent of the tree morphology. To see the effects of the phylogeny, one needs to consider fluctuations, i.e., the variance $`\mathrm{\Delta }\rho =\rho ^2`$

$`\mathrm{\Delta }\rho `$ $`=`$ $`\left({\displaystyle \frac{1}{N}}{\displaystyle \underset{i}{}}\sigma _i^a\sigma _i^b\right)^2\left({\displaystyle \frac{1}{N}}{\displaystyle \underset{i}{}}\sigma _i^a\right)^2\left({\displaystyle \frac{1}{N}}{\displaystyle \underset{i}{}}\sigma _i^b\right)^2`$ (45)   
$`=`$ $`{\displaystyle \frac{1}{N^2}}{\displaystyle \underset{ij}{}}\sigma _i^a\sigma _j^a\sigma _i^b\sigma _j^b{\displaystyle \frac{1}{N^4}}{\displaystyle \underset{ij}{}}\sigma _i^a\sigma _j^a{\displaystyle \underset{kl}{}}\sigma _k^b\sigma _l^b`$ (46)   
$`=`$ $`{\displaystyle \frac{1}{N^2}}\left[N+{\displaystyle \underset{ij}{}}\tau ^{2d_{i,j}}\right]{\displaystyle \frac{1}{N^4}}\left[N+{\displaystyle \underset{ij}{}}\tau ^{d_{i,j}}\right]^2.`$ (47)   
The first equality in the above equation was obtained by rewriting Eq. (44) as $`\rho =\rho _1\rho _2`$ and noting that $`\rho _1\rho _2=\rho _2^2`$. The final expression can be simplified using $`_{ij}\tau ^{2d_{i,j}}=N(N1)G_2(\alpha ^2,k)`$ with $`G_2(\alpha ^2,k)`$ the pair correlation of Eq. (4), considered as a function of $`\alpha ^2`$. The following expression for the variance is obtained

$$\mathrm{\Delta }\rho =\left[\frac{1}{N}+\left(1\frac{1}{N}\right)G_2(\alpha ^2,k)\right]\left[\frac{1}{N}+\left(1\frac{1}{N}\right)G_2(\alpha ,k)\right]^2.$$
(48)   
For the star morphology the leading order of the fluctuations is independent of the mutation rate and it scales as the familiar $`N^1`$. For the binary tree morphology, there are again two regimes, characterized by $`p>p_c`$ or $`p<p_c`$ where $`p_c`$ is now defined by $`2\alpha ^4=1`$, i.e.,

$$p_c=\frac{1}{2}\left(1\left(\frac{1}{2}\right)^{\frac{1}{4}}\right).$$
(49)   
When $`p<p_c`$, the phylogeny plays a significant role and the variance is exponentially enhanced $`\mathrm{\Delta }\rho \alpha ^{4k}`$, while when $`p>p_c`$, the variance is still statistical in nature $`\mathrm{\Delta }\rho A\mathrm{\Delta }^{}\rho `$ with $`A>1`$. Hence, it is more likely to observe large values of $`\rho `$ in the tree morphology than it is in the star morphology, even when the sites evolve independently. Since correlations and variance play opposite roles, they are influenced in different ways by the phylogeny.

## VII Summary

In summary, we have studied the influence of the phylogeny on correlations between the tree’s nodes. In general, for sufficiently small mutation rates, the morphology plays a minor role. For sufficiently high mutation rates large correlations that can be attributed to the phylogeny may occur. The transition between the two regimes of behavior is sharp and is marked by a critical mutation rate. Below this critical point all correlations are well described by the average, while above it, correlations decay much slower than the average. Underlying this transition is the competition between the multiplicity and the the degree of correlations between genetically close and distant leafs. This competition also leads to larger fluctuations in the correlation between different sites, even when these evolve independently.

We have also seen that this behavior is robust and appears to be independent of many details of the model. While the overall behavior generally holds, specific details such as the location of the critical point and the decay rate in the regime $`p>p_c`$ depend on a specific tree dependent parameter: the average number of children.

The above results can be extended in several directions. It will be interesting to see whether the recursive methods can be generalized to stochastic tree morphologies and in particular to the continuous time case. This methods should still be applicable even when the mutation rates are time dependent or disordered. In such cases it will be interesting to determine which parameters determine the critical point, the decay constants, etc.

Correlations can serve as useful measure of the diversity of a system since small correlations indicate large diversity and vice versa. If the diversity can be measured in an experiment where the phylogeny is controlled, its time dependence can be used to infer the mutation probability. Similarly, if the mutation probability can be controlled, than the degree of correlation/diversity can be used to infer characteristics of the phylogeny. Thus, our results may be useful for inferring statistical properties of actual biological systems.

This research is supported by the Department of Energy under contract W-7405-ENG-36.