On the establishment of a mutant.

How long does it take for an initially advantageous mutant to establish itself in a resident population, and what does the population composition look like then? We approach these questions in the framework of the so called Bare Bones evolution model (Klebaner et al. in J Biol Dyn 5(2):147-162, 2011. https://doi.org/10.1080/17513758.2010.506041) that provides a simplified approach to the adaptive population dynamics of binary splitting cells. As the mutant population grows, cell division becomes less probable, and it may in fact turn less likely than that of residents. Our analysis rests on the assumption of the process starting from resident populations, with sizes proportional to a large carrying capacity K. Actually, we assume carrying capacities to be [Formula: see text] and [Formula: see text] for the resident and the mutant populations, respectively, and study the dynamics for [Formula: see text]. We find conditions for the mutant to be successful in establishing itself alongside the resident. The time it takes turns out to be proportional to [Formula: see text]. We introduce the time of establishment through the asymptotic behaviour of the stochastic nonlinear dynamics describing the evolution, and show that it is indeed [Formula: see text], where [Formula: see text] is twice the probability of successful division of the mutant at its appearance. Looking at the composition of the population, at times [Formula: see text], we find that the densities (i.e. sizes relative to carrying capacities) of both populations follow closely the corresponding two dimensional nonlinear deterministic dynamics that starts at a random point. We characterise this random initial condition in terms of the scaling limit of the corresponding dynamics, and the limit of the properly scaled initial binary splitting process of the mutant. The deterministic approximation with random initial condition is in fact valid asymptotically at all times [Formula: see text] with [Formula: see text].

Introduction

There has been much work in stochastic adaptive dynamics and evolutionary branching, see Dieckmann and Law (1996), Metz et al. (1996), Champagnat et al. (2002), Champagnat and Méléard (2011) and Sagitov et al. (2013), to mention just a few. Here we confine ourselves to a simple mathematical model for evolution, where an established resident population is invaded by a mutant. From that moment on, the two populations compete for resources. At the moment of invasion the resident, wild-type, population is assumed to have the size near its carrying capacity . Here K should be thought of as large, and is fixed. The size of the mutant population is initially negligible as compared to K, since it starts from one individual. It has a reproductive advantage over the resident, but as its progeny grows this advantage diminishes.

We want to answer the question of how long it takes for a mutant to become established, i.e. to grow to a size comparable to the host population. And what is the population composition then? Already the simplified model of two competing populations we consider, will require new mathematical techniques and lead to insightful results. We show that the deterministic approximation with a random initial condition is valid for times with any fixed and a large K. However, unlike in the classical case on deterministic approximation (Kurtz 1970; Barbour 1980), some stochasticity remains and enters as a random initial condition.

The Bare Bones evolutionary model

This simple but basic model of species reproducing under interaction with their environment was introduced in Klebaner et al. (2011). It builds upon asexual binary splitting and evolves in discrete time. Thus, each individual either gets two children in the next generation or none. However, interaction with environment and population size is allowed—in contrast to classical stochastic approaches—but drastically condensed. Following the idea of Malthus, populations reach sizes proportional to available resources, and we assume that the the habitat is characterised by a carrying capacity, , thought of as large. Given the population size, individuals reproduce independently. Initially only the resident, wild-type, population is present and, at population size z, the individual probability of successful splitting is taken to be . Here is a constant, which determines the population size at its macroscopic (quasi-)equilibrium: when , the probability of splitting is 1/2. On the average, thus, a population of this size produces one child per individual. As a result, the population size fluctuates around this (quasi-)steady state for what is presumably a very long time, cf. Jagers and Klebaner (2011).

In that stage, the population will experience its first mutation giving rise to a new population. The new, mutant population starts from a single individual, its ancestor. The basis of adaptive dynamics can then be said to be furnished by the branching mechanism, which forces the new population to either die out or else embark on exponential growth, in which case the old resident dies out, or the two populations will coexist for a time span that turns out to be exponential in the carrying capacity.

Mathematically, this dynamics can be described as follows. The branching process starts from a pair of positive integers , the first component denoting the size of the resident and the second that of the mutant population, at time 0, when the mutation appears. We assume that the established original population is at equilibrium at the moment of invasion , , and . Each population develops by binary splitting with probabilities dependent on the numbers of cells, with transitions from generation n to described by the recursionThe random variables are independent, given the preceding, and only depend upon the last generation , with probabilitieswhere is the parameter, which controls the mutant equilibrium population size, and is the interaction coefficient, assumed to satisfy . The biological meaning of is that cells of one type encroach less upon the reproduction of the other cell type than do cells of the same type. That is the same in both probabilities means that influence is symmetric between the cell types.

In the absence of mutants, the established population thus has a critical reproduction, whereas the mutant population starts supercritically, provided , as is assumed throughout the paper, see (C) below.

Stochastic nonlinear dynamics for the evolution of the density

Important insights into the behaviour of populations with state dependent reproduction and large carrying capacity is provided by their density process (Klebaner 1984, 1993). It allows representation of the process as having stochastic nonlinear dynamics, which can be separated into a deterministic part and a random perturbation. This is useful not only for the mathematical analysis but also for the biological interpretation.

The density process is the population sizes relative to KNote that the splitting probabilities (and hence the offspring distributions) in (2) are in fact functions of the density; denoting the density state by we see thatAccordingly, the offspring mean at is also a function of the densityThe underlying deterministic dynamicsis determined by the function ,This can be easily seen from (1) by writing the density process asThe first term on the r.h.s. of (5) gives the deterministic dynamics (3), and the second term acts as the random perturbation,withThese random variables have zero mean and variance , where are the splitting probabilities. Therefore the random noise term in (6) is of order and the density process can indeed be viewed as generated by a nonlinear dynamical system, perturbed by a small random disturbance.

Note that in the view of the above discussion, the trajectory of the deterministic system (3) depends on K through the initial condition . Similarly, the process generated by the stochastic dynamics (6), depends on K through and the noise term. Whenever appropriate, we will leave this dependence implicit, omitting it from the notation.

Deterministic dynamics

If we neglect the small random noise in (6), we obtain the deterministic dynamics (3). Fixed points (solutions to ) play an important role in the behaviour of such systems. The trajectories are repelled from the unstable fixed points and attracted to the stable ones. Our system, generated by the function in (4), has four fixed points,Since we are concerned with both populations, the relevant case is when both coordinates of are nonnegative. This is true if the following co-existence condition holdsIt is easy to see by examining the Jacobian matrix , see (22) below, that the point is stable, and unstable. The points and are saddle points, that is, stable in one direction and unstable in another. In our theory the point plays a special role due to proximity of the initial condition . In the absence of a mutant, is the stable equilibrium for the resident population, and 0 is unstable for the mutant population.

The large capacity limit of the stochastic dynamics

A rigorous treatment for neglecting small noise is given by the classical results in perturbation theory of dynamical systems, see e.g. Kurtz (1970), Barbour (1980), Freidlin and Wentzell (1984) and Kifer (1988). They assert that as the noise converges to zero, that is, when , the trajectory of the stochastic system (6) converges on any bounded time interval to that of the deterministic dynamics (3), started from the initial condition . Namely, for an arbitrary but fixed integer N,In our setup, the initial condition turns out to be the fixed point ,Therefore the corresponding limit trajectory is constant, for all Consequently, the limit (8) fails to provide any information on the transition to a new coexistence equilibrium. We shall see that if such a transition occurs, it becomes visible much later, at a time increasing with K, in fact, of order .

Recently, limit theorems, capable of capturing this transition, were obtained in Barbour et al. (2015, 2016), Chigansky et al. (2018) and Baker et al. (2018). They involve a time shift which grows logarithmically in K. In Barbour et al. (2015) this shift is random and the process is approximated by the trajectory of the deterministic system (3) with a random shift. We have learnt from a referee that a precursor to random shift theory in Barbour et al. (2015) in the context of epidemic models can be found in Metz (1978), where precise conjectures were stated and later proved in an unpublished manuscript for the simple SIR epidemic model, (Altmann 1993; Mollison 1995).

In Barbour et al. (2016), Chigansky et al. (2018) and Baker et al. (2018), the shift is deterministic, and converges to a trajectory of (3), started from a random initial condition. While the two approaches, the random shift and the random initial condition, are related, they are not equivalent. The main building block in the random initial condition theory is a certain scaling limit of the deterministic flow, which does not appear in the random shift theory. Existence of this limit was so far established only in the one dimensional case.

This work is the first such result in two dimensions. Having established it, we can complement the “random shift” picture in Barbour et al. (2015) with that of “random initial condition” for the Bare Bones model. Recently heuristics for similar random initial conditions for selective sweeps in large populations in one dimension were given in Martin and Lambert (2015). Other stochastic approaches involving carrying capacity can be found in Lambert (2005, 2006).

Main results

In what follows we consider the stochastic process generated by (6) or, equivalently, by (1). As mentioned in Introduction, the resident population initially has a critical reproduction, and is at equilibrium, when a single mutant appears, so that . Even though the probability of a mutant present at any time n is positive, , we do not say that it established itself until its numbers are proportional to its carrying capacity, in other words proportional to K. This can be formalized asFor example, as we have seen above for any fixed n as . This conveys that the mutant is not established by any fixed time n. We show however, that it may establish itself at a time, which grows logarithmically with K. More precisely, we prove that at time with a certain constant b,whereas for ,and, therefore, in particular.

The logarithmic order of time of the mutant’s establishment can be roughly explained as follows. As the process starts near , the state dependent splitting probabilities can be approximated, at least initially, by their values at , giving probabilities of division 1/2 and for the resident and the mutant populations respectively. Note that due to coexistence condition (C), the mutant process is supercritical with meanHence it grows at the rate , and it takes timefor it to grow to the size proportional to K, as . In fact, this heuristics is correct, and made precise in the following result, which implies both (9) and (10). We denote the fractional part of by .

Theorem 1

There exist a non-degenerate scalar random variable and a function , whose entries are positive on the open half-plane , such thatIn particular, along the subsequence of exact powers , ,

Let us now detail about the random variable W and the function appearing in this theorem. The approximate mutant process, mentioned in the heuristic explanation above, has the same splitting probability as the mutant component of at . More precisely, it is a supercritical Galton-Watson binary splitting, started with a single ancestor, , and for defined iteratively bywhere the offsprings are i.i.d. random variables with the constant splitting probability .

It is well known that is a non-negative martingale. As such it converges almost surely to a limit,which is the random variable appearing in (11).

The function in Theorem 1 is the limit of the n-fold iterated map along the unstable manifold of the dynamics in (3).

Theorem 2

Under the basic assumptions stated, the limitexists, and the convergence is uniform on compacts.

Remark 1

(1) It is easy to see that solves the Abel functional equation , subject to . While much is known of such equations in one dimension, in higher dimensions the theory is more involved.

(2) Numerical calculations indicate that is constant with respect to the resident population component , see Fig. 1. This is consistent with the criticality of that population at the density : the global stability of the monomorphic dynamics makes those perturbations shrink to 0 when is iterated.

Fig. 1

Numerical approximation of the function

The next result describes the density process after establishment of the mutant, at times , ; it shows that the population composition is governed by the deterministic nonlinear dynamics started at the random point, as illustrated on Fig. 2. Furthermore, it equally holds when n is a negative integer. Denote the random vector appearing in Theorem 1 by

Fig. 2

The density process versus the approximation with the random initial condition

Corollary 1

For any ,

The next corollary to Theorem 1 answers the question what is the probability of the successful establishment of the mutant? Since the argument is short we present it here. It is known that is exactly the extinction probability of the Galton-Watson process Y(n). It is easily calculated to be . But on the event , . Since on the complimentary event , and for , we have the following corollary of (11).

Corollary 2

With probability the mutant establishes itself alongside the large original population.

Proofs

A preview

The proof is inspired by the observation that supercritical populations, which start from a small number of individuals and develop on a habitat with a large but bounded capacity, grow initially as the Galton-Watson supercritical branching and then follow closely a deterministic curve, determined by the underlying nonlinear dynamics. This heuristics dates back at least to 50’s, e.g., Kendall (1956) and Whittle (1955), and the already mentioned (Metz 1978). A rigorous proof for epidemics is given in Mollison (1995), and in a wider context the rigorous implementations are relatively recent, see Barbour et al. (2015, 2016) and Chigansky et al. (2018).

Let us briefly sketch the ideas. The main ingredient is the Galton-Watson branching process , whose components mimic the behaviour of those of at the moment of mutant’s appearance, that is, around the equilibrium point . Thus its first, resident component is critical and starts at and its second, mutant component is supercritical with the offspring mean , and it starts from a single individual. The two processes and are constructed on the same probability space and are coupled in such a way, that they remain close at least until time with some fixed constant .

Following the above heuristics the density is approximated by gluing the linearized stochastic process with the deterministic nonlinear dynamics,where is the density of the Galton-Watson branching. The assertion (11) of Theorem 1 follows if we show that The main technical difficulty in proving (14) is to control the difference on the time interval , which itself grows with K. It turns out that the usual technique, based on straightforward linearization of the dynamics, does not provide bounds sharp enough in this case. Instead we construct a suitable coupling in Sect. 3.3, which involves several additional auxiliary Galton-Watson processes.

the process does indeed approximate the target density at time ,

and the approximation behaves asymptotically as claimed in Theorem 1,

The key to the limiting expression in (15) is the representationIt shows that the convergence in Theorem 1 follows once we prove the limit of Theorem 2 and check thatThe random variable W is the martingale limit of the supercritical branching , cf. (12). The most challenging element of the proof of this part is convergence (13), see Sect. 3.2 below. Previously, it has been proved in Chigansky et al. (2018) in dimension one, and analysis in higher dimensions, in our case two, requires a completely different approach. When convergence (11) is proved, the assertion of Corollary 1 follows by continuity of .

The limit in equation (10) can be proved in a similar way: note that in this case, cf. (16)where, in view of (17),We omit the proof of this part, which closely follows that of Theorem 1 with obvious adjustments.

The limit

In this section we construct limit (13), by means of a convergent telescopic series.

An auxiliary recursion in dimension one

Let us start with an auxiliary one dimensional quadratic recursionsubject to initial condition with , where and are constant coefficients. In what follows we will need the following estimate on its solution.

Lemma 1

There exists a finite function , such that

Proof

Let us first note that no generality will be lost if is assumed. Indeed, if (18) is multiplied by C, the recursionwith is obtained for the rescaled sequence . Hence if (19) holds for with some , then it holds for with . From here on we set .

Since proving the desired bound amounts to showingTo this end, consider the Schröder functional equationwhere and is the inverse of the parabola on . The function f(x) satisfies the following conditions Under these conditions (Seneta 1969) shows that the limitexists, solves (21) and satisfies the following properties Changing the variable in (21) to , we getand, inverting, we obtain the conjugacyHenceIn particular, (20) and, therefore also (19), hold.

f is continuous and strictly increasing on

and for

as

f(x) is concave (and therefore f(x)/x is decreasing on )

for all

on (nontrivial limit),

is monotone on

is invertible1

Properties of

Let us summarize some relevant properties of the function , which governs the deterministic dynamics in (3). Since and change signs across the lines and respectively, as shown at the phase portrait (Fig. 3), the coexistence equilibrium is globally stable. The local behaviour around the unstable fixed point is determined the Jacobian of at ,

Fig. 3

The phase-portrait of the deterministic system (3); the trajectory from a small vicinity of the resident equilibrium to the coexistence equilibrium is depicted in red (colour figure online)

To study perturbations around it will be convenient to consider the translationwith and Jacobian . In particular, existence of the limit (13) follows from that of .

Note that for and , formulas (4) for the entries of the function and the configuration of the fixed points (7) implyandHence the subset is forward invariant under , namely, . Then by (23) the subsetis forward invariant under .

In what follows, stands for the norm for vectors and the corresponding operator norm for matrices. In particular, the matrix A defined in (22) satisfies . The linear subspace is invariant under A andBelow C, , etc. stand for constants, which depend only on , and and whose values may change from line to line.

The first coordinate of can be written asand, similarly,Hence has the formwhere matrix satisfies the boundwith a constant C. Similar calculation also shows that for where matrix satisfiesThese formulas and the bound from Lemma 1 give the following growth estimate.

Lemma 2

For any and all n large enough,with a finite function , .

For any and all n large enough and, since E is invariant, for all m. Hence by (25), the sequence satisfiesBy induction , where solves (18) subject to , and the claim follows from Lemma 1.

Proof of Theorem 2

We will argue that the increments of are absolutely summable, uniformly over compacts in . Let n be large enough so that and therefore, by invariance, for all . Consider the arrayFor , due to (25),where and, in view of (26), is a sequence of vectors, whose norm is bounded uniformly in n. Both vectors and depend continuously on , which is omitted from the notations. For , (27) impliesand, letting , we getSince , by virtue of (29) and (28)where constants ’s depend continuously on .

Let us now bound the remaining term in (30). To this end, observe thatSince and is invariant under A, we have due to (24). Therefore for all where all ’s depend continuously on . ConsequentlyPlugging this and (31) into (30) yieldswhere depends continuously on . This implies that converges as , uniformly on compacts. Existence of the limit in (13) now follows from (23), and Theorem 2 is proved.

The main approximation

In this section we construct the random variable W and prove convergence (11). To this end, let U(n, j) and V(n, j) be i.i.d. random variables distributed uniformly over the unit interval [0, 1] and define and in (1) asDefine Galton–Watson branching process with componentssubject to and , and the corresponding density process . Note that coincides in distribution with the process defined in (12) andFinally, for a fixed constant defineIn particular, , cf. Theorem 1.

As explained in Sect. 3.1, the limit (11) follows from (14) and (15).

Proof of (15)

Since and we haveand hence for any ,This along with (33) implies (17), and in view of representation (16), the limit in (15) follows by the continuous mapping theorem and the uniform convergence in (13).

Proof of (14)

Sinceit suffices to prove thatandLet us first prove (36). Recall that the density process solves (6),and hence the difference satisfiessubject to . A direct calculation shows that the Jacobian of is boundedHence is -Lipschitz on with respect to norm andLet , thenwhere convergence holds if c is chosen close enough to 1. This proves (36).

To check (37), writeSince, by (33) and (35), the sequence converges to (0, W) in probability and, by Theorem 2, the sequence converges uniformly on compacts to , it suffices to show thatthat is,where has been already fixed in the previous calculations. We prove (38) for , omitting the similar proof for the case .

To this end, choose arbitrary constantsand, using the same random variables U(n, j) and V(n, j) as in (32), define two additional auxiliary Galton–Watson branching processes and with the entrieswhereandwhereDefine the exit timesThe random variable is a coupling time for the above processes, in the sense thatandHenceConvergence (38) for holds, if we show thatandsince for any .

The limit in (40) holds, sinceIt is left to prove (41). Sinceit suffices to check that each of the exit times, c.f. (39),satisfy for any. To this end, we can reuse the auxiliary processes and , defined above, with appropriately chosen parameters , and . Define exit timesTo prove (42c), we can choose . By construction,and thereforewhere holds by Doob’s inequality (Shiryaev 1996, Theorem VII.3.3, p. 493, eq. (11)), applied to the martingale .

To prove (42b), let us choose . By construction,and since the process is a submartingale,where is another variant of Doob’s inequality (Shiryaev 1996, Theorem VII.3.1, p. 492, eq. (1)) and holds by the Jensen inequality. The first term satisfieswhere the convergence holds by the choice . The second term satisfiesFinally, to prove (42a), let us choose , thenwhere holds since and we used Doob’s inequality as before. This verifies (41) and, in turn, (38) for . The proof for is done similarly and (37) follows. This completes the proof of (14).