Existence of traveling wave solutions for a nonlocal reaction-diffusion model of influenza A drift

In this paper we discuss the existence of traveling wave solutions for a nonlocal reaction-diffusion model of Influenza A proposed in Lin et. al. (2003). The proof for the existence of the traveling wave takes advantage of the different time scales between the evolution of the disease and the progress of the disease in the population. Under this framework we are able to use the techniques from geometric singular perturbation theory to prove the existence of the traveling wave.


1.
Introduction. Recently work has been developed in ecological and epidemiological models using integro-differential equations to model spatial processes. Kendall [9], studied the traveling wave solutions for an integro-differential SIR model, where the integral term models the nonlocal contacts in among the susceptible and infectious populations. In his model Kendall assume that the rate of infection is given by: where I(x, t) represent the density of the infected individuals at time t, and K(x) is a kernel function of mass one that measure the contribution of the infected individuals at location y to the infection of susceptible individuals at location x. The integral models the fact that two individuals can influence each other even if their separation in subtypes is large. For more about these epidemiological models the reader should The variables in the model are represented by: S(t) represent the number of susceptible host at time t, while i(x, t) and r(x, t) are the distribution of infected individuals with variant x and recovered individuals from variant x at time t, respectively. Observe that the nonlinear terms in (2b,2c) are given by the nonlocal terms. The nonlocal term in (2b) represent the per capita growth rate of infection by hosts carrying variant x, and the diffusion term in (2b) is the drift of the virus. In [13], the authors carried out a preliminary analysis that suggested the existence of traveling wave solutions for the system. In our paper, we want to formalize those findings in a more mathematical way, at least for certain class of kernels. In [13], the authors used the Monoid kernel K(v) = v/(v + a), while in this paper we proposed a more standard exponential kernel that satisfied the conditions on [13].
The proof for the existence of traveling waves solutions for (2) uses geometrical singular perturbation theory [3,8] to study the dynamics of the system. The analysis of the proof for the existence of the traveling wave solutions takes the advantage of the smallness of the diffusion coefficient. It is realistic to assume that the kernel, because recovered individuals are protected against close by variants of the virus than from those further away has the form: for some a > 0. The parameter a measures the nonlocality of the disease transmission, so that for a large this kernel gives a strong nonlocality, i.e. if a is large strains of the value that are far away in the variants spaces become more important. The main result of this paper is the existence of traveling wave solution for (2). The first part of the argument considers the extended system, that is (2) and the equations given by the new variables (5) below. Thus first, by studying the slow-fast dynamics of the system for the case when ǫ = 0 we establish the existence of the traveling wave in the critical manifold.
Theorem 1.1. The unperturbed system when ǫ = 0 has a heteroclinic solution connecting the quasi steady-states (0, 0, 0, For the perturbed system (2), we construct the so-called slow manifold in an ǫ-neighborhood of the critical manifold. By a careful study of the asymptotic behavior of the solution we were able to transform the problem of the existence of the solution to a corresponding fixed point problem of a nonlinear map. We apply the Contraction Mapping Principle in some weighted Banach space to show that the existence of the traveling wave solution is also true for ǫ > 0, but small. Theorem 1.2. System (51) has a traveling wave solution, and furthermore the solution is in a neighborhood of the unperturbed solution.
The rest of the paper is organized in the following manner: in Section 2 we discuss the structure of the traveling wave, and some preliminaries analysis for the existence of the traveling wave solutions for the whole system. The construction of the traveling wave solutions on the critical manifold. In Section 3 we prove the robustness of the traveling wave for the perturbed system.

2.
Structure of the traveling waves. In traveling wave form, with z = x − ct, the original system can be reformulated as: with the boundary conditions (i(±∞), r(±∞)) = (0, 0), such traveling wave solutions are sometime referred to as traveling pulses. Recall that we are interested in finding a homoclinic orbit connecting the origin at z → ±∞. We assume that, in this framework, S * represents a quasi steady-state for the susceptible population, and we only need to focus on the equations (4b) and (4c). Defining the new variables: Observe that K ′ (x) = a(1 − K(x)), thus by computing the derivatives of R and I we obtain: Thus, by computing the second derivative we have that R(u) satisfies the following second order ODE: Similarly, we obtain that I(u) satisfies the second ODE: Adding these new equations to (4b-c), and replacing the integral terms by the new variables (5) we obtain the following system of equations.
The trade off of doing this is that the number of equations is increased from two to four. Analytical experiments suggest that ǫ is small while a is large when compare to ǫ. Thus, we assume that the cross-protection parameter a is inversely proportional to the mutation rate ǫ, that is a = m ǫ . Converting system (8) to a first order system of equations, we have Note that the problem of finding a homoclinic orbit for the original system (2) or (4) connecting the origin at ±∞, becomes a problem of finding a heteroclinic orbit for the system (9) that connects the steady-states solutions (0, 0, 0, 0, 0, 0, I −∞ ) to (0, 0, 0, 0, 0, R +∞ , 0), for some I −∞ > 0, and R +∞ > 0 So, Thus, integrating (13) from z to 0, or To help constructing the weighted Banach spaces in Section 3 let us compute the eigenvalues of the asymptotic systems at (0, 0, 0, 0, 0, 0, I −∞ ) and at (0, 0, 0, 0, 0, R +∞ , 0). At −∞ the linear system is given by: And the set of eigenvalues for the matrix A −∞ is given by, Observe that there are two negative eigenvalues, three positive eigenvalues, and two zero eigenvalues. This implies that at (0, 0, 0, 0, 0, 0, I −∞ ) there is a 3-dimensional unstable manifold, a 2-dimensional stable manifold, and there is also a 2-dimensional center manifold. Similarly, for the asymptotic system at +∞ the linearized matrix is given by

JOAQUIN RIVERA AND YI LI
The eigenvalues of A +∞ are given in the following set Thus, there are three positive eigenvalues, two negatives and two zero eigenvalues.
From (17), and (18) we obtain that: and the two eigenvalues that we will use to estimate the weights in Section 3.

Remark 1.
We can think of the heteroclinic orbit connecting (0, 0, 0, 0, 0, 0, I −∞ ) and (0, 0, 0, 0, 0, R +∞ , 0) as the intersection of the unstable manifold of the asymptotic system at −∞ with stable manifold of the asymptotic system at +∞. Observe that the critical speed of the wave c = c(ǫ) depends on the parameter ǫ, as it is illustrated in (19), and (20) 2.2. Construction of the traveling wave solution. Note that the system (9) can be formulated in slow-fast variables by setting the new variables x = (i 1 , R 1 , I 1 ) T and y = (i, r, R, I) T , which are the fast and slow variable, respectively. Thus, rewriting system (9) as: where the functions f , and g are given by the right hand sides of system (9). The motivation for this separation is to take advantage of the geometric singular perturbation methods to study the dynamics for the fast and slow system. Setting ǫ = 0 we take the critical manifold to be the set given by: for some c > 0. Thus, we obtain the equations for the slow flow on the critical manifold M 0 by substituting the relations on M 0 into the remaining equations: The only condition that we need check on M 0 is that it is normally hyperbolic.
Lemma 2.1. The critical Manifold M 0 is normally hyperbolic relative to the fast system.
Proof. Note that we only need to show that the linearization for the fast system evaluated on M 0 has four eigenvalues with zero real part and three with nonzero real part. Thus, it suffices to show that D x f (x,ŷ, 0) has three eigenvalues with Re(λ) = 0 for any (x,ŷ) ∈ M 0 The Jacobian matrix for f (x,ŷ, 0) is given by: which has nonzero eigenvalues λ = −c, ±m. Hence M 0 is normally hyperbolic.
From the Lemma we have that the critical manifold has a 2 dimensional stable manifold and a 1 dimensional unstable manifold. Thus by applying all of the Fenichel's Theorem there is, for some ǫ > 0 but small a manifold M ǫ which is diffeomorphic and close to M 0 . Before we deal with the perturbed system in M ǫ , we will discuss the dynamics on the critical manifold, and then we will connect it to the dynamics on the slow manifold.
3. Dynamics of the model.

3.1.
Dynamics on the critical manifold. In this section we show the existence of a traveling wave solution of (24) connecting the equilibrium points (0, 0, 0, I −∞ ) and (0, 0, R +∞ , 0). Note that this is a four dimensional system of equations, so the typical techniques to establish the existence of heteroclinic curves, like phase-plane analysis, are not applicable here. In order to tackle the problem we will rewrite the system as a single equation of fourth order, then construct a solution for that equation.
Note that by substituting the equations I ′ = −i, and R ′ = r into the equations for i ′ and r ′ we obtain the new second order equations: Now, observe that we can solve for R in the first equation obtaining the relation: thus, by computing the firsts and second order derivative of R, by using (27), i.e. and Thus by substituting this into the the second equation we have the following fourth order differential equations on I(z): where I (j) denotes the j-th derivative of I with respect to z. Observe that now we only need to find a heteroclinic solution of (30) that connects the points 0 to I * , where I * = I −∞ . This may seems like a more difficult problem than the original system (24), since it is of higher order nonlinear equation. But by a change of variables we are able to reduce the equation (30) to a second order singular equation.
Let v = I(z) and dI dz = g(v). Computing all the derivatives in (30), and combining all the terms we have that g(v) is the solution of the following singular second order ordinary differential equation: dv + C 1 , for some constants C 1 and C 2 . Notice that singular nature of equation is not due to limiting behavior as ǫ → 0, but to the fact the equation changes its order at any point v 0 ∈ (0, I * ] such that g(v 0 ) = 0. By using the asymptotic behavior of I(z) at ±∞ we can obtain the boundary conditions for (31) in [0, I * ]. Thus, as z → −∞, it follows that v → I * , so g(v) → g(I * ) = 0. Similarly, as z → +∞ we have that g(v) → 0. Thus, the boundary conditions for (31) are g(0) = 0 = g(I * ).
From equation (24d) we obtain that I(z) is a decreasing function and recall that g(v) = dI dz < 0, and since g(0) = 0 = g(I * ) it must be the case that there is a v 0 so that min v∈[0,I * ] g(v) = g(v 0 ). By (31) we can obtain the following approximations for the function near 0 and I * : where M , and N are positive constants, and the exponents α, and β are to be determined. We now substitute (32a) and (32b) into (31) to estimate the exponents α, and β. Thus, for (32a) we obtain: Multiplying (33) by v 2−2α we obtain, where h.o.t denotes the remaining higher order terms. So, if α > 1 we have that Now, if α < 1 we have from (34) that Therefore, α = 1. Similarly, we can deduce that β = 1. Considering the test solution: for some M > 0. Substituting (37) into (31) and solving for the constants terms C 2 , M , and I * to get a solution. Thus, we obtain the constants: as a test solution, where M , and I * are given by (38). By using this function and by reversing the change of variables we obtain the following: for some choice of C 1 , and f = I(z). Thus after solving (39) for f we obtain that Thus substituting the constants (38) into (40), and letting ρ = 1 − 2e c , we can simplify (40) as From (41) we can obtain R(z) just by substituting the first and second order derivative of (41) into (27), so Since we are looking for such a solution R(z) with the property lim z→−∞ R(z) = 0, we must have ρc − 1 + R 0 S * = 0. And by solving for S * we obtain S * = 2e R 0 . After simplifying (42) we have Finally from (24c) and (24d) we can get the functions r(z) and i(z) It is easy to check that the set {i(z), r(z), R(z), I(z)} are solutions of the system (24) with the correct asymptotics. Hence, we have just constructed heteroclinic orbit on the critical manifold M 0 that connects the points (0, 0, 0, I −∞ ) and (0, 0, R +∞ , 0).
The fast system is given by Note that if we let ǫ = 0 the new reduced system is given by: with (i, r, R, I) all constants. Thus, we can argue that the only dynamics of interest for the case when ǫ = 0 is given by the slow system. For that reason the problem of interest when ǫ > 0, but small, is to show that the traveling wave solution that we constructed in the previous sections persists for the perturbed problem.

3.2.
Dynamics on the slow manifold. The purpose of this section is to study the dynamics of the perturbed system (8) for ǫ > 0, but sufficiently small. By Lemma 2.1 we have that the critical manifold M 0 is normally hyperbolic, so the existence of a slow manifold M ǫ follows by the invariant manifold theorem. The following characterization gives the relation between the slow manifold and the critical manifold: Note that by invariant manifold theorem we have that the functions f i are smooth. Thus, since ǫ is sufficiently small, and by taking the Taylor expansion of the f i 's about ǫ = 0, we can rewrite the perturbations terms as: for i = 1, 2, 3. Notice that f i1 (i, r, I, R, 0) = 0. We also disregard the higher order terms for the moment. Then, we approximate the perturbed terms in (47) by only the the linear terms in (48). Thus, the first order approximation for the perturbation is: The goal is to use this perturbations to study the flow on M ǫ . Now, we want to find explicit expression for the f i on (49).
Substituting the perturbed terms (49) into (9a) and letting ǫ = 0 we obtain the functional expressions for f i . Thus, the first order perturbations are given by and the flow is given by Note that system (51) is a regular perturbation of the slow system (24).

3.3.
Construction of the singular solution. Now, we will look for solutions of (51) near the unperturbed solutions (41), (43), and (44). We will also show that the solutions on the slow manifold M ǫ are perturbations of the solutions of the unperturbed problem. In order to simplify the expressions we let ǫ = δc ρ , where δ ≪ 1. Thus, we obtain the following system of equations: where 0 < δ * is fixed.
First, we consider a linear system associated to (52) where the functions ĩ (u),r(u),R(u),Ĩ(u) represent known functions near the unperturbed solutions. Notice that by integrating the equations in (53) we obtain the following integral operator: 2R 0Ĩ (t) + 1 − ρc, and x = (i, r, I, R) ′ . Let X be the space of continuous functions C(R). Define the weighted Banach space to be the space X equipped with the following norm: where the exponents γ 1 < 0 and γ 2 > 0 are the asymptotic values of the growth function α(t) at −∞, and +∞, respectively. We denote the Banach space (X, ) by X. Let x c be the unperturbed solution. If for x ∈ X 4 such that ||x − x c || ≤ δ for any δ ≤ δ * , then the asymptotic behavior of the function α(t) is given by where the terms H i (t) decay exponentially fast as t → ±∞. Let B δ0 = (i, r, I, R) ∈ X 4 : ||(i, r, I, R) − (i c , c , R c )|| ≤ δ 0 be the ball of radius δ 0 center at the unperturbed solution with 0 < δ < δ 0 .
It is easy to check that: The argument for the existence of the singular solution will follow directly from the Contraction Mapping Principle. By examining the definition of T we can observe that all the terms in each of its components have δ or δ * as a coefficient. This would be fundamental in the proof that T is contraction for 0 < δ < δ * << 1. The following two results summarize the complete argument that we are using to establish the existence of the fixed point solution.
Proposition 2. The T operator is a contraction on B δ0 for all sufficiently small δ > 0.
In order to prove Proposition 2 we need the following Lemma.
Lemma 3.1. The growth function α(t) is a Lipshitz function.
for the i function. Note that α only depends on R, so we will rewrite it as: for some constants A, B so that A > R 0 c (1 − 2δ), and B > δR 2 0 c 2 ρ . Also note that Therefore, the growth function is Lipshitz with constant M .

Remark 2.
As a consequence of Lemma 3.1 and to the fact that the weight function W (u) = e γ1u + e γ2u grows exponentially at both u = ±∞, we have that Proof of Proposition 2 We will prove the required condition on T by proving that each component of T is a contraction. The first component is α(s)ds i(t)r(t)dt.

Conclusion.
Despite the simplicity of the ideas of the model considered in this paper, we can observe that very interesting mathematics arise from that model due to various factors: (1) The difference in time scales between the mutation process and the epidemiological process, (2) The nonlocal nonlinearities, and (3) The dependence of the traveling wave on the parameters. This analysis suggests that the pulse approaches 0 as z → ±∞, but the speed of the wave depends on ǫ for ǫ > 0, but small. The analysis for the solution is the result of a singular perturbation argument, and therefore no conclusion can be drawn for larger ǫ.
In conclusion, we can state that the pulse solutions are robust, in the sense that they persist under small perturbation, and they look qualitatively the same as the unperturbed solution.