Stochastic dynamics for adaptation and evolution of microorganisms
We present a model for the dynamics of a population of bacteria with a continuum of traits, who compete for resources and exchange horizontally (transfer) an otherwise vertically inherited trait with possible mutations. Competition influences individual demographics, affecting population size, which feeds back on the dynamics of transfer. We consider a stochastic individual-based pure jump process taking values in the space of point measures, and whose jump events describe the individual reproduction, transfer and death mechanisms. In a large population scale, the stochastic process is proved to converge to the solution of a nonlinear integro-differential equation. When there are only two different traits and no mutation, this equation reduces to a non-standard two-dimensional dynamical system. We show how crucial the forms of the transfer rates are for the long-term behavior of its solutions. We describe the dynamics of invasion and fixation when one of the two traits is initially rare, and compute the invasion probabilities. Then, we study the process under the assumption of rare mutations. We prove that the stochastic process at the mutation time scale converges to a jump process which describes the successive invasions of successful mutants. We show that the horizontal transfer can have a major impact on the distribution of the successive mutational fixations, leading to dramatically different behaviors, from expected evolution scenarios to evolutionary suicide. Simulations are given to illustrate these phenomena.
Keywords: horizontal gene transfer, bacterial conjugation, stochastic individual-based models, long time behavior, large population approximation, interactions, fixation probability, trait substitution sequence, adaptive dynamics, canonical equation.
MSC 2000 subject classification: 92D25, 92D15, 92D30, 60J80, 60K35, 60F99.
1 Introduction and biological context
A distinctive signature of living systems is Darwinian evolution, that is, a propensity to generate as well as self-select individual diversity. To capture this essential feature of life while describing the dynamics of populations, mathematical models must be rooted in the microscopic, stochastic description of discrete individuals characterized by one or several adaptive traits and interacting with each other. In this paper, we focus on the mathematical modeling of bacteria evolution, whose understanding is fundamental in biology, medicine and industry. The ability of a bacteria to survive and reproduce depends on its genes, and evolution mainly results from the following basic mechanisms: heredity, i.e. transmission of the ancestral traits to offspring (also called vertical transmission); mutation which occurs during vertical transmission and generates variability of the traits; selection which results from the interaction between individuals and their environment; exchange of genetic information between non-parental individuals during their lifetimes (also called horizontal gene transfer (HGT)). In many biological situations, competition between individuals and vertical and horizontal transfers are involved. The combined resulting effects may have a key role in the transmission of an epidemic, in the development of antibiotic resistances, in epigenetics or for the bacterial degradation of novel compounds such as human-created pesticides. There are several mechanisms for horizontal gene transfer: transformation, where some DNA filaments directly enter the cell from the surrounding environment; transduction, where DNA is carried and introduced into the cell by viruses (phages); and conjugation, when circular DNA (plasmids) replicates into cells and is transmitted from a cell to another one, independently of the chromosome. Conjugation plays a main role for infectious diseases since the genes responsible for virulence or antibiotic resistance are usually carried by plasmids. In this paper, we focus on conjugation modeling in order to understand transmission of pathogens and the evolution of antibiotic resistances.
We propose a general stochastic eco-evolutionary model of population dynamics with horizontal and vertical genetic transmissions. The stochastic process describes a finite population of discrete interacting individuals characterized by one or several adaptive phenotypic traits, in the vein of the models developed in . Other models for HGT have been proposed in the literature, based on the seminal contribution of Anderson and May on host-pathogen deterministic population dynamics  (see also [18, 24]) or on a population genetics framework without ecological concern (see [3, 22, 25]). Additionally, the previous models assume unilateral transfer, dividing the population into two classes: donors and recipients. In the present paper, we relax most of the previous limitations.
The stochastic model is a continuous time pure jump process with values in the space of point measures. The jump events are
births with or without mutation, horizontal transfers and deaths, with dynamics depending on the trait values of each individual and on the global population.
Our model covers both cases of frequency- and density-dependent horizontal transfer rates; these dependencies appear as special cases of a more general form of transfer rate, that we call Beddington-DeAngelis by analogy with a similar model used to describe predator-prey contacts ([4, 9]).
In a large population limit, using ideas developed in Fournier and Méléard , the stochastic process is shown to converge to the solution of a nonlinear integro-differential equation whose existence and uniqueness are proved. (See also Billiard et al.  for different evolutionary behaviors depending on the order of magnitude of population size, mutation probability and mutation step size). In the case where the trait support is composed of two values, the equation reduces to a non-standard two-dimensional dynamical system whose long time behavior is studied. This study highlights the impact of HGT on the maintenance of polymorphism and the invasion or elimination of pathogen strains. When a trait is initially rare in the population (e.g. a mutation of the common trait), its subpopulation is purely stochastic and we explain how HGT can drastically influence its probability of invasion and time to fixation. To do so, we combine the stochastic behavior of the mutant population size with the deterministic approximation of the resident population size. Then we assume that mutations are rare enough to imply a separation between the competition and mutation time scales, following ideas of Champagnat et al.  in a case without HGT. Here, under an Invasion-Implies-Fixation assumption, a pure jump process is derived from the population size process at the mutation time scale, for which the jump measure is strongly affected by the horizontal transfer. In the last section we present simulations in a case of unilateral transfer, which highlight the effect of HGT on evolution. In particular, we show that HGT can completely change the evolutionary outcomes. Depending on the transfer rate, we can obtain dramatically different behaviors, from expected evolution scenarios to evolutionary suicide.
2 A general stochastic individual-based model for vertical and horizontal trait transmission
2.1 The model
We model a bacteria population with a stochastic system of interacting individuals (cf. Fournier-Méléard , Champagnat-Ferrière-Méléard [7, 8]). Each individual is characterized by a quantitative parameter , called trait, which belongs to a compact subset of and summarizes the phenotype or genotype of the individual. The trait determines the demographic rates. It is inherited from parent to offspring (reproduction is asexual), except when a mutation occurs, in which case the trait of the offspring takes a new value. It can also be transmitted by horizontal transfer from an individual to another one. The demographic and ecological rates are scaled by which is taken as a measure of the "system size" (resource limitation, living area, carrying capacity, initial number of individuals). We will derive macroscopic models from the individual process by letting tend to infinity with the appropriate renormalization for individuals’ weight.
At each time , the population is described by the point measure
is the size of the population at time and the trait of the -th individual living at , individuals being ranked by lexicographical trait values.
Let us now describe the transitions of the measure-valued Markov process .
An individual with trait gives birth to a new individual with rate . With probability , the new individual carries the trait and with probability , there is a mutation on the trait. The trait of the new individual is chosen in the probability distribution .
An individual with trait dies with intrinsic death rate or from the competition with any other individual alive at the same time. If the competitor has the trait , the additional death rate is . Then in the population , the individual death rate due to competition is .
In addition, individuals can exchange genetic information. Horizontal transfers can occur in both directions: from individuals to or the reverse, possibly at different rates. In a population , an individual with trait chooses a partner with trait at rate . After transfer, the couple becomes . In the specific case of bacterial conjugation, the recipient acquires the trait of the donor (i.e. ). This occurs for instance when the donor transmits a copy of its plasmid to individuals devoid of plasmid (in that case, transfer is unilateral). We refer to the paper of Hinow et al.  for other examples.
We denote by the set of point measures on weighted by and by the set of finite measures on . The generator of the process is given for measurable bounded functions on and by
In particular, if we consider the function for a bounded and measurable function on and with the notation , we have
Assuming that for any , the functions , , and are bounded, it is standard to construct the measure valued process as the solution of a stochastic differential equation driven by Poisson point measures and to derive the following moment and martingale properties (see for example  or Bansaye-Méléard ).
Under the previous assumptions and assuming that for some ,
, we have the following properties.
For all measurable functions from into such that for some constant and for all , , the process
is a càdlàg -martingale starting from .
For a bounded and measurable function on ,
where is a càdlàg square integrable martingale starting from with quadratic variation
3 Large population limit and rare mutation in the ecological time-scale
3.1 A deterministic approximation
We derive some macroscopic approximation by letting the scaling parameter tend to infinity with the additional assumption of rare mutation, i.e.
The timescale is unchanged. It is called the ‘ecological’ timescale of births, interactions (competition and transfer), and deaths.
The next hypotheses describe the different scalings considered in the paper.
(i) When tends to infinity, the stochastic initial point measures converge in probability (and for the weak topology) to the deterministic measure and
(ii) When , the functions and (respectively ) converge uniformly on (respectively on ) to the continuous functions and (respectively to ).
(iii) We assume that for any ,
This means that in absence of competition, the subpopulation with trait is super-critical and that the regulation of the population size comes from the competition. We denote by
the intrinsic growth rate of the subpopulation of trait .
(iv) When , the functions converge uniformly on to a continuous function . This function depends on the mechanism of transfer. More precisely, we assume
where is continuous on .
The form (3.2) is derived from the so-called “Beddington-DeAngelis” functional response in the ecological literature ([4, 9]). This function covers different interesting cases regarding HGT. The horizontal transfer rate for an individual with trait in the population is . Assuming or very small gives a density-dependent HGT rate (denoted DD): the individual transfer rate is proportional to the density of the recipients in the population. Assuming or very large gives a frequency-dependent HGT rate (denoted FD): the individual transfer rate is proportional to the frequency of the recipients. Finally, assuming and gives a mixed HGT rate between frequency and density-dependent HGT rates (denoted BDA). This general case models some experimental observations for plasmids, for which a correlation between the form (density- versus frequency-dependent) of the transfer rate and the size of the population (low size versus close to carrying capacity) was suggested (Raul Fernandez-Lopez, pers. com.). We will show that the choice of has important consequences on the dynamics.
Assume and (3.1). Let . When , the sequence converges in probability in to the deterministic function defined for any continuous function as unique solution of
Let us note (by choosing ) that the total size of the population satisfies the equation
This equation is not closed and cannot be easily resolved except when is a constant function. Nevertheless by Assumptions , the functions and are bounded above and below by positive constants on . Then the process is bounded above and below by the solutions of two logistic equations which converge to strictly positive limits when . For example, where
with the notation and .
The proof is standard and consists in a tightness and uniqueness argument. The reader will follow the steps detailed in  or in : uniform moment estimates on finite time interval, tightness of the sequence of laws, continuity of the limiting values, identification of the limiting values as solutions of (3.3), uniqueness of the solution of (3.3). The last point deserves attention. Let us consider and two continuous solutions of (3.3) with the same initial condition . From the remark after (3.4), we have that
Let be a real bounded measurable function on such that . We obtain
By an elementary computation using Assumptions and (3.5), we obtain that for any ,
where is a positive constant. So taking the supremum over all functions such that and applying Gronwall’s Lemma we conclude that for all
Therefore uniqueness holds for (3.3). ∎
3.2 Trait replacement and the bacteria conjugation subcase
We now emphasize on the case where horizontal transmission results from the replacement of the recepient’s trait by the donor’s trait, i.e. and .
In this case, (3.3) becomes:
We note that the behavior of the deterministic dynamical system is influenced by HGT only through the ‘horizontal flux’ rate
In Section 5.1, we will show that in contrast, the fully stochastic population process depends not only on the flux but also on itself.
The horizontal flux rate quantifies the asymmetry between transfers in either directions and can be positive as well as negative (or zero in the case of perfectly symmetrical transfer). Note that bacteria conjugation is a subcase: a plasmid is transferred from the plasmid bearer to the empty individual , while the reverse is not possible (emptiness can not be transferred). This corresponds to the case where and and .
Assume that the initial measure is absolutely continuous with respect to the Lebesgue measure, then this property propagates in time and for any , , where is a weak solution of the integro-differential equation
Let us mention that at our knowledge, the long time behavior of a solution of this equation is unknown, except in the case without transfer where it has been studied by Desvillettes et al. . Some close equations with transfer have also been considered and studied in the long time by Hinow et al.  and by Magal-Raoul .
4 The two traits case
4.1 The dynamical system
Let us now assume that the population is dimorphic and composed of only two subpopulations characterized by the traits and . We set and define Let us assume that converges in probability to the deterministic vector . Then Proposition 3.2 is stated as follows.
When , the stochastic process converges in probability to the solution of the following system of ordinary differential equations (ODEs):
When , we get the classical competitive Lotka-Volterra system. Point is an unstable equilibrium and there are 3 stable equilibria: a co-existence equilibrium and two monomorphic equilibria and , where
is the unique stable equilibrium of the standard logistic equation
It is well known that the sign of the invasion fitness function, defined as
governs the stability. If and , the system converges to while if and , the system converges to and if and , the system converges to a non trivial co-existence equilibrium. In the case where the competition kernel is constant and is a monotonous function, the fitness function is equal to , which prevents co-existence in the limit.
When , the behavior of the system is drastically different as it can be seen in the phase diagrams of Figure 4.1.
Figure 4.1 shows eight possible phase diagrams for the dynamical system (4.1), where the circles and stars indicate stable and unstable fixed points, respectively. Figures (1)-(4) are possible for all forms of HGT rates, Figures (5)-(6) can happen in frequency-dependent or Bedington-deAngelis cases, while Figures (7)-(8) can be observed only in Bedington-deAngelis case. Compared to the classical two-species Lotka-Volterra system, at least 4 new phase diagrams are possible: Figures (5)-(8).
Now, define the invasion fitness of individuals with trait in the -resident population by
4.2 Properties of the dynamical system (4.1)
Let us now analyze the behavior of the system (4.1).
We first exclude the possibility of cycles contained in the positive quadrant. Recall that a Dulac function on is a smooth non vanishing function such that
has the same sign on the domain .
Assume that and . Then the function is a Dulac function in . As a consequence, the system (4.1) has no cycle in .
From this result and the Poincaré-Bendixson theorem ([12, Section 1.7] or [16, Th.1.8.1, p.44]) we conclude that any accumulation point of any trajectory starting inside the positive quadrant is either a fixed point or is on the boundary.
Expressing (4.1) in terms of the size of the population and proportion of trait , , we obtain:
These equations are generalizations of the classical equations of population genetics with two alleles under selection , in which we have made the influence of demography explicit. Eq. (4.5) is useful to investigate the dynamics on the boundary of the positive quadrant which is an invariant set.
Let us recall that
The points , and are the only stationary points on the boundary. The origin is unstable and the two other points are stable for the dynamics on the boundary. Their transverse stability/instability is given by the sign of the fitness function given in (4.1).
The proof is left to the reader. This implies that any accumulation point of any trajectory starting inside the positive quadrant is a fixed point. We now investigate the fixed points inside the positive quadrant.
Besides the fixed points in the boundary, there is
in the BDA case, , there are at most stationary points,
in the FD case ( ), there are at most stationary points,
in the DD case (), there is at most stationary point,
or a line of fixed points inside .
It is easier to consider the system in its form (4.5). The stationary points are denoted by for convenience. They satisfy
If and , we deduce from the first equation that
for . Replacing by this quantity, we write the second equation as
When and (BDA case), the term between the large brackets is a priori a polynomial in of degree . But explicit computation shows that the term of order vanishes. Then this polynomial is of degree and there are at most stationary points inside the domain. In FD cases, the expression simplifies as times a polynomial of degree and there are at most two stationary points. The DD case reduces to a Lotka-Volterra system. ∎
Let us first remark that the trace of the Jacobian matrix of any fixed point inside , is equal to
This implies that any fixed point inside the positive quadrant is either a sink (index ), a saddle (index ) or a non-hyperbolic point of index with a negative eigenvalue of the Jacobian matrix (because the vector field is analytic, see [12, Th.6.34]). We use the circuit with anticlockwise orientation drawn in Fig. 4.2. The largest radius is chosen large enough such that there are no fixed points outside the loop. The fixed points and on the boundaries are denoted by and on Fig. 4.2.
The arrows represent the directions of the vector field along the different parts of the circuit. It can be shown in all cases that for a radius large enough, the large arc contributes to the index.
Assume all fixed points are hyperbolic. The only possibilities are as follows:
- if and are unstable points, the index of the circuit is and there is either one stable point inside the domain or fixed points: stable nodes and one saddle point.
- if and are stable points, the index is -1 and there is either one saddle point inside or fixed points: saddle points and one stable point.
- if one of the points or is an unstable node and the other one a saddle point, then the index is and we have either fixed point or two fixed points: one saddle point and one stable point.
This proposition follows from the Poincaré-Hopf theorem: the index of the curve is equal to the sum of the indices of the fixed points inside the domain (see [12, Prop.6.26, p.175] or [16, Prop.1.8.4, p.51]). Combining this result with Proposition 4.3, one can decide between the different possibilities depending on the parameters.
The diagrams in Figure 4.1 realize the different situations described above. However, there may exist other diagrams in accordance with Proposition 4.5 that we have never observed numerically. We are yet unable to prove or disprove the existence of such other diagrams. One can nevertheless show that in the case where and are sufficiently similar, the phase diagrams of Figure 4.1 are the only possible ones (cf. ). In the case of non hyperbolic fixed points inside the positive quadrant (with index as mentioned previously), an analogue of Proposition 4.5 can be established. This situation is however exceptional since it implies a nonlinear (polynomial) relation between the coefficients.
4.3 The case of constant competition
Assume that the competition kernel is constant for all . Eq.(4.5) gives:
Let us consider separately the cases of FD transfer rate and DD or BDA transfer rates.
Frequency-dependent horizontal transfer rate. With and , (4.6) shows that there are only two equilibria for the second equation: or (Figures 4.1 (1)-(2)). Therefore there is no polymorphic fixed point and we get a very simple “Invasion-implies-Fixation” criterion: trait will invade a resident population of trait and get fixed if and only if
Thus, compared to a system without HGT, horizontal transfer can revert the direction of selection (i.e. and have opposite signs) provided that
This implies that HGT can drive a deleterious allele to fixation.
Density-dependent or BDA horizontal transfer rate. When , there exists a polymorphic fixed point when
If and are both positive, the above expression is negative and there is fixation of . If and are both negative, which never happens since the left hand side is positive and the right hand side is negative. So there is fixation of in this case. When and have opposite signs, there may exist a non-trivial fixed point which is stable if
In contrast to the classical Lotka-Volterra competition model in which constant competition prevents stable coexistence, HGT with DD or BDA transfer rates allows the maintenance of a deleterious trait () in a stable polymorphic state; this requires that the flux rate be positive and large enough in favor of to .
5 Rare mutation probability in the evolutionary time-scale
As seen in Section 3, it is not possible to capture the effect of rare mutations () at the ecological time scale. We have to consider a much longer time scale to observe this effect. The mutation time scale is of order and we will assume in the following that when is large enough,
A separation of time scales between competition phases and mutation arrivals results from this assumption. Indeed, mutations being rare enough, the selection will have time to eliminate deleterious traits or to fix advantageous traits before the arrival of a new mutant.
Let us now give a rigorous approach of the mechanism governing the successive invasions of successful mutants.
5.1 Probability and time of invasion and fixation under competition with horizontal transfer
As an intermediate step, we investigate the fate of a newly mutated individual with trait in a resident population in which trait is common. We assume that the invasion fitness of trait defined in (4.1) is positive, . This includes both cases of an advantageous trait (), or a deleterious trait () provided that the horizontal transfer rate from to is high enough. Figure 5.1 gives illustration of the determining effect of the transfer in the evolution. It shows the different stochastic dynamics one can obtain under frequency or density-dependent HGT in the simple case of unilateral transfer and that could not be realized without transfer. Figure 5.1 shows that a deleterious trait can invade a resident population and go to fixation (Fig.5.1(b) and (c)), or can stably coexist with the resident one (Fig.5.1(a) and (d)). Fig.5.1(a) especially shows that both traits stably coexist even though competition is constant, which is made possible by density-dependent HGT (we have recalled in Subsection 4.1 that it cannot occur for a usual Lotka-Volterra system).
Consider an individual with trait introduced in a resident population of individuals with trait , whose size is very close to the equilibrium . During the first phase, is very small with respect to . It can be approximated by a linear birth and death branching stochastic process, at least until it reaches the threshold , for a given . In this birth and death process, the transfer acts as a birth term and the transfer as a death term. When tends to infinity, the probability for the process to reach is approximatively the survival probability of the process (e.g. [6, 8]) and is given by
5.2 Times of invasion and fixation
We refer here to , where the results are rigorously proved. As the selectively advantageous trait increases from rare, the first phase of the -population growth has a duration of order . If reaches the threshold , then the second phase begins, where the processes stay close to the dynamical system (4.1). The deterministic trajectory, which has a duration of order 1, can reach one of two final states: either both types of individuals stably coexist, or individuals with trait invade the population and the -population density reaches the threshold (i.e. ). Should the latter happens, the third phase begins and can be approximated by a subcritical linear birth and death branching process, until is fixed and is lost. In this birth and death process, the transfer acts as a birth term and the transfer as a death term. The third phase has an expected duration given by (see [20, Section 5.5.3, p.190])
When , , which means that the third phase is of order in duration. Summing up, the fixation time of an initially rare trait going to fixation is of order
where the expressions for and are given in (4.1) and is a negligible term.
5.3 The Trait Substitution Sequence
The limiting population process at the mutation time scale will describe the evolutionary dynamics of invasions of successful mutants.
Let us assume in what follows that the ecological coefficients impede the coexistence of two traits. This is known as the Invasion Implies fixation (IIF) assumption.
Assumption (IIF): Given any and Lebesgue almost any : either is a stable steady state of (4.1), or we have that and are respectively unstable and stable steady states, and that any solution of (4.1) with initial state in converges to when .
From Section 4 we know that invasion does not necessarily imply fixation, even when the invasion fitnesses of the two types have opposite signs, as shown by Fig. 4.1 (5) and (6). In these cases, fixation depends on initial conditions and is usually not achieved when the invading type starts from a small density. Considering the special case of constant competition, however, invasion does imply fixation (cf. Section 4.3) if HGT rates are FD or when condition (4.8) is not satisfied if HGT rates are DD or BDA.
Assumption (5.1) together with Assumption (IIF) imply that for a monomorphic ancestral population, the dynamics at the time scale can be approximated by a jump process over singleton measures on whose mass at any time is at equilibrium. More precisely, we have
We work under Assumptions , (5.1) and (IIF). The initial conditions are with , and .
Then, the sequence of processes converges in law to the -valued process