Bahram Houchmandzadeh

A general method for the estimation of moments of substitution number in molecular evolution

The random variable substitution number Nt, i.e. the number of mutations that have accumulated in a sequence under neutral evolution during a time t, is a cornerstone of the eld of molecular evolution. We show here that a complete solution of the moments of this random variable can be obtained explicitly by recurrence, using a simple matrix shift method. This result leads to simple expressions for both the short and long time limits of the moments that can be computed directly from the substitution matrix used to model the neutral evolution. The method developed here is also used to compute the moments of the complementary variable, Tn, the time it takes to accumulate n mutations. The method we develop here necessitates only elementary operations on the substitution matrix and does not involve spectral decomposition.

I. INTRODUCTION.

Substitution models are the cornerstone of the eld of molecular evolution and describe the Markovian stochastic process of genetic variation through xed mutations.

One of their most important application is to evaluate the length of the branches in a phylogeny tree, and hence the time of divergence from an ancestor, trough the estimation of the number of substitutions N that has occurred between two sequences. A substitution model describes the rate of mutation between symbols of an alphabet, that can describe for example base pairs in DNA (4 symbols), codons (61 symbols, excluding the STOP), amino acids (20 symbols), or more complicated alphabets taking into account various intricacy of the evolutionary process under consideration. Since their inception by Jukes and Cantor[1] and Kimura and Ohta [START_REF] Kimura | On the stochastic model for estimation of mutational distance between homologous proteins[END_REF][START_REF] Kimura | A simple method for estimating evolutionary rates of base substitutions through comparative studies of nucleotide sequences[END_REF], many substitution models have been proposed, documented in textbooks [START_REF] Graur | Fundamentals of Molecular Evolution[END_REF][START_REF] Yang | Molecular Evolution: A Statistical Approach[END_REF] and incorporated into phylogenetic software such as IQ-Tree [START_REF] Nguyen | IQ-TREE: A Fast and Effective Stochastic Algorithm for Estimating Maximum-Likelihood Phylogenies[END_REF]. A review of these models and advances in this eld can be found in [START_REF] Arenas | Trends in substitution models of molecular evolution[END_REF].

Once a substitution model Q has been chosen, it can be used to characterize the random variable N t of the number of substitutions that have occurred during a time t. It is straightforward to deduce the mean number of substitutions ⟨N t ⟩ from the substitution matrix Q. However, the mean is only the rst step to characterize a random variable and one need its higher moments (variance, skewness, kurtosis, ...) to set the condence limits on this estimation. The task of evaluating the higher moments has proved extremely dicult and for example has led to the much discussed controversy of the overdispersion of the molecular clock [811].

Zheng [START_REF] Zheng | On the dispersion index of a Markovian molecular clock[END_REF] was the rst to use Markov chains to investigate the variance of substitution number as a solution of a set of dierential equations. His investigation was further developed by Bloom, Raval and Wilke [START_REF] Bloom | Thermodynamics of Neutral Protein Evolution[END_REF] who gave the general solution in terms of the spectral decomposition of the substitution matrix; this solution was extended by Raval [START_REF] Raval | Molecular Clock on a Neutral Network[END_REF] for a specic class of matrices used for random walk on neutral graphs. Minin and Suchard [START_REF] Vladimir | Counting labeled transitions in continuous-time Markov models of evolution[END_REF] used the same spectral method to derive an analytical form for the generating function of a binary process. Note that the substitution matrix is singular and asymmetric, its eigenvalues are in general complex and the spectral decomposition cannot in general be obtained exactly. We recently published [START_REF] Houchmandzadeh | A Simple, General Result for the Variance of Substitution Number in Molecular Evolution[END_REF] a simple, exact result for the variance of the substitution number and showed, as a consequence, that for the large class of Generalized Time Reversible matrices, the dispersion index, i.e. the variance to mean ratio is always larger than one. Although the result for the variance was simple, the method to obtain it was intricate and could not be generalized to compute the higher moments of the substitution number.

In the present article, we show that there exists a simple transformation of the substitution matrix that decomposes the problem into a series of simple equations that can be solved by elementary matrix algebra and allows for the computation of moments of arbitrary order.

The same method allows us also to compute directly the moments of the complementary random variable T n , the time it takes to reach n substitutions. This provides for a more direct approach to estimate the divergence time from an ancestor. This article is organized as follows. In the next section, we highlight the main results of this article and their comparison to numerical simulations of the stochastic process. Section III then details the methods used to compute the moments of N t . Section IV address the complementary method of computing the variance of T n . Section V is devoted to discussions and conclusions. The appendices contain details of some of the computations.

II. MAIN RESULTS

A. Denitions and notations.

The problem of substitution number N can be formulated as follows. A discrete nite Markov process X has K possible values and is described by substitutions rates Each time a transition occurs, the random variables (N, T ) are incremented by (1, ti) where ti is an exponential time depending on the state that is left. (b) A realization of this process can be represented by a random path in the (N, T ) space (solid lines for 200 such random paths, simulated by a Gillespie algorithm). Given a time t, one could then study the probability P (n; t) of observing n substitutions at time t and the moments of the random variable Nt (vertical red slices). On the other hand, given a substitution number n, one could study the probability density p(t; n) of the time it takes to reach n substitutions, and the moments of the random variable Tn (horizontal, blue slices). of transition from one state α into another β (Fig. 1, inset). Each time a transition occurs, the random variable N is incremented by one unit, while the random variable time T is incremented by t α , an exponentially distributed time with a parameter depending on the state α that the system leaves (Fig 1a). We can represent each realization of this process by a random path (N, T ) (Fig 1b).

m α β α, β ∈ {1 . . . K} (
The problem of the substitution number is to evaluate the moments (variance and higher order moments) of N t , the number of substitutions observed at time t. A complementary and valuable information is the moments of the random variable T n , the time it takes to reach n substitution numbers. To our knowledge, this last problem has not been addressed in the literature. The substitution number is used to evaluate the length of the branches in a phylogeny tree [START_REF] Yang | Molecular Evolution: A Statistical Approach[END_REF]. Under the assumption of molecular clock [START_REF] Bromham | The modern molecular clock[END_REF][START_REF] Simon | Molecularclock methods for estimating evolutionary rates and timescales[END_REF], the time of divergence can be estimated from the substitution number. Its intrinsic uncertainty could be in principle deduced from the moments of the substitution number. Computing the moments, and in particular the variance, of T n is a more direct approach to the problem of divergence time.

The main result of this article is that a simple transformation, a matrix shift, allows one to compute the quantities described above, both for N t and T n .

First, let us recall the basic denitions used in the eld of molecular evolution. The substitution matrix

Q = M -D (2)
contains all the information about the rates : The matrix M contains the transition rates m α β (α ̸ = β) between states while the diagonal matrix D = diag(m α ) contains the rate of leaving state α

m α = K β=1 m α β (3)
The matrix Q is not symmetric and we must make a clear distinction between column (right) and row (left) vectors (i.e. linear forms). The Dirac notations [START_REF] Dirac | A new notation for quantum mechanics[END_REF] are extremely useful to manage this distinction and will be used through this article. A column vector is noted |u⟩, a row vector ⟨u|, left (right) multiplication of a matrix by a row (column) vector are ⟨u| R and R |u⟩; ⟨u|v⟩ is a scalar while |u⟩ ⟨v| is a matrix.

By denition, each column of Q sums to zero, therefore it has a zero eigenvalue :

⟨1| Q = ⟨0| (4)
where ⟨1| = (1, 1, ..., 1) is the left (row) eigenvector. The right (column) eigenvector associated to the zero eigenvalue is denoted |π⟩:

Q |π⟩ = |0⟩ (5)
and is called the equilibrium probability of the states.

The average rate m, i.e. the weighted average of the leaving rates is dened as :

m = K α=1 m α π α = ⟨1|D|π⟩ = ⟨1|M |π⟩ (6)
The rst moment (average) of N t and T n are straightforward to compute :

Nt = mt (7) Tn = 1 m n (8)

B. Moments of Substitution number Nt

As we discussed above, the average quantities are not sucient to evaluate a random variable and one must possess the higher central moments such as variance, third and fourth moments. An exact estimation of the higher moments had proved extremely dicult. The source of this diculty resides precisely in the existence of the zero eigenvalues of Q (see section III). We show here that this diculty can be circumvented by making a simple transformation of the substitution matrix : using the matrix |π⟩ ⟨1| whose columns are all equal and formed of the eigenvector |π⟩ and dening the invertible matrix

Q = m |π⟩ ⟨1| -Q (9)
one can obtain the higher moments of N t in terms of simple manipulations of the matrix Q, its inverse Q-1 and the matrix M = M -mI.

In the long time limit the second (variance, µ 2 ), third (µ 3 ) and fourth (µ 4 ) central moments of N t are as follows. For the variance, we obtain

µ 2 Nt = 1 + 2 δ 2 m + O(1/t) (10)
The quantity 2δ 2 / m is is called the overdispersion of the molecular clock and δ 2 is found to be

δ 2 = ⟨1| M Q-1 M |π⟩ (11)
Expression ( 10) is equivalent to the expression published in our earlier work [START_REF] Houchmandzadeh | A Simple, General Result for the Variance of Substitution Number in Molecular Evolution[END_REF].

For the third central-moment, the result is found to be

µ 3 Nt = 1 + 6δ 3 m + O(1/t) (12)
where

δ 3 = ⟨1| M Q-1 M ( Q-1 M + I) |π⟩ (13)
Note that for a Poisson process, µ 3 / N = 1.

For the fourth central-moment, the result is found to be

µ 4 Nt 2 = 3 1 + 4 δ 2 m + 2 δ 2 m 2 + O(1/t) (14)
For a Poisson process, µ 4 / N 2 = 3 + O(1/ N ).

Comparison to numerical simulations for µ 3 is given in gure 2 and detailed in section III. For a two state system with only one free parameter (beyond time scaling), a simple formula for µ 3 can be derived (equation 38 and gure 7).

Similar expression for the short time limit can also be found by the same method (see section III).

C. Variance of Tn, time to reach n substitution.

The matrix

M = M D -1
contains the jump probabilities from state α to β,knowing that a transition has occurred (see section IV). This matrix has one eigenvalue equal to unity, and its left and The long term variance of T n is then given as

1 n Var(T n ) = 2 1|D -1 (I -F ) -1 D -1 |ψ -1|D -1 |ψ 2 + O(1/n)
where the term in O(1/n) can be exactly evaluated (equation 57) and the residual is of the form F n that converges exponentially to zero. Figures 3 shows the comparison of this result to numerical simulations. For a two state system with only one free parameter, a simple formula for Var(T n ) can be derived (equation 59 and gure 7).

The next sections are devoted to the derivation of the above results.

III. COMPUTING MOMENTS OF Nt

A. Background.

Consider a random process X that can take K values. The rate of transition from a state β to a state α is noted m β α . The probability P α (t) of observing, at time t, the system in the state α is given by the master equation

d dt P α (t) = K β=1 m β α P β (t) -m α P α (t) (15)
where m α = β m α β is the rate of leaving the state α. Collecting the components P α into the vector |P ⟩ we can rewrite equation [START_REF] Vladimir | Counting labeled transitions in continuous-time Markov models of evolution[END_REF] in matricial notations as

d dt |P ⟩ = M |P ⟩ -D |P ⟩ (16)
where the elements of matrix M are given by the rates m β α (with m α α = 0) and D is the diagonal matrix of the elements m α . The matrix Q = M -D, called the substitution matrix (relation2), has one zero eigenvalue with the associated left and right eigenvectors ⟨1| and |π⟩(relations 4-5).

We can now set a counter N called the substitution number (Figure 4) : each time a transition occurs, N is incremented by one unit. The joined probability P n α (t) of observing at time t the system in state α after n substitutions is given by the Master equations

d dt P 0 α (t) = - K β=1 m α β P 0 α (t) (17)
d dt P n α (t) = K β=1 m β α P n-1 β (t) -m α P n α (t) ; n > 0 (18)
Figure 4. The stochastic process X as a function of time t.

The probability density of jumping from state α to state β is m α β (equ. 1). Each time a transition occurs, the substitution number Nt is incremented by one unit.

or, in matricial notations introduced above :

d dt |P n ⟩ = M |P n-1 ⟩ -D |P n ⟩ (19)
where, without loss of generality, we suppose that the system at t = 0 is in equilibrium. :

|P 0 (t = 0)⟩ = |π⟩
From the expression of |P n ⟩ the average of the random variable N t is computed by

N = ⟨N t ⟩ = ∞ n=0 ⟨1|nP n (t)⟩
and is trivial (see appendix A) to show that :

N = mt (20)
The focus of this article is the higher moments of N t . Let us dene the central moment of order k for the state α as

µ k α (t) = n n -N k P n α (t)
We group these moments into the vectorial moments

|µ k (t)⟩ = µ k 1 (t), . . . µ k K (t)
T

The true quantity of interest, the (scalar) central moment for the substitution number is therefore

µ k = K α=1 µ k α (t) = 1|µ k (t) (21)
B. Deriving the moments equation.

From the recurrence equation ( 19) we can directly derive a set of hierarchical ordinary dierential equations (see appendix A) for the central moments :

d dt |µ k ⟩ = Q |µ k ⟩ + k M |µ k-1 ⟩ + M k ℓ=2 k ℓ |µ k-ℓ ⟩ (22)
where

M = M -mI (23)
and therefore, by denition,

1| M |π = 0
Time derivative of a moment of order k depends only on the moments of equal and lower order. The zero eigenvalue of matrix Q causes the moments to grow as a power of time. This zero eigenvalue however prevents us to solve the equations in terms of convolutions of exp(-Qt) with some functions, as the matrix Q -1 does not exist. One could theoretically use the spectral decomposition of Q to solve equation ( 22) ; Q being asymmetric in general, some of the eigenvalues are complex and this is not a practical solution, specially when the number of states K > 2. Surprisingly, this diculty can be circumvented by a simple matrix shift :

Q = m |π⟩ ⟨1| -Q (24)
The matrix |π⟩ ⟨1| is just a matrix whose columns are the repetition of the eigenvector |π⟩. The matrix Q has the same eigenvectors as the matrix Q. The zero eigenvalue is shifted to m:

Q |π⟩ = m |π⟩ (25)
while the sign of other eigenvalues is reversed : for |e i ⟩ associated to the eigenvalue λ i (< 0 ) of Q, noting that ⟨1|e i ⟩ = 0,

Q |e i ⟩ = -λ i |e i ⟩
The matrix Q is invertible, as all its eigenvalues have their real part greater than 0. Using relation (24), equation ( 22) transforms into :

d dt |µ k ⟩ = -Q |µ k ⟩ + mµ k |π⟩ + k M |µ k-1 ⟩ + M k ℓ=2 k ℓ |µ k-ℓ ⟩ (26)
while the main quantity of interest, the scalar moment µ k = 1|µ k , is obtained from ( 22) by projecting the equation on the ⟨1| linear form:

d dt µ k = k 1| M |µ k-1 + k ℓ=2 k ℓ 1|M |µ k-ℓ (27)
The set of equations (26,27) allows us now to obtain, by the usual techniques of ordinary dierential equations, all the moments. The algorithm is as follows :

Knowing the vectorial moment up to |µ k-1 ⟩, compute the scalar moment µ k from equation ( 27).

This necessitates only to compute the quantities 1|M |µ k-ℓ , i.e. left and right multiplication of the Matrix M by known vectors.

Knowing the scalar moment µ k , compute the vectorial moment |µ k ⟩from the relation (26). This necessitates only the knowledge of |µ k-j ⟩ and µ k . By denition, |µ 0 ⟩ = |π⟩, so the recurrence in the above process can be processed step by step. This two step resolution allows for the isolation of the eect of the zero eigenvalue from the other ones, as detailed below.

C. Solving the moments equations.

Equations (26,27) are best studied in the Laplace do-

main |μ k (s)⟩ = ˆ∞ 0 e -st |µ k (t)⟩ dt
where linear dierential equations are transformed into algebraic ones. The long time limit t → ∞ can be obtained by the limit s → 0 in the Laplace domain. On the other hand, the short range limit t → 0 is obtained by the limit s → ∞ of the Laplace transform.

Using the basic rules of Laplace transforms, equations (26,27) are transformed into

μk = k s 1| M |μ k-1 + 1 s k ℓ=2 k ℓ 1|M |μ k-ℓ (28)
and

|μ k ⟩ = m s + m μk |π⟩ + kU s M |μ k-1 ⟩ + U s M k ℓ=2 k ℓ |μ k-ℓ ⟩ ( 29
)
where we have used the shorthand notation

U s = sI + Q -1 As (sI + Q) |π⟩ = (s + m) |π⟩, the simplication U s |π⟩ = 1 s + m |π⟩
has been used to derive expression (29). Moreover, development of U s in power of s is straightforward

U s = Q-1 -s Q-2 + . . .
and is used below to obtain the long time limit of the moments.

Let us know apply the resolution algorithm to solve for the rst moments of N t . The resolution process, as we described above, uses the successive steps

μ1 ⇒ |μ 1 ⟩ ⇒ μ2 ⇒ |μ 2 ⟩ ⇒ μ3 . . .
We already know that μ1 = 0 by denition. As |µ 0 ⟩ = |π⟩, from eq.( 29), we obtain

|μ 1 ⟩ = 1 s U s M |π⟩ (30)
From the above expression, we compute the variance of N t and its long time behavior:

μ2 = 2 s 2 1| M U s M |π + m s 2 (31) = 2 s 2 1| M Q-1 M |π + m s 2 + O(1/s) (32)
The dispersion index R, i.e. the variance to mean ratio,

converges to R = 1 + 2 δ 2 m ( 33
)
where

δ 2 = 1| M Q-1 M |π
The above expression is equivalent to the expression obtained, by a more intricate method, in [START_REF] Houchmandzadeh | A Simple, General Result for the Variance of Substitution Number in Molecular Evolution[END_REF] that has been shown to be exact by numerical simulations.

Continuing the resolution, from expression (31), we compute |μ 2 ⟩ :

|μ 2 ⟩ = m s 2 (s + m) 1| M U s M |π + m |π⟩ + 2 s U s M 2 |π⟩ + 1 s U s M |π⟩
The expression for μ3 is trivially obtained if we recall that 1| M |π = 0 :

μ3 = 6 s 2 1| M U s M 2 |π + 6 s 2 1| M U s M |π + m s 2 (34)
In the long time limit (s → 0), keeping the leading terms in the above expression, we have

µ 3 Nt = 1 + 6 m 1| M Q-1 M Q-1 M + I |π + O(1/t) (35)
The explicit computation of µ 4 (relation A) is given in the appendix A. ison between numerical simulation and theory for two states system. For each value of q, the 2 × 2 matrix M (relation 36) is generated and is normalized to have m = 1, so all numerical simulations take the same time. For each matrix, 10 7 (A series, blue circles) or 10 6 (B series, red squares) paths are generated for a total time of T = 1000 using a standard

Gillespie algorithm [START_REF] Daniel T Gillespie | Exact stochastic simulation of coupled chemical reactions[END_REF] ; each simulation for a given matrix takes approximately 50 minutes. The value of µ 3 at each time step and its limit value is computed from these paths. Solid line is the theoretical expression (38).

D. Applications.

As an example, consider a two states system. Noting the two transition rates mq and m(1 -q), we have

M = m 0 1 -q q 0 ; |π⟩ = 1 -q q (36)
and m = 2mq(1 -q). All the matricial computations can be performed symbolically and we have, for large times,

µ 2 / Nt = 1 + (1 -2q) 2 (37)
µ 3 / Nt = 1 + 3(1 -2q) 4 (38)
gure 5 shows the excellent agreement between the above theoretical results and numerical simulations of the stochastic process.

For higher dimensional matrices, the symbolic expressions are too cumbersome to compute symbolically. Fig-

ure 2 shows the comparison of the numerical simulations to the theory for random 4 × 4 substitution matrices.

E. Short term limit

The short term limit of the expression for the moments can also be simply retrieved by taking the s → ∞ limit of equations (31,34). Note that in this limit, the operator

U s = sI + Q -1 = 1 s I - 1 s 2 Q + . . . Therefore, expression (31) is developed into μ2 = m s 2 + 2 s 3 1| M 2 |π + O(1/s 4 ) µ 2 = N + 1| M 2 |π m2 N 2 + O(t 3 ) (39)
Note that the coecient ν m = 1| M 2 |π is a variance like expression of the leaving rates weighted by equilibrium probabilities :

ν m = α (m α ) 2 π α - α m α π α 2
By the same token, expression (34) is developed into

μ3 = m s 2 + 6 s 3 ν m + O(1/s 4 ) = N + 3 ν m m2 N 2 + O(t 3 ) (40)

IV. DERIVING THE VARIANCE EQUATION FOR Tn

A. Evolution of the probability density of Tn.

Consider again the stochastic process X, but this time as a function of the substitution number and not time (Fig; 6). We note P α (n) the probability of observing the system in state α, knowing that n substitutions have taken place. Knowing that the system was in state α and that a transition has occurred, the probability of nding the system in state β is

mα β = m α β /m α . (41)
Collecting the coecients P α into the vector |P (n)⟩, and the coecients mα β into the matrix M (not to be confused with the matrix M of preceding section), we can write the evolution of the system as

|P (n)⟩ = M |P (n -1)⟩ (42)
By denition, the columns of the matrix

M = M D -1 (43) sum to 1 : ⟨1| M = ⟨1| (44)
therefore ⟨1| is the left (row) eigenvector of M associated to the eigenvalue 1. We denote by |ψ⟩ the right (column) normalized eigenvector of M associated to the same eigenvalue 1 :

M |ψ⟩ = |ψ⟩ ; ⟨1|ψ⟩ = 1 (45)
Figure 6. The stochastic process X as a function of the number of transitions n. Knowing that a transition has occurred, the probability of jumping from state α to state β is mα β (equ. 41). The residency time τα in the state α is exponentially distributed with parameter m α . The total time tn that the system takes to reach n transition is the sum of these residency times : tn = n i=1 τn.

|ψ⟩ is the equilibrium probability of the X process, considered as a function of n. Without loss of generality, we suppose that initially, the system is in equilibrium

|P (0)⟩ = |ψ⟩ (46)
The vector |ψ⟩ is related to the zero eigenvector of the matrix Q (relation 2) through the relation

D |π⟩ = m |ψ⟩ (47)
The probability density f (τ ; α) that the system resides a time τ in the state α is exponentially distributed with parameter m α :

f (τ ; α) = m α e -m α τ
We can now derive the expression for the random variable T n . The probability density p α (t; n) of the total time t to reach n substitutions and to be in state α is therefore

p α (t; n + 1) = ˆt τ =0 β p β (t -τ ; n) mβ α f (τ ; β)dτ. ( 48
)
In plain words, the probability density that the total time at n + 1 is t and the system in state α equals the probability density that (i) at n the total time was t -τ and (ii) the system was in β and made a transition β → α and (iii) the system resided a time τ in state β, summed over all β and u.

Collecting the coecients p α (t; n) into the vector |p(t; n)⟩, equation ( 48

⟨1| L n (0) = ⟨1| M n = ⟨1| (52) L n (0) |ψ⟩ = M n |ψ⟩ = |ψ⟩ (53)
the mean of the random variable T n can be computed by

⟨T n ⟩ = 1| d ds |p(s; n)⟩ s=0 = ⟨1|L ′ n (s = 0)|ψ⟩ (54)
The matrix K′ (s) = M (D + sI) -2 and K(s) don't commute in general. Therefore, the derivative of the matrix L n (s) has to be computed explicitly :

L ′ n (s) = n-1 m=0 L m (s) K′ (s)L n-m-1 (s)
Using the simplications (44,52,53) and the relation (47), it is straightforward to compute the average

⟨T n ⟩ = n 1|D -1 |ψ = n m
where m is the average leaving rate dened in equation [START_REF] Nguyen | IQ-TREE: A Fast and Effective Stochastic Algorithm for Estimating Maximum-Likelihood Phylogenies[END_REF].

For the second moment,

T 2 n = 1| d 2 ds 2 |p(s; n)⟩ s=0
We note that

L ′′ n (s) = n-1 m=0 L ′ m (s) K′ (s)L n-m-1 (s) + n-1 m=0 L m (s) K′′ (s)L n-m-1 (s) + n-1 m=0 L m (s) K′ (s)L ′ n-m-1 (s)
and therefore (see appendix C)

T 2 n = 2 n-1 m=1 (n -m) 1|D -1 M m D -1 |ψ + 2n 1|D -2 |ψ (55)
As the matrix M has one eigenvalue equal to 1, we cannot explicitly compute the geometric series (55) as it is.

However, we can use the shifted matrix method used in the preceding section and write

M = |ψ⟩ ⟨1| + F

By the above denition, we have

⟨1| F = ⟨0| ; F |ψ⟩ = |0⟩ (56)
and therefore

M m = |ψ⟩ ⟨1| + F m
The Matrix (F -I) is therefore invertible and

n-1 m=1 (n -m)F m = F (I -F ) -2 [I -F n -n(I -F )]
All the eigenvalues of F are ∈ [-1, 1[ and in the limit of large n, F n → 0 exponentially. Note that the special case when F has a -1 eigenvalue does not change the result (see for example, the two states system below). Finally, grouping all terms, it is found that

1 n Var(T n ) = 2 1|D -1 (I -F ) -1 D -1 |ψ -1|D -1 |ψ 2 - 2 n 1|D -1 F (I -F ) -2 D -1 |ψ (57)
where Var(

T n ) = T 2 n -⟨T n ⟩ 2 .
We observe that Var(T n )/n converges as 1/n toward its limit. We stress that expression (55) involves only simple matricial operations: inversion, and left and right multiplication by row and column vectors. ues of Var(Tn)/n for a two states system (eq. 58) for four different values of q. Each matrix is normalized in order to have m = 1. A standard Gillespie algorithm [START_REF] Daniel T Gillespie | Exact stochastic simulation of coupled chemical reactions[END_REF] is used to simulate 10 6 path for each matrix for a duration of nmax = 100.

The Variance at each index n is computed using these paths.

For each random matrix, its numerical variance (symbols) is compared to the theoretical prediction of eq. 59).

C. Applications.

As a rst example, consider the two state system (equation 36). From the denition of the matrix M , it is trivial to show that

M = 0 1 1 0 ; |ψ⟩ = 1 2 1 1 (58)
This is a particular case where M possesses a -1 eigenvalue with the left and right associated eigenvector ⟨i| = (1, -1) and |ψ 1 ⟩ = (1/2)(1, -1) T . The matrix M can therefore be written as

M = |ψ⟩ ⟨1| -|ψ 1 ⟩ ⟨i| M m = |ψ⟩ ⟨1| + (-1) m |ψ 1 ⟩ ⟨i|
For this simple system, we don't have to use the large n approximation and the result is found to be

m2 n Var(T n ) = 1 + (1 -2q) 2 1 + a n n (59)
where a n = 1 2 (1 + (-1) n )

gure 7 shows the excellent agreement between the above expressions and numerical simulations of the stochastic process.

For higher dimensional matrices, the symbolic expressions are too cumbersome to compute explicitly. Figure 3 shows the comparison of the numerical simulations to the theory for random 4 × 4 substitution matrices, and the excellent agreement between them.

A nal remark is in order. We saw in section III that the variance to mean ratio for the random variable N t is always superior to 1 and is overdispersed compared to a Poisson process. For the random variable T n , the situation is reversed. We can write

T n = n i=1 τ i
where the random variable τ i is the time the system takes to make one transition. These variables are exponentially distributed and in particular,

Var(τ i ) = 2 1|D -2 |ψ -1|D -1 |ψ 2
If τ i were independent, we would have Var(T n ) = nVar(τ )

Our numerical simulations show that, although Var(T n ) is linear in n,

Var(T n ) < nVar(τ )

i.e. correlations between τ i decreases the variance for T n .

V. DISCUSSION AND CONCLUSIONS.

The characterization of the random variable substitution number N t is a problem that is stated simply, but its solution has proved, until now, to be extremely dicult. In this article, we have developed a simple matrix shift method, that allows for the exact computation of the moments of this random variable and its associated variable T n : the results obey a moment closure scheme, i.e. results from step k, which are obtained by simple matrix algebra, are used to obtain results for step k + 1.

Even though we have restricted the computation to the rst four moments, a formal algebra software can be programmed in principle to derive the result for arbitrary moments.

Note that in molecular evolution, the main observable is the probability p d (t) that two dierent sequences are dierent at a given site [20] :

p d (t) = 1 -tr(e Qt Π) ( 60
)
where Π is the diagonal matrix formed by the elements of the vector |π⟩. One can estimate p d (t) from the fraction of observed dierences between two sequences p.

For sequences of length L, p is given by a binomial distribution B(L, p) and the variance of the distance estimator

Figure 1 .

 1 Figure 1. The random process of substitution number. Inset : a random process X can have K dierent values (here the four bases of DNA) with rates of transition m α β between states. (a)

Figure 2 .

 2 Figure 2. The third central moment µ 3 (relation 35) : compar-ison between numerical simulation and theory for four states system. 80 random 4×4 substitution matrix M are generated and normalized to have m = 1, so all numerical simulations take the same time. For each matrix, 6 × 10 7 paths are generated for a total time of T = 1000 using a standard Gillespie algorithm[START_REF] Daniel T Gillespie | Exact stochastic simulation of coupled chemical reactions[END_REF] ; each simulation for a given matrix takes approximately 5 hours. (a) For each matrix, the theoretical value of µ 3 is plotted against its numerical value. Solid lines represents the diagonal. (b) The relative absolute dierence between these two values for each matrix.

Figure 3 .

 3 Figure 3. Comparison between theoretical and numerical limit of Var(Tn)/n for 200 random 4 × 4 substitution matrices. Each matrix is normalized in order to have m = 1. A standard Gillespie algorithm[19] is used to simulate 10 6 path for each matrix for a duration of nmax = 200. The Variance at each index n is computed using these paths. The numerical limiting value is extracted from the last 100 indexes. (a) for each random matrix, its theoretical variance is plotted against its numerical value. (b) The relative absolute dierence between these two values for each matrix. The relative dierence is inversely proportional to the number of paths used for the numerical simulation.

Figure 5 .

 5 Figure 5. The third central moment µ 3 (relation 35) : compar-

  ) can be written as |p(t; n + 1)⟩ = ˆt τ =0 M exp(-Dτ ) |p(t -τ ; n1)⟩ dτ (49) where we have used the relation M D = M . The above expression is a convolution product, therefore using the Laplace transform |p(s; n)⟩ = ˆ∞ t=0 e -st |p(t; n)⟩ dt equation (49) transforms into a simple recurrence relation : |p(s; n + 1)⟩ = K(s) |p(s; n)⟩ (50) where K(s) = M (D + sI) -1 is the Laplace transform of the kernel K(τ ) = M exp(-Dτ ) in relation (49). Formally, the solution of the recurrence relation (50) is |p(s; n)⟩ = K(s) the variance of Tn Various moments of the random variable T n can be computed by the derivations of the expression (51) as a function of s for s = 0. Let us dene L n (s) = K(s) n and note that, from equations (44,45)

Figure 7 .

 7 Figure 7. Comparison between theoretical and numerical val-

  By eliminating time in relations (7,60), it is then possible to relate the estimators d (of ⟨n⟩) and p d = f (p)

d can be deduced from relation (61). This quantity however is very dierent from the intrinsic variance of the substitution number that we have computed here.

The above method is very general and can be used for any stochastic process that has a unique equilibrium probability.

Appendix A: Deriving the moments equation for Nt Consider the system of equations

with the initial condition

In order to simplify the computations involving index shift and without loss of generality, we extend the domain of N to negative integers :

as the left hand side is the probability of nding the system in a given state, whatever the number of substitution.

We dene the moment-vector of order k as

The moments of the random variable N are obtained simply by summing over all the components of the corresponding vector-moment :

Consider the mean -vector

Applying n n to both sides of equations (A1) we obtain

where Q = M -D and we have used the change of index n → n + 1 to simplify the expression

Summing over the component of |N ⟩ by applying the left form ⟨1| to expression (A4), we have

m α π α is simply the weighted average of the leaving rates. Integrating equation (A5), we have the linear growth law of the average :

We apply the same method to obtain higher orders moments. We dene the central moment-vector of order

Deriving the above expression as a function of time, we have

We can now use equation (A1) to develop the last term of the above expression. To compute the term

we make the usual change of index n → n + 1 and use the binomial expression

Grouping all the terms and recalling that Q = M -D, we obtain

The above expression is an exact moment closure : the dierential equation of moments-vector of order k depends only in moments-vector of order ≤ k. Projecting on the ⟨1| form and recalling that ⟨1| Q = ⟨0|, we get nally the expression for the centered moments

where M = M -mI.

where

and its matrix element ⟨1|L ′′ (0)|ψ⟩. Using the relation (C5), the rst line (C1) of the above relation is :

Using the simplications

where the second line has been obtained by re-indexing and changing the order of summation.

Appendix D: Numerical simulations.

Numerical simulations are performed by a classical

Gillespie algorithm [START_REF] Daniel T Gillespie | Exact stochastic simulation of coupled chemical reactions[END_REF]. The algorithm is implemented using the Julia language [START_REF] Bezanson | Julia: A fresh approach to numerical computing[END_REF]. The process is simulated as described in gure 1. To generate one path for a given substitution matrix M , we begin with a state α 0 drawn from the equilibrium state |π⟩or |ψ⟩ depending on the problem. The time t 1 to the the next transition is drawn from an exponential distribution with parameter m α0 ; the next state α 1 is drawn at random using the α 0 column of the matrix M of transition probabilities. A loop until the nal index n is used to continue the process and the path is recorded as the vector (t 1 , . . . t n ).

[