## Abstract

In this paper we use a system of non-local reaction–diffusion equations to study the effect of host heterogeneity on the phenotypic evolution of a pathogen population. The evolving phenotype is taken to be the transmission rate of the pathogen on the different hosts, and in our system there are two host populations present.

The central feature of our model is a trade-off relationship between the transmission rates on these hosts, which means that an increase in the pathogen transmission on one host will lead to a decrease in the pathogen transmission on the other. The purpose of the paper is to develop a classification of phenotypic diversity as a function of the shape of the trade-off relationship and this is achieved by determining the maximum number of phenotypes a pathogen population can support in the long term, for a given form of the trade-off. Our findings are then compared with results obtained by applying classical theory from evolutionary ecology and the more recent *adaptive dynamics method* to the same host–pathogen system. We find our work to be in good agreement with these two approaches.

## 1. Introduction

In Gudelj *et al*. (2004), motivated by the existence of sibling species among the fungal pathogens of agricultural crops, the authors studied the evolution of pathogen diversity in the following evolutionary host–pathogen system. Two sympatric crops act as hosts to a pathogen which is categorized into a number of possible types according to its transmission rate on the two hosts. Only one of these is used as an explicit measure of diversity, namely the transmission rate on host one, denoted *x* and taking on values in a discrete set contained within [0,1]. This value is termed the phenotype of the pathogen.

The underlying biological feature of the system is the following trade-off relationship: pathogens with higher transmission rates on host one have lower transmission rates on host two. In order to represent this relationship, a monotonically decreasing function *f*(*x*) is used to provide the transmission rate of a pathogen on host two directly from the value of *x*.

Gudelj *et al*. found that the number of pathogen types present and the values of their corresponding phenotypes depend on the shape of the trade-off function *f*. However, many of these results have been ascertained from numerical simulations conducted for particular trade-offs. In this paper we build on this work in order to obtain a qualitative description of what form of trade-off function will ensure the long-term presence of one, two or more phenotypes, where infinitely many are possible in principle. This ultimately provides a classification of pathogen diversity as a function of the trade-off relationship, something that would not be possible using the methods of Gudelj *et al*. (2004).

This is an important step because the experimental literature has for a long time acknowledged the importance of trade-offs in the evolution of host–parasite systems. Indeed, there are an increasing number of studies conducted to identify those components that are traded during the evolutionary process (e.g. Boots & Begon 1993; Kraaijeveld & Godfray 1997; Messenger *et al*. 1999). However, the shapes of these trade-offs are extremely difficult to obtain experimentally.

Let us now concentrate on the particular host–pathogen problem that we wish to study and consider the system of differential equations proposed in Gudelj *et al*. (2004):(1.1a)(1.1b)(1.1c)where *i* is an index which counts the number of resident pathogen phenotypes present in the system. The state variables *H*_{1} and *H*_{2} represent densities of healthy tissue of the two host populations, while *P*_{i} represents the density of host tissue infected with a pathogen of phenotype *x*_{i}. In the absence of pathogens, the host populations are governed by a simple *immigration and birth minus death* process of the form *r*−*μH*, which leads to the long-term persistence of each host. The parameter *γ* denotes the disease-induced mortality rate of the host tissue and it is assumed to be sufficiently small so that a non-trivial, positive steady state of (1.1) exists and is stable.

Equation (1.1) possesses a discrete number of phenotypes and the rules which describe the phenotypic mutation process and the conditions that define the fate of each mutant can be found in Geritz *et al*. (1998), this is known as the *adaptive dynamics* method. There are many other ways of representing mutations in phenotypic space and some examples can be found in Bonhoeffer & Nowak (1994), Sasaki & Godfray (1999), Bowers & Hodgkinson (2001) and more recently in Calsina & Cuadrado (2004). Calsina & Cuadrado (2004) consider the evolving trait to be a continuous variable and mutations are included by the use of an convolution operator with a mutation kernel giving the probability distribution of the offspring phenotype, for a given parental phenotype. This yields a system of integro-differential equations and the authors relate their findings to the discrete trait approach used in the adaptive dynamics method.

In Calsina *et al*. (1994), Calsina & Perelló (1995), Tsmiring & Levine (1996) and Saldaña *et al*. (2003), the authors use an approach where a diffusion operator represents phenotypic mutations. At an informal level this corresponds to applying a continuum limit to (1.1) in the case where mutations are given by nearest-phenotypic-neighbour connections. We also adopt this strategy and show that, in cases where a comparison is appropriate, our results agree with the findings of Gudelj *et al*. (2004) where the adaptive dynamics method is used, as well as many of the conclusions of classical evolutionary ecology (see Levins 1968; Lawlor & Maynard-Smith 1976). Note that while classical theory assumes that populations evolve to maximize their fitness, neither the model (1.1) nor the one presented in this paper possess any global notion of fitness. Moreover, the methodology used in Gudelj *et al*. that is applied to (1.1) refers to a local notion of fitness of a phenotype, while there no concept of fitness is used in the remainder of this paper.

## 2. A non-local PDE

Let us begin this section with the formal definition of a trade-off function: if satisfiesthen we call it a trade-off. From the mean-value theorem there is a *θ*∈(0,1) such that *f*′(*θ*)=−1 and we also make the following restriction on the oscillatory properties of *f*

If *Θ*=1, then *f* is concave or convex. If *Θ*=2, then we shall say that *f* is *sigmoidal* and if *Θ*≥3, then *f* is said to be a *staircase*. A sigmoidal trade-off is said to be *convex–concave* if there are intervals *I*_{1}=[0,*a*) and *I*_{2}=(*a*,1] such that *f* is convex on *I*_{1} and concave on *I*_{2}. A concave–convex sigmoidal trade-off is similarly defined.

We now propose a non-local PDE system as our model of evolution by natural selection. Mutations in phenotypic space are represented by a diffusion process with a constant diffusion coefficient *ϵ*, which is motivated directly by the work of Calsina *et al*. (1994) and Calsina & Perelló (1995). The system represents the equations of motion for the probability density of the phenotype resident in the pathogen population as a function of evolutionary time.

Let us denote the state variables *H*_{i}(*t*) for the density of host *i*, for *i*∈{1,2}, that is susceptible to infection at time *t*. Also, *P*(*x*,*t*) is the total density of hosts (both 1 and 2) infected with a pathogen with phenotypic trait *x*. Consequently, represents total biomass at time *t*. With this format we propose the following model:together with Neumann boundary conditions *P*_{x}(0,*t*)=*P*_{x}(1,*t*)=0 such that *H*_{1}(0)≥0, *H*_{2}(0)≥0, where and *P*_{0}(*x*)≥0 for all *x*∈[0,1]. The remaining parameters in (*HP*) are defined as in (1.1).

Although the host dynamics of (*HP*) are very simple, namely when *P*≡0, the results of this paper generalize to cover other types of interaction. We have used the particular form of (*HP*) in order to keep the exposition as straightforward as possible.

### (a) Species

So as not to confuse terminology surrounding the term *species* with any currently in use we shall use the following, less ambitious terminology.

By the term modal phenotype at time t for a solution of (HP) we understand a point *x*∈[0,1] such that *P*_{x}(*x*,*t*)=0 and *P*_{xx}(*x*,*t*)<0. A modal phenotype *x*∈(0,1) is said to be a *generalist* and a modal phenotype *x*∈{0,1} is said to be a *specialist*.

We consider the term *strategy* to be a synonym for phenotype.

### (b) Critical pathogenicity

Let us define some terminology. If *X* is any space of functions on [0,1] containing all continuous functions, a subscript zero, as in *X*_{0} denotes the fact that functions in *X* vanish at zero and one. Now define the Banach spaces , *X*≔*C*^{0}[0,1] and that will be needed, moreover *Z*_{+} will denote the cone of non-negative functions in *Z*. The Sobolev spaces for integer *k*, *p*≥1 will be used on occasion and we shall write *H*^{k}(0,1) for *W*^{k,2}(0,1). Here, *u*^{(k)} denotes the *k*th weak derivative and the operation of differentiating twice will sometimes be represented by the operator , as in and .

By *γ*_{0}, we denote the unique eigenvalue of the problem of finding a *ϕ*∈*H*^{2}(0,1) such that and(2.1)

Clearly, *γ*_{0} is positive and by well-known properties of eigenvalues of elliptic problems, it is an algebraically simple eigenvalue. The value *γ*_{0} will be called the *critical pathogenicity*. Let *ϕ*_{0} denote the eigenfunction corresponding to *γ*_{0} such that . Finally, let *ϕ*_{1}∈*H*^{2}(0,1) be the unique positive function that satisfies , and , for some eigenvalue *γ*_{1}>0. Of course, the value of *γ*_{0} depends on *f* and we may write *γ*_{0}(*f*) to emphasize this. It is also true that *γ*_{0}(*f*)>*γ*_{1} if *f*>0 and that *γ*_{0}(0)=*γ*_{1}.

The existence and positivity of solutions of (*HP*) is a simple consequence of the maximum principle and the theory of analytic semigroups, and the following lemma gives some simple boundedness properties which show that the hosts are persistent.

*Suppose that* (*P*_{0},*h*_{1},*h*_{2})∈*Z*_{+} *represents an initial datum for* (*HP*)*. There is a unique non-negative solution of* (*HP*), *with this initial datum. This solution satisfies*(2.2)

*If we define* *(which is finite* a priori*), then there is an upper bound*(2.3)*for all t*>0*, so that total biomass is bounded above. If* *(which is finite a priori), then there is a lower bound*(2.4)

*Consequently, the total biomass is eventually bounded below*

The existence, uniqueness and positivity property is a simple exercise in well-known theory of analytic semigroups (see Henry 1981, theorem 3.3.3, p. 54 and exercise 8, p. 61).

To prove (2.2), simply note that *H*_{i} satisfies the inequality , which integrates to give the desired bound. To prove (2.3), write and note thatand this inequality can be integrated to show that

To prove (2.4) for *i*=1, note that for some *M*>0 depending on initial conditions and therefore

This can be integrated to give

The bound on *H*_{2} is similarly obtained. ▪

The following lemma shows that if the pathogenicity has a greater value than the critical pathogenicity, the fate of all solutions of (*HP*) is to converge to the trivial steady state whereby the pathogen distribution is the zero function.

*The set* *is invariant for* (*HP*) *and any orbit in* *converges to* (0, *r*_{1}/*μ*, *r*_{2}/*μ*) *as t*→∞. *If γ*>*γ*_{0} *then any solution of* (*HP*) *converges eventually exponentially to* (0, *r*_{1}/*μ*, *r*_{2}/*μ*) *as t*→∞ *in the norm on Z given by* .

The first part of the lemma is obvious, so let us suppose that *γ*>*γ*_{0} and define the functional , evaluated along a solution of (*HP*). Now *V* is differentiable for *t*>0 andand by lemma 2.2, for every sufficiently small *η*>0, there exists a *T*>0 such that and , for all *t*>*T*. Consequently, there exists an *α*<0 such that , for all *t*≥*T* and therefore eventually. The result now follows. ▪

## 3. Classifying equilibrium states

Consider the equilibrium problem of (*HP*)(3.1)subject to the restrictions *P*>0, *P*′(0)=*P*′(1)=0. Throughout the paper we shall make use of the function(3.2)defined whenever *P* is a positive solution of (3.1).

*If (F1) holds and f*″(*x*)>0 *for all x*∈[0,1] *(f is strictly convex), then* *, r is concave and any zeros of r, of which there is at least one in* (0,1) *and at most two in* [0,1]*, are necessarily transverse*

*Moreover, r*′ *has at most one zero in* [0,1] *and it is transverse, so that r is either monotonic on* [0,1] *or has exactly one local extremum in* [0,1] *which is a global maximum*.

Since , we see that *r* is strictly concave and so has at most two zeros. However, since *P*′(0)=*P*′(1)=0, we find *P*″(*x*_{0})=0 for some *x*_{0}∈(0,1), whence *r*(*x*_{0})=0. Also, . The remainder of the proof is elementary and uses the mean-value theorem and the convexity of *f*. ▪

Similarly, we have the following.

*If (F1) holds and f*″(*x*)<0 *for all x*∈[0,1] *(f is strictly concave), then* *, r is convex and any zeros of r, of which there is at least one in* (0,1) *and at most two in* [0,1]*, are necessarily transverse. Moreover, r*′ *has at most one zero in* [0,1] *and it is transverse, so that r is either monotonic on* [0,1] *or has exactly one local extremum in* [0,1] *which is a global minimum*.

Denote the right-hand side of (3.1) by the smooth operator *F*(*P*,*γ*), where , then we can find a *δ*>0 such that *F* is an analytic mapping on the open set , where *P*≫−*δ* means that *P*(*x*)>−*δ* for all *x*∈[0,1]. Let us also define a solution seta set of positive solutions, , and a set of non-trivial solutions .

Note that if *f*∈*C*^{2} as per assumption (F1) and *P* is a solution of (3.1), then *P*∈*C*^{4} immediately follows. According to lemma 2.3, there can be no positive solution of (3.1) for *γ*>*γ*_{0}; we now sketch a proof that in fact *γ*=*γ*_{0} is a global bifurcation point to positive solutions of (3.1).

Let us first state a preliminary lemma.

*Suppose that* (*F1*) *holds,* (*P*,*γ*)∈*N and P*≥0 *then* (*P*,*γ*)∈*S*_{+}. *More generally, if* (*P*,*γ*)∈*N, then the zeros of P are transverse and therefore finite in number*. *Consequently, if C*⊂*N is a connected set that satisfies C*∩*S*_{+}≠∅*, then C*⊂*S*_{+}.

Suppose that (*P*,*γ*)∈*N* satisfies *P*(*x*_{0})=*P*′(*x*_{0})=0 for some *x*_{0}∈[0,1], then define numbers and . It follows that *P* is a solution of an ordinary differential equation with initial condition *P*(*x*_{0})=*P*′(*x*_{0})=0, but then *P*(*x*)≡0. This is a contradiction since (*P*,*γ*)∈*N* and therefore no such *x*_{0} exists. This immediately implies that non-negative solutions of (3.1) are strictly positive.

Now suppose that (*P*,*γ*)∈*C*∩*S*_{+} and define a mapping by . Then *i* is continuous on *C* by the transversality of zeros of *P* and integer valued. It is locally constant and therefore constant on the connected set *C*, so that *i*(*C* )=0 as *C* contains an element also in *S*_{+}. Hence, *C*⫅*S*_{+}. ▪

*For each γ*∈(0,*γ*_{0}), *the Neumann problem* *(3.1)* *has at least one solution. If we denote a locus of such solutions by P*_{γ}*, then*

*If C*(*γ*_{0})⊂*N is the maximal connected subset of N which contains the point* (*P*,*γ*)=(0,*γ*_{0}) *in its closure, then C*(*γ*_{0})⊂*S*_{+}.

(sketch) First, note that *F*(0,*γ*)≡0. Now, because where *L* is defined in (2.1) and , by the defining property of *γ*_{0} one may apply the theorem on bifurcation from a simple eigenvalue to obtain the existence of a local bifurcation. The Crandall and Rabinowitz global alternative (see Zeidler 1986) can then be used to demonstrate the existence of *C*(*γ*_{0}) in *S*_{+}. We omit the details as they are straightforward.

Integrating (3.1) over [0,1] shows that *γ*>0 must hold for any non-negative solution which is not identically zero. The second part then follows from a continuity argument because if there is a positive sequence (*γ*_{n}) such that *γ*_{n}→0 and (*P*_{n}) is a sequence of positive solutions of (3.1) which is *L*^{1}-bounded, the equation(where *α* is a constant chosen to ensure the invertibility of the linear operator −*ϵΔ*+*αI* independently of *n* in various spaces) yields a *W*^{2,1}-bound, hence a *C*^{1} bound and using (F1), a *C*^{3}-bound follows for (*P*_{n}) that is independent of *n*. As a result, there is a subsequence of (*P*_{n}) that converges in *C*^{2} to a non-negative solution of (3.1) when *γ*=0, which is a contradiction as none exist. ▪

### (a) Steep trade-offs

The following result shows that if the trade-off function is sufficiently *steep* in the sense of having most of its mass concentrated around *x*=0, then there is only one solution of (3.1) and it is close (in a *C*^{1} sense) to a strictly monotonic increasing function on (0,1) and so its mass is largely concentrated near *x*=1. This corresponds to solutions of (*HP*) that are close to one that supports only one modal phenotype at the specialist point *x*=1. Such behaviour was identified by Gudelj *et al*. (2004, fig. 3(b), p. 2190) using the adaptive dynamics approach and so we ought to expect to observe this result in (*HP*).

*Suppose that γ*<*γ*_{1}. *There is a ρ*>0 *such that if f is any trade-off function satisfying* *, then there is exactly one positive solution P*(=*P*(*f*)) *of* *(3.1)*. *Moreover, P*(*f*)→*p in H*^{2}(0,1) *(and hence in C*^{1}[0,1]*) as* *where p is the strictly monotonic increasing function on* (0,1) *that satisfies* *(3.1)* *with f*=0.

Fix *γ*<*γ*_{1}(≤*γ*_{0}(*f*)) and first consider an auxiliary problem.

Under the assumptions of the theorem, the problem(3.3)obtained by setting *f*=0 in (3.1) has a unique positive solution and it is strictly monotonic increasing on (0,1).

(of claim) The *existence* part of this claim follows from lemma 3.4. This shows that (3.3) has a global bifurcation to a branch of positive solutions emanating from the trivial solution at *γ*=*γ*_{1} and that 0<*γ*<*γ*_{1} holds along this branch.

The *uniqueness* part of the claim follows because if (3.3) has two positive solutions *p*,*q*>0 and we consider their Wronskian *w*≔*p*′*q*−*pq*′, then

Now *w*(0)=*w*(1)=0 and so *w*′(*x*_{1})=0 follows for some *x*_{1}∈(0,1) and therefore 〈*x*,*p*−*q*〉=0, so that *w*(*x*)≡0. From this it follows easily that *p*=*q*.

The *monotonicity* part of the claim is proven as follows. Since *p*′(0)=*p*′(1)=0 there is an *x*_{0}∈(0,1) such that *p*″(*x*_{0})=0 and therefore

This implies must be the unique point in (0,1) where *p*″ can be zero. If were true for some , then there would be two points in (0,1) where *p*″ was zero, one in and one in . This contradiction ensures that no such exists. Since *ϵp*″(0)=*γp*(0)>0 it follows that *x*=0 is a local minimum for *p*, consequently *p* must be strictly monotonic increasing on (0,1). ▪

Now for the proof of the theorem. Let us denote the solution of (3.3) by *p* and define the operator bywhere is a closed subspace of *H*^{2}(0,1). Now *N* is a *C*^{1} mapping and from the above claim *N*(0,*p*)=0, but *p*∈*H* satisfies *p*∈*Y* by a simple regularity argument and so provides the only solution of *N*(0,*P*)=0. Moreover, the Fréchet derivative of *N* with respect to *P* at (*f*,*P*)=(0,*p*) iswhence

We may write , where *g* is the smooth positive function and *A* is a self-adjoint, second-order linear differential operator such that *Ap*=0. Since *A* is a Fredholm operator of index-0, so too is *d*_{P}*N*(0,*p*) as it is a rank-one perturbation of *A* and to show that this derivative is an isomorphism we only need to prove injectivity. Now, *d*_{P}*N*(0,*p*)[*Q*]=0 yields

As a result, 〈*x*,*Q*〉=0 and therefore *AQ*=0 follows, from where *Q*=*αp* for some real constant *α* and the simplicity of the null space of *A*. The same integral identity yields *α*=0 and we deduce that *d*_{P}*N*(0,*p*) is an isomorphism from *H* to *L*^{2}, so that we can apply the implicit function theorem to solve *N*(*f*,*P*)=0 locally to (0,*p*) for *P*=*P*(*f*).

Local uniqueness is thus established, but for global uniqueness suppose that there are two distinct solutions of (3.1) and say, that are found for *f*=*f*_{n} where as *n*→∞. Note that holds by (F1) and so (*f*_{n}) tends to zero in each *L*^{p}-space for *p*≥1.

Let *P* be a solution of (3.1). By integrating over [0,1] the *f*-independent *L*^{1}-bound holds. Integrating from 0 to *z*∈(0,1) we obtainwhich is bounded above by *r*_{1}+*r*_{2}. The elementary inequality now givesas a result and since , we find . Integrating this givesso that . Choosing an such that is an isomorphism, a solution *P* of (3.1) satisfieswhere we have used 0≤*f*(*x*)≤1 that comes from (F1).

As *H* is a Hilbert space we can extract weakly convergent subsequences so that and in *H* as *n*→∞ for some *P*^{1}, *P*^{2}∈*H* and by compact embedding results, we can choose a new subsequence if necessary to ensure that strongly in *H*^{1}(0,1) as *n*→∞. Since *f*_{n}→0 in *L*^{2}, it follows from (3.1) thatin *L*^{2} for each *i*=1, 2. As a result, strongly in *H*^{2}(0,1) and because of the continuity of , *P*^{1} and *P*^{2} must both be solutions of *N*(*f*,*P*)=0 when *f*=0. But then *P*^{1}=*P*^{2} and it immediately follows from the previously proven local uniqueness result that for all sufficiently large *n* as both sequences must eventually lie in *H*-neighbourhoods of *p*. ▪

One could restrict the trade-off further in the statement of theorem 3.5 if we wished in such a way that the resulting operator *P*(*f*) would yield monotonic solutions of (3.1), rather than solutions that are *C*^{1}-close to being monotonic. However, we have not done this for brevity.

## 4. Convex and concave trade-offs

Let *f*∈*C*^{2}[0,1] and *P*∈*Y*∩*C*^{4}[0,1] be a corresponding positive solution of (3.1). In Gudelj *et al*. (2004), the authors deduced that the convexity or concavity properties of particular trade-offs resulted in the presence of one or two pathogen *strains*. The purpose of §§4 and 5 is to use the properties of the function *r* defined in (3.2) to deduce how many modal phenotypes a solution of (3.1) can support. Indeed, we deduce a similar result.

*Assuming* (*F1*) *and f*″(*x*)≠0 *for all x*∈[0,1] *(so that* (*F2*) *holds), then any solution of* *(3.1)* *has at most three distinct local extrema in* [0,1]. *Consequently, there is at most one local extremum of such a solution in* (0,1).

Suppose that a solution *P* of (3.1) has four, or more, local extrema , so that *P*′(*x*) vanishes at all of these points. This yields three zeros of *P*″ in (0,1) and therefore *r* too has three zeros. Now *r*′(*x*) has at most one zero since(4.1)and *f*′ is monotonic, a contradiction. Since the boundary conditions on *P* are Neumann, the boundary points provide two local extrema. ▪

It immediately follows from lemma 4.1 that any steady-state solution of (*HP*) has either one modal phenotype, which we call monomorphism, or two modal phenotypes, termed dimorphism. The manner in which monomorphism and dimorphism occurs when *f*″(*x*)≠0 is quantified in the following two propositions.

*If* (*F1*) *holds and f*″(*x*)>0 *for all x*∈[0,1] *then any solution of* *(3.1)* *is either monotonic or has a unique global minimum and two local maxima at x*=0 *and x*=1.

By lemma 4.1, either there exists an *x*_{0}∈(0,1) such that *P*′(*x*_{0})=0 or there does not. In the latter case *P* is monotonic, so let us suppose that such an *x*_{0} exists and note that it must be unique by lemma 4.1. From lemma 3.1 there are four cases to consider:

*r*has two zeros in (0,1),*r*has one zero in (0,1) and*r*(0)=0,*r*(1)≠0,*r*has one zero in (0,1) and*r*(0)≠0,*r*(1)=0 and*r*has one zero in (0,1) and*r*(0)≠0,*r*(1)≠0.

In case (1), since *f* is convex *r* is concave and therefore has *sign pattern* −+− on [0,1], as shown in the left-hand, topmost diagram of figure 1. Suppose that *r*(*x*)>0 on an interval *I*⊂(0,1) and *r* is zero on ∂*I*, if *x*_{0}∉*I*, then *P*″ must vanish outside , and therefore so too must *r*. This contradiction ensures that *x*_{0}∈*I*. Since , we see that *P*′(*x*_{0})=0 and *P*″(*x*_{0})>0, so that *x*=*x*_{0} is a local minimum for *P*. Since *P* is monotonic on (0,*x*_{0}) and (*x*_{0},1), *P* must have a global minimum at *x*=*x*_{0}. The boundary conditions ensure that *P* has two local maxima at the boundary of [0,1].

In cases (2–4), the mean-value theorem ensures that *P*′ cannot have a zero in (0,1) and *P* is therefore monotonic on (0,1). ▪

*If (F1) holds and f*″(*x*)<0 *for all x*∈[0,1] *then any solution of* *(3.1)* *is either monotonic or has a unique global maximum in* (0,1) *and two local minima at x*=0 *and x*=1.

Using lemma 3.2, the proof is entirely analogous to that of proposition 4.2. The only possible non-monotonic solution is illustrated in figure 2. ▪

### (a) Phenotypic structure as ϵ→0

In this section we perform an asymptotic analysis to locate the modal phenotype of (3.1) in the case where the trade-off function *f* is concave. We concentrate on a specific form of trade-off but, provided *f* is smooth, our analysis can be carried out for more general trade-off functions. We therefore have in mind equation (3.1) where

It is helpful to define new constants *Γ*_{1} and *Γ*_{2} byso that (3.1) becomes(4.2)

We are taking *r*_{1}, *μ* and *r*_{2} as given, but *Γ*_{1} and *Γ*_{2} are not; equation (4.2) is essentially a linear problem, but the nonlinearity will arise when we have to fix the values of *Γ*_{1} and *Γ*_{2}.

From the results obtained in §3 in the case of concave trade-offs, we anticipate that the solution of (4.2) is characterized by a bump centred at some point inside [0,1]. Therefore, we start by defining the coordinate *X* given bywe are expecting the solution to be concentrated in a layer of width (*ϵ*^{1/4}), but we do not yet know where, so *x*_{0} is unknown and is to be found in the course of the following analysis.

Let us expand formally(4.3)(4.4)(4.5)as we are expecting to obtain the elements of (4.2) in powers of *ϵ*^{1/4}. Suppose also that *f*(*x*) has a well-behaved Taylor series expansion about *x*_{0}, so we writewhere *f*_{0}=*f*(*x*_{0}), and .

An essential point is to recognize that the double derivative term in (4.2) is now formally (*ϵ*^{1/2}) (when the multiplying *ϵ* is taken into account) and we can then start to fix various constants appropriately. It follows that the expression *γ*−*Γ*_{1}*x*−*Γ*_{2}*f*(*x*) must vanish at (1) and (*ϵ*^{1/4}) and for this we require(4.6a)(4.6b)

Note that the second relation in (4.6*b*) is a higher-order effect which we shall not need again.

Taking the (*ϵ*^{1/2}) terms in equation (4.2) we obtain the expression

Let us definewhich is positive owing to the concavity of *f*(*x*) and note that we may translate *X* by defining , where is a constant chosen so that the linear term in the latter equation is removed. Details of the constant shift are not important for the leading-order analysis, but what is left is the determining equation(4.7)where is some constant that depends on the coefficients in the expansions of *Γ*_{1} and *Γ*_{2} in (4.4) and (4.5). Now we would like solutions of this equation that decay away from *x*_{0} that is as . This is perfectly feasible as the equation is just a scaled parabolic cylinder problem, the first solution of which isfor an appropriate and any constant . In determining the value of , a further restriction must be imposed on the coefficients in (4.4) and (4.5) because must hold for to be a solution of (4.7).

Having identified the leading-order eigenfunction, the stage is now set for us to tie down the actual location of the bump. Although we have not been concerned about the scaling of the eigenfunction yet, it is at this point that we observe that it must be of size (*ϵ*^{−1/4}) as it occupies a zone of width (*ϵ*^{1/4}) and we suspect that the overall integral of *P* is (1). Suppose thatthen as *P* is concentrated around *x*_{0}, we find at leading order that

Consequently, we deduce that and . Equation (4.6*b*) now gives(4.8)and (4.6*a*) tells us that(4.9)

Thus, we have found that for a given concave function the eigensolution is concentrated in a thin layer around *x*=*x*_{0}. This point is determined as the solution of (4.8) and (4.9) (for given *r*_{1}, *r*_{2} and *f*(*x*), the only unknowns are *x*_{0} and *A*). The eigensolution itself is a Gaussian profile whose decay rate is determined as a function of the underlying trade-off *f*(*x*).

## 5. Sigmoidal and staircase cases

*If* (*F1*) *and* (*F2*) *hold with Θ*=2*, then any solution of* *(3.1)* *has at most four distinct local extrema in* [0,1]. *Consequently, there are at most two local extrema of such a solution in* (0,1).

Suppose that a solution *P* of (3.1) has five, or more local extrema , so that *P*′(*x*) vanishes at all of these points. This yields four zeros of *P*″ in (0,1) and therefore *r* too has four zeros from which it follows that *r*′ has three zeros in (0,1). From (4.1) it follows that *f*″ has two zeros in (0,1), a contradiction. Since the boundary conditions on *P* are Neumann, the boundary points provide two local extrema. ▪

*If* (*F1*) *and* (*F2*) *hold with Θ*=2*, then r has at least one zero in* (0,1) *and at most three zeros in* [0,1] *counted according to multiplicity. Moreover, r*′ *has at most two zeros in* [0,1] *counted according to multiplicity, so that r is either monotonic, has exactly one local extremum in* [0,1] *or has two local extrema in* [0,1].

Since *P*′(0)=*P*′(1)=0, *P*″(*x*_{0}) for some *x*_{0}∈(0,1), whence *r*(*x*_{0})=0. Now suppose that *r* has four zeros in [0,1], . This yields three zeros of *r*′ in (0,1) and consequently two zeros of *r*″ in (0,1). Therefore, it follows that *f*″ has two zeros in (0,1), a contradiction.

The same contradiction is obtained if *r* has two non-transverse zeros, so that *r* has at most three zeros counted according to multiplicity and similar reasoning can be used to deduce that *r*′ has at most two zeros, either two transverse zeros or one double zero. ▪

*If* (*F1*) *and* (*F2*) *hold with Θ*=2 *and the trade-off f is sigmoidal convex–concave, then any solution of* *(3.1)* *is either*

*monotonic,**has a unique global minimum and two local maxima at x*=0*and x*=1,*has a unique global maximum and two local minima at x*=0*and x*=1*or**has one local maximum and one local minimum in*(0,1)*, with a local maximum at x*=0*and a local minimum at x*=1*(see**figure 3**a)*.

By lemma 5.1, either there exists an *x*_{0}∈(0,1) such that *P*′(*x*_{0})=0, there exist such that , or *P*′ has no zeros in (0,1). In the latter case *P* is monotonic, so let us suppose that *P*′ has one or two zeros in (0,1). From lemma 5.2 there are eight cases to consider, initially supposing that all the zeros of *r* are transverse:

*r*has one zero in (0,1) and*r*(0)=0,*r*(1)≠0,*r*has one zero in (0,1) and*r*(0)≠0,*r*(1)=0,*r*has one zero in (0,1) and*r*(0)≠0,*r*(1)≠0,*r*has one zero in (0,1) and*r*(0)=0,*r*(1)=0,*r*has two zeros in (0,1) and*r*(0)≠0,*r*(1)=0,*r*has two zeros in (0,1) and*r*(0)=0,*r*(1)≠0,*r*has two zeros in (0,1) and*r*(0)≠0,*r*(1)≠0 or*r*has three zeros in (0,1).

In cases (1–4), the mean-value theorem ensures that *P*′ cannot have a zero in (0,1) and *P* is therefore monotonic on (0,1). In cases (5–7), using the arguments from propositions 4.2 and 4.3, it follows that if *P* is not monotone, then *P* either has a unique global minimum and two local maxima at *x*=0 and *x*=1 or has a unique global maximum and two local minima at *x*=0 and *x*=1.

In case (8), it is possible for *P*′ not to have a zero in (0,1) and *P* is therefore monotone. However, if *P*′ has one zero, then it must have another one. Since *f* is sigmoidal convex–concave, *r* is sigmoidal concave–convex and therefore has *sign pattern* −+−+ on [0,1] as shown in figure 3*a* (see left-hand diagram). Therefore, there are *x*_{1}, *x*_{2} and *x*_{3} such that and *r*(*x*)>0 on (*x*_{1},*x*_{2}) and *r*(*x*)<0 on (*x*_{2},*x*_{3}), while . Suppose that the two zeros of *P*′ in (0,1) are *x*_{0} and . The mean-value theorem ensures that *x*_{0}∈(*x*_{1},*x*_{2}) and . Since , *x*=*x*_{0} is a local minimum and is a local maximum. The boundary conditions ensure that *x*=0 is a local maximum and *x*=1 is a local minimum, as described in case (iv) in the statement of the theorem.

If *r* has a non-transverse zero of multiplicity 3, then this situation is analogous to case (3) above. The case whereby *r* has one transverse and one double zero is inadmissible as then *P* would have two consecutive non-degenerate local minima. This covers all possible non-transverse cases by lemma 5.2. ▪

*If* (*F1*) *and* (*F2*) *hold with Θ*=2 *and f is sigmoidal concave–convex, then any solution of* *(3.1)* *is either*

*monotonic,**has a unique global minimum and two local maxima at x*=0*and x*=1,*has a unique global maximum and two local minima at x*=0*and x*=1*or**has one local maximum and one local minimum in*(0,1)*, with a local minimum at x*=0*and a local maximum at x*=1*(see**figure 3**b)*.

The proof is entirely analogous to proposition 5.3. ▪

From the preceding sequence of results the following can be deduced directly.

## 6. Concluding remarks

The asymptotic analysis carried out in §4*a* gives a leading-order approximation to the location of the modal phenotype for smooth, concave trade-offs, which are also evolutionary stable strategies (ESSs) in this context (for a definition of an ESS see Lawlor & Maynard-Smith 1976). We can compare this with the location predicted by Gudelj *et al*. (2004) using the adaptive dynamics method and we now give the details of such a calculation.

Let us assume that *f* is a concave trade-off, so that from Gudelj *et al*. (2004), we know that the adaptive dynamics of (1.1) possesses a single ESS at some phenotypic location in [0,1]. Let us consider the case where *i*∈{1,2}, so that *P*_{1} represents the resident and *P*_{2} the mutant pathogen densities. The resulting four-dimensional system of ODEs has the mutant-free invariant set , and when restricted to this flow has the globally attractive steady state , whereand *P*_{1}(*x*) is given by *P* in the equation(6.1)

Now, the introduction of a mutant corresponds to an examination of the flow of (1.1) when perturbed away from . However, is unstable as an invariant manifold if * e* is unstable as an equilibrium, so we consider the right-most eigenvalue of the linearization of (1.1) about

*. This is the quantity given by (*

**e***x*

_{1},

*x*

_{2}), where is often called the local fitness of the mutant.

From eqn (3) of Lawlor & Maynard-Smith (1976), the phenotypic location of the ESS occurs where the rate of change of the fitness with respect to variations in the mutant phenotype is zero, namely at that value *x* which satisfies

This is the equation(6.2)where *P* is obtained from (6.1).

In §4*a*, the location of the *unique modal phenotype* was shown to be determined to leading order by (4.8) and (4.9) when the trade-off is concave, but these equations are identical to (6.2) and (6.1), respectively. Thus, we have deduced that the location of the ESS obtained in Gudelj *et al*. which uses the fitness maximization procedure of Lawlor & Maynard-Smith (1976), agrees to lowest order with our asymptotic findings, in the context of concave trade-off functions. Moreover, our classification also accords with the findings of Gudelj *et al*. even in the cases where *f* is not concave, but convex and even sigmoidal, although our results are more general because Gudelj *et al*. often rely on specific choices of trade-off function.

Although the classification of phenotypic diversity given in this paper relies heavily on the second derivative of the trade-off, it only provides an upper bound on the number of modal phenotypes that exist for a given trade-off. Changes in other physical parameters could lead to a diminishing of phenotypic diversity. Consider, e.g. figure 5*a*, where for small *γ* the pathogen has a bimodal distribution with a specialist phenotype situated on the left boundary at *x*=0. As *γ* increases the modal phenotype at this phenotype vanishes and the pathogen population becomes unimodal (see figure 5*b*) and infective to both hosts. In this case the trade-off is sigmoidal, but the same comment applies to other shapes too.

It would be a laudable, if intractable, goal to completely describe the set of modal phenotypes as a function of *t* for a given initial phenotypic distribution. To formalize a discussion in this direction, let us denote by *σ*(*t*)⊂[0,1] the set of modal phenotypes at time *t*, i.e.

If there is a *T*>0 and a sequence *ϵ*_{n}↓0 such that for all *n* (where a hash (#) denotes set cardinality), then *diversification* is said to have taken place at *t*=*T*. On the other hand, if , then *extinction* is said to have occurred.

Now, a function such that *s*(*t*)∈*σ*(*t*) for all *t*>0 is called a selection function. If diversification takes place at *t*=*T* and there exists a selection of *σ*(.) which is discontinuous, then an evolutionary *leap* is said to have taken place at *t*=*T*. On the other hand, if every selection function of *σ*(.) is continuous at *t*=*T*, then evolutionary *branching* is said to have taken place.

Different diversification and extinction events can be characterized according to bifurcations of the smooth bifurcation problem *P*_{x}(*x*,*t*)=0, with *t* playing the role of parameter. It is clear that both evolutionary branching *and* leaping events can occur in smooth solutions of (*HP*) and are caused by well-known codimension-1 singularities. For instance, evolutionary branching occurs when *P*_{x}(*x*,*t*)=0 has a pitchfork bifurcation, whereas an evolutionary leap occurs when the same equation has a saddle-node bifurcation point (see figure 6 for pictorial representations of both cases). From this observation, it is apparent that evolutionary leaping ought to be the generic diversification mechanism and robust to perturbations in the initial state of the phenotypic distribution, whereas evolutionary branching should be a rare, non-generic event. A related discussion can be found in Stewart (2003) where the author reaches a similar conclusion but starting from different modelling assumptions. We also note that an evolutionary leap has some resonance with Gould's concept of punctuated equilibrium.

## Acknowledgments

I.G. was supported by a NERC-EMS fellowship, grant no. NE/B501998/1. C.D.C. was supported by a Nuffield Newly Appointed Lecturers grant, number NAL/00688/G.

- Received October 5, 2004.
- Accepted August 1, 2005.

- © 2005 The Royal Society