{'fraction_non_alphanumeric': 0.05252532449796335, 'fraction_numerical': 0.010794827033382023, 'mean_word_length': 3.9859671975842264, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 16, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 14, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Mutation, selection, and ancestry in branching models: a variational approach

We consider the evolution of populations under the joint action of mutation and differential reproduction, or selection. The population is modelled as a finite-type Markov branching process in continuous time, and the associated genealogical tree is viewed both in the forward and the backward direction of time. The stationary type distribution of the reversed process, the so-called ancestral distribution, turns out as a key for the study of mutation-selection balance. This balance can be expressed in the form of a variational principle that quantifies the respective roles of reproduction and mutation for any possible type distribution. It shows that the mean growth rate of the population results from a competition for a maximal long-term growth rate, as given by the difference between the current mean reproduction rate, and an asymptotic decay rate related to the mutation process; this tradeoff is won by the ancestral distribution. Our main application is the quasispecies model of sequence evolution with mutation coupled to reproduction but independent across sites, and a fitness function that is invariant under permutation of sites. Here, the variational principle is worked out in detail and yields a simple, explicit result.

Abstract. We consider the evolution of populations under the joint action of mutation and differential reproduction, or selection. The population is modelled as a finite-type Markov branching process in continuous time, and the associated genealogical tree is viewed both in the forward and the backward direction of time. The stationary type distribution of the reversed process, the so-called ancestral distribution, turns out as a key for the study of mutation-selection balance. This balance can be expressed in the form of a variational principle that quantifies the respective roles of reproduction and mutation for any possible type distribution. It shows that the mean growth rate of the population results from a competition for a maximal long-term growth rate, as given by the difference between the current mean reproduction rate, and an asymptotic decay rate related to the mutation process; this tradeoff is won by the ancestral distribution.
We then focus on the case when the type is determined by a sequence of letters (like nucleotides or matches/mismatches relative to a reference sequence), and we ask how much of the above competition can still be seen by observing only the letter composition (as given by the frequencies of the various letters within the sequence).
If mutation and reproduction rates can be approximated in a smooth way, the fitness of letter compositions resulting from the interplay of reproduction and mutation is determined in the limit as the number of sequence sites tends to infinity.
Our main application is the quasispecies model of sequence evolution with mutation coupled to reproduction but independent across sites, and a fitness function that is invariant under permutation of sites. In this model, the fitness of letter compositions is worked out explicitly. In certain cases, their competition leads to a phase transition.

Introduction
Evolution is often understood as an optimization process of some kind, and there is a long tradition to consider evolutionary models, particularly those from population genetics, from a variational perspective. The most popular result in this context is known as the Fundamental Theorem of Natural Selection (FTNS). In its simplest form it states that, in the deterministic selection equation for a single locus in continuous time, mean fitness can only increase along trajectories (i.e., it is a Lyapunov function), and the rate of this increase equals the variance in fitness, cf. [5,Ch. I.10.3]. More sophisticated versions in the context of quantitative genetics and multiple loci, along with a general discussion of optimality principles for the selection equation, are discussed in [5, Ch. II.6.3-II. 6.6] and [11,Ch. 2.9,7.4.5 and 7.4.6]; see also [8].
If, rather than selection alone, the joint dynamics of selection and mutation is considered, results become sparse. The FTNS may be generalized to house-of-cards mutation (i.e., mutation rates are independent of the parent type), see [1] and [19]. If mutation is reversible, a Lyapunov function is available for a certain L 2 -renormalized version of the dynamics, but not for the original mutation-selection equation [33].
The above approaches refer to the genetic (or, more generally, type) composition at the population level. In contrast, this article is concerned with a variational principle in mutation-selection models (and closely related branching processes) from the point of view of individual lineages through time, their ancestry and genealogy. This principle is related to the (stochastic) processes that take place along such lines of descent, with a special emphasis on the relation between the present and the past. We will, however, not include genetic drift (i.e., resampling) into our models; therefore, our backward point of view differs from that of the coalescent process (see [17] for a recent review of this area).
The paper is organized as follows. In Section 2, we will set up our model(s) and recapitulate a few fundamental facts. Section 3 provides an informal preview of the results that will be detailed (and proved) in the remainder of the article. Section 4 will develop the lineage aspect that will be required furtheron. Looking at the mutation process along individual lines, we will obtain a fairly general variational principle (Section 5, Thm. 1), which quantifies the tradeoff between the mean reproduction rate along a line and the asymptotic rate at which it is lost; it further implies a connection between the type processes that emerge in the forward and backward directions of time. In Sections 6 and 7, we will specialize on the case where types are sequences over a finite alphabet. If mutation is independent and fitness is additive across sites, the original high-dimensional variational principle may be reduced to a simpler, low-dimensional one (Section 6.3, Thm. 2). The same holds asymptotically if mutation rates and fitness function allow for a suitable smooth approximation when the number of sites gets large (Section 7, Thms. 3, 4 and 5). The corresponding approximate maximum principle will be derived explicitly for the quasispecies model of sequence evolution (Section 8, Thm. 6).
The paper ties together, unifies and generalizes various aspects that have appeared in previous publications. Special cases of the low-dimensional maximum principle were first described in [18], and applied to concrete examples. An extension appeared in [3]; it relies on methods from linear algebra and asymptotic analysis, but makes no connection to the stochastic processes on individual lines, nor does it include worked examples. The connection to the backward point of view relies on earlier work on branching processes [22,23] and was investigated in [18] and [15]. These results will reappear here as parts of a larger picture.

Models
Consider a finite set of types S (with |S| > 1) and a population of individuals, each of which carries one of these types. (We think of individuals as haploid, and of types as alleles.) 2.1. 1 The parallel mutation-reproduction model. Let us start with the most basic mutation-reproduction model in which mutation and reproduction occur in parallel, that is, independently. As depicted in Fig. 1 individual of type i ∈ S may, at every instant in continuous time, do either of three things: It may split, i.e., produce a copy of itself (this happens at birth rate B i ≥ 0), it may die (at rate D i ≥ 0), or it may mutate to type j (j = i) (at rate U ij ≥ 0). Different meanings may be associated with this verbal description. Probabilists will take it to mean a multi-type Markov branching process in continuous time (see [2,Ch. V.7], or [25,Ch. 8] for a general overview). That is, an i-individual waits for an exponential time with parameter A i = B i + D i + j:j =i U ij , and then dies, splits or mutates to type j = i with probabilities B i /A i , D i /A i , and U ij /A i , respectively. The number of individuals of type j at time t, Z j (t) ∈ Z ≥0 := {0, 1, 2 . . . }, is a random variable; the collection Z(t) = Z j (t) j∈S is a random vector. The corresponding expectation is described by the first-moment generator A = U + R. Here, U is the Markov generator U = (U ij ) i,j∈S , where the mutation rates U ij for j = i are complemented by U ii := − j:j =i U ij for all i ∈ S. Further, R := diag{R i | i ∈ S}, where R i := B i − D i is the net reproduction rate (or Malthusian fitness). More precisely, we have E i (Z j (t)) = (e tA ) ij , where E i (Z j (t)) is the expected number of j individuals at time t in a population started by a single i-individual at time 0.

Deterministic aspects.
Ignoring stochastic effects and focussing on the mean behaviour of the population, one often considers the deterministic mutation-reproduction model where y(t) = (y i (t)) i∈S is the row vector associating to each type i its abundance y i (t) ∈ R ≥0 (i.e., the size of the subpopulation of type i). As y(t) = y 0 e tA , the deterministic model describes the expectation of the corresponding branching process, provided the initial condition is chosen accordingly. However, the independent reproduction of individuals as implied so far is unrealistic for large populations. They usually experience density regulation; in the simplest case, this is modelled by an additional death term γ(t) ≥ 0, that is, D i is replaced by D i + γ(t) (for all i ∈ S), where γ(t) may depend on time (maybe through total population size), but not on the type. Then, of course, (1) generalizes tȯ where I is the identity matrix. In theoretical ecology, a wide variety of models is in use that specify γ for the many biological situations that may arise. In population genetics, however, one is usually more interested in the relative frequencies q i (t) := y i (t)/ j y j (t). Differentiating this and inserting (2) leads tȯ independently of γ. Here we think of the row vector q = (q i ) i∈S as a probability measure, of the column vector R = (R i ) i∈S as a function on S (known as the fitness function), and of the scalar product q(t), R = i∈S q i (t)R i as the associated expectation, namely, the mean fitness of the population at time t. Eq. (3) is the well-known parallel (or decoupled ) mutation-selection model, which goes back to [6, p. 265]. Although we have derived it here for haploid populations (and will adhere to this picture), it is well known, and easily verified, that the same equation describes diploids without dominance (in an approximation using Hardy-Weinberg proportions). For a comprehensive review of the model and its properties, see [5,Ch. III]. Rather than considering deterministic and stochastic models separately, we aim at a unifying picture and note that the branching process is particularly versatile: Its expectation fulfills (1), and the solution of (1), in turn, implies that of (3) (via normalization). Properties of the branching process will, therefore, immediately translate into properties of the mutationselection equation (but not, necessarily, vice versa). For this reason, we will consider the branching process as our primary model throughout this paper. Let us, therefore, return to branching populations and look at alternatives to the parallel model.
2.1. 3 The coupled mutation-reproduction model. In this model one assumes that mutations occur on the occasion of reproduction events (see Fig. 2): An i-individual again dies at rate D i and gives birth at rate B i , but every time it gives birth, the offspring is possibly mutated (of type j with probability P ij ), while the parent itself survives unchanged. The corresponding first-moment generator A has elements An example of the coupled model will be studied in Sec. 8.

Di
BiPij i i j Fig. 2. The coupled mutation-reproduction model.

General splitting rules.
Both the parallel and the coupled models are special cases of the general Markov branching model as depicted in Fig. 3: An i-individual lives for an exponential time τ i with prescribed parameter A i and then produces a random offspring N i = (N ij ) j∈S with distribution p i on Z S + and finite means E(N ij ) for all i, j ∈ S. More precisely, N ij ∈ Z ≥0 is the number of children of type j, and p i (κ) = P(N ij = κ j , ∀j ∈ S). The first-moment generator A has elements For the coupled and the general branching rules, the first-moment generator may again be written in the 'parallel' form A = U + R where U is a Markov generator and R is a diagonal matrix; this decomposition is uniquely given by U ij = A ij for i = j, U ii = − j:j =i U ij , and R i = j∈S A ij for all i ∈ S. At the time being, this is a formal decomposition, but will receive its branching process interpretation later in Sec. 4.2. The corresponding deterministic models then all take the form (1) and (3), provided the parameters are interpreted in the above way.

Fundamental facts 2.2.1 Forward view and long-time characteristics.
We will assume throughout that A (or, equivalently, U ) is irreducible. Perron-Frobenius theory then tells us that A has a principal eigenvalue λ (namely a real eigenvalue exceeding the real parts of all other eigenvalues) and associated positive left and right eigenvectors π and h which will be normalized so that π, 1 = 1 = π, h , where 1 = (1) i∈S is the vector with all coordinates equal to 1. We will further assume that λ > 0, i.e., the branching process is supercritical. This implies that the population will, in expectation, grow in the long run, as is obvious from (1); in individual realizations, it will survive with positive probability, and then grow to infinite size with probability one, see (6) below.
The asymptotic properties of our models forward in time are, to a large extent, determined by λ, π, and h, and provide further connections between the stochastic and the deterministic pictures. The left eigenvector π holds the stationary composition of the population, in the sense that lim t→∞ q(t) = π for the differential equation (3), and, for the branching process, lim t→∞ Z(t) Z(t) 1 = π with probability one, conditionally on survival, (5) where Z(t) 1 := j∈S Z j (t) is the total population size. This is due to the famous Kesten-Stigum theorem, see [27] is the asymptotic growth rate (or equilibrium mean fitness) of the population. Here the first equality follows from the identity λ = πA, 1 = π, A1 = π, R ; the second one is an immediate consequence of (1) and Perron-Frobenius theory, and the third is from [15] and holds with probability one in the case of survival. Finally, the i-th coordinate h i of the right eigenvector h measures the asymptotic mean offspring size of an i individual, relative to the total size of the population: For more details concerning this quantitity, see [18] and [15] (for the deterministic and stochastic pictures, respectively).

Backward view and ancestral distribution.
In the above, we have adopted the traditional view on branching processes, which is forward in time. It is less customary, but equally rewarding, to look at branching populations backward in time. To this end, consider picking individuals randomly (with equal weight) from the current population and tracing their lines of descent backward in time (see Fig. 4). If we pick an individual at time t and ask for the probability that the type of its ancestor is i at an earlier time t − τ , the answer will be α i = π i h i in the limit when first t → ∞ and then τ → ∞. Thus the distribution α = (α i ) i∈S describes the population average of the ancestral types and is termed the ancestral distribution, see [15,Thm. 3.1] for details. Likewise, the time average along ancestral lines also converges to α in the long run, see [15,Thm. 3.2]. If we pick individuals from the population at a very late time (so that its composition is given by the stationary vector π), then the type process in the backward direction is the Markov chain with generatorḠ = (Ḡ ij ) i,j∈S , G ij = π j (A ji − λδ ij )π −1 i , as first identified by Jagers [22,23]. The corresponding time-reversed process has generator G = (G ij ) i,j∈S , where it has been considered in [15], has been termed the retrospective process, and may be understood as the forward type process along the ancestral lines leading to typical individuals of the present population. By definition, G andḠ both have stationary distribution α.

Preview of results
In this Section, we will give an informal preview of the results that will be obtained in the remainder of the article. This overview will not aim at full generality, nor will it dwell on specific technical conditions that are required to make things precise. Rather, we will try and motivate the concepts and explain the results in the context of the model. The details will be worked out in the later Sections.
We will work our way from the more general to the more specific. We will start with a general variational principle, valid for all model variants of the previous Section, irrespective of the type space and of the parameters. Next, we will specialize on the case where types are sequences, and mutation and reproduction rates are invariant under permutation of sites. This will allow to dissect the variational problem into two simpler problems, which are easier to solve. Finally, we will treat one specific example, namely, the quasispecies model of sequence evolution, in full detail.

The general variational principle
A main object of this paper is to show that the asymptotic growth rate λ of the population can be understood as the result of a competition between the mutation and reproduction processes along a typical ancestral line. In this informal Section, however, we will avoid the family tree picture, and rather imagine we are observing just one line. To start with, we even ignore reproduction, and consider only the simple Markov process {M (t)} t≥0 on S with generator U ; i.e., the type process which associates with t the type at time t under the mutation model U . A crucial quantity in what follows will be the corresponding empirical measure i.e., the random vector with components L i (t) := 1 t t 0 I{M (τ ) = i}dτ , where I{.} denotes the indicator function. This quantity measures the fraction of time the process spends in the various states, and hence is also known as occupation time measure. Clearly, L(t) is a random element of P(S), the set of all probability measures on S. It is well-known by the ergodic theorem for Markov chains that, for t → ∞, one has L(t) → ρ with probability one, where ρ is the stationary distribution of U . It is, perhaps, less wellknown that the rate of convergence may be characterized asymptotically by a so-called large deviation principle, which may be informally put as that is, the probability that L(t) is close to some measure ν decays exponentially, for large time, with a decay rate (or rate function) I U (ν) which can be written down explicitly (see (25) and (30) below). I U is nonnegative, and I U (ν) = 0 precisely for ν = ρ, in line with the above fact that, in the long run, only the stationary measure ρ will survive. Let us now add reproduction, i.e., turn to the branching process. As a consequence of the above large deviation principle, we will obtain, in Thm. 1, a link between the forward-time stationary distribution π of the branching process, the reproduction rate R, the asymptotic growth rate λ, the mutation process U and the ancestral distribution α, namely, the equation This variational principle may be understood in terms of a competition between all possible distributions for a maximal long-term growth rate, as given by the difference between the current mean reproduction rate ν, R , and the asymptotic decay rate I U (ν). The first quantity is maximized by those measures that put mass only on the fittest type(s); the second one is minimized by ρ; the tradeoff is won by α. Furthermore, (11) connects the forward and the retrospective point of view in that the maximum equals the mean fitness π, R of the stationary population. Note that the mean fitness of the ancestral population exceeds the mean fitness of the stationary one by I U (α), which is positive unless α = ρ (which implies R i = const., i.e., there is no selection). This reflects the fact that the present population carries with it a tail of (mainly unfavourable) mutants that are present at any time, but do not survive in the long run.
We will see in Sec. 5 that this 'competition of distributions' can be made more concrete, namely, in terms of a competition of lines of descent, by considering the empirical distributions L ω (t) of types along distinct lines ω. But before we can embark on this, we must first develop a way of constructing trees, lines, and processes on lines, in a consistent way; this will be taken up in the next Section.
It is interesting to note that the above variational principle resembles the thermodynamic maximum principles in statistical physics. Indeed, our reproduction rates may be identified with an energy, and the rate function with an entropy; in fact, the rate function for the continuous-time Markov chain M (t) can be naturally derived from the usual entropy governing the so-called pair-empirical measure of a discrete-time Markov chain, cf. [20, Ch. IV].

Sequence space models
The variational principle (11), valuable as it is conceptually, is not very useful if one aims at an explicit solution; this is because maximization is over a large space (the set of probability measures on S). However, it turns out that, in certain models of sequence evolution, this task boils down to a much simpler one if the original problem is dissected into two, one of which can be solved explicitly. Let us first describe this 'divide and conquer' strategy.
Assume that the type of an individual is characterized by a sequence of nucleotides, amino acids, matches/mismatches with respect to a reference sequence of nucleotides, or, in general, letters from some alphabet Σ. Thus N is typically large. However, if the mutation and reproduction rates are invariant under permutations of sequence sites, all relevant information on a sequence σ = (σ k ) 1≤k≤N ∈ Σ N is already contained in its letter histogram (or letter composition) which indicates how often each letter ℓ shows up in σ. In other words, it is sufficient to look at the reduced type space with d = |Σ|, which consists of all possible letter compositions. This lumping procedure induces a model on S that is again a Markov branching process; its reproduction rates R = (R i ) i∈S and mutation generator U = (U ij ) i,j∈S are uniquely determined by the corresponding rates of the original process on Σ N . Many models of sequence evolution allow for such a lumped representation; as a particularly realistic example, let us mention the mutationselection model for regulatory DNA motifs [16], which also involves analysis of sequence data.
To get back to the variational problem, we will classify the possible distributions ν ∈ P(S) according to the value of their mean ν, id ∈ R d ; here id denotes the identity function on S defined by id i = i for all i ∈ S, and, in line with previous usage, the scalar product gives the expectation of this vector-valued function under the measure ν. Keeping in mind that S arises from lumping a sequence space Σ N as in Fig. 5, we think of ν, id as the expected value of a random letter composition with distribution ν, i.e., the mean histogram if histograms have distribution ν.
Let us now foliate the variational problem (11) according to these mean letter frequencies. That is, we write with Λ : conv S → R given by Here we write conv S for the convex hull of S, that is, the set of convex combinations of elements of S, or, in other words, the set of all possible mean letter compositions. The (unique) maximizer of Λ isẑ := α, id , i.e., the mean ancestral letter composition. The function Λ(z) describes the growth rate resulting from the competition between all distributions with mean letter composition z; we will therefore call it the constrained mean fitness of z. In analogy with the interpretation of the unconstrained variational principle (11), the competing distributions may be identified with empirical letter compositions along lines of descent, and Λ(z) will turn out as the asymptotic growth rate of the lines with empirical letter histogram (close to) z; this will be shown in Prop. 2. It follows that the growth rate of the total population coincides with the growth rate Λ(ẑ) of the subpopulation consisting of all lines with empirical letter histogram close to the mean ancestral one. Now, the main point is that Λ(z) can be calculated explicitly in two interesting situations, namely: (1) All sites of the sequence mutate independently and according to the same Markov process in continuous time, and fitness is additive across sites. Thus R i = R(i) and U ij = U j−i (i) are linear functions on S (that will be extended to conv S). In Thm. 2, we then obtain Λ(z) explicitly, and exactly, as where the sum is over all possible mutational steps, and ν (z) = Mult N,z/N is the multinomial distribution with mean z.
(2) The reproduction and mutation rates have a continuous approximation of the form with functions r and u k that are smooth enough. Under further technical conditions, an analogue of (16) will be obtained in Thm. 4, namely, where Strictly speaking, the approximation (18) is only true when e(z) is concave; otherwise e(z) has to be replaced by its concave envelope, and the distribution attaining the constrained maximum Λ(z) will show distinct peaks. This behaviour, which indicates some kind of phase transition, will be the subject of Thm. 5. For z =ẑ, this phenomenon means that the total growth rate λ is determined by two or more coexisting subpopulations with distinct empirical letter histograms.

The quasispecies model
We will finally consider the coupled sequence space model on {0, 1} N , known as the quasispecies model ; more precisely, we will use a slightly adapted version of the original in [9]. It will be assumed that births and deaths occur at rates that are invariant under permutation of sites, and mutations occur on the occasion of birth events, independent across sites, and at probabilities v = µ/N and w = ν/N from 0 to 1 and vice versa, where µ and ν are positive and independent of N . Then, lumping may be performed into S := {0, 1, . . . , N } by counting the number of 1's in a sequence. If the birth and death rates of the resulting model on S have a continuous approximation analogous to that of (17), namely, takes the role of e(z) in (19).

Trees, lines, and processes on lines
To understand the probabilistic significance of the variational principle previewed above, it is necessary to develop a detailed picture of the branching process that includes the full family tree. However, to keep technicalities at a minimum we confine ourselves, in the first subsection, to the parallel model; in this case, a particularly simple construction is available which is sufficient for our needs. A more versatile procedure for general splitting rules will be sketched in Subsection 4.2.

The parallel model
Let us explain the construction for the parallel model, as illustrated in Fig. 6. The population is started by a single individual (the root ) of type i. In a first step, we ignore all death events and consider only the splitting events. Then all lines are infinite and can be labeled by a sequence ω ∈ {0, 1} Z ≥1 =: Ω, where ω n tells us whether the n-th offspring corresponds to the upper (0) or lower (1) branch in the graphical representation of the tree, or, equivalently, whether it is counted as 'first' or 'second' at birth. Next, individuals are defined as (finite) initial segments of the infinite lines, i.e., x = (ω 1 , . . . , ω n ) is an n-th generation individual. The empty initial string ∅ of length 0 corresponds to the root and is counted as generation 0. The set X := {∅} ∪ ( n≥1 {0, 1} n ) then comprises all individuals that may possibly occur (and do occur as long as death events are ignored). A realization of the Markov branching process described informally in Sec. 2 may then be specified by associating with every line ω the times at which it splits, its type (as a function of time), and the time it dies (by a death event). For convenience, the construction proceeds in two steps: we first grow a tree by splitting and mutation alone (with the appropriate exponential waiting times); the death events are then superimposed in a second step to determine which lines are still alive. This way, lines that have already died live on virtually and may continue to divide and mutate. However, this does not influence the lines that are alive; only these constitute the realization of the branching process. In particular, we denote by X(t) ∈ X the set of individuals alive at time t; note that this is a mixture of various generations. (We remain a bit informal here; for one of the various possible ways of a rigorous construction, see [15].) The branching process with mutation and binary splitting. Bullets mark death events; line segments that are alive are shown in black, virtual ones in grey, types are indicated by various line styles. The (randomly chosen) representative line is marked fat; its initial segment shown here is the first child of the second child of the first child of the root, i.e., the individual x = (0, 1, 0). Since it has experienced three splitting events, it is a third generation individual; but it is virtual, as in fact already its mother (0, 1) died. The 'black' tree is a realization of the branching process. The individuals alive at time t constitute the population X(t) (a mixture of various generations); here, X(t) = {(0, 0), (1, 1, 0), (1, 1, 1)}; the other individuals at time t are virtual.
For each line ω, we consider now the following families of random variables: {M ω (t)} t≥0 , the type process, which associates with t the type of ω at time t; {β ω (t)} t≥0 , the number of birth events along ω before t; and T ω , the time line ω dies (if ω survives forever, this time is infinite). Both the birth and the death process depend on the type process, but not vice versa. The crucial information on {M ω (t)} t≥0 is contained in its empirical measure cf. (9). For an individual x at time t, the empirical measure only depends on the initial segment of ω that describes x. With this in mind, we will sometimes also write L x (t) rather than L ω (t).
The above families of random variables are not independent between lines (they are dependent through common ancestry), but, by symmetry between the two offspring at every splitting event, they share the same marginal laws for all ω ∈ Ω. In particular, since mutation is not influenced by the reproduction events, the type process on any given line (regardless of the others) is a copy of the mutation process generated by U . Let us choose one particular such line ω * , for example, by setting ω * = (000 . . .), or by tossing a coin. The line ω * may or may not survive, but it will always be present at least virtually. We will call it the representative line for reasons to become clear in a moment, and set β(t) We will now see that, once we know the laws of these quantities, they can tell us a lot about the entire tree.
The basic observation is that, in generation n, there are 2 n possible (real or virtual) individuals, all with the same marginal laws for the random variables just discussed. This allows us to express the expected population size of a population started by a single i individual at time 0 as follows: Now, conditionally on L(t), the random variables {T > t} and β(t) are independent, having probability exp(−t L(t), D ) resp. the Poisson distribution Poi t L(t),B with parameter t L(t), B . Therefore, (both independently of the type of the root), where the former relies on the fact that, for a random variable Y with distribution Poi λ , one has E(2 Y ) = e λ . Therefore, (22) turns into Note that the remaining expectation (and the outer one where expectations are nested) is with respect to L(t). We also remark that the underlying tree construction lurks behind the above derivation, but in the simple case at hand it need not be made more explicit.

General splitting rules
We have, so far, restricted ourselves to the decoupled model with parallel mutation, reproduction and death. The crucial simplifiying feature here is the fact that, forward in time on every line, we have a copy of the mutation process generated by U . Therefore, we could consider any line as representative.
Outside the parallel model, the decomposition A = U + R is formal to start with, and the generator U has no immediate interpretation. But with the help of a more advanced tree construction, one can again obtain a representative line with its type process M (t) generated by U . We will only give a rough sketch here; for the full picture we refer the reader to [15].
The construction relies on a so-called size-biased tree with random spine (or trunk ). The general concept was introduced in [21,30] and [28]; the particular (continuous-time) version required here can be found in [15,Remark 4.2]. Informally, one constructs a modified tree with a randomly selected, distinguished line (called the trunk or spine), along which time runs at a different rate and offspring are weighted according to their size; in particular, there is always at least one offspring along the trunk so that the trunk survives forever. The children off the trunk get ordinary (unbiased) descendant trees; see Fig. 7. Fig. 7. A realization of a size-biased tree with its trunk (the fat line). An individual of type j, off the trunk, has offspring Nj with distribution p j after an exponential waiting time τj with mean 1/Aj; an individual of type i along the trunk bears offspring e Ni with biased distributionp i after an exponential waiting timeτi with mean 1/AiE( Ni 1).
More precisely, for each type i ∈ S, we introduce the size-biased offspring distributionp Starting at the root, an individual of type i on the trunk waits for an exponential time with parameter A i E( N i 1 ) and then produces offspring N i according top i ; one of these offspring is chosen randomly (with equal weight) as the successor on the trunk. It is easily verified that the type process on the trunk is a Markov chain generated by U . The trunk takes the role of the representative line, and the considerations of the previous Subsection carry over. We do not spell this out here explicitly; for the complete picture and many details, in particular on how the trunk may be used to extract further information about the tree, see [15]. To avoid misunderstandings, we would like to emphasize that the size-biased tree as applied to the parallel model does not reduce to the simple special construction of the previous Subsection. In particular, unlike the representative line of this construction, the trunk of the size-biased tree is certain to survive forever. However, both constructions share the essential property that the mutation process along the trunk or representative line, respectively, is generated by U , and the fact that many properties of the entire tree may be extracted from this distinguished line.

Variational characterization of the asymptotic growth rate
We are now in a position to derive the variational characterization (11) of the asymptotic growth rate λ. The idea is to observe both the mutation process and the reproduction rate along the representative line of the tree. The appropriate tool for analyzing the tradeoff between these processes is the large deviation principle for the mutation process.

Using the large deviation principle
Let us, for the moment, restrict ourselves to the parallel model; we will see later that our results hold automatically for general splitting rules. For the parallel model, we can combine (6) and (23) to obtain that is, the growth rate can be determined by observing the types and the associated reproduction rates along the representative line. The competition between reproduction and mutation will lead to a variational formula for λ, which can immediately be derived from the variational formulas of large deviation theory. The basic fact is the following large deviation principle for where the supremum is taken over all v ∈ R S >0 , and the fraction is to be understood component-wise, i.e., U v/v is the vector with components for any closed set C ⊂ P(S), and lim inf for any open set O ⊂ P(S). Furthermore, I U is continuous, strictly convex and nonnegative, and I U (ν) = 0 precisely for ν = ρ, the stationary distribution of U .
For an informal statement of the LDP recall (10). (Although we have stated the LDP here only for the special case we need, it is indeed quite a general principle that applies to many common types of random variables. We refer the interested reader to the monographs [7] or [20].) Returning to (24), we now see that, on the right-hand side, the exponential factor e t L(t),R is integrated over a probability measure that behaves essentially like e −tIU . It may thus be evaluated by Varadhan's lemma on the asymptotics of exponential integrals, which is a far-reaching generalization of Laplace's method; see [20,Thm. III.13] which may be understood as a 'largest exponent wins' principle. Let us continue with a series of comments. 5.1.1 Relation to the retrospective process. The maximum principle (26), though derived by considering the branching process forward in time, is directly connected to the retrospective process of (8). In analogy with (25), the rate function for the empirical measure of the retrospective process (generated by G of (8)) reads I G (ν) = sup w>0 [− ν, (Gw)/w ]. This, however, is closely related to I U (ν). Indeed, Again, I G (ν) is nonnegative, strictly convex, and vanishes if and only if ν = α, the stationary distribution of G. It follows that the ancestral distribution α is the unique maximizer in (26). We may thus summarize our findings in the following theorem (recall (6) for the first identity). Theorem 1. The forward-time stationary distribution π, the reproduction rate R, the asymptotic growth rate λ, the mutation process U and the ancestral distribution α are linked via the equation

The mutation rate function at the ancestral distribution. Thm. 1 yields the additional relation
i.e., the value of the mutational rate function at the optimum equals the long-term loss of offspring due to mutation, wherefore it was previously termed mutational loss function; see [18, Sec. 5 and Appendix A] for the biological implications.

Balance of mutation and reproduction.
On every line ω, the mutation process runs randomly through a sequence of histories, and hence determines an evolution of empirical measures L ω (t) ∈ P(S). As t → ∞, the empirical measures ν = L ω (t) that differ from the stationary distribution ρ of U become exponentially less probable at asymptotic rate I U (ν). In particular, ρ is the (almost-sure) long-term time average on the line ω in the forward direction of time. In spite of this, the long-term population average π of (5) differs from ρ, in general. This is because mutation is counterbalanced by reproduction, at rate R M ω (t) at instant t, and at mean rate L ω (t), R for the entire line segment up to time t. We note that in realistic biological models the largest reproduction rates typically belong to types that are improbable under the stationary mutation distribution ρ ('good' types are rare under mutation alone, otherwise it would not require selection to establish them!). Hence, empirical measures with a large mean reproduction rate tend to differ markedly from ρ. The resulting tradeoff between the mean reproduction rate of a line and its asymptotic rate of decay is won by those lines ω for which L ω (t) = ν maximizes the difference, ν, R − I U (ν). According to Thm. 1, these are precisely the lines having the ancestral distribution α as their time average. It is therefore this α that is successful in the long run and that we see when looking back into the past.
5.1.4 Extension to general splitting rules. In our proof of Thm. 1 above, we used a probabilistic argument that relied on the parallel model and the associated tree construction. So it might seem that this theorem is limited to this particular model. Note, however, that all quantities appearing in Thm. 1 are solely determined by the first-moment generator A of the process, so that it is a property of A rather than the underlying process. For an arbitrary Markov branching process, we can use the formal decomposition A = U +R of its first-moment generator to build a parallel model with the same A. Since the theorem holds for the latter process, it also holds for the former; this is some kind of "invariance principle". All that is lost is the probabilistic interpretation given in the previous comment; such an interpretation may be regained with the help of the size-biased tree construction of Subsec. 4.2, but is then more involved.

Reversible mutation rates, and symmetrization
We will now discuss the important special case that U is reversible, in that ρ i U ij = ρ j U ji for all i, j ∈ S. This is assumed in most models of nucleotide evolution, see, e.g., [12,Ch. 13]. The interest in this case comes from the following facts.
5.2.1 Explicit form of the rate function. For reversible U , the maximization in (25) can be carried out explicitly, so that the rate function takes the closed form [20, p. 50, Ex. IV.24] here both the square root and the fraction are to be read componentwise, and u, v ρ denotes the Dirichlet form i u i v i ρ i for vectors u, v, and ρ. (It is an interesting fact that no such simplification exists for reversible Markov chains in discrete time.) Noting that ρ i > 0 for all i ∈ S by irreducibility, using the reversibility in the form ρ i /ρ j U ij = U ij U ji , and recalling that U ii = − j:j =i U ij , one readily finds that Eq. (29) is equivalent to

Estimation of the reproduction rate from the ancestral distribution.
The reversibility of U immediately implies that the vector ρh := (ρ i h i ) i∈S is a left eigenvector of A = U + R for the principal eigenvalue λ, cf. [3]. Hence π = ρh up to a normalization factor, and therefore α = ρh 2 , or h = α/ρ, again up to a normalization factor. (As before, the square root and the fraction are to be read componentwise.) This in turn means that α, together with ρ, determines the reproduction rate R up to an additive constant. Indeed, suppose that R and R ′ are two reproduction rates (for the same mutation matrix U ) having the same ancestral distribution α = α ′ . Then h = h ′ , whence (R − R ′ )h = (λ − λ ′ )h. As h is strictly positive, it follows that all components of R − R ′ agree.

Symmetrized mutation rates.
For reversible U , one can introduce the matrix A := (Ã ij ) i,j∈S byÃ ij = √ ρ i A ij / √ ρ j , which is symmetric and has the same spectrum as A = U + R. The maximum principle of Thm. 1 can therefore also be derived from the Rayleigh-Ritz (or Courant-Fisher) variational principle for the leading eigenvalue of A; see [3, Sec. 2]. We emphasize, however, that the large deviation approach to (26) is not tied to reversible matrices and, as we have shown above, admits a natural interpretation in terms of the underlying family tree. Nevertheless, we will take advantage of the symmetrization A in Sect. 7 below. In particular, we will use the (unique) decomposition A = F + E into a symmetric Markov generator F = (F ij ) i,j∈S , defined through and a diagonal matrix E := diag(E i | i ∈ S) with elements 2 6. Unfolding the variational principle As we have seen, the maximum principle of Thm. 1 provides some general insight into the competition, and resulting tradeoff, between mutation and reproduction. In general, however, it can not be solved explicitly. This is because both the maximization over the space P(S), and the eigenvalue equations determining π, h and thus α, are |S|-dimensional, and S is typically large. It is thus natural to ask whether one can obtain a low-dimensional variational principle in a specific setting. In the rest of this paper we will therefore confine ourselves to genetic models of sequence type where each type is specified by a sequence of letters from a finite alphabet. The variational problem can then be split into two simpler ones, a constrained variational principle with fixed mean letter composition, and a maximization over all possible constraints. In some cases, each of these two subproblems may be treated explicitly or, at least, approximately.

Lumping of sequence types, or: Choice of a type space
As previewed in Subsection 3.2, we will now assume that the type of an individual is characterized by a sequence of letters from some finite alphabet Σ, which leads to the type space Σ N . If we assume that the mutation and reproduction rates are invariant under permutations of sequence sites, as we will do in what follows, this sequence space can be lumped into the smaller space recall Fig. 5. For example, this is possible for sequence space models with parallel mutation and reproduction, in which (L1) all sites mutate independently and according to the same (Markov) process (a natural first assumption made in many models of sequence evolution) and (L2) the fitness function is invariant under permutation of sites (a less natural, but still common assumption that applies, for example, if fitness only depends on the sequence through the number of mutated positions (i.e., the Hamming distance) relative to a reference sequence, often termed the 'wildtype'); see, e.g., [18,14] or [16] for previous work on this case. (As an alternative to the choice (33), one can use the constraint ℓ∈Σ i ℓ = N to remove an element a ∈ Σ by setting Σ * = Σ \ {a} and work instead with where now d = |Σ * | = |Σ| − 1.) Specifically, if the reproduction and mutation rates on Σ N are given by R σ and U στ (σ, τ ∈ Σ N ), then, by permutation invariance, there is a vector R = (R i ) i∈S and a Markov generator U = (U ij ) i,j∈S so that R σ = R H(σ) and τ :H(τ )=j U στ = U H(σ),j for all σ ∈ Σ N ; here H is as in (12). R and U then define a Markovian branching process with type space S. 3 In fact, assumption (L1) even implies that the mutation rates U ij of the lumped model are linear in i ∈ S (or affine in i ∈ S * ). This is seen as follows: If w ℓm is the mutation rate (at every site) from letter ℓ to letter m, then the corresponding transition in the lumped model (based on Σ) is i → i − e ℓ + e m (where e j is the the unit vector in R d having a 1 at coordinate j), and occur at rates i ℓ w ℓm , due to independence of the sites. If, instead, one removes one dimension by setting i a = N − ℓ∈Σ * i ℓ and then works with Σ * , one obtains the additional transitions i → i − e ℓ at rate w ℓa i ℓ , and i → i + e m at the (affine) rate w am (N − ℓ∈Σ * i ℓ ).
Assumption (L2) is less specific than (L1); the fitness function in the lumped model will, in general, be nonlinear due to interactions between sites. It will, however, turn linear (or affine) if fitness contributions are additive across sites, as is usually assumed in, e.g., models of codon bias (where Σ is the set of possible codons). Additivity reflects independent fitness contributions of the sites and means that, for i ∈ S, one has R i = ℓ∈Σ r ℓ i ℓ (based on Σ), or R i = r a N + d ℓ∈Σ * (r ℓ − r a )i ℓ (if Σ * is used), where r ℓ ∈ R for ℓ ∈ Σ. We will examine such linear models in Subsec. 6.3.

Fixing the empirical mean
The only property of the special choices (33) or (34) of the type space S we need at the moment is that S ⊂ R d . This provides S with the structure of an abelian group (elements of S can be added and subtracted), and allows us to classify the possible empirical distributions ν ∈ P(S) according to the value of their mean ν, id ∈ R d . In particular, for the random measure L ω (t) of (21), L ω (t), id is a random vector in R d , namely the empirical mean, or empirical mean letter composition along the line ω up to time t. If S is obtained through lumping a sequence space Σ N as in Fig. 5, the ℓ'th coordinate of L ω (t) indicates the total fraction of time up to t for which some site in the sequence characterizing an individual on the line ω shows letter ℓ ∈ Σ. Note that this involves a twofold averaging, namely an average over time and a (non-normalized) average over sequence sites.
As indicated in (14) and (15), we will now foliate the variational formula (26) by prescribing the mean of the underlying type distribution. That is, where is the constrained mean fitness of z ∈ conv S. As before, the maxima are attained by continuity, and the maximizer in (36) is unique by the strict convexity of I U . The function Λ is strictly concave; this follows again from the strict convexity of I U , together with the linearity of · , R and · , id . In particular, Λ is continuous on the relative interior of conv S [32, p. 82]. In general, the relative interior rint D of a set D ⊂ R d is defined as the interior of D relative to the smallest affine subspace containing D. 4 Moreover, since α is the unique maximizer in (28), there exists a uniqueẑ ∈ conv S that maximizes Λ, namelŷ i.e., the unique maximizerẑ in (35) is the ancestral type average. If U is reversible, we may restrict the maximization in (35) to those z that are strict convex combinations of the elements of S. This is obvious from the explicit form of I U in (30): If at least one component of ν vanishes, one has (∂/∂ν i )I U (ν) = +∞ for some i. Therefore, the maximum will be located in rint conv S, so that Eq. (35) can be replaced by If the function Λ(z) were known explicitly, the variational problem of Thm. 1 would boil down to a maximization over a subset of R d ; for small d one could aim at explicit solutions. Such low-dimensional variational principles for λ were recently derived for several examples, by methods from linear algebra and asymptotic analysis [3,13,14,18]. However, a plausible understanding for the resulting function to be maximized has been lacking so far. The next Proposition reveals the probabilistic meaning of Λ(z): It is nothing but the asymptotic growth rate of the lines having empirical type average (close to) z. Together with (37), this shows that the growth rate of the total population coincides with the growth rate Λ(ẑ) of the subpopulation consisting of all individuals with empirical type average close to the ancestral one.
Proposition 2. For all z ∈ rint conv S, the solution of the constrained variational problem (36) satisfies Proof. Consider first the parallel model. By the reasoning leading to (23), the growth rate of the subpopulation consisting of all individuals with empirical mean close to z, up to some maximal deviation ε > 0, is equal to Here we have used the conventions ∞·1 = ∞ and ∞·0 = 0 in the third step, and Varadhan's lemma in the fourth, in analogy with (26); the maximum over ν is attained since the condition ν, id − z 1 ≤ ε defines a compact subset of P(S). As Λ is continuous on rint conv S, the last expression converges to Λ(z) as ε → 0, as asserted. For a general splitting rule, the argument is the same except that the particular tree construction of Subsec. 4.1 has to be replaced by the sizebiased tree described in Subsec. 4.2. In fact, one simply has to omit the second line of (39) above and instead invoke Eq. (4.4) of [15] which shows that the first line of (39) coincides with the third; the random measure L(t) in the third line is then again the empirical measure of a Markov chain with generator U , namely the mutation process along the spine of the size-biased tree.
⊓ ⊔ Like the unconstrained variational problem (26) leading to λ, the constrained problem (36) defining Λ provides insight into the mutation-reproduction process, but does not, in general, lead to an explicit solution if S is large. ¿From the point of view of explicit calculations, it rather expresses one difficult problem (the leading eigenvalue of a large matrix) in terms of another difficult problem (the maximization over a large space). But if U is reversible, there are two cases in which (36) may be solved explicitly or, at least, asymptotically. These are the cases when fitness and mutation are linear (already hinted at in Sec. 6.1), or when they allow a continuous approximation in the limit as the number N of sequence sites grows large. These will be discussed in the next Subsection and in Section 7.

Exact results for linear reversible models
In this Subsection we have a closer look at the sequence space models of Sec. 6.1 that describe the independent evolution of N sites with a finite alphabet Σ and lead, after lumping, to models with state space S as in (33), with linear fitness and mutation, and mutational transitions i → i + k restricted to those with k ∈ S := {e m − e ℓ | m, ℓ ∈ Σ, m = ℓ}. In line with standard assumptions on sequence evolution (see, e.g., [12, Ch. 14]), we posit that the mutation process acting at the sites is reversible, that is, the mutation rates (w ℓm ) ℓ,m∈Σ define an irreducible and reversible Markov generator with a reversible distribution γ on Σ. After lumping, the associated mutation process on S then has rates U ij = U j−i (i) , i, j ∈ S, given by The reversibility of (w ℓm ) ℓ,m∈Σ readily implies that the mutation generator U = (U ij ) i,j∈S is also reversible; its reversible distribution is ρ := Mult N,γ , the multinomial distribution for N samples from the distribution γ on Σ. As motivated in Sec. 6.1, we will also assume here that the reproduction rates are linear, in that R i = R(i) for all i ∈ S, for a linear function R of the form R(z) = r · z, r, z ∈ R d . Here and below we write '·' for the scalar product of vectors in R d , in contrast to ., . , which we have reserved for scalar products of vectors in R S . In this setting, the constrained variational problem (36) admits an explicit solution as follows. Due to (38), we may -and will -restrict ourselves to considering means in Theorem 2. In the situation described above, for every z ∈ rint conv S the restrained maximum of (36) is given by where ν (z) = Mult N,z/N is the multinomial distribution with mean z.
Proof. Let z ∈ rint convS be given, and consider any ν ∈ P(S) with ν, id = z. It is then clear that ν, R = R(z) by linearity. Let us rewrite Eqn. (30) in the form I U (ν) = 1 2 k∈S x k − y k 2 2 , where x k = (x k,i ) i∈S and y k = (y k,i ) i∈S are the vectors with components i ∈ S, k ∈ S. In the boundary case when i ∈ S but i + k / ∈ S, we have U k (i) = 0 by definition, and likewise U −k (i + k) = 0 when i / ∈ S but i + k ∈ S. Hence x k,i = y k,i = 0 unless i, i + k ∈ S. By linearity of the U k , it follows that x k 2 2 = U k (z) and y k 2 2 = U −k (z) for all k ∈ S. As the distance between any two vectors is minimized when the vectors are parallel, we conclude further that with equality if and only if there is a positive constant C k so that This, however, is the case when ν = ν (z) because, for each i ∈ S, ν (z) i = e β·i ρ i for β = log(z/N γ) (where the fraction and the logarithm are taken componentwise), and ρ is reversible; in fact we have C k = e −β·k . Combining the preceding observations we get the result.

⊓ ⊔
If we turn from the linear model to the affine one, by removing one coordinate as indicated at the close of Sec. 6.1, Theorem 2 clearly remains true, with the middle expression in (40) expressed in terms of the reduced coordinates. Eq. (40) has been derived previously for certain specific choices for the mutation rates [14,18]; remarkably, the above result provides both an extension (to arbitrary reversible models), and a simplification of the proof.

Partial convex conjugation
On our way to the second case of an explicit version of Λ(z), we need a general intermediate step: a relation between Λ(z) and the mean growth rate λ for a suitably modified reproduction rate R. This relationship is based on partial convex conjugation, a standard procedure of convex analysis which will be spelt out here for our purposes. In Sec. 7, this will allow us to determine the asymptotic behaviour of Λ(z) when the number N of sequence sites gets large.
Let us rewrite equation (26) in the form indicating the dependence on R; U will be considered as fixed. The following proposition asserts that the function z → −Λ(z) of constrained extrema is a partial convex conjugate of the function R → λ(R).
Here, β z ∈ R d is the negative slope vector of any tangent plane to Λ at z, and α z is the unique ancestral distribution corresponding to the reproduction rate R + β z · id for any such β z . In particular, the function β → λ(R + β · id) is differentiable on R d , and ∇ β λ(R + β · id) | β=β z = α z , id = z.
Proof. For any z ∈ rint conv S and β ∈ R d we have, writing ν z for the maximizer in (36) and using Thm. 1, Taking the infimum over β we arrive at the inequality To show equality we recall that Λ is strictly concave and finite on a (relative) neigbourhood of z and therefore admits a tangent plane at z. That is, there exists some β ∈ R d such that with strict inequality for y = z. Denoting by α β the ancestral distribution for the reproduction rate R + β · id and letting y = α β , id we find Together with (41) it follows that equality holds everywhere in (42). Hence y = z, (41) holds with equality, and the infimum is attained for any β determining a tangent to Λ at z. In general, there may be several such tangents, e.g., if S is contained in a hyperplane of R d . However, the associated ancestral distribution is uniquely determined. For, suppose there exist β 1 = β 2 both determining a tangent to Λ at z, and let α 1 and α 2 be the ancestral distributions for the reproduction rates R + β 1 · id and R + β 2 · id, respectively. The preceding argument then holds for every β in the segment [β 1 , β 2 ], whence (42) holds with equality everywhere for all these β. We can thus conclude that the function β → λ(R + β · id) is affine on [β 1 , β 2 ]. In particular, using Thm. 1 and the shorthand f i (ν) = ν, R + β i · id − I U (ν) we find Since f 1 and f 2 are strictly concave, this is only possible if they have the same maximizer. That is, α 1 = α 2 . Finally, using the equality in (41) and the convex duality lemma [7, Lemma 4.5.8] we find that the function β → λ(R + β · id) is the convex conjugate of the strictly convex function −Λ, and thus differentiable; see [32,Thm. 26.3,p.253]. Its gradient at β z necessarily coincides with z.

⊓ ⊔
In the case of a reversible mutation matrix U , the preceding proposition can be complemented as follows. We write T = span(S − S) ⊂ R d for the linear space generated by the set of differences of elements of S. Corollary 1. For reversible U the following additional statements hold. (a) The function Λ defined in (36) is differentiable 5 on rint conv S, and its conjugate function β → λ(R + β · id) is strictly convex on T . Moreover, for z ∈ rint conv S and β ∈ T we have β = −∇Λ(z) if and only if z = ∇ β λ(R + β · id). (b) The function Λ on rint conv S remains unchanged under symmetrization, i.e., by replacing U with the matrix F of (31), and R with the function E defined in (32).
Proof. (a) Let z ∈ rint conv S and β, β ′ be two negative slope vectors of Λ at z. In view of the uniqueness of α z and the remarks in 5.2.2, the scalar product (β − β ′ ) · i is then independent of i ∈ S. This means that β − β ′ is orthogonal to T , so that there is a unique negative slope vector β z ∈ T . By concavity, the uniqueness of the tangent plane is equivalent to differentiability; cf. [32,Thm. 25.1,p. 242]. By the proof of Prop. 3, this is also equivalent to strict convexity of λ(R + β · id) on T . The final statement comes from the observation that both assertions are equivalent to the iden- , so that their principal eigenvalues agree. The result thus follows from Prop. 3 by minimization over β. ⊓ ⊔

Smooth approximations
While still adhering to a lumped sequence model, we will now turn to a situation complementary to that of Thm. 2: we consider nonlinear reproduction and mutation rates that allow for a continuous approximation if the number of sequence sites becomes large; this approximation is only required locally, which provides much more freedom, and, in particular, removes constraints imposed by the boundary (recall the boundary conditions U k (i) = 0 for i ∈ S, i + k / ∈ S in Thm. 2). For a large family of models with reversible U , an asymptotic low-dimensional maximum principle for λ is available then [3], but no connection to the constrained mean fitness (36) has been made there, and the ensuing probabilistic interpretation was still lacking. On the basis of Prop. 3, this can now be provided.
In view of Corollary 1(b), the case of a reversible mutation matrix U can be reduced to the case of a symmetric mutation matrix F . That is, instead of the first moment generator A = U + R we can and will consider the symmetrized version A = E + F defined in (31) and (32). In Subsection 7.1 we will present a slight refinement of an asymptotic maximum principle derived in [3,Thm. 1]. In Subsection 7.2 we will derive an approximation of Λ(z) in two particularly interesting situations. An application to the quasispecies model follows in Sec. 8.

Approximation of the asymptotic growth rate λ
Consider the following setup. For each N let -S = S(N ) ⊂ Z d be a state space as in (33) or (34).
The rescaled set 1 N S is then contained in a simplex D ⊂ R d , viz. either D = conv {e 1 , . . . , e d } or D = conv {0, e 1 , . . . , e d }, with e 1 , . . . , e d the unit vectors of Z d . (In the first case, D is contained in a hyperplane, whence in the following we will always consider the relative interior of D rather than simply its interior.) In the limit as N → ∞, 1 N S becomes dense in D. For each N let also -F be a symmetric Markov generator on S, and E := diag(E i | i ∈ S) a diagonal matrix.
We assume that F and E admit a continuous approximation as follows: There exist real functions e and f k on D, and an "approximation domain" A ⊂ rint D such that the following conditions hold.
(A1) e is C 2 on A and, as N → ∞,  We postpone the proof until Subsection 7.3, discussing first the significance of the assumptions and the result.
7.1.1 Formal comments. The above approximation for the principal eigenvalue of A clearly also holds for the similar matrix A = R + U . Note also that only the function e remains relevant in the limit; the f k play no role. This means that A = E +F provides the 'right' decomposition into the 'relevant' E-term, and an F -term whose contribution to the leading eigenvalue vanishes in the limit.
It is also interesting to observe that the approximation assumption (A1) is only required in a neigbourhood A of a single maximizer z * ∈ rint D of e; further maxima may appear, even on the boundary, but these do not matter. This locality of the approximation domain is the main difference to Thm. 1 of [3] which requires a globally uniform approximation. (As the example of linear mutation in Thm. 2 shows, it often happens that the derivatives diverge at the relative boundary of conv S, so that a global approximation is not feasible. This is also the case for the quasispecies model considered in Sec. 8.) As a global requirement we need only the bounds in (A3).

Significance of the assumptions for the model.
Our setup implies that replacing i ∈ S by i/N ∈ S/N will yield a continuous type variable z ∈ D in the limit. Accordingly, the matrix elements are required to become smooth functions of z as N → ∞ -at least locally, in line with (A1).
Condition (A2) says that the mutation rates must decay fast enough with distance to the target type -again, at least locally. This assumption may appear to be rather special at first sight, but actually it is very natural: As we have seen in Sec. 6.1, independent mutation at the sites of a sequence leads to nearest-neighbour mutation on S, hence (A2) is trivially fulfilled.
For the corresponding quasispecies model (to be described below), still with independent mutation at the sites, the decay of f k with k is exponential, rather than only cubic as required in (A2); this will be shown in Sec. 8.
In many concrete examples, the reproduction and mutation rates have their own continuous approximations each, i.e., with C 2 (D, R) functions r and u k . Moreover, the range of all mutational steps is finite (on S, and independently of N ); that is, there is a finite symmetric (i.e., S = −S) set S ⊂ Z d with the property that, for all N , U ij = 0 whenever j − i / ∈ S. Then (A1) is automatically satisfied for any A on which u k (z) is C 2 for all k ∈ S; inspecting the matrix elements of E in (32) and noting that k∈S u k (z) = 0 one finds that It is interesting to observe that the expression above is formally identical with Λ(z) of Thm. 2, although we are considering quite a different situation here. Special cases of (43) have appeared in [3,18] in the context of parallel sequence space models, and the resulting maximum principle turned out as a key to determine the mutation load, genetic variance, and the existence of error thresholds. 7.1.3 Locality of the ancestral distribution. Under the additional (but generic) assumptions that the function e admits a unique maximizer z * ∈ rint D and the Hessian of e at z * , restricted to T = span(D − D), is (strictly) negative definite, one can also characterize the ancestral distribution, which is connected to λ through the general variational principle of Thm. 1. Namely, by Thm. 2 of [3], this distribution is concentrated in a neighbourhood of z * whose width decreases with 1/ √ N . More precisely: For every 0 < ε ≤ 1, there is a constant c > 0, independent of N , so that, for N large enough, By Cor. 3 of [3], it follows that α, id = z * + O(N −1/3 ), i.e., the ancestral type average coincides with the unique maximizer of e up to a small error term. The constant in the error term here depends on those in the assumptions and some bounds separating the spectrum of the Hessian (restricted to T ) of e at z * from −∞ and 0. The proofs given in [3] are solely based on a local approximation and thus remain valid under our weaker assumptions.

Approximations of the constrained mean fitness Λ
Our next goal is an approximation for the partial maximum Λ(z) of (36). In fact, the similarity of the expression (43) for e(z) and the expression for Λ(z) in Thm. 2 leads one to ask whether the asymptotic identity of the global maxima of Λ and e, as asserted by Theorem 3, can be extended to an asymptotic relation between these functions as a whole. On the basis of Prop. 3 such an approximation can indeed be given. We consider first the most salient points of e, i.e., the points where e coincides with its concave envelope. Let us say z ∈ rint D is an exposed smoothness point of e if e(y) < t z (y) := e(z) + ∇e(z) · (y − z) for all y = z, i.e., z is the unique point where e hits its tangent plane t z at z. -e is C 2 on a neighbourhood of z, and the Hessian of e at z, as a bilinear form on T = span(D − D), is negative definite.
If e is strictly concave on D, the first condition is trivially satisfied. If e is also C 2 , the second condition just covers the generic case of strict concavity. In other words, for a generic strictly concave C 2 -function e, every z ∈ rint D is an exposed smoothness point.
To state the hypotheses of the next theorem we recall that the assumptions (A1) and (A2) only involve an approximation on a local set A, while (A3) imposes some global bounds, including an upper bound on E i in terms of sup e. We now replace the constant sup e by suitable tangent planes of e, thereby turning (A3) into the hypothesis (A3 ′ ) For all z ∈ A and suitable constants C ′ , C ′′ < ∞ we have  The proofs of this and the subsequent theorem follow in the next subsection.
Theorem 4 raises the question of what happens if e touches a tangent plane at two or more distinct points of int D. Let z be a strict convex combination of these points, and β the negative slope of this plane. The ancestral distributionᾱ for the reproduction rateĒ = E +β ·id/N is then expected to split into distinct peaks located at the competing maximum points of the associatedē; its mean type ᾱ, id will remain close to z, but the reproduction rate Λ(z) will be the corresponding convex combination of the values of e at the maximum points ofē. So it may be conjectured that, in general, Λ(z) is approximated by the concave envelope of e at z. The next theorem shows that this is indeed the case. We note that this kind of behaviour is related to the phenomenon of error thresholds and phase transitions described in detail in [18].
For a given function e on D we let e(z) := inf a − β · z a ∈ R, β ∈ R d , a − β · y ≥ e(y) ∀ y ∈ D , z ∈ D, be the concave envelope of e. (For an example see Fig. 8.) We consider the situation when e deviates from strict concavity, so thatê is affine on a nontrivial set B. Let us say that B is a basin of e if B has nonempty relative interior and for suitable a ∈ R, β ∈ R d . Note that a basin B is necessarily convex and compact. We write ex B for the set of its extremal points. Let us say that a basin B of e is determined by smooth hills of e if there exists a relatively open neigbourhood H of ex B in D such that -H \ B consists of exposed points of e, and e is C 2 on H, and its Hessian (restricted to T ) is negative definite with a spectrum which is bounded from below and bounded away from zero uniformly on H.
This is the situation one typically encounters when a smooth e deviates from strict concavity. The theorem below provides an approximation of the restrained maximum Λ(z) defined in (36). uniformly for all z ∈ B.

Proofs
We now turn to the proofs of the three Theorems of this section.
In the proof of Prop. 4 of the original article, starting from the second display (p. 97), we split the sum into We now note that the display in the middle of p. 97 implies that, for |i/N − By the assertion on the error terms in Thm. 3 and in paragraph 7.1.3, the error term here is locally uniform in z.
Next we note that the (N -dependent) mapping φ : z →z from A into D is a homeomorphism. For, φ is the composition of −∇e and β → ∇ β λ(E + β · id/N ). Now, ∇e is a diffeomorphism from A into T = span(D − D) because, by assumption, the Hessian of e (restricted to T ) is nondegenerate everywhere on the convex set A, so that ∇e(x) = ∇e(y) only if x = y by the mean value theorem. On the other hand, Corollary 1(a) shows that ∇ β λ(E + β · id/N ), as a function of β ∈ T , has the inverse y → −∇Λ(N y); these gradients are continuous by Corollary 25.5.1 of [32]. Now let C ⊂ A be compact and C ′ ⊂ A, say, a convex polytope containing C in its relative interior. φ moves the faces of C ′ by at most a distance of κN −1/3 , for some constant κ < ∞. Hence φ(C ′ ) ⊃ C for large N . For these N we can invert φ on C to get φ −1 (y) = y + O(N −1/3 ) uniformly for all y ∈ C. Since |∇e| is bounded on C, it follows that uniformly for all y ∈ C.
Next we fix some ε > 0. By hypothesis, for each k = 1, . . . , ℓ we can find a point y k ∈ H \ B and a relatively open convex neighbourhood A k of y k such that |y k − z k | < ε and A k consists of exposed smoothness points of e. Theorem 4 thus asserts that Λ(y k ) = e(y k ) + O(N −1/3 ). In view of the assumed uniform bounds on the spectrum of the Hessians, the error term here is independent of k and the choice of y k . Letting ε → 0, we thus can conclude that Λ(z k ) = e(z k ) + O(N −1/3 ), and therefore by concavity On the other hand, since the upper estimate on E i in (A3 ′ ) also holds for z = z k , assumptions (A1)-(A3) hold withĒ i = E i +β ·i/N andē := e+β ·id in place of E i and e, respectively; here β is as in (45). Prop. 3 and Thm. 3 therefore imply that, for each k, where λ(Ē) stands for the principal eigenvalue of the matrixĒ + F with reproduction rateĒ. Taking the average over k we find The proof is therefore complete. ⊓ ⊔ 8. Application to the quasispecies model

The model and its large-N asymptotics
We will now illustrate and apply the results of the preceding Section to the coupled counterpart of the parallel sequence space model of Subsec. 6.1.
The coupled sequence space model, known as the quasispecies model, was introduced in [9] and has, since then, been the subject of numerous investigations. It assumes that mutations occur on the occasion of reproduction events, that is, they represent replication errors. Let us assume that mutation is, again, independent across sites and occurs at probabilities v = µ/N and w = ν/N from 0 to 1 and vice versa, where µ and ν are positive and independent of N . This is a slight generalization of the original model [9] with symmetric mutation, and the factor 1/N in the mutation rate is introduced to obtain a suitable limit 6 . The matrix of mutation probabilities, P = (P στ ) σ,τ ∈Σ N , is then given by where the tensor product reflects the independence across sites. The quasispecies model is complete if we further specify birth rates B σ and death rates D σ for all σ ∈ Σ. When a birth event occurs to a σ individual, it survives unchanged and produces an offspring of type τ with probability P στ ; at a death event, a σ individual dies (as in Fig. 2 with i, j replaced by σ, τ ). We will assume that, for all σ ∈ Σ, B σ and D σ are invariant under permutation of sites. Since the same holds, by construction, for the mutation probabilities (47), we have a situation analogous to (L1) and (L2) for the parallel model, and may perform lumping into S := {0, 1, . . . , N } by the mapping σ → H(σ) ∈ S, where H(σ) is the number of sites occupied by letter 1 (see Sec. 6.1). The resulting model on S has birth rates B i , death rates D i , and mutation probabilities P ij , where B σ = B i , D σ = D i , and τ :H(τ )=j for any σ with H(σ) = i. In the lumped model, given the current type i, the distribution of jumps is obviously given by the convolution where Bin n,p denotes the binomial distribution with parameters n, p > 0, and Bin n,p its image under the reflection of Z at the origin; we further identify Bin 0,p with the point measure located at 0. Explicitly, The Markov chain so defined is reversible with respect to ϕ = (ϕ i ) i∈S = Bin N,v/(v+w) ; this is most easily seen by noting that P (on sequence space) is reversible with respect to the Bernoulli measure on {0, The lumped Markov branching process has first-moment generator A with elements A ij = B i P ij − D i (cf. Eq. (4)), and has been much studied, see [10] for a review of early work, and [24] for a review of recent theoretical developments, and their connection to experimental results on virus evolution. In particular, the error thresholds displayed by this model have attracted a lot of attention.
The function e(z) that would simplify the model's analysis does not seem to have appeared so far; it is far less obvious than its parallel counterpart (43), and will be established in what follows. We start by decomposing A of (4) into a Markov generator U and a diagonal matrix R, which gives U ij = B i P ij for i = j, U ii = −B i (1 − P ii ), and R i = j∈S A ij = B i − D i . Since P = (P ij ) i,j∈S is reversible, U is also reversible; its reversible distribution ρ is given by ρ i = cϕ i /B i for a normalizing constant c > 0. The elements of the symmetrized matrices E and F of (32) and (31) therefore emerge as and After these preparations, let us identify conditions under which Thms. 3, 4, and 5 are applicable. We will consider the approximation of the birth and death rates as given; we will then show that the Poisson approximation to the distribution P i,i+• , namely p • (i/N ) = Poi µ(N −i)/N * Poi νi/N , will also lead to the 'right' approximation to the matrix elements (51)-(53). In line with previous notation, Poi λ is the Poisson distribution with parameter λ > 0, Poi λ its reflected version, and Poi 0 is identified with the point measure at 0. This will give us the following result. Postponing the proof for a moment, let us first look at an example.

An example.
For the purpose of illustration, let us consider the quasispecies model with a 'smoothed' version of truncation selection (where a gene tolerates a certain number of mutations and then deteriorates rapidly). Let the birth and death rate functions be given by (i.e., we assume a mixture of fecundity and viability selection). Fig. 8 shows the fitness function, and the function e together with its concave envelope.

8.1.2
Connection to the parallel model. The quasispecies model is closely related to the lumped parallel sequence space model with birth rates B i , death rates D i , and mutation rates U i,i+1 = µ(1 − i/N ), U i,i−1 = ν i/N and U ij = 0 for |j − i| > 1 (where µ and ν are now mutation rates per site rather than probabilities). In fact, the latter may be considered as the former's weak-selection weak-mutation limit (cf. [5, Ch.II.1.2], and [19]). It leads to the simpler expression cf. (43), and [3,18]. Indeed, this function is easily identified as the weakselection weak-mutation limit of (55) by replacing b by 1 + δb, d by 1 + δd, µ by µδ, ν by νδ, and e by e/δ; the last replacement means that time is measured in units of δ. e(z) of (57) then emerges from (55) in the limit δ → 0.

Proof
The proof of Thm. 6 consists in verifying the assumptions of Thms. 3, 4 and 5 for the matrices E and F in (51)-(53). The main difficulty will be to establish the O(1/N ) approximation as required in (A1). Besides the Proof. If a = 0 or b = 0, p k (a, b)p −k (a, b) = 0 except for k = 0, so that the assertion is trivial. For a, b > 0 we can write, using Lemma 1 and Markov's inequality: By symmetry, the last sum is the variance of Poi √ ab * Poi √ ab and thus equal to 2 √ ab.

⊓ ⊔
We need a similar tail estimate for the geometric symmetrization of the matrix P defined in (50). Note that P depends on N .
for a constant C depending on µ and ν but not on N .

⊓ ⊔
The crucial step is the following Poisson approximation of the geometric symmetrization of P . Proof. Consider an arbitrary i ∈ S and suppose z := i/N is bounded away from 0 and 1. (We will generally suppress the i-dependence of all abbreviations to be introduced below.) The main difficulty of the proof comes from the fact that the parameters of the probabilities P ji depend on the variable j rather than i. Fortunately, P is reversible, and Lemmas 2 and 3 allow us to confine ourselves to the j's with |j − i| ≤ √ N . We proceed by a comparison of upper and lower bounds on r k := P i,i+k P i+k,i , k ∈ S − i .
Step 1: A lower estimate. Since P is reversible w.r.t. Bin N,v/(v+w) , we have for all k ∈ S − i. Since also m!/n! ≥ n m−n for all m, n, it follows that r k ≥ r k := P i,i+k s k for all k, where s := (iw/(N − i)v) 1/2 = (zν/(1 − z)µ) 1/2 . We now take the sum over k. Using (50) and the binomial formula we can write The inequality follows from the fact that for any n ≥ 1 and |x| ≤ n; see [31, 3.6.2., p. 266]. In the last step we used that z is bounded away from 0 and 1.
Step 2. An upper bound. Arguing as for (67) we find where ψ(x) := (1+x/z) 1/2 (1−x/(1−z)) −1/2 . Next, for each k an expansion of ψ gives for some ϑ k between 0 and k/N . As long as z is bounded away from 0 and 1 and |k| ≤ √ N , the ψand ψ ′ -expressions on the right-hand side are bounded from above, so that ψ(k/N ) k ≤ 1 + k 2 c k N for some c < ∞. Also, using (64) we find Collecting all estimates we arrive at the upper bound r k ≤r k := e va+wb p k (a, b)s k 1 + k 2 c k N .
Next, a summation over k gives Step 3: Conclusion. Consider now q k := p k (z)p −k (z). By (61), q k = p k (z)s k . It is also immediate that p k (z) ≤ e va+wb p k (a, b) with a, b as in (70). Hence q k ≤r k for all k. Combining this with Lemmas 2 and 3 and the bounds (67) and (71) we find Now, (68) and (72) show that the first sum exceeds the second only by a term of order 1/N , and (72) together with (63) and Lemma 2 imply that the first sum exceeds the third by at most a term of order 1/N . This completes the proof.

⊓ ⊔
Besides the preceding key approximation in the interior of D, we also need a uniform bound which will be used close to the boundary of D. Here, an error bound of order 1/ √ N is sufficient.
where the O(1/ √ N ) term depends on µ and ν but not on i.
Proof. In view of Lemma 3, we only need to estimate the sum over all j = i + k ∈ S with |k| ≤ √ N . For these j we find, using the Poisson bound (64) and writing again z = i/N , k (a, b).
so that eã +b = 1 + O(1/ √ N ); the error terms do not depend on i. The claim thus follows from Lemmas 3, 2, 1 and the fact that ( After these preparations we are now ready to read off the approximating functions e and f k for the lumped quasispecies model, that is, we can proceed to the Proof (of Thm. 6). The main point of the proof is to establish condition (A1) for any compact interval A ⊂ ]0, 1[. The (asymptotic) boundedness of the B j 's and Lemma 3 imply that, for each i ∈ S, (73) The asymptotics of the B j 's implies further that the sum on the right-hand side is equal to where z := i/N and β(x) := (log b(x))/2. By hypothesis, β ∈ C 2 b ([0, 1]). Hence β(z + k/N ) − β(z) = β ′ (z)k/N + O(1/N ) for |k| ≤ √ N , so that the last expression takes the form where δ := β ′ (z). Next we can omit the exponential e kδ/N , making an error of order 1/N only. Indeed, using inequality (65) The second inequality is obtained by taking formula (60) for p k ( √ ab, √ ab), using e |k||δ|/N ≤ e ℓ|δ|/N e m|δ|/N , and summing up. If z = i/N is bounded away from 0 and 1, we can finally apply Prop. 4, Lemma 2 and the identity (63) to obtain Taking this together with the assumed approximation of the D i , we arrive at the approximation (A1) of E, and the diagonal elements of F , by the functions e and f k defined in (55) and (58). But the approximation of the nondiagonal elements of F is also guaranteed, since the estimates leading to (77) all hold term by term. This completes the proof of (A1). Next, condition (A2) follows directly from Lemma 1 and the fact that Poisson distributions and their convolutions have a finite third moment. We note further that the upper bound on F ij and the lower bound on E i in (A3) and (A3 ′ ) are obvious.
Before turning to the upper bound on E i let us discuss the particular context of Theorems 3, 4 and 5. We observe that the function g is continuous on [0, 1] and smooth on ]0, 1[ with g ′ (0) = −g ′ (1) = −∞, while b > 0 and d are C 2 functions on [0, 1]. This entails that e is C 2 on ]0, 1[ and attains its absolute maximum at a point z * ∈ ]0, 1[; in particular, z * is contained in an interval A satisfying (A1) and (A2), as is required for Thm. 3. In addition, e ′′ is negative in a neigbourhood of 0 and 1, and has only finitely many zeroes by assumption. This implies that each basin of e is determined by smooth hills, as is necessary for applying Thm. 5. Now let A ′ be any set of exposed smoothness points of e which is bounded away from 0 and 1. If δ > 0 is sufficiently small, the set A ′ δ of all y satisfying e(y) ≥ t z (y) − δ for all z ∈ A ′ is still bounded away from 0 and 1. For all i with i/N ∈ A ′ δ , the upper bound on E i in (A3 ′ ) follows directly from (77). For all other i's, this bound follows from (73)-(76) as soon as N is so large that the O(1/ √ N )-term in Prop. 5 is less than δ/ max b. This completes the proof of (A3 ′ ) under the conditions of Thms. 4 and 5. Since ]0, 1[ splits into finitely many intervals forming basins and smooth hills of e, the stated approximation result for Λ follows. Finally, as z * is also an exposed point, the choice A ′ = {z * } gives us the upper bound on E i in (A3). Thm. 3 can therefore be applied, proving the approximation of λ as stated.