Stochasticity in the adaptive dynamics of evolution: the bare bones.

First a population model with one single type of individuals is considered. Individuals reproduce asexually by splitting into two, with a population-size-dependent probability. Population extinction, growth and persistence are studied. Subsequently the results are extended to such a population with two competing morphs and are applied to a simple model, where morphs arise through mutation. The movement in the trait space of a monomorphic population and its possible branching into polymorphism are discussed. This is a first report. It purports to display the basic conceptual structure of a simple exact probabilistic formulation of adaptive dynamics.

1. Driving simplicity to its extreme

One aspect of the fascination with biology is the contrast between the overwhelming detail of nature and its order. Deus creavit, Linnœus disposuit, was Linnæus’ famous, not so humble summary of the two sides and who cared for which. The friendly and modest country gentleman Darwin had little of Linnæus's self-image, but shared his obsession with the nitty-gritty of flora and fauna, since early childhood. The idea of evolution must have come forth as gratifyingly simple to him and to any biologist immersed in the richness of nature. To a mathematician, it invites further reduction into a more or less axiomatic theory, where everything can be derived from first principles.

Adaptive dynamics is an attempt at this. In many regards, it is successful albeit (of course!) criticized for simplistic genetics as well as ecology. A critique more on grounds of principle is that it treats populations as infinite and even as a continuous matter. The renowned canonical equation [5] is an example of this.

Recently, Méléard, Champagnat, Lambert and co-authors have published a sequel of seminal articles [2-4], and others, reformulating classical adaptive dynamics [7,20,21] in terms of population-size-dependent birth-and-death processes, then passing to various limits by letting first the carrying capacity grow, so that populations become infinite, then mutational steps shrink to zero. The present paper presents another such oversimplified model, going back to the ideas of [1,26,27], and indeed in a certain sense the simplest possible, still retaining indispensable properties of a population under evolution. In our view, they are as follows.

Any such population is finite and consists of reproducing individuals.

Individuals are characterized by their traits.

Together with population size (of the own species as well as competitors) traits may influence reproduction.

But they do not determine it: there is a variation in the reproductive behaviour between individuals of the same kind and, seemingly, in the same circumstances.

At reproduction, there is a probability for mutation, resulting in a child with a trait different from its mother.

The resulting population survives or not, depending upon chance and fitness.

The finiteness requirement is essential: natural populations are finite. Still, passage to limits of infinite populations may help to uncover underlying structure, or exhibit its aspects more clearly. But this is a tool only and it should occur explicitly, so as not to hide difficulties, avoid mistakes and also render investigations into the accuracy of approximations feasible. There is no all-embracing continuity of properties from the finite to the infinite; infinite populations need not die out, even though their finite counterparts necessarily do, if bounded and allowing variation in reproduction between individuals. Usually, in models, furthermore, not only is the carrying capacity of the habitat large (tending to infinity) but also mutation steps small (infinitesimal), and two limits thus involved, the order between them possibly significant.

Variation is another intrinsic feature of nature. A deterministic description of individual life remains an illusion: there is randomness in life. The simplest stochastic reproduction conceivable assumed here is asexual binary splitting. Following the adaptive dynamics tradition, we disregard mating. This will lead to underestimation of extinction time, but simplifies theory radically. Thus, we assume that each individual either gets two children in the next generation or none.

Similarly, interaction with environment and population is drastically condensed. The habitat concerned has a carrying capacity K, and if there is only one population present, of size n, the probability of an individual successfully splitting is taken to be K/(K + n). Given the population size, individuals reproduce independently. If there are several morphs present, the number n is replaced by a linear combination of the various population numbers and coefficients expressing the strength of competition. Again, a Gaussian form of these coefficients will be assumed in the last section. But first we consider just one or two morphs, with given coefficients.

The carrying capacity can be different for different species; K can depend upon the trait considered. Then we write K(t), t being a real number, the trait, which can be thought of as abstract, like a reproductivity measure, or in a very concrete manner, like the length of a beak. In the last section, we shall work with a specific Gaussian form of dependence of carrying capacity upon trait.

At reproduction, mutation may occur with some small probability which we will assume to be neither trait nor population-dependent; with a probability μ, a newborn with a mother with trait t turns out to be a mutant, then with a random trait normally distributed around t.

Time structure is stripped to its bare bones: only the successive (non-overlapping!) generation size are counted, or alternatively everybody is supposed to have the same life-span of one season.

Clearly, when compared with the many-faceted phenomenon of evolution, such a mathematical sketch is nothing but a caricature. Still, though reproduction and time structure are the simplest possible, they retain a fundamental affinity to reality. Other properties are more arbitrary, like the form of dependence on the environment or the normal distribution of mutational steps or in the form of certain coefficients. The former has some historical precedence, in the papers cited and earlier, but nothing else to speak for it than its ability to yield a reproductive behaviour changing with the relation of population size to carrying capacity and prevalence of competing species. The assumption of normal distributions is also ad hoc though a central limit style argument from many small effects might have passed through the minds of the first proponents of this distribution for mutational steps. (In other contexts, it has been argued that mutational steps should follow extreme value distributions, since only the most advantageous mutations result in viable children [23].) But actually, we shall not use more of the normal distribution of mutational steps than the fact that deviations are of the order of the standard deviation, a fact that is much more generally true. And in many regards, this specific model is just a representative, and many of its properties, like those building on second moment properties, large deviation theory, quasi-stationarity and early growth and ultimate extinction, extend to more general structures.

In the next two sections, neither the mutation probability nor any mutation step distribution will occur, though in Section 3 the fact that we consider two populations, one of which called the resident starts from its carrying capacity, the other called the mutant population starting from one individual, certainly reflects having a mutation probability in mind, which is taken to be small enough for the first population to have grown to its carrying capacity before the mutation occurs. In Section 4, these mechanisms will be made explicit. The difference from the series of papers by Champagnat, Ferrière, Méléard and Lambert lies in our taking discrete time, non-overlapping generations, i.e. Galton–Watson-style processes as the simple ‘pure’ structure for first study.

2. A single trait

The model population is, thus, a binary, population-size-dependent Galton–Watson branching process. It starts from a positive integer number Z0 = z and is then recursively given by

The random variable ξ is the number of children of individual i in generation n. These are independent and identically distributed, given the population size Z, or indeed the whole past population history, Z0, Z1,…, Z. We refer to the process value Z as population size (at time or generation n). Since reproduction is identically distributed for all individuals in the same generation and the distribution, given Z, is the same for all generations n, we shall often delete the suffices, at least when not referring to several individuals in one context.

The trait dependence of the carrying capacity K is suppressed, since we consider only one trait, and no mutations yet. K should be thought of as a large number, and later we shall let it tend to infinity.

It is important to grasp this single trait process not only for its own sake but also for the analysis of the complete model with mutation and several traits.

Obviously, it behaves like a subcritical branching process whenever the population is greater than K, like a supercritical process for all population sizes smaller than K, and it is critical if the size is precisely K (then necessarily an even integer).

Since thus the population has a bounded expectation, it follows that the extinction probability is P(Zn → 0, as n ∞) = 1 (cf. [8, pp. 108 ff.] or [13]). Ultimately, it dies out.

Under the level K, the population is, however, supercritical. Therefore, it tends to increase, while this is the case, and it is prone to reach large values (around K) before ultimate extinction. In fact, we shall show that the time to extinction is very large, with an overwhelming probability. The population will settle at a quasi-stationary regime, in the sense that it will fluctuate around the value K for a time period that is exponentially long in K, until it finally dies out.

For a ≥ 0, write Ta for the first time the population reaches, or passes, a, from below or above, depending upon the starting position. For short, let T = T0 be the time of extinction. What is the relation between these two random variables for large a? Indeed, it is only rarely that the process dies out before reaching a level dK, if 0 < d < 1 is small or the starting number z large.

Let 0 < d < 1. Then for any 1 ≤ z ≤ dK,

In this, and elsewhere, the probability or expectation indexed by z, P, E means the probabilities of such a population starting from size Z0 = z. No index either indicates z = 1, or else some completely different probability measure.

Such assertions are proved by comparison with suitably chosen (not population size dependent) simple Galton–Watson branching processes, about which much is known. In the present case, consider such a binary splitting process with the probability of begetting zero children being d/(d + 1). Call it . Since x/(1 + x) is an increasing function of x, any k < dK yields

Hence, as long as population size stays below dK, the probability of producing no offspring is smaller than the corresponding probability pertaining to this classical Galton–Watson process. Therefore, clearly the probability that our process becomes extinct by time n, without crossing dK, is smaller than the corresponding probability for the binary Galton–Watson process . The latter must be smaller than the Galton–Watson probability of ultimate extinction, P( → 0) = . It is well known (cf. [8, p. 113]) that is the smallest root of the quadratic equation which is simply d. Hence,

Thus, with a positive probability, the population will not die out but reach sizes at the order of the carrying capacity. Sooner or later, it will then die out. But when?

Whatever the starting number z, carrying capacity K and time (generation) n, and [T] ≤ e.

Write

The elementary inequality yields

Hence, and the upper bounds on the probabilities follow by induction and another elementary inequality, 0 < 1 – u < e− for 0 < u < 1. The second assertion follows by summation:

A corresponding lower estimate can be deduced but it is irrelevantly small. Instead, note that we already know that with a positive probability the process will move from any starting point up to the vicinity of the carrying capacity. By large deviation theory, it will stay there a very long time until by a random perturbation it leaves the area (if it were deterministic, it would get stuck) and by chance becomes extinct. This is why the upper exponential bound derived indeed gives a survival time of the right order; persistence time is exponential in the carrying capacity. We shall have a closer look at this now.

Introduce the density process X = Z/K, and note that the offspring distributions are in fact functions of the density,

We have that where the ξ are the offspring numbers in the nth generation, explicitly indexed, and are indeed independent and distributed like twice a Bernoulli random variable with the parameter p(X) = 1/(1 + X). Writing m(X) = 2p(X) for the offspring mean when the density is X = Z/K, and centring the variables of the sum in the usual central limit fashion, we obtain

With and we can write this

Clearly, E[η] = 0 and

Hence, a central limit theorem argument shows that the noise term in Equation (2) is of order 1/√K.

It follows that the process X, which does of course depend upon K, converges (in various senses), as K → ∞, to the deterministic sequence {x}, defined by an initial value x0 and the recurrence x = f(x) [14]. This is a sort of law of large numbers; a corresponding central limit theorem is available in [16].

However, this convergence also pinpoints the difference between stochastic and deterministic systems. Indeed, if x0 = 0, then for all n, x = 0 since 0 is a fixed point of f. In its awkward way, this is how the deterministic system tells us that extinction is certain. For other outcomes, we need Z0/K → x0 > 0, as K → ∞. In other words, the population should become comparable to K in order to get attracted to the stable fixed point 1. However, by Proposition 1, we know that this will occur with a high probability, and after that X will be approximated, for large K, by x with a positive starting value. The process will be near 1 for late times n and this means that the original population size process Z will be near K. The deterministic dynamics being globally attractive, Z will aim at moving towards K; only the event of extinction can prevent this.

Now, due to the random nature of the density process, exhibiting small but persistent perturbations from x, ultimately there will be an exit from the levels [d K, d'K], 0 < d < 1 < d', due to large deviations of random terms.

The main theory for such systems is analogous to the Freidlin–Wentzell theory and was developed in [12], but the present special setting, the problem of exit from a domain, was solved in [17]. Under some assumptions, it was shown that if the deterministic system has a stable fixed point, then the perturbed system will stay in the domain of attraction of this point for an exponentially long time, e for some c > 0. Since in our case the domain of attraction is the whole positive line (but not zero), this result applies to the extinction time. For any δ > 0 and a starting point z ≥ d K,

In the special case of our model, where summands are (twice) binomial random variables, there is a direct proof, using a large deviation result due to Janson [11], stating c explicitly. Indeed, see also [9]. Since, in the present context, our interest in this process is rather as a representative of a broad class of population models, proof and results specific for that model will be published elsewhere [10].

If the process starts from a small number z of individuals, then the time it takes to reach the level dK (if it does not die out) is at least of an order approximated by the binary Galton–Watson process with mean 2/(1 + d) and at most of an order corresponding to such a process with mean reproduction 2K/(K + 1), i.e. of order log dK, since for any Galton–Watson process we will have

By exponential growth, the conclusion will follow that the time to reach d K is of the order log d K; for details cf. [10].

Close to K, our process is near critical and approach will be slower; it is well known that a non-extinct critical branching process tends to grow linearly [8], and our process being on the supercritical side will tend to reach the carrying capacity quicker than in linear time (1 – d)K. Actually, it can be proved that linear time is correct [10]. In conclusion, if the population does not die out before reaching high levels, the time until extinction is of the order e, only a minor part (of the order K) is spent growing up to a quasi-stationary size. That one also has a limiting distribution.

Consider the distribution of Xn, given that Xn > 0, for fixed K. As n → ∞, this converges to a proper distribution function called the quasi-stationary distribution. When in this setting K → ∞, all probability concentrates at the point 1.

This was shown in [15] and is also contained in the more general setting of [18,25]. Indeed, the existence of quasi-stationary distributions is a consequence of the Krein–Rutman theory of positive operators.

Simulations exhibit first the quick growth and ensuing quasi-stabilization of five populations starting from one ancestor (Figure 1) and then the pseudo-stable distribution (Figure 2).

Figure 1.

Five population developments with K = 50.

Figure 2.

Histogram of a population size for the last 500 of 10,000 generations with K = 50.

3. Two morphs

After some, probably a long, time, the population will experience its first mutation giving rise to a new population. By then, the original resident population will be in its quasi-stationary stage; mutation probabilities are thus taken to be small. The new, mutant population has just one member, its ancestor. The basis of adaptive dynamics can then be said to be furnished by the branching process fact that the new population either dies out or else embarks on exponential growth, in which case the old resident dies out, or the two populations will coexist for a time span that is exponential in the carrying capacity. We proceed to make this precise in terms of our model.

The process starts from a pair (Z(1)0, Z(2)0) of positive integers, the first component denoting the size of the resident and the second that of the mutant population, when the mutation appears.

Thus, Z(1)0 is in the vicinity of the carrying capacity and Z(2)0 = 1. Transitions from generation n to n + 1 are described by where ξ(, k = 1, 2, …, are independent given the preceding generation and satisfy and

Now, this process has two competing types, but it is not genuinely multi-type, in the sense that none of them can produce children with the other trait. We also disregard, presently, further mutation to new traits. The interaction coefficient γ is supposed to be positive but less than one, and the same in both directions. This means that interaction is taken as symmetric, an assumption that is not at all necessary. Indeed, taking them different can lead to several interesting phenomena, like that of evolutionary suicide [9,24]. It may be noted that though the specific form of dependence is certainly ad hoc, it has some historical precedence, going back at least to [19].

We assume that the resident population starts exactly at its carrying capacity, which is a1K. Then the mean reproduction of the mutant in a generation with z mutants present will be

It follows that the mutant is subcritical throughout, if a2 < γa1, and indeed dominated by a subcritical binary Galton–Watson process with the splitting probability a2/(a2 + γa1). By [22], the probability that the maximum M of the latter reaches j is of the order

Hence, for large carrying capacities, the mutant population will die out directly, never coming close to K with overwhelming probability.

If instead and K large enough, the mutant population will instead start supercritically, and either die out as above or else grow at a geometric rate for some time, much as in the preceding section. This is the case of possible invasion. It is illustrated in the following figure, which shows five runs of a population system with γ = 0.7, a1K = 40 and a2K = 70. In three cases the resident prevailed, the invader dying out very quickly, and in two runs, the invader took over (Figure 3).

Figure 3.

Five competitive population evolutions, K = 100, γ = 0.7, a1 = 0.4 and a2 = 0.7. In two of them (red and orange lines), the invader takes over.

We proceed to have a look at the qualitative behaviour of the process and the domains of super- and subcriticality (Figure 4). They will be given in terms of the (x1 = z1/K, x2 = z2/K)-plane in which the process has a supercritical, critical and subcritical reproduction, respectively. Solving the inequalities for the mean offspring number greater than one yields that these regions are

Figure 4.

Regions of super- and subcriticality.

defined by:

Note that the boundary is where one of the populations is critical; it consists of two segments of straight lines. The point at which both processes are critical is

Since 0 < γ < 1, there are three possibilities for this point

The remaining case x∗1 < 0, x∗2 < 0 can be ruled out, since it implies both that γ > a1/a2 and that γ > a2/a1, and hence γ > 1.

As an example, we give the following picture for x∗1 < 0, x∗2 < 0 made with

and where we take K as the scale unit of the graph.

The domain A above the two lines corresponds to subcriticality of both populations, the left intermediate domain B is where the mutant is subcritical and the resident supercritical and the right intermediate domain C corresponds to supercriticality of {Z(2)} and subcriticality of {Z(1)}, whereas the domain below the two lines corresponds to their being supercritical. This is precisely the situation where coexistence may be realized.

Like in the single trait setup, the process {(Z(1), Z(2)n)} scaled by K, the density process, can be approximated by a two-dimensional dynamical system for large values of K. The latter has a fixed point at (x∗1, x∗2) and is attracted to it. Combine this with the fact that the density process is non-negative to see that the co-existence case is when the fixed point lies in the positive quadrant x∗1 > 0, x∗2 > 0 and corresponds to the condition

Like in the single trait case, the two-trait process must die out. But the process being supercritical for small population sizes, both populations will settle to a quasi-stationary regime, in the sense that their sizes will fluctuate around the values K x∗1 and K x∗2 for an exponentially long time (in K), provided the new trait gets established, in the sense that its size becomes comparable to K.

Thus, the probability of this happening is needed. Estimates from below and from above the probability of more or less direct extinction are established in the following proposition, pointing at γa1/a2 as a rough indicator. Superscripts refer to the corresponding populations. We omit the suffixing starting point (which is (a1 k, 1)).

The probability of the mutant dying out while the resident remains above the level dK, 0 < d < a

Similarly, the probability of the mutant population dying out before it reaches the level b K, b < a

In this b, u can be so chosen that γu + b < a2.

Proof On the set {T(2) < T(1)} ∩ {T(2) < T(2)} ∩ {n < T(2)}, Z(1) < u K and Z(2) < b K. Hence, the probability of the new trait to reproduce successfully is

Provided γu + b < a2, the splitting probability is larger than one-half. Selecting u in the range a1 < u < a2/γ allows us to choose 0 < b < a2 – γu < 1. Consequently the trait extinction probability does not exceed that of a supercritical Galton–Watson process with the splitting probability as above. This leads to the upper bound. Bounding the splitting probability from above on the set {T(2) < T(1)} yields

and the lower bound for the extinction probability of the mutant.

We proceed to bounding the probability of the mutant dying out, while the resident remains close to the carrying capacity.

Let Z satisfies

We note that for large K, we will have that say.

For 1 ≤ k ≤ n, let

For i ≥ 1 and j ≥ d K, since the expression in the middle decreases with i. Hence,

Thus, for any k ≤ n which, evidently, implies proving Equation (8).

A similar argument can be used to find an exponential bound for mutant survival probabilities, given that the resident does not surpass high levels [10].

Like the single trait process, considered in Section 2, the present process behaves as a deterministic dynamical system perturbed by small noise, which is, however, not to be forgotten.

Write the density process X(1) = Z(1)/K, X(2) = Z(2)/K and note that the offspring distributions, as given in the beginning of this section, are again functions of the density. As before,

Denote the offspring mean by m(x) when the density is x, m = (m1, m2), and write

The deterministic dynamics is given by and

It is also easy to see that the noise term when the system is at x is of order 1/√K. Indeed, has mean zero and variance one. Hence, we can write Equations (9) and (10) as

The function f in Equation (11) that generates this system has two fixed points 0 and

This is exactly the point at which both reproductions are critical. The determinant Δ of partial derivatives shows that 0 is repelling, since Δ = 4 there. At the point x∗, the determinant of partial derivatives is which allows deriving conditions for stability Δ < 1 of that fixed point in terms of the parameters γ, a1, a2.

Indeed, the deterministic system has four equilibrium points (0, 0), (a1, 0), (0, a2) and (x∗1, x∗2). As pointed out, the origin is always unstable. The fixed point corresponding to (x∗1, x∗2) lies in the positive quadrant only if γ is less than the maximum of a2/a1 and a1/a2. In that case, it is also attracting, and the other two points are unstable. If γ > a1, the point (a1, 0) is stable, so the second morph goes extinct. If γ > a1/a2, (0, a2) is stable, and the former resident goes extinct.

4. A simple stochastic model with mutation and adaptive dynamics

Having investigated one lone species and the joint behaviour of a resident and a mutant, we return to the well-established special [6] model sketched in Section 1, determined by Gaussian functions. It is controlled by five parameters: K, μ, σ1, σ2 and σ3, with the following interpretation.

K is the maximal carrying capacity of a monomorphic population (large).

μ is the mutation rate per gene per generation (and very small).

σ21 describes how the carrying capacity varies.

σ22 is often referred to as the competition width.

σ2 is the variance of mutational steps (again small).

Important compound parameters are

Note that −1 < δ < 1, and

A morph with trait t ∊ ℝ has the carrying capacity explaining the role of σ21, in what sense it describes how carrying capacity varies for different traits, with the maximum K(0) = K for the trait value t = 0.

In a polymorphic population, the competition between two morphs t1 and t2 is quantified by the coefficient depending only on the difference |t1 – t2| and reaching its maximum value γ(t, t) = 1 for individuals with the same trait value.

Mutation in a monomorphic population occurs at a rate of the order μ mutations per individual and time unit. From the preceding analysis, the mutation rate should be taken to satisfy e− ≪ μ ≪ 1/K, the constant c defined in Equation (3), so as to avoid repeated mutations during the time while one mutant establishes itself, while guaranteeing mutations during the residence. In adaptive dynamics, this would be formulated as conditions on mutations not occurring in ecological time spans, but (almost) certainly in an evolutionary time scale. Given a mutation, the new trait value t2 then has the normal conditional density function implying symmetry in the step distribution and that larger changes in trait value are less probable.

The process starts from a finite number of individuals of a finite number of different traits. Traits in the various generations are denoted by vectors (t1, t2, …, t), where p can vary from generation to generation. Denote the number of individuals with trait t in the nth generation by Z(t). Given the vector (Z(t1), Z(t2), …, Z(t)), individuals of the nth generation reproduce independently, splitting into two with the probability otherwise begetting no child. Each child is a mutant with probability μ, the mutant's trait being normally distributed around its mother's, and with the variance σ2.

This is an informal (but precise!) description of a discrete-time stochastic process, whose states are counting measures on ℝ. If it starts from one individual at say t1, it will develop according to Section 2 (with K(t1) replacing K), until the first mutation occurs. Then we find ourselves in the circumstances of Section 3 with γ as defined by Equation (15), a1 = exp(–t21/(2σ21)), and a2 = exp(–t22/(2σ21)), K = K(0), according to Equation (14).

After some algebra from the inequality a2 > γa1 (Equation (4)), we arrive at a simple condition for the mutants possibly being able to establish themselves: If t2 < t1, then also t2 < δt1. If t1 > 0, the second condition is empty, whereas if t1 < 0, t2 must be at a considerable distance left of t1, which is not feasible due to the smallness of σ. We conclude that residents with traits on the right half-axis can only establish mutants to the left of them, and vice versa. Indeed, the situation is summarized by the condition being necessary for establishment, whatever be the position of t1, and equivalent to t2/t1 ∊ (−∞, 1) ∪ (1/δ, ∞).

This means that in the monomorphic phase, the trait value will move towards zero until it reaches a σ-neighbourhood of zero. Then the trait value may jump in both directions and moreover the first evolutionary branching becomes plausible. Indeed, the reverse invasion condition is equivalent to t2/t1 ∊ (−∞, δ) ∪ (1, ∞). When both conditions are fulfiled, the pair of trait values (t1, t2) may coexist.

In a slightly more precise wording, we start the process from a monomorphic population, say with trait T0 = t > 0. The population will either die out or (quasi-)establish itself at its carrying capacity. In the latter case, sooner or later at a random time v1, a successful mutant will replace it, and the trait will have moved to TV1, where it will stay until the next successful mutation v2, until finally it would have come so close to zero, so as to render branching possible, the trait being of the order of σ. Then {T} will take the value of two traits. Thus, we obtain a stochastic process with values in ∪n ℝ, about which many interesting questions can be posed.

What is the time to the first successful mutation? In other words, what is the distribution of ν1?

When will the first branching occur, or rather, what is the distribution of the waiting time until a monomorphic population turns dimorphic, the distribution to the next change in number of morphs, etc.

How can the stochastic movement {Tn} in trait space be described?

And how the degree of polymorphism (the development of the number of coexisting species)?

What limiting behaviours can be discerned, when K → ∞ and μ → 0?

For birth-and-death processes, many of these matters have, as pointed out, been discussed by Méléard, Champagnat et al. (see in particular [3]). Our more mathematical findings on the present simple model will appear in [10]. As one of the referees pointed out, the basic method used here, domination by not-size-dependent branching processes during time periods when populations sizes satisfy certain inequalities, invites usage also in more realistic, age-dependent or general branching processes.

Acknowledgements

The paper has profited from a knowledgeable and dedicated refereeing. We thank the three referees. F.C. Klebaner was supported by the Australian Research Council Grant DP0881011. S. Sagitov was supported by the Bank of Sweden Tercentenary Foundation. V.A. Vatutin was supported by the RAS programme ‘Mathematical Control Theory’ and INTAS 03-51-5018. P. Haccou was supported by the Netherlands Organisation for Scientific Research (NWO). P. Jagers was supported by the Swedish Research Council.