Existence and qualitative properties of travelling waves for an epidemiological model with mutations

# Existence and qualitative properties of travelling waves for an epidemiological model with mutations

Quentin Griette 111Centre d’Écologie Fonctionnelle et Évolutive, UMR 5175, CNRS, 1919 Route de Mende, 34293 Montpellier, France. E-mail: quentin.griette@cefe.cnrs.fr  222I3M, Université Montpellier 2, place Eugène Bataillon, 34095 Montpellier, France. E-mail quentin.griette@univ-montp2.fr Gaël Raoul 333Centre d’Écologie Fonctionnelle et Évolutive, UMR 5175, CNRS, 1919 Route de Mende, 34293 Montpellier, France. E-mail: raoul@cefe.cnrs.fr.
###### Abstract

In this article, we are interested in a non-monotone system of logistic reaction-diffusion equations. This system of equations models an epidemics where two types of pathogens are competing, and a mutation can change one type into the other with a certain rate. We show the existence of minimal speed travelling waves, that are usually non monotonic. We then provide a description of the shape of those constructed travelling waves, and relate them to some Fisher-KPP fronts with non-minimal speed.

## 1 Introduction

Epidemics of newly emerged pathogen can have catastrophic consequences. Among those who have infected humans, we can name the black plague, the Spanish flu, or more recently SARS, AIDS, bird flu or Ebola. Predicting the propagation of such epidemics is a great concern in public health. Evolutionary phenomena play an important role in the emergence of new epidemics: such epidemics typically start when the pathogen acquires the ability to reproduce in a new host, and to be transmitted within this new hosts population. Another phenotype that can often vary rapidly is the virulence of the pathogen, that is how much the parasite is affecting its host; Field data show that the virulence of newly emerged pathogens changes rapidly, which moreover seems related to unusual spatial dynamics observed in such populations ([24, 33], see also [29, 25]). It is unfortunately difficult to set up experiments with a controlled environment to study evolutionary epidemiology phenomena with a spatial structure, we refer to [4, 27] for current developments in this direction. Developing the theoretical approach for this type of problems is thus especially interesting. Notice finally that many current problems in evolutionary biology and ecology combine evolutionary phenomena and spatial dynamics: the effect of global changes on populations [32, 13], biological invasions [36, 26], cancers or infections [21, 18].

In the framework of evolutionary ecology, the virulence of a pathogen can be seen as a life-history trait of the pathogen [34, 17]. To explain and predict the evolution of virulence in a population of pathogens, many of the recent theories introduce a trade-off hypothesis, namely a link between the parasite’s virulence and its ability to transmit from one host to another, see e.g. [3]. The basic idea behind this hypothesis is that the more a pathogen reproduces (in order to transmit some descendants to other hosts), the more it ”exhausts” its host. A high virulence can indeed even lead to the premature death of the host, which the parasite within this host rarely survives. In other words, by increasing its transmission rate, a pathogen reduces its own life expectancy. There exists then an optimal virulence trade-off, that may depend on the ecological environment. An environment that changes in time (e.g. if the number of susceptible hosts is heterogeneous in time and/or space) can then lead to a Darwinian evolution of the pathogen population. For instance, in [6], an experiment shows how the composition of a viral population (composed of the phage and its virulent mutant cl857, which differs from by a single locus mutation only) evolves in the early stages of the infection of an E. Coli culture.

The Fisher-KPP equation is a classical model for epidemics, and more generally for biological invasions, when no evolutionary phenomenon is considered. It describes the time evolution of the density of a population, where is the time variable, and is a space variable. The model writes as follows:

 ∂tn(t,x)−σΔn(t,x)=rn(t,x)(1−n(t,x)K). (1)

It this model, the term models the random motion of the individuals in space, while the right part of the equation models the logistic growth of the population (see [39]): when the density of the population is low, there is little competition between individuals and the number of offsprings is then roughly proportional to the number of individuals, with a growth rate  ; when the density of the population increases, the individuals compete for e.g. food, or in our case for susceptible hosts, and the growth rate of the population decreases, and becomes negative once the population’s density exceeds the so-called carrying capacity . The model (1) was introduced in [16, 28], and the existence of travelling waves for this model, that is special solutions that describe the spatial propagation of the population, was proven in [28]. Since then, travelling waves have had important implications in biology and physics, and raise many challenging problems. We refer to [42] for an overview of this field of research.

In this study, we want to model an epidemics, but also take into account the possible diversity of the pathogen population. It has been recently noticed that models based on (1) can be used to study this type of problems (see [7, 2, 9]). Following the experiment [6] described above, we will consider two populations: a wild type population , and a mutant population . For each time , and are the densities of the respective populations over a one dimensional habitat . The two populations differ by their growth rate in the absence of competition (denoted by in (1)) and their carrying capacity (denoted by in (1)). We will assume that the mutant type is more virulent than the wild type, in the sense that it will have an increased growth rate in the absence of competition (larger ), at the expense of a reduced carrying capacity (smaller ). We assume that the dispersal rate of the pathogen (denoted by in (1)) is not affected by the mutations, and is then the same for the two types. Finally, when a parent gives birth to an offspring, a mutation occurs with a rate , and the offspring will then be of a different type. Up to a rescaling, the model is then:

 ⎧⎨⎩∂tw(t,x)−Δxw(t,x)=w(t,x)(1−(w(t,x)+m(t,x)))+μ(m(t,x)−w(t,x)),∂tm(t,x)−Δxm(t,x)=rm(t,x)(1−(w(t,x)+m(t,x)K))+μ(w(t,x)−m(t,x)), (2)

where is the time variable, is a spatial variable, , and are constant coefficients. In (2), represents the fact that the mutant population reproduces faster than the wild type population if many susceptible hosts are available, while represents the fact that the wild type tends to out-compete the mutant if many hosts are infected. Our goal is to study the travelling wave solutions of (2), that is solutions with the following form :

 w(t,x)=w(x−ct),m(t,x)=m(x−ct),

with . (2) can then be re-written as follows, with :

 (3)

The existence of planar fronts in higher dimension () is actually equivalent to the case (), our analysis would then also be the first step towards the understanding of propagation phenomena for (2) in higher dimension.

There exists a large literature on travelling waves for systems of several interacting species. In some cases, the systems are monotonic (or can be transformed into a monotonic system). Then, sliding methods and comparison principles can be used, leading to methods close to the scalar case [40, 41, 35]. The combination of the inter-specific competition and the mutations prevents the use of this type of methods here. Other methods that have been used to study systems of interacting populations include phase plane methods (see e.g. [38, 15]) and singular perturbations (see [20, 19]). More recently, a different approach, based on a topological degree argument, has been developed for reaction-diffusion equations with non-local terms [5, 2]. The method we use here to prove the existence of travelling wave for (3) will indeed be derived from these methods. Notice finally that we consider here that dispersion, mutations and reproduction occur on the same time scale. This is an assumption that is important from a biological point of view (and which is satisfied in the particular phage epidemics that guides our study, see [6]). In particular, we will not use the Hamilton-Jacobi methods that have proven useful to study this kind of phenomena when different time scales are considered (see [30, 7, 9]).

This mathematical study has been done jointly with a biology work, see [23]. We refer to this other article for a deeper analysis on the biological aspects of this work, as well as a discussion of the impact of stochasticity for a related individual-based model (based on simulations and formal arguments).

We will make the following assumption,

###### Assumption 1.1.

, and .

This assumption ensures the existence of a unique stationary solution of (2) of the form (see Appendix A.2). It does not seem very restrictive for biological applications, and we believe the first result of this study (Existence of travelling waves, Theorem 2.1) could be obtained under a weaker assumption, namely:

 r∈(1,∞),K∈(0,1),μ∈(0,K).

Throughout this document we will denote by and the terms on the left hand side of (3):

 fw(w,m):=w(1−(w+m))+μ(m−w),fm(w,m):=rm(1−(w+mK))+μ(w−m). (4)

We structure our paper as follows : in Section 2, we will present the main results of this article, which are three fold: Theorem 2.1 shows the existence of travelling waves for (3), Theorem 2.2 describes the profile of the fronts previously constructed, and Theorem 2.3 relates the travelling waves for (3) to travelling waves of (1), when and are small. sections 3, 4 and 5 are devoted to the proof of the three theorems stated in Section 2.

## 2 Main results

The first result is the existence of travelling waves of minimal speed for the model (2), and an explicit formula for this minimal speed. We recall that the minimal speed travelling waves are often the biologically relevant propagation fronts, for a population initially present in a bounded region only ([10]), and it seems to be the one that is relevant when small stochastic perturbations are added to the model ([31]). Although we expect the existence of travelling waves for any speed higher than the minimal speed, we will not investigate this problem here - we refer to [5, 2] for the construction of such higher speed travelling waves for related models. Notice also that the convergence of the solutions to the parabolic model (2) towards travelling waves, and even the uniqueness of the travelling waves, remain open problems.

###### Theorem 2.1.

Let satisfy Assumption 1.1. There exists a solution of (3), such that

 ∀x∈R,w(x)∈(0,1),m(x)∈(0,K),
 liminfx→−∞(w(x)+m(x))>0,limx→∞(w(x)+m(x))=0,
 c=c∗,

where

 c∗:=√2(1+r−2μ+√(r−1)2+4μ2) (5)

is the minimal speed for which such a travelling wave exists.

The difficulty of the proof of Theorem 2.1 has several origins:

• The system cannot be modified into a monotone system (see [38, 8]), which prevents the use of sliding methods to show the existence of traveling waves.

• The competition term has a negative sign, which means that comparison principles often cannot be used directly.

As mentioned in the introduction, new methods have been developed recently to show the existence of travelling wave in models with negative nonlocal terms (see [5, 2]). To prove Theorem 2.1, we take advantage of those recent progress by considering the competition term as a nonlocal term (over a set composed of only two elements : the wild and the virulent type viruses). The method of [5, 2] are however based on the Harnack inequality (or related arguments), that are not as simple for systems of equations (see [12]). We have thus introduced a different localized problem (see (13)), which allowed us to prove our result without any Harnack-type argument.

Our second result describes the shape of the travelling waves that we have constructed above. We show that three different shapes at most are possible, depending on the parameters. In the most biologically relevant case, where the mutation rate is small, we show that the travelling wave we have constructed in Theorem 2.1 is as follows: the wild type density is decreasing, while the mutant type density has a unique global maximum, and is monotone away from this maximum. In numerical simulations of (2), we have always observed this situation (represented in Figure 1), even for large . This result also allows us to show that behind the epidemic front, the densities and of the two pathogens stabilize to , , which is the long-term equilibrium of the system if no spatial structure is considered. For some results on the monotony of solutions of the non-local Fisher-KPP equation, we refer to [14, 1]. For models closer to (2) (see e.g. [2, 7]), we do not believe any qualitative result describing the shape of the travelling waves exists.

###### Theorem 2.2.

Let satisfy Assumption 1.1. There exists a solution of (3) such that

 limx→−∞(w(x),m(x))=(w∗,m∗),limx→∞(w(x),m(x))=(0,0),

where is the only solution of .

The solution satisfies one of the three following properties:

(a)

is decreasing on , while is increasing on and decreasing on for some ,

(b)

is decreasing on , while is increasing on and decreasing on for some ,

(c)

and are decreasing on .

Moreover, there exists such that if , then there exists a solution as above which satisfies .

Finally, we consider the special case where the mutant population is small (due to a small carrying capacity of the mutant, and a mutation rate satisfying ). If we neglect the mutants completely, the dynamics of the wild type would be described by the Fisher-KPP equation (1) (with ), and they would then propagate at the minimal propagation speed of the Fisher-KPP equation, that is . Thanks to Theorem 2.1, we know already that the mutant population will indeed have a major impact on the minimal speed of the population which becomes , and thus shouldn’t be neglected. In the next theorem, we show that the profile of is indeed close to the travelling wave of the Fisher-KPP equation with the non-minimal speed , provided the conditions mentioned above are satisfied (see Figure 2). The effect of the mutant is then essentially to speed up the epidemics.

###### Theorem 2.3.

Let , , and (see Theorem 2.1 for the definition of ), , a solution of (3) such that

 liminfx→−∞(w(x)+m(x))>0,limx→∞(w(x)+m(x))=0.

There exists , and such that if , then

 ∥w−u∥L∞⩽CKβ,

where is a traveling wave of the Fisher-KPP equation, that is a solution (unique up to a translation) of

 (6)

with speed

The Theorem 2.3 is interesting from an epidemiological point of view: it describes a situation where the spatial dynamics of a population would be driven by the characteristics of the mutants, even though the population of these mutants pathogens is very small, and thus difficult to sample in the field.

## 3 Proof of Theorem 2.1

We will prove Theorem 2.1 in several steps. We refer to Remark 2 for the conclusion of the proof.

### 3.1 A priori estimates on a localized problem

We consider first a restriction of the problem (3) to a compact interval , for . More precisely, we consider, for ,

 ⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩w,m∈C0([−a,a]),−cw′−w′′=fw(w,m)χw≥0χm≥0,−cm′−m′′=fm(w,m)χw≥0χm≥0,w(−a)=w∗,m(−a)=m∗,w(a)=m(a)=0, (7)

where we have used the notation (4), and are defined in the Appendix, see Subsection A.2.

#### 3.1.1 Regularity estimates on solutions of (7)

The following result shows the regularity of the solutions of (7).

###### Proposition 3.1.

Let satisfy Assumption 1.1 and . If satisfies

 {−cw′−w′′=fw(w,m),−cm′−m′′=fm(w,m), (8)

on , where are defined by (4), and , then

###### Proof of Proposition 3.1.

Since for any , the classical theory ([22], theorem 9.15) predicts that the solutions of the Dirichlet problem associated with (8) lies in This shows that for any But then for any (thanks to Sobolev embeddings). It follows that is a function of the variable (see (4) for the definition of ). Let us choose one such . Now we can apply classical theory ([22], theorem 6.14) to deduce that . But then and verify some uniformly elliptic equation of the type

 −c(w′′)′−(w′′)′′=g,
 −c(m′′)′−(m′′)′′=h,

with , and we can apply again ([22], theorem 6.14). This argument can be used recursively to show that for any , so that finally, . ∎

#### 3.1.2 Positivity and L∞ bounds for solutions of (7)

In this subsection, we prove the positivity of the solutions of (7), as well as some bounds.

###### Proposition 3.2.

Let satisfy Assumption 1.1, , and . If is a solution of (7), then and satisfy are positive, that is and for all .

###### Proof of Proposition 3.2.

We observe that

 fw(w,m)=w(1−(w+m))+μ(m−w)=w(1−μ−w)+m(μ−w),

so that if then Let such that and the connex compound of the set that contains Since over and the weak minimum principle imposes , and thus But then reaches its global minimum at so the strong maximum principle imposes that or else would be constant. We deduce then from our hypothesis that That shows that in

To show that , we notice that

 fm(w,m) = rm(1−w+mK)+μ(w−m) = m(r−μ−rKm)+w(μ−rKm),

so that if then . The end of the argument to show the positivity of can the n be reproduced to show that .

###### Proposition 3.3.

Let satisfy Assumption 1.1, , and . If is a positive solution of (7), then and satisfy

 ∀x∈(−a,a),w(x)<1,
 ∀x∈(−a,a),m(x)
###### Proof of Proposition 3.3.

Let a positive solution of (7).

• We assume that there exists such that Let then the connex compound of the set that contains Then in we have

 −cw′−w′′ = w(1−μ−w−m)+μm⩽w(−μ−m)+μm = m(μ−w)−μw⩽0,

along with so that the weak maximum principle states in which is absurd because Therefore, for all

• We assume that there exists such that Let then the connex compound of the set that contains Then in we have

 −cm′−m′′ = m(r−μ−rK(w+m))+μw⩽m(−μ−rwK)+μw = w(μ−rKm)−μm⩽0,

Thanks to Assumption 1.1. Since , the weak maximum principle states in which is absurd because Therefore, for all .

• Now if , we still have the estimate

 −cw′(x)−w′′(x)⩽m(x)(μ−w(x))+w(x)(1−μ−w(x))⩽0,

so that if there exists such that , then is locally equal to thanks to the strong maximum principle. But in that case

 0=(−cw′−w′′)(x0)=−m(x0)+μ(m(x0)−1)<0,

which is absurd. Hence, . Similarly, if , we get

 0=(−cm′−m′′)(x0)=−Kμ+w(x0)(μ−r)<0,

which is absurd, and thus .

#### 3.1.3 Estimates on solutions of (7) when c≥c∗ or c=0

The next result shows that the solutions of (7) degenerate when if the speed is larger than a minimal speed (see Theorem 2.1 for the definition of ).

###### Proposition 3.4 (Upper bound on c).

Let satisfy Assumption 1.1. There exists such that for and , any solution of (7) satisfies

 ∀x∈[−a,a],max(w(x),m(x))≤Ce−c−√c2−c2∗2(x+a).
###### Proof of Proposition 3.4.

Let , and

 M:=(1−μμμr−μ).

Since is a positive matrix, the Perron-Frobenius theorem implies that has a principal eigenvalue and a positive principal eigenvector (that is for ), given by

 (9)

The function with and is then a solution of the equation

 −cψ′η−ψ′′η=Mψη=h+ψη.

We can define , which is a closed subset of . is non-empty since and are bounded while for .

Consider now . Then , , and there exists such that either or . We first consider the case where . Then

 −c(w−(ψη)1)′(x0)−(w−(ψη)1)′′(x0)⩽−w(x0)(w(x0)+m(x0))⩽0

over . The weak maximum principle ([22], theorem 8.1) implies that

 sup[−a,a](w−(ψη)1)=max((w−(ψη)1)(−a),(w−(ψη)1)(a)),

and then, thanks to the definition of , . Since this means that , and thus

 η0=b−w1−r+√(1−r)2+4μ2eλ−a.

The argument is similar if , which concludes the proof. ∎

The following Proposition will be used to show that .

###### Proposition 3.5.

Let satisfy Assumption 1.1, and . Every positive solution of (7) with satisfies the estimate

 max[−a0,a0](w+m)⩾K2(1−μ). (10)
###### Proof of Proposition 3.5.

We assume that , , and that (10) does not hold. We want to show that those assumptions lead to a contradiction. For the function defined by

 ψA(x)=Acos⎛⎝√1−μ2x⎞⎠,

is a solution of the equation over . Since over and are bounded, the set is a closed bounded nonempty set in . Let now We still have over , and then, since (10) does not hold and ,

 −(w−ψA0)′′ ⩾ (1−max[−a0,a0](w+m)−μ)w−1−μ2ψA0 ⩾ 1−μ2(w−ψA0)⩾0.

 −(m−ψA0)′′⩾1−μ2(m−ψA0)⩾0.

The weak minimum principle ([22], theorem 8.1) then imposes

 min(inf[−a0,a0](w−ψA0),inf[−a0,a0](m−ψA0)) =min((w−ψA0)(−a0),(w−ψA0)(a0),(m−ψA0)(−a0),(m−ψA0)(a0)).

But the left side of the equation is by definition of while the right side is strictly positive since . This contradiction shows the result. ∎

###### Remark 1.

Notice that Propositions 3.1, 3.2,3.3, 3.4 and 3.5 also holds if is a solution of

 (12)

where .

### 3.2 Existence of solutions to a localized problem

To show the existence of travelling waves solutions of (3), we will follow the approach of [2]. The first step is to show the existence of solutions of (7) satisfying the additional normalization property , that is the existence of a solution to

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩(c,w,m)∈R×C0([−a,a])2,−cw′−w′′=fw(w,m)χw⩾0χm⩾0,−cm′−m′′=fm(w,m)χw⩾0χm⩾0,w(−a)=w∗,m(−a)=m∗,w(a)=m(a)=0,max[−a0,a0](w+m)=ν0, (13)

where are defined by (4), and are defined in Appendix A.2.

We introduce next the Banach space , with and . We also define the operator

 Kσ:X⟶X,(c,w,m)⟼(c+max[−a0,a0](~w+~m)−ν0,~w,~m) (14)

where is the unique solution of

The solutions of (13) with are then the fixed points of in the domain .

We define

 Ω:={(c,w,m)∈R+×C0([−a,a])2;c∈(0,c∗),∀x∈[−a,a],−1

where is defined by (5).

###### Lemma 3.6.

Let satisfy Assumption 1.1, and . Then, , defined by (14), is a family of compact operators on , that is continuous with respect to .

###### Proof of Lemma 3.6.

We can write where is defined by

 (LD)−1(c,g,h)=(~c,~w,~m),

where is the unique solution of

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩−c~w′−~w′′=g on (−a,a),−c~m′−~m′′=h on (−a,a),~w(−a)=w∗,~m(−a)=m∗,~w(a)=~m(a)=0,~c=c+max[−a0,a0](~w+~m)−ν0,

and is the mapping

 Fσ(c,w,m)=(c,w(1−(w+σm))+μ(σm−w),rm(1−σw+mK)+μ(σw−m)).

is a continuous mapping from to , and is a continuous application from into itself (see Lemma A.2), it then follows that is a continuous mapping from to . Finally, the operator is compact (see Lemma A.2), which implies that is compact for any fixed . ∎

We now introduce the following operator, for :

 Fσ:=Id−Kσ. (15)

Similarly, we introduce the operator

 Kτ:X⟶X,(c,w,m)⟼(c+max[−a0,a0](~w+~m)−ν0,~w,~m) (16)

where is the unique solution of

 (17)

The argument of Lemma 3.6 can be be reproduced to prove that is also a continuous family of compact operators on , and we can define, for , the operator

 Fτ:=Id−Kτ. (18)

Finally, we introduce, for some that we will define later on,

 ~Ω:={(c,w,m)∈R+×C0([−a,a])2;c∈(¯c,c∗),∀x∈[−a,a],−1

In the next Lemma, we will show that the Leray-Schauder degree of in the domain is non-zero as soon as is large enough. We refer to chapter 12 of [37] or to chapter 10-11 of [11] for more on the Leray-Schauder degree.

###### Lemma 3.7.

Let satisfy Assumption 1.1. There exists such that the Leray-Schauder degree of in the domain is non-zero as soon as .

###### Proof of Lemma 3.7.

We first notice that for , the solution of (17) is independent of , and then,

 F0(c,w,m)=(ν0−max[−a0,a0](wc+mc),w−wc,m−mc),

where is the solution of (17) with , that is

 (wc,mc)(x):=(w∗(e−cx−e−caeca−e−ca),m∗(e−cx−e−caeca−e−ca)),

for , and for . The solutions of then satisfy and . In particular, the solutions of satisfy and on , and then,

 (c,w,m) ∉ {(~c,~w,~m)∈R×C0([−a,a])2;∃x∈[−a,a],~w(x)∈{−1,1}} ∪{(~c,~w,~m)∈R×C0([−a,a])2;∃x∈[−a,a],~m(x)∈{−K,K}}).

The solutions of also satisfy

 max[−a0,a0](wc∗+mc∗)≤2ec∗a0ec∗a−1,

so that if for some . It follows that has no solution in , provided . Finally, for , the solutions of satisfy , so that

 max[−a0,a0](wc+mc)>max[−a0,a0](w0+