THE INFLUENCE OF INFECTIOUS DISEASES ON POPULATION GENETICS

Malaria is the vector-transmitted disease that causes the highest morbidity and mortality in humans. Motivated by the known influence of sickle-cell anemia on the morbidity and mortality of malaria-infected humans, we study the effect of malaria on the genetic composition of a host (human) population where sickle-cell anemia is prevalent and malaria is endemic. The host subpopulations are therefore classified according to three genotypes, AA, AS, and SS. It is known that AA malaria-infected individuals experience higher malaria-induced mortality than AS or SS individuals. However, individuals carrying the S gene are known to experience a higher mortality rate in a malaria-free environment than those who lack such a gene. The tradeoffs between increased fitness for some types in the presence of disease (a population level process) and reduced fitness in a disease-free environment are explored in this manuscript. We start from the published results of an earlier model and proceed to remove some model restrictions in order to better understand the impact on the natural hosts’ genetics in an environment where malaria is endemic.

1. Introduction.Efforts to understand a host's evolutionary dynamics (often slow) in the context of disease dynamics in "chemically" treated environments (e.g., chemotherapy) or the influence of diseases on the host's genetic variability or both have been carried out recently (see, for example, Andreasen [1], Beck [2], Castillo-Chavez and Feng [3], Feng et al. [5,6], Galvani et al. [8,9], Hsu-Schmitz [13], Kribs et al., May and Anderson [10], McKenzie [11,12]) in the context of human diseases such as HIV, tuberculosis, malaria and others.Feng et al. [5] studied the influence of malaria dynamics (assumed to be fast) on hosts' fitness in their efforts to assess the effect of vector-borne diseases on the genetic composition of a host population.A typical setting was adopted by assuming that the human population is characterized by a two-allele single loci system.Here, we remove some of their model simplifications and show that the qualitative dynamics do not change.However, we briefly illustrate the potentially significant effect that the incorporation of additional host characteristics (age) may have on disease dynamics and on 468 ZHILAN FENG, CARLOS CASTILLO-CHAVEZ the genetic structure of the host population.In other words, incorporating the lifehistory dynamics of the host can have a dramatic effect on the qualitative dynamics of the disease and the gene distribution in the host population.Here, we have only scratched the surface.
It was assumed (see [5]) that the human population could be properly characterized by a two-allele single loci system in the study of a malaria-human system.It is assumed that only three genotypes are possible: AA, AS, and SS, where the letters A and S denote the two alleles.Under a simplifying assumption that the SS individuals do not survive, the following model was introduced and analyzed in [5]: ( Here, u i (t) denotes the host population of uninfected i-individuals (i = 1, 2 with AA = 1, AS = 2), v i (t) denotes the host population of infected i-individuals (i = 1, 2), and N (t) = 2 i=1 (u i (t) + v i (t)) denotes the total host population at time t.Pi (t) (i = 1, 2) denotes the fraction of each genotype born into the host population and it is assumed that these proportions satisfy the following relationships (Hardy-Weinberger proportions): where q(t) = r(t)/2 is the frequency of the S-gene, p(t) = 1 − q(t) is the frequency of the A-gene, and Following the approach of the classical Ross McDonald model for the spread of malaria, we use the fraction of infected mosquitoes, z(t), instead of the total number of infected mosquitoes in our model.The rationale behind this assumption comes from the fact that what typically determines "observed" vector densities is the number of vector breeding sites where there is fierce competition and an oversupply of eggs.Susceptible human hosts of type i become infected at the rate β hi zu i (i = 1, 2), while the fraction of susceptible mosquitoes becomes infected (from biting an infected human host of type i) at the rate (1 − z)β vi v i /N .Here, β hi and β vi (i = 1, 2) denote the transmission coefficients from humans to mosquitoes and mosquitoes to humans, respectively.The malaria disease-induced death rate is α i while the average time before a victim of malaria recovers is denoted by 1/γ i .The documented resistance of AS individuals to malaria infection is modeled by a reduction in susceptibility to disease invasion.Hence, it is assumed that The hosts' demography is important, so b(N ) denotes the human per capita birth rate, possibly density dependent, and m i = m+ν i , where m is the assumed constant per capita natural mortality of humans and ν i (i = 1, 2) the excess per capita death rate ascribed to S-gene carriers.Here, we take ν 1 = 0, ν 2 > 0. Definitions of all variables and parameters can be found in Table 1.
Was the omission of SS individuals justified in model (1)?Here, we remove such a restriction, making our work relevant to vector-transmitted infectious diseases for which all three genotypes of individuals need to be considered.Although there are obvious quantitative differences, we show that the qualitative results are equivalent to those generated by the simpler model.In other words, on the fast time scale of malarial dynamics, the disease level reaches an equilibrium and, consequently, the host's gene frequency distribution is influenced only by the equilibrium level of malaria.However, the generality of this conclusion comes into question when the host's population structure is incorporated.The rest of the manuscript is organized as follows: Section 2 describes the extended model; Sections 3 and 4 justify the simplifications that result from the identification of two processes whose dynamics take place at highly distinct temporal scales; Section 5 looks at the fitness of the S-gene, and Section 6 provides a critical discussion of the results while suggesting that a key missing component, the host's population structure, may indeed make the robustness of the above results questionable.
(2) will be changed to where and In this case the system (1) can be extended to become the following seven-dimensional system: ui All parameters have the same meaning as in (1), with some additional assumptions concerning the SS individuals; i.e., ν 2 ≤ ν 3 , and As in [5], we take advantage of the fact that the time scales of evolutionary processes compared to those of the disease are very different, which allows us to use a singular perturbation approach to separate the fast and slow dynamics of the system.For this purpose we introduce new variables, x i = u i /N and y i = v i /N .These are the fractions of corresponding subpopulations of human hosts.Using the new notation, we can rewrite the frequencies of the AS and SS individuals as r = x 2 + y 2 and s = x 3 + y 3 , respectively.We note for clarification that x 1 + y 1 + x 2 + y 2 + x 3 + y 3 = 1 and x 1 + y 1 = 1 − r − s.We can replace the set of variables in the system (4) with the new set of variables: y 1 , y 2 , y 3 , z, r, s, and N , which obviously describe important epidemiological, demographic, and population genetic quantities.System (4) is equivalent to the following system: Following the approach in [5] our mathematical analysis is for the specific density dependent per capita birth function, b(N ) = b(1 − N/K), where b is a constant (the maximum birth rate when the host population size is small) and K is approximately the density dependent reduction in birth rate.We remark that the full model in [5] is a five-dimensional system, and the slow system consists of two equations.Here, the full system ( 4) is seven-dimensional, and we will show next that the slow system consists of three equations.

3.
Fast dynamics of epidemics.The relevant parameters vary across many orders of magnitude.For example, the demographic parameters (b and m i ) and the genetic parameters (α i ) are on the order of 1/decades, and the malaria disease parameters (β hi , γ i , β vi , and δ) are on the order of 1/days.Hence, although the malaria disease dynamics and the changes in genetic composition are two coupled processes, the former occurs on a much faster time scale than the latter.Let m i = mi , α i = αi , and b = b with > 0 being small.We can use this fact to simplify the mathematical analysis of the full model with the use of singular perturbation techniques, which allows us to separate the time scales of the different processes.By letting = 0 we obtain the following system for the fast dynamics: where r and s are considered parameters which represent genotype frequencies of AS and SS individuals.Thus, the fast system (7) describes the epidemics of malaria for a given distribution of genotypes determined by r and s.Let w = r + s.Then, q ≤ w ≤ 2q.Note that q(t) → 0 if and only if w(t) → 0 as t → ∞.Therefore, the variable w provides a good measure of the S-gene frequency when studying persistence properties.
The basic reproductive number for the malaria epidemics can be calculated using the next generation matrix (see [4]).We keep in mind that, on the fast time scale, r and s are treated as parameters.Let E 0 = (0, 0, 0, 0) be the disease-free equilibrium (DF E) of the fast system (7).Then from the Jacobian matrix at E 0 we obtain the next generation matrix: whose leading eigenvalue is This quantity gives the basic reproductive number.However, to simplify the notation we define our R 0 to be this number squared: which involves parameters associated with malaria transmission between mosquitoes and humans of genotype i.In fact, R i (or √ R i ) is the basic reproductive number when the host population consists of entirely humans of genotype i.When R 0 < 1, all eigenvalues of G are bounded by 1, and at least one eigenvalue of G will exceed 1 when R 0 > 1.This implies that the Jacobian at E 0 has all eigenvalues with negative real parts if R 0 < 1 and has at least one eigenvalue with positive real part if if R 0 > 1.We have proved the following result.
Theorem 1.The DF E E 0 of the fast system (7 The following result also holds as in most endemic models: Theorem 2. The fast system (7) has no endemic equilibrium when R 0 < 1.When R 0 > 1, a unique endemic equilibrium exists and is l.a.s.
Proof.To simplify notations we introduce the variable w = r + s, which represents the frequency of individuals carrying the S-gene.Using (5) we have Then we can rewrite R 0 as ) be a nontrivial equilibrium of the fast system (7); i.e., all components of E * are positive.Then where and z * is a positive solution of a quadratic equation with Here we have used the fact that T h3 = T h2 which is due to (5).
To show that the equation ( 9) has no positive solution, it suffices to show that k 1 > 0. Let T 0 = max{T h1 , T h2 }.Then T h1 + T h2 − T 0 R 0 ≥ 0, and Since k 0 is always positive, it follows that both solutions of (9) are negative.Therefore, E * is not biologically feasible.
Notice that M ≥ 0; i.e., all elements of M are nonnegative (recall that 1−r−s−y and D is a diagonal matrix with positive diagonal elements.It is known that all eigenvalues of H have negative real parts if and only if the dominant eigenvalue of the matrix M D −1 is less than one.M D −1 has a double zero eigenvalues and two other eigenvalues given by where Obviously, λ − < 1.Next we show that λ + < 1.Using the following equalities (which are obtained by setting the right-hand side of the y 1 and z equations in (7) equal to zero), , or and noticing that β hi z * + γ i > γ i , i = 1, 2, 3, we obtain Similarly we can obtain the following inequalities: Substitution of these inequalities into the expression of λ + yields Since 0 < z * < 1, it follows that λ + < 1 and that E * is l.a.s.This finishes the proof.
4. Slow dynamics of population genetics.By using the rescaled time τ = t, we can rewrite the full system (6) as where This system has a three-dimensional slow manifold: M = (y 1 , y 2 , y 3 , z, r, s, N ) : which is normally hyperbolically stable, as it consists of a set of such equilibria E * of the fast system (7).The functions y * i , i = 1, 2, 3, and z * are given in ( 8) and ( 9).The slow dynamics on M are described by the equations .Phase portraits of system (12) in the (q, N ) plane (q = r/2 + s).For (a) and (b) all parameters are equal except ν 2 and ν 3 .In (a) ν 2 = ν 3 = 5 × 10 −5 for which F > 0. It shows that there is a stable interior equilibria (represented by a "•").In (b) ν 2 = ν 3 = 6 × 10 −5 for which F < 0. It shows that there is no interior equilibrium and that the only stable equilibrium is on the N -axis.Other parameter values are given in the text. where and again y * 1 , y * 2 , y * 3 and z * are given in ( 8) and ( 9).Since M is normally hyperbolically stable, geometric theory of singular perturbations due to Fenichel [7] allows us to study the system (11) by studying the reduced slow system (12) (see [6] for more details).In other words, if the dynamics of system (12) can be characterized with bifurcations, then the bifurcating dynamics on the slow manifold M are structurally stable and hence robust, subject to perturbations.Therefore, results from the slow system will provide bifurcation properties of the system (11) as well as the full system (6).Some interesting scenarios are illustrated in Figure 1, which shows numerical solutions of the slow system for two sets of parameter values.The frequency of the S-gene is q = r/2 + s.The two sets of parameters used for the figures are chosen to give positive (Figure 1(a)) and negative (Figure 1(b)) values, respectively, of the S-gene fitness F (whose definition is given in the next section).Figure 1(b) shows that the boundary equilibrium on the N -axis of the slow system is a global attractor for negative fitness.When the fitness is positive, Figure 1(a) shows that it is possible for the S-gene to be maintained provided that the initial density of the S-gene is not too low.One of the important features of this coupled system is the possibility of multiple interior equilibria (see Figure 1(a)) and the bistability of two nontrivial equilibria (represented by the solid circle).The parameter values used are β h1 = 0.04, , and K = 10000.These outcomes are usually absent in population genetic models without a dynamic disease process.

Fitness of the S-gene.
Recall that q = r/2 + s is the S-gene frequency.The invasion ability of the S-gene can be described by 1 q dq dτ q=0 .
We will use this quantity to define the fitness of the S-gene, F. Consider the case when ν 3 = ν 2 > 0; i.e., the extra death rate due to the S-gene in AS and SS individuals is equal.Then m 3 = m 2 , and using the r and s equations in (12) we get Since 0 ≤ r, s ≤ q, we have r = s = 0 if q = 0. Using ( 8) and ( 9), and noticing that α3 = α2 and T h3 = T h2 , we have 1 q Here, we have used the fact that : We can also show that and 1 q It follows from ( 13)-( 16) that the fitness of the S-gene is given by where Clearly, F is determined by the difference of weighted death rates between the nonsickled (AA) and sickled (AS and SS) individuals, and the weights, W 1 and W 2 , depend only on epidemiological parameters.This allows us to explore the effect of malaria epidemics on the distribution of genotypes.
The following result confirms that the fitness F defined above indeed determines the invasion ability of the S-gene.Theorem 3. The slow system (12) has a trivial equilibrium U 0 = (0, 0, N 0 ), where Proof.We first notice that (see ( 14), (15), and (18)) Then the Jacobian matrix at U 0 can be written where a " * " denotes a number that does not affect the eigenvalues of J. Since α3 = α2 and m3 = m2 we know that the three eigenvalues of J are Therefore, J has three negative eigenvalues if F < 0 and has a positive eigenvalue if F > 0. This completes the proof.
Theorem 3 implies that the S-gene cannot invade if the fitness is negative and that the invasion is possible if the fitness is positive.This result is confirmed by numerical computations of the slow system (see Figure 1).Figure 1(a) is for the case F > 0, which shows that there is a locally asymptotically stable interior equilibrium, and Figure 1(b) is for the case F < 0, which shows that there is no interior equilibrium and that the boundary equilibrium on the N -axis is a global attractor.These numerical results also suggest that the fitness measure given by F provides a criterion not only for invasion but also for the possible maintenance of the S-gene.Figure 1(a) shows that the S-gene can establish itself if its initial frequency is not too low.Notice from (17) and (18) that the fitness F is influenced by several of the parameters associated with malaria epidemiology.Figure 2(a) illustrates the joint effect of β h1 (malaria infection rate of type 1 humans) and γ 1 (recovery rate from malaria of type 1 humans) on the fitness.A contour plot of Figure 2(a) is shown in Figure 2(b), which can also be viewed as a bifurcation diagram in the parameter plane (γ 1 , β h1 ) indicating the region in which U 0 is stable (F < 0) or unstable (F > 0).This bifurcation is confirmed by our numerical simulations of the system (see Figures 1 and 3).
The slow system (12) can also have multiple interior equilibria and bistable nontrivial equilibria when the fitness is positive.Figure 3(a) is a phase portrait showing that there are two interior equilibria with one of them locally asymptotically stable (represented by a solid circle) while the other one unstable (a saddle, represented by a solid triangle).There is another stable equilibrium on the the N -axis.Figure 3(b) is a time plot showing the two locally stable equilibria.Parameter values for this figure are the same as those used for Figure 1 except that β v2 = β v3 = 0.12.6. Discussion.We extended the model in [5] by including all three genotypes of individuals.More complexity is introduced by adding the additional subpopulation of SS individuals, which increased the dimension of the full system from 5 to 7 and the dimension of the slow system from 2 to 3.However, the qualitative behaviors of the two models seem to be the same.By coupling the malaria epidemiology and the sickle-cell genetics we investigate how the epidemiological and demographic parameters may affect the fitness of the sickle-cell gene.We derived threshold conditions which allow us to explore the impact of malaria disease dynamics on the genetic composition of the human population.From these threshold conditions, we can draw conclusions similar to those obtained in [5].For example, whether the rare gene will go to extinction or persist in a population (on the slow time scale) is determined by the fitness coefficient of the gene, and this fitness coefficient depends not only on parameters related to the genetics of the S-gene but also on the parameters associated with malaria epidemiology which determine the disease prevalence.Our analytic results consider the ability of a rare gene to invade a host population composed of mainly wild-type individuals.This invasion ability provides a measure for the fitness of the S-gene F. Our numerical studies indicate that when F > 0 the S-gene not only can invade but also will be maintained in the host population.Moreover, the slow system (host population genetics) can have multiple interior equilibria and multiple stable nontrivial equilibria.This occurs when F is negative but close to zero.In this case, although the S-gene cannot invade when its initial frequency is low, it may be able to establish itself if its value becomes large due to some stochastic perturbations (e.g., environmental changes).
The main findings from this model are 1) the S-gene may be selected for if the selection force from malaria is strong; and 2) multiple stable equilibria are possible.These types of population dynamics are not observed in most models of population genetics models without coupling malaria epidemiology.One would expect this malaria-sickle cell system to exhibit oscillatory behaviors under certain conditions for coexistence.Intuitively, the frequency of the sickle-cell gene will decrease in the absence of malaria due to a higher S-gene-related death rate in the SS and AS individuals than in the homozygote wild-type individuals (AA).On the other hand, the S-gene may be selected for if the endemic level of malaria is sufficiently high, and consequently polymorphism in the host population may be maintained in an oscillatory manner.However, as is seen in [5], the ODE system (12) or (4) does not produce periodic solutions.To observe such oscillatory dynamics we looked at a further extension of the model to include a time delay in the density dependence of the birth function.For example, if we choose the birth function in system (12) to be of the form b[1 − N (t − T )/K], where T denotes a constant time delay, then the modified system is capable of producing oscillations (see Figure 4).From ( 13) and ( 16) we see that the introduction of delay will not affect the formula for the fitness F given in (17).Notice that m 2 − m 1 = ν, where ν is the extra death rate of individuals with the S-gene.Rewrite this formula (17) as F = W 1 α1 − W 2 α2 − ν.
The first two terms involve only malaria-associated parameters with the property W 1 α1 − W 2 α2 > 0. Therefore, the sign of F is determined by the genetic parameter ν if all the malaria related parameters are fixed.
Our numerical simulations of the slow system (12) show that there exists a critical value ν c such that solutions converge to a positive equilibrium for ν > ν c (top panel of Figure 4) and that oscillations occur for ν < ν c (middle panel of Figure 4) for a given value of the delay T = 0.25.We can also fix the value of ν and vary T to observe the switching from a stable equilibrium to periodic oscillations.For example, the top and bottom figures in Figure 4 5 shows that the full system (4) also exhibits similar dynamics as that of the slow system (12).From the oscillatory solutions (see Figure 6) we observe that for most of the time the change in the gene frequency q (= r/2 + s) follows the change in the endemic level of malaria (i.e., the fraction of total infected people y =: y 1 + y 2 + y 3 ).For example, q increases with increasing y until y reaches its maximum.Soon after y starts decreasing, q also starts decreasing, until y reaches its minimum.
Other modifications of the model may also lead to periodic solutions.We have only presented a simple example.In any case, by explicitly coupling the malaria disease dynamics with changes in the frequency of the S-gene, our model allows

Figure 1
Figure1.Phase portraits of system(12) in the (q, N ) plane (q = r/2 + s).For (a) and (b) all parameters are equal except ν 2 and ν 3 .In (a) ν 2 = ν 3 = 5 × 10 −5 for which F > 0. It shows that there is a stable interior equilibria (represented by a "•").In (b) ν 2 = ν 3 = 6 × 10 −5 for which F < 0. It shows that there is no interior equilibrium and that the only stable equilibrium is on the N -axis.Other parameter values are given in the text.

Table 1 .
Definition of variables and parameters i malaria-induced death rate of humans of type i α 1 ≥ α 2 = α 3 νi S-gene related death rate of humans of type i ν1 = 0 < ν2 ≤ ν3 b(N ) per capita birth rate of humans P1 fraction of total births of genotype AA (1 − r/2 + s) 2 P 2 fraction of total births of genotype AS (1 − r/2 + s)(r + 2s) P3 fraction of total births of genotype SS (r/2 + s) 2