# Traveling Waves in a Hybrid Reaction-Diffusion-Difference SEIR Epidemic Model With Nonlocal Transmission and Protection Phase With Delay

Hanene Meghelli, Abdennasser Chekroun, Mostafa Adimy

► To cite this version:

Hanene Meghelli, Abdennasser Chekroun, Mostafa Adimy. Traveling Waves in a Hybrid Reaction-Diffusion-Difference SEIR Epidemic Model With Nonlocal Transmission and Protection Phase With Delay. Mathematical Methods in the Applied Sciences, 2025, 48 (14), pp.13380-13398. 10.1002/mma.11110. hal-05447865

HAL Id: hal-05447865

https://inria.hal.science/hal-05447865v1

Submitted on 8 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Traveling waves in an hybrid reaction-diffusion-difference SEIR epidemic model with nonlocal transmission and protection phase with delay

Hanene Meghelli<sup>a</sup>, Abdennasser Chekroun<sup>a*</sup> and Mostafa Adimy<sup>b</sup>

March 26, 2025

## Abstract

In this paper, we study the existence and non-existence of traveling waves for a SEIR epidemic model with diffusion. This system is coupled to an age-structured protection equation that takes into account a proportion of susceptible individuals who protect themselves from the disease through medical treatment or vaccination with temporary immunity, for example. The model includes a compartment of exposed individuals to represent diseases with an incubation period. In addition, the model contains a nonlinear saturated incidence rate to describe the interaction between susceptible and infected individuals. By solving the equation associated with the exposed compartment, we obtain a nonlocal integral representation of the exposed population as a function of the nonlinear incidence rate. This reformulation reduces the model to an hybrid system of reaction-diffusion equations coupled to a continuous difference equation that includes a delay and a nonlocal term, reflecting the movement of individuals during the protection and exposure phases. For wave speeds above a certain threshold, we construct lower and upper solutions for the resulting hybrid system. Finally, we apply Schauder's fixed point theorem in an appropriate function space to establish the existence of progressive wave solutions.

MSC: 35C07, 35K57, 35Q92, 92D30

**Keywords:** SEIR epidemic model, age-structured model, reaction-diffusion system with delay, traveling waves.

## 1 Introduction

Mathematical modeling in epidemiology is an important tool for understanding the dynamics of infectious diseases and predicting their spread within populations. It is also useful for testing hypotheses, evaluating the effectiveness of public health interventions and planning appropriate responses to epidemics [10]. Mathematical models in epidemiology are generally based on differential equations describing interactions between different compartments of individuals in

<footnote><sup>a</sup> Laboratoire d'Analyse Nonlinéaire et Mathématiques Appliquées, University of Tlemcen, Tlemcen 13000, Algeria</footnote>

<footnote><sup>b</sup> Inria, CNRS, Ecole Centrale de Lyon, INSA Lyon, Universite Claude Bernard Lyon 1, Université Jean Monnet, ICJ UMR5208, 69603 Villeurbanne, France</footnote>
a population. Most of the systems used are based on compartment models. These include the well-known SIR (Susceptible-Infectious-Recovered) model, or variants such as SIS, SEIR (with exposed individuals) and other models with additional compartments to capture more detailed dynamics, see the books [6, 10, 20, 22] for several other epidemiological models.

For the record, the SIR epidemic model was first introduced by Kermack-McKendrick in 1927 [19]. This is one of the most fundamental mathematical models in epidemiology, used to predict the spread of an epidemic within a population. It should be noted that the basic Kermack-McKendrick model does not take into account population demographics, nor many other factors such as natural mortality, mortality due to disease, vaccination, protection, etc. In this case, the total population $N(t) = S(t) + I(t) + R(t)$ is assumed to be constant. The basic reproduction number $R_0$ is a threshold that is generally introduced to determine the dynamics of disease transmission within a population. More precisely, it represents the number of infected individuals produced from a single infected individual during its lifetime [11]. For the SIR model, the epidemic potential is obtained as a function of this number: if $R_0 > 1$, the number of infected individuals first increases to a maximum, then decreases to zero, giving rise to an epidemic. However, when $R_0 < 1$, the number of infected individuals decreases to zero, so there is no epidemic. It should be noted that such a situation and this threshold result are common to most epidemiological models.

The Kermack-McKendrick SIR model described above is based on the assumption that the population is spatially homogeneous and does not move. In practice, however, individuals may come into contact with infected people in different locations and be exposed to the disease. Many authors introduce the spatial variable to study the spread of an epidemic, see for example [25, 26], for more details. Reaction-diffusion models are the systems most frequently used to describe population movements during an epidemic. The Kermack-McKendrick model [19] has been extensively studied, taking into account the effects of spatial diffusion of individuals. For example, the following reaction-diffusion epidemiological model was considered by Hosono and Ilyas in 1995, [17], for $t > 0$ and $x \in \mathbb{R}$,

$$ \left\{ \begin{aligned} \frac{\partial}{\partial t} S(t, x) &= \Delta S(t, x) - \beta S(t, x) I(t, x), \\ \frac{\partial}{\partial t} I(t, x) &= d \Delta I(t, x) + \beta S(t, x) I(t, x) - \gamma I(t, x). \end{aligned} \right. \qquad (1) $$

These authors proved that if $S(0, x) = S_0$ is a given constant initial condition such that $\beta S_0 / \gamma > 1$, then for each

$$ c \geq c^* := 2\sqrt{\beta S_0 d \left(1 - \frac{\gamma}{\beta S_0}\right)}, $$

the system (1) has a traveling wave solution $(S(x + ct), I(x + ct))$, where $c$ is the wave speed. This solution satisfies the following limit conditions $S(+\infty) = S_0$, $S(-\infty) = \epsilon < S_0$ and $I(\pm\infty) = 0$. The spread of epidemics in time and space is associated with the existence of this type of solution, called traveling waves, which represent the propagation of a profile in a specific direction at a constant speed [23, 29]. On this subject, many authors have obtained relevant information on the existence and non-existence of traveling wave solutions. In particular, most of them have provided minimum wave speeds. For more detailed results on this subject, we refer you to the following references, [4, 7, 12, 13, 16, 21, 28, 31, 32, 36, 37, 39, 38].

Vaccination is one of the most effective ways of stopping the transmission of infectious diseases, by boosting the population's immunity. Vaccination strategies are adapted to the
type of infectious pathogen, the transmission mechanism and the most appropriate age to
optimize disease control. Recently, for example, a vaccination campaign against infection by
the COVID-19 virus was organized to contain the epidemic. To maintain the immunity provided
by certain vaccines, it is sometimes necessary to update them regularly. Our model takes into
account the regular renewal of protection against an infectious disease. Many researchers have
taken the duration of protection into account in their epidemiological models. In general, to
take account of the period of protection, authors use terms with a time delay (duration of
protection, whether by vaccination, medication or even confinement), see [2, 14, 15] for some
studies on the subject.

In this paper, we consider a general diffusive Kermack-McKendrick epidemic model in-
troduced in [4] that takes temporary protection into account. As for many diseases, see
[18, 27, 34, 35, 37], we take into account the incubation period by introducing an additional
compartment of individuals exposed to the disease but not yet contagious. We use the sat-
uration incidence rate, which is more realistic, rather than the bilinear incidence rate. The
originality of our model is that it combines these four aspects in a simple model: a temporary
protection phase, an exposure phase that can be more or less long depending on the disease,
an incidence with saturation, and space. Hepatitis B, which has an incubation period of 1 to
6 months, and Tuberculosis, that has an incubation period of several weeks to several years
and can sometimes remain inactive for decades, are two well-known examples of this type of
disease. During the incubation period, infected individuals are not contagious because the virus
(or bacteria) present in their organism is not at a high enough level to allow transmission. It is
therefore important to take into account the incubation period in order to properly study the
propagation of these diseases. This is particularly important when it comes to protecting sus-
ceptible individuals (through vaccination, for example) and/or taking into account the spatial
diffusion of individuals. These infectious diseases with their specific characteristics motivated
us to develop this type of model. However, our system remains generic and can be adapted
to several other infectious diseases. This work also follows our previous study [1] on a hybrid
differential-difference system describing a SEIR model but without taking space into account.
The main results of [1] concern the local and global stability of disease-free and endemic steady
states.

This paper is organized as follows. In Section 2, we first formulate the age- and space-
structured model and use the characteristic method to reduce the main model to a specific
hybrid differential-difference system with delay. We also transform the equation for exposed
individuals into an integral equation with a nonlocal term. The result is a system composed of
two reaction-diffusion equations and a difference equation, with nonlocal terms that take into
account the movement of exposed and protected individuals. In Section 3, we introduce some
preliminary results on traveling wave solutions for our reduced hybrid system. In section 4,
using an iterative method for monotone operators, we construct a lower and an upper solution
that will be used later for the existence of traveling waves. In Section, 5, we reformulate the
problem of the existence of traveling waves as a fixed-point problem. This allows us to show the
existence of a suitable solution and thus to provide the existence of traveling wave solutions.
We also establish the non-existence conditions for traveling waves. Finally, we end this paper
with a discussion in Section 6.
## 2 Model formulation

We first present the basic elements of our epidemiological model. The total population is divided into five classes: Susceptible $S(t, x)$, Exposed $E(t, x)$, Infected $I(t, x)$, Recovered $R(t, x)$ and Protected individuals $P(t, x)$, where $t \ge 0$ and $x \in \mathbb{R}^n$, $n \ge 1$, represent time and location, respectively. Our model focuses on studying the propagation of a short-lived infectious disease, so birth and natural mortality are neglected. The evolution of $S, E, I, R$ and $P$ is represented by a reaction-diffusion system. We suppose that $D_S, D_E, D_I, D_R$ and $D_p$ describe the diffusion coefficients in each compartment. Susceptible individuals protect themselves temporarily with a rate $h$, e.g. by getting vaccinated or taking medication, for example. The incidence rate is represented by the function $I \mapsto g(I)$, which satisfies certain conditions that will be presented later. The parameter $\mu$ describes the infection rate of exposed individuals, which means that $1/\mu$ describes the average incubation time of the disease. It is assumed that individuals will be immunized at a rate $r$ (recovery rate). This means that the average duration of infection is $1/r$. The duration of the protection phase is supposed to be equal to $\tau$. We suppose that a part of the protected individuals renew their protection at the end of this phase, at a rate $\eta \in [0, 1)$. We denote by $p(t, x, a)$ the age distribution of protected individuals, where $a \in [0, \tau]$ is the time spent by an individual in this temporary protection phase. The evolution of density $p$ is described by the following equation,

$$ \frac{\partial p}{\partial t} + \frac{\partial p}{\partial a} = D_p \Delta p, \quad t > 0, \ x \in \mathbb{R}^n \text{ and } 0 < a < \tau, $$

satisfying the following boundary condition,

$$ p(t, x, 0) = hS(t, x) + \eta p(t, x, \tau), \quad t > 0, \ x \in \mathbb{R}^n. $$

This condition expresses the fact that not only new individuals in the susceptible compartment protect themselves, but also already protected individuals that renew their protection at the end of this phase $a = \tau$. We have the following formulation of the total population of protected individuals, for $t > 0$ and $x \in \mathbb{R}^n$,

$$ P(t, x) = \int_{0}^{\tau} p(t, x, a) da. $$

All the parameters used to describe the interactions between the compartments are positive and are summarized in the following table.

<table><thead><tr><th>Parameters</th><th>Description</th></tr></thead><tbody><tr><td>$D_S, D_E, D_I, D_R, D_p$</td><td>Diffusion coefficients</td></tr><tr><td>$h$</td><td>Temporary protection rate (vaccination, medication, etc.)</td></tr><tr><td>$0 \le \eta < 1$</td><td>Renewal protection rate</td></tr><tr><td>$\mu$</td><td>Infection rate of exposed individuals</td></tr><tr><td>$r$</td><td>Infection recovery rate</td></tr><tr><td>$\tau$</td><td>Duration of protection phase</td></tr></tbody></table>

Using the heat kernel and the method of characteristics, we solve the equation of $p$ for a large time ($t > \tau$), to obtain

$$ p(t, x, \tau) = \int_{\mathbb{R}^n} \Gamma_p(\tau, x - y)p(t - \tau, y, 0)dy, \quad t > 0, \ x \in \mathbb{R}^n, \qquad (2) $$
where the fundamental function $\Gamma_p$ is given by

$$ \Gamma_p(t, x) = \frac{e^{-|x|^2/(4D_p t)}}{(4\pi D_p t)^{n/2}}, \quad t > 0, \ x \in \mathbb{R}^n. $$

The other components of the SEIR epidemiological model are represented by the following system of equations, for $t > 0$ and $x \in \mathbb{R}^n$,

$$ \begin{aligned} \frac{\partial S}{\partial t} &= D_S \Delta S - hS(t, x) - S(t, x)g(I(t, x)) + (1 - \eta)p(t, x, \tau), \\ \frac{\partial E}{\partial t} &= D_E \Delta E + S(t, x)g(I(t, x)) - \mu E(t, x), \\ \frac{\partial I}{\partial t} &= D_I \Delta I + \mu E(t, x) - rI(t, x), \\ \frac{\partial R}{\partial t} &= D_R \Delta R + rI(t, x). \end{aligned} \quad (3) $$

Using the same approach as in [2, 4, 5] and putting $u(t, x) := p(t, x, 0)$, the equation (2) becomes, for $t > \tau > 0$ and $x \in \mathbb{R}^n$,

$$ p(t, x, \tau) = \int_{\mathbb{R}^n} \Gamma_p(\tau, x - y)u(t - \tau, y)dy. $$

This last expression and a translation in time give the following system, for $t > 0$ and $x \in \mathbb{R}^n$,

$$ \left\{ \begin{aligned} \frac{\partial S}{\partial t} &= D_S \Delta S - hS(t, x) - S(t, x)g(I(t, x)) + (1 - \eta) \int_{\mathbb{R}^n} \Gamma_p(\tau, x - y)u(t - \tau, y)dy, \\ \frac{\partial E}{\partial t} &= D_E \Delta E + S(t, x)g(I(t, x)) - \mu E(t, x), \\ \frac{\partial I}{\partial t} &= D_I \Delta I + \mu E(t, x) - rI(t, x), \\ \frac{\partial R}{\partial t} &= D_R \Delta R + rI(t, x), \\ u(t, x) &= hS(t, x) + \eta \int_{\mathbb{R}^n} \Gamma_p(\tau, x - y)u(t - \tau, y)dy. \end{aligned} \right. $$

The initial conditions are

$$ S(0, x) = S_0(x), \quad E(0, x) = E_0(x), \quad I(0, x) = I_0(x), \quad R(0, x) = R_0(x), $$

and

$$ u(t, x) = u_0(t, x), \quad \text{for} \quad -\tau \le t \le 0. $$

The above system is a coupled system between four reaction-diffusion equations including non-local terms in time and space and a continuous difference equation. The difference equation describes the number of new individuals entering the protection phase; on the one hand, susceptible individuals who protect themselves from the disease at a rate $h$, and on the other, already protected individuals who renew their protection at the end of this phase at a rate $\eta$. Note that, the nonlocal term

$$ \int_{\mathbb{R}^n} \Gamma_p(\tau, x - y)u(t - \tau, y)dy $$
includes the movement of individuals during the protection period. The case where $\eta = h = 0$ leads to the classic SEIR model with diffusion already discussed in several previous works, see for example [27, 34, 35, 37].

Throughout this paper, the function $g$, which represents the nonlinear incidence rate, is assumed to be continuously differentiable and to satisfy the following assumptions,

$$g'(y) \ge 0 \text{ and } 0 \le g(y) \le g'(0)y, y \in \mathbb{R}_+,$$ $$g(y) \ge g'(0)(y - ay^2), y \in [0, 1/a], \text{ for some } a \in \mathbb{R}_+^*,$$

and

$$\sup_{z \ge 0} g(z) + \sup_{z \ge 0} g'(z) < +\infty.$$

It is clear that $g(0) = 0$. The above conditions are satisfied by the classical function

$$g(I) = \frac{\beta I}{1 + aI}, \quad \beta, a \in \mathbb{R}_+^*,$$

known as the saturation incidence rate, see [7, 37].

The equations of $S$, $E$, $I$ and $u$ are independent of the equation of $R$. We will therefore only consider the system of $S$, $E$, $I$ and $u$. That is, for $t > 0$ and $x \in \mathbb{R}^n$,

$$\left\{
\begin{aligned}
\frac{\partial S}{\partial t} &= D_S \Delta S - hS(t, x) - S(t, x)g(I(t, x)) + (1 - \eta) \int_{\mathbb{R}^n} \Gamma_p(\tau, x - y)u(t - \tau, y)dy, \\
\frac{\partial E}{\partial t} &= D_E \Delta E + S(t, x)g(I(t, x)) - \mu E(t, x), \\
\frac{\partial I}{\partial t} &= D_I \Delta I + \mu E(t, x) - rI(t, x), \\
u(t, x) &= hS(t, x) + \eta \int_{\mathbb{R}^n} \Gamma_p(\tau, x - y)u(t - \tau, y)dy.
\end{aligned}
\right.
\quad (4)$$

Independently of the initial conditions, we have the following integral representation of $E$, see [27], for $t > 0$ and $x \in \mathbb{R}^n$,

$$E(t, x) = \int_{0}^{+\infty} \int_{\mathbb{R}^n} k(s, y)S(t - s, x - y)g(I(t - s, x - y))dyds, \quad (5)$$

where the modified heat kernel $k$, associated to the operator $D_E\Delta - \mu$, is given by

$$k(s, y) := \frac{e^{-\mu s - |y|^2/(4D_E s)}}{(4\pi D_E s)^{n/2}} = e^{-\mu s} \Gamma_E(s, y), \quad s > 0, \ y \in \mathbb{R}^n.$$

We note here that in the expression (5), we calculate $E(t, x)$ using the entire past of populations $S$ and $I$ between $-\infty$ and $t$. We introduce the following notation, which considers a term that is twice nonlocal in time and space,

$$(k * * f)(t, x) := \int_{0}^{+\infty} \int_{\mathbb{R}^n} k(s, y)f(t - s, x - y)dyds.$$
The integral in (5) can then be rewritten as follows,

$$E(t, x) = [k * * Sg(I)](t, x).$$

We put $G(s, y) = \mu k(s, y)$. Then, we have, for $s > 0$ and $y \in \mathbb{R}^n$,

$$G(s, y) \geq 0, \quad G(s, y) = G(s, -y) \quad \text{and} \quad \int_{0}^{+\infty} \int_{\mathbb{R}^n} G(s, y) dy ds = 1.$$

In fact, we have the following remark.

**Remark 1.** For $y = (y_1, y_2, \dots, y_n) \in \mathbb{R}^n$ and $s > 0$, $G$ satisfies

$$G(s, y) = \frac{\mu e^{-\mu s - |y|^2/(4D_E s)}}{(4\pi D_E s)^{n/2}} = \mu e^{-\mu s} \prod_{k=1}^{n} \frac{e^{-y_k^2/(4D_E s)}}{(4\pi D_E s)^{1/2}} = \mu e^{-\mu s} \prod_{k=1}^{n} \Gamma_E^1(s, y_k),$$

where

$$\Gamma_E^1(s, x) := \frac{e^{-x^2/(4D_E s)}}{(4\pi D_E s)^{1/2}} \quad \text{and} \quad \int_{-\infty}^{+\infty} \Gamma_E^1(s, x) dx = 1.$$

Substituting the obtained expression of $E$ into the third equation of System (4), we get, for $t > 0$ and $x \in \mathbb{R}^n$,

$$\left\{
\begin{aligned}
\frac{\partial S}{\partial t} &= D_S \Delta S - hS(t, x) - S(t, x)g(I(t, x)) + (1 - \eta) \int_{\mathbb{R}^n} \Gamma_p(\tau, x - y)u(t - \tau, y)dy, \\
\frac{\partial I}{\partial t} &= D_I \Delta I + [G * * Sg(I)](t, x) - rI(t, x), \\
u(t, x) &= hS(t, x) + \eta \int_{\mathbb{R}^n} \Gamma_p(\tau, x - y)u(t - \tau, y)dy.
\end{aligned}
\right.
\quad (6)$$

Our study is based on the system (6) and we aim to establish the existence and nonexistence of nontrivial traveling wave solutions.

### 3 Preliminary results on traveling wave solutions

In this section, we establish the basis for studying the existence of traveling wave solutions for the system (6). Let $S(0, x) = S_0$ and $u(0, x) = u_0$, $x \in \mathbb{R}^n$, be two constants representing susceptible and protected individuals at the beginning of the epidemic (before these individuals become infected). Note that $(S_0, 0, u_0)$ is a free steady state of the system (6) if and only if

$$S_0 = \frac{(1 - \eta)u_0}{h}.$$

A traveling wave solution is a positive solution defined as follows,

$$(S(t, x), I(t, x), u(t, x)) = (\phi_1(x.e + ct), \phi_2(x.e + ct), \phi_3(x.e + ct)), \quad t > 0 \text{ and } x \in \mathbb{R}^n,$$
where $\phi_i \in C^2(\mathbb{R}, \mathbb{R}^+)$, $i = \{1, 2, 3\}$, $C^k$ indicates that the $k$-th derivative of the function is continuous, and $c > 0$ describes the wave speed, see [29, 33] for a detailed description of such wave solutions. It is assumed that $e = (1, 0, ..., 0)$ which is the first unit vector of $\mathbb{R}^n$. Thus, $x.e = x_1 \in \mathbb{R}$. Setting $z = x_1 + ct$ and taking into account Remark 1, we arrive at the wave system, for $z \in \mathbb{R}$,

$$ egin{cases} c\phi_1'(z) = D_S\phi_1''(z) - h\phi_1(z) - \phi_1(z)g(\phi_2(z)) + (1-\eta)\int_{-\infty}^{+\infty} \Gamma_p^1(\tau, y)\phi_3(z-y-c\tau)dy, \\ c\phi_2'(z) = D_I\phi_2''(z) + [G_E^1 \star \phi_1 g(\phi_2)](z) - r\phi_2(z), \\ \phi_3(z) = h\phi_1(z) + \eta \int_{-\infty}^{+\infty} \Gamma_p^1(\tau, y)\phi_3(z-y-c\tau)dy, \end{cases} \quad (7) $$

with $G_E^1(s, y) = \mu e^{-\mu s} \Gamma_E^1(s, y)$ and the convolution is defined by

$$ (G_E^1 \star w)(z) = \int_0^{+\infty} \int_{-\infty}^{+\infty} G_E^1(s, y)w(z-y-cs) dy ds, \quad z \in \mathbb{R}. $$

Here, it is natural to look for bounded solutions such that $0 \le \phi_1(z) \le S_0$, $0 \le \phi_2(z)$ and $0 \le \phi_3(z) \le u_0$, for all $z \in \mathbb{R}$. Furthermore, we impose the following asymptotic boundary conditions,

$$ \begin{aligned} \phi_1(-\infty) &= S_0, & \phi_2(-\infty) &= 0, & \phi_3(-\infty) &= u_0, \\ \phi_1(+\infty) &= S_\infty < S_0, & \phi_2(+\infty) &= 0, & \phi_3(+\infty) &= u_\infty < u_0. \end{aligned} \quad (8) $$

In order to study the existence of solutions to the system (7)-(8), we start by uncoupling this problem from the third equation. Let $X$ be the Banach space of all bounded and uniformly continuous functions from $\mathbb{R}$ to $\mathbb{R}$ with the usual supremum norm $|\cdot|_X$. Define the following bounded linear operator $B$,

$$ (Bv)(z) = \eta \int_{-\infty}^{+\infty} \Gamma_p^1(\tau, y)v(z-y-c\tau)dy, \quad z \in \mathbb{R}, v \in X. $$

Since $\eta \in [0, 1)$, then $B$ is a contracting operator i.e. $|B|_X < 1$. This means that the operator $Id - B$ defined from $X$ to $X$ is invertible and its inverse is given by

$$ (Id - B)^{-1}(v) = \Psi * v, \quad \text{with} \quad \Psi(z) = \sum_{n=0}^{+\infty} \Psi_n(z), \quad (9) $$

$$ \Psi_n(z) = \eta^n \Gamma_p^1(n\tau, z - nc\tau) = \frac{\eta^n}{2\sqrt{nD_p\pi\tau}} \exp\left(-\frac{(z - nc\tau)^2}{4nD_p\tau}\right), \quad z \in \mathbb{R}, n \in \mathbb{N}, $$

and

$$ (\Psi * v)(z) = \int_{-\infty}^{+\infty} \Psi(z-y)v(y)dy. $$
In fact, we have

$$((Id - B)^{-1}v)(z) = \sum_{n=0}^{+\infty} (B^n v)(z), \\ = \int_{-\infty}^{+\infty} \sum_{n=0}^{+\infty} \frac{\eta^n}{2\sqrt{nD_p\pi\tau}} \exp\left(-\frac{(z-y-nc\tau)^2}{4nD_p\tau}\right) v(y)dy, \\ = \int_{-\infty}^{+\infty} \Psi(z-y)v(y)dy, \\ = (\Psi * v)(z).$$

Using the definition of the operator B, we can write the last equation of System (7) as follows,

$$\phi_3(z) = h\phi_1(z) + (B\phi_3)(z), \quad z \in \mathbb{R}.$$

This is equivalent to,

$$((Id - B)\phi_3)(z) = h\phi_1(z).$$

So, according to (9), we get

$$\phi_3(z) = h((Id - B)^{-1}\phi_1)(z) = h(\Psi * \phi_1)(z), \quad z \in \mathbb{R}. \qquad (10)$$

By replacing the expression (10) in the system (7), we obtain an equivalent system, $z \in \mathbb{R}$,

$$\left\{ \begin{aligned} c\phi_1'(z) &= D_S\phi_1''(z) - h\phi_1(z) - \phi_1(z)g(\phi_2(z)) + h(1-\eta)\int_{-\infty}^{+\infty} \Gamma_p^1(\tau, y)(\Psi * \phi_1)(z-y-c\tau)dy, \\ c\phi_2'(z) &= D_I\phi_2''(z) + [G_E^1 * \phi_1 g(\phi_2)](z) - r\phi_2(z). \end{aligned} \right. \qquad (11)$$

We set, for $z \in \mathbb{R}$,

$$\Upsilon_{\tau}(\Psi * \phi_1)(z) = \int_{-\infty}^{+\infty} \Gamma_p^1(\tau, y)(\Psi * \phi_1)(z-y)dy.$$

Then, we have

$$\Upsilon_{\tau}(\Psi * \phi_1) = \Gamma_p^1(\tau, \cdot) * (\Psi * \phi_1) = (\Gamma_p^1(\tau, \cdot) * \Psi) * \phi_1.$$

Let $\kappa_a, a \in \mathbb{R}$, be the translation operator defined on $X$ by $(\kappa_a f)(z) := f(z - a)$. As a consequence, we obtain

$$\Upsilon_{\tau}(\Psi * \phi_1)(z - c\tau) = (\kappa_{c\tau}\Upsilon_{\tau}(\Psi * \phi_1))(z) = [(\kappa_{c\tau}(\Gamma_p^1(\tau, \cdot) * \Psi)) * \phi_1](z). \qquad (12)$$

Finally, we consider $\psi: \mathbb{R} \to \mathbb{R}^+$ defined by $\psi = (1 - \eta)\kappa_{c\tau}(\Gamma_p^1(\tau, \cdot) * \Psi)$. Then, from (12), we get

$$\Upsilon_{\tau}(\Psi * \phi_1))(z - c\tau) = \frac{1}{1 - \eta}(\psi * \phi_1)(z), \quad z \in \mathbb{R}.$$
Clearly, we have

$$ \int_{-\infty}^{+\infty} \psi(x)dx = (1 - \eta) \int_{-\infty}^{+\infty} \Psi(x)dx = 1. $$

The system (11) becomes, for $z \in \mathbb{R}$,

$$ \left\{ \begin{array}{l} c\phi'_{1}(z) = D_{S}\phi''_{1}(z) - h\phi_{1}(z) - \phi_{1}(z)g(\phi_{2}(z)) + h(\psi * \phi_{1})(z), \\ c\phi'_{2}(z) = D_{I}\phi''_{2}(z) + [G_{E}^{1} * \phi_{1}g(\phi_{2})](z) - r\phi_{2}(z). \end{array} \right. \qquad (13) $$

In the sequel, we base our study on the wave problem (13) instead of (7). We add the conditions

$$ \begin{aligned} \phi_{1}(-\infty) &= S_{0}, & \phi_{2}(-\infty) &= 0, \\ \phi_{1}(+\infty) &= S_{\infty} < S_{0}, & \phi_{2}(+\infty) &= 0. \end{aligned} \qquad (14) $$

We need the following result.

**Proposition 1.** Let $\lambda^{-}(c)$ and $\lambda^{+}(c)$ be given, for $c > 0$, by

$$ \lambda^{-}(c) := \frac{c}{2D_{p}} \left( 1 - \sqrt{1 + \frac{4D_{p}}{\tau c^{2}} \ln\left(\frac{1}{\eta}\right)} \right) < 0, $$

$$ \lambda^{+}(c) := \frac{c}{2D_{p}} \left( 1 + \sqrt{1 + \frac{4D_{p}}{\tau c^{2}} \ln\left(\frac{1}{\eta}\right)} \right) > 0. $$

Then, for $\lambda \in (\lambda^{-}(c), \lambda^{+}(c))$, we have

$$ 1 - \eta e^{D_p \tau \lambda^2 - c \tau \lambda} > 0 \quad \text{and} \quad (\psi * e^{\lambda})(z) = \frac{(1 - \eta) e^{D_p \tau \lambda^2 + (z - c \tau) \lambda}}{1 - \eta e^{D_p \tau \lambda^2 - c \tau \lambda}}. $$

*Proof.* We use the following fundamental property, which can be obtained by a simple calculation,

$$ (\Gamma_{p}^{1}(s, \cdot) * e^{\lambda^{\cdot}})(z) = e^{D_{p}s\lambda^{2} + \lambda z}, \quad s > 0, \quad z \in \mathbb{R}. $$

Hence,

$$ (\psi * e^{\lambda^{\cdot}})(z) = (1 - \eta)e^{D_p \tau \lambda^2 + (z - c\tau)\lambda} \sum_{n=0}^{+\infty} (\eta e^{D_p \tau \lambda^2 - c\tau\lambda})^n. $$

By choosing $\lambda \in (\lambda^{-}(c), \lambda^{+}(c))$, we get $\eta e^{D_p \tau \lambda^2 - c \tau \lambda} < 1$. This proves the proposition. $\square$

To achieve our goal, we now need to define the transcendental characteristic function $\Delta_c$ for the linearization of the second equation of the reduced system (13) around the disease free point. That is

$$ \Delta_c(\lambda) := c\lambda - D_I\lambda^2 + r - S_0 g'(0) \int_0^{+\infty} \int_{-\infty}^{+\infty} G_E^1(s, y) e^{-cs\lambda - y\lambda} dy ds. \qquad (15) $$
Note that the integral on the right-hand side is not well defined for all $\lambda \in \mathbb{R}^+$. In fact, as stated before, we have

$$ \int_{\mathbb{R}} \Gamma_E^1(s, y)e^{-\lambda y}dy = e^{D_E s \lambda^2}, \quad s > 0. $$

Moreover, we have $\mu + c\lambda - D_E\lambda^2 > 0$, for $\lambda \in [0, \lambda_{\max}(c))$, with

$$ \lambda_{\max}(c) := \frac{c + \sqrt{c^2 + 4\mu D_E}}{2D_E}. $$

Then, for all $\lambda \in [0, \lambda_{\max}(c))$, (15) is well defined and becomes

$$ \Delta_c(\lambda) := c\lambda - D_I\lambda^2 + r - \frac{\mu S_0 g'(0)}{\mu + c\lambda - D_E\lambda^2} \quad \text{and} \quad \lim_{\lambda \to \lambda_{\max}(c)} \Delta_c(\lambda) = -\infty. \quad (16) $$

The following lemma gives some results on the nonnegative real roots of the function $\Delta_c$, see Figure 1. The proof is standard and therefore omitted.

**Lemma 1.** Assume that

$$ \frac{S_0 g'(0)}{r} > 1. \quad (17) $$

Then, there exists $c^* > 0$ and $0 < \lambda^* < \lambda_{\max}(c)$ such that,

$$ \frac{\partial}{\partial \lambda} \Delta_{c^*}(\lambda^*) = 0 \quad \text{and} \quad \Delta_{c^*}(\lambda^*) = 0. $$

In addition, we have the properties:

* If $c > c^*$, then $\Delta_c(\lambda) = 0$ has two real roots $\lambda_1(c)$ and $\lambda_2(c)$ such that $0 < \lambda_1(c) < \lambda_2(c) < \lambda_{\max}(c)$ and $\Delta_c(\lambda) > 0$ for $\lambda \in (\lambda_1(c), \lambda_2(c))$.

* If $0 < c < c^*$, then $\Delta_c(\lambda) < 0$, for all $\lambda \in (0, \lambda_{\max}(c))$.

**Remark 2.** The case where $D_I = D_E$ and $r = \mu$ can be discussed in more detail. In fact, the two roots $\lambda_1(c)$ and $\lambda_2(c)$ of $\Delta_c(\lambda) = 0$ are given, for $c > c^*$, by

$$ \lambda_1(c) = \frac{c}{2D_I} \left( 1 - \sqrt{1 - \frac{4rD_I}{c^2} \left( \sqrt{\frac{S_0 g'(0)}{r}} - 1 \right)} \right) > 0, $$

and

$$ \lambda_2(c) := \frac{c}{2D_I} \left( 1 + \sqrt{1 - \frac{4rD_I}{c^2} \left( \sqrt{\frac{S_0 g'(0)}{r}} - 1 \right)} \right) > 0. $$

When, $\lambda \in (\lambda_{\max}(c), +\infty)$, other roots are possible, but are not useful for our study. In this case, the minimal wave speed is given explicitly by

$$ c^* = 2\sqrt{rD_I \left( \sqrt{\frac{S_0 g'(0)}{r}} - 1 \right)}. $$
To achieve our study, it is necessary to ensure that $\lambda_1(c) < \lambda_{\max}(c)/2$. Note that by a direct verification, we always have $\lambda_1(c) < \lambda_{\max}(c)/2$ when $D_I = D_E$ and $r = \mu$. We can expect that even without these equalities but for close values, the condition $\lambda_1(c) < \lambda_{\max}(c)/2$ remains valid. For the general case, since

$$ \frac{\partial}{\partial c} \Delta_c(\lambda) > 0 \quad \text{and} \quad \lim_{c \to +\infty} \Delta_c(\lambda) = +\infty, $$

the graph of the concave function $\lambda \mapsto \Delta_c(\lambda)$ increases with $c$, see Figure 1. Consequently, there exists a speed $c$, denoted again $c^*$, such that $\lambda_1(c) < \lambda_{\max}(c)/2$. However, with this reasoning, it is not certain that this speed threshold expresses the minimum speed.

According to Proposition 1, Lemma 1 and the above discussion, we keep in mind the following thresholds,

$$ \lambda^+(c) := \frac{c}{2D_p} \left( 1 + \sqrt{1 + \frac{4D_p}{\tau c^2} \ln\left(\frac{1}{\eta}\right)} \right) \quad \text{and} \quad \tilde{\lambda}(c) := \frac{\lambda_{\max}(c)}{2} = \frac{c + \sqrt{c^2 + 4\mu D_E}}{4D_E}. $$

<figure><img src="image_2.png" /><figcaption>Figure 1: Some curves of the function $\Delta_c(\lambda)$, for several values of $c$, are plotted. The graph of the concave function $\lambda \mapsto \Delta_c(\lambda)$ increases as $c$ increases. Parameters are set arbitrarily: $D_I = 0.1$, $r = 0.02$, $\mu = 0.15$, $S_0 = 20$, $g'(0) = 0.5$ and $D_E = 0.1$.</figcaption></figure>

Through the rest of the paper, we assume that (17) is satisfied and we fix $c > c^*$. For the sake of simplicity, we often remove the dependency on $c$.

# 4 Existence of lower and upper solutions

In this section, motivated by the iteration method for monotone operators, we construct a lower and an upper solution that will be used to show the existence of traveling waves. We assume that

$$ \frac{S_0 g'(0)}{r} > 1 \quad \text{and} \quad c > c^*, $$
where the minimum wave speed $c^*$ is given by Lemma 1.

**Definition 2.** An upper solution $(\overline{\phi_1}, \overline{\phi_2}) \in C(\mathbb{R}, \mathbb{R}^2)$ and a lower solution $(\underline{\phi_1}, \underline{\phi_2}) \in C(\mathbb{R}, \mathbb{R}^2)$ of (13) are two pairs satisfying the following conditions,

$$ egin{cases} c\overline{\phi_1}'(z) \geq D_S\overline{\phi_1}''(z) - h\overline{\phi_1}(z) - \overline{\phi_1}(z)g(\underline{\phi_2}(z)) + h(\psi * \overline{\phi_1})(z), & \text{a.e. in } \mathbb{R}, \\ c\underline{\phi_2}'(z) \geq D_I\underline{\phi_2}''(z) + [G_E^1 * \overline{\phi_1}g(\underline{\phi_2})](z) - r\underline{\phi_2}(z), & \text{a.e. in } \mathbb{R}, \end{cases} $$

and the reversed inequality

$$ egin{cases} c\underline{\phi_1}'(z) \leq D_S\underline{\phi_1}''(z) - h\underline{\phi_1}(z) - \underline{\phi_1}(z)g(\overline{\phi_2}(z)) + h(\psi * \underline{\phi_1})(z), & \text{a.e. in } \mathbb{R}, \\ c\overline{\phi_2}'(z) \leq D_I\overline{\phi_2}''(z) + [G_E^1 * \underline{\phi_1}g(\overline{\phi_2})](z) - r\overline{\phi_2}(z), & \text{a.e. in } \mathbb{R}, \end{cases} $$

where the first and the second derivative of the functions $\underline{\phi_1}, \underline{\phi_2}, \overline{\phi_1}$ and $\overline{\phi_2}$ exist almost everywhere (a.e.).

In the first inequality of Definition 2, we use $g(\underline{\phi_2})$ instead of $g(\overline{\phi_2})$; this is motivated by the fact that the nonlinearities in our model are mixed quasimonotone (see, Definition 1.2, Chapter 8 of the book by [24]). Based on the roots of the characteristic function (15) (or (16)), we can show the following lemma.

**Lemma 3.** The pair of functions $(\overline{\phi_1}, \overline{\phi_2})$, given by,

$$ (\overline{\phi_1}, \overline{\phi_2}) : \mathbb{R} \to (\mathbb{R}^+)^2, \\ z \mapsto (\underline{\phi_1}(z), \underline{\phi_2}(z)) = (S_0, e^{\lambda_1 z}), $$

is an upper solution of (13). Furthermore, if we choose $\epsilon_1 \in (0, \min\{\lambda^+, \lambda_1\})$, $\epsilon_2 \in (0, \min\{\epsilon_1, \lambda_2 - \lambda_1\})$, $\mu_1 > S_0$ and a constant $\mu_2 > 1$ which will be clarified later, then the pair of functions

$$ (\underline{\phi_1}, \underline{\phi_2}) : \mathbb{R} \to (\mathbb{R}^+)^2, \\ z \mapsto (\underline{\phi_1}(z), \underline{\phi_2}(z)) = (\max\{0, S_0 - \mu_1 e^{\epsilon_1 z}\}, \max\{0, e^{\lambda_1 z} - \mu_2 e^{(\lambda_1 + \epsilon_2)z}\}), $$

is a lower solution of (13). In particular, we have

$$ \underline{\phi_1}(z) \leq \overline{\phi_1}(z) \quad \text{and} \quad \underline{\phi_2}(z) \leq \overline{\phi_2}(z), \quad z \in \mathbb{R}. $$

**Proof.** We start to show that $(\overline{\phi_1}, \overline{\phi_2})$ is an upper solution. We have, for $z \in \mathbb{R}$,

$$ \overline{\phi_1}'(z) = \overline{\phi_1}''(z) = 0. $$

Given the above calculation and the fact that $\underline{\phi_2}$ is nonnegative, we get

$$ c\overline{\phi_1}'(z) - D_S\overline{\phi_1}''(z) + h\overline{\phi_1}(z) + \overline{\phi_1}(z)g(\underline{\phi_2}(z)) - h(\psi * \overline{\phi_1})(z) = S_0g(\underline{\phi_2}(z)) \geq 0. $$

Moreover, as $g(y) \leq g'(0)y$, $y \in \mathbb{R}^+$, we obtain

$$ c\overline{\phi_2}'(z) - D_I\overline{\phi_2}''(z) - [G_E^1 * \overline{\phi_1}g(\overline{\phi_2})](z) + r\overline{\phi_2}(z) \\ \geq c\lambda_1 e^{\lambda_1 z} - D_I\lambda_1^2 e^{\lambda_1 z} - S_0 g'(0)[G_E^1 * \overline{\phi_2}](z) + re^{\lambda_1 z}, \\ = e^{\lambda_1 z} \left( c\lambda_1 - D_I\lambda_1^2 + r - S_0 g'(0) \int_0^{+\infty} \int_{-\infty}^{+\infty} G_E^1(s, y)e^{-cs\lambda_1 - y\lambda_1} dy ds \right), \\ = e^{\lambda_1 z} \Delta_c(\lambda_1) = 0. $$
Therefore, from the previous inequalities, we proved that $(\overline{\phi_1}, \overline{\phi_2})$ is an upper solution of (13).

Next, we consider the pair $(\underline{\phi_1}, \underline{\phi_2})$. Recall that we have to show that, for $z \in \mathbb{R}$, the have the following inequality,

$$c\underline{\phi_1}'(z) \leq D_S\underline{\phi_1}''(z) - h\underline{\phi_1}(z) - \underline{\phi_1}(z)g(\overline{\phi_2}(z)) + h(\psi * \underline{\phi_1})(z).$$

Note that

$$\underline{\phi_1}(z) = \begin{cases} 0, & z \geq z_0, \\ S_0 - \mu_1 e^{\epsilon_1 z}, & z < z_0, \end{cases}$$

where $z_0 = \frac{1}{\epsilon_1} \ln \left( \frac{S_0}{\mu_1} \right) < 0$. On the one hand, for $z \geq z_0$, we have

$$c\underline{\phi_1}'(z) - D_S\underline{\phi_1}''(z) + h\underline{\phi_1}(z) + \underline{\phi_1}(z)g(\overline{\phi_2}(z)) - h(\psi * \underline{\phi_1})(z) = -h \int_{-\infty}^{+\infty} \psi(z-y)\underline{\phi_1}(y)dy \leq 0.$$

By using Proposition (1) and since $g(y) \leq g'(0)y$, $y \in \mathbb{R}^+$, we have, for $z < z_0$,

$$\begin{aligned} c\underline{\phi_1}'(z) - D_S\underline{\phi_1}''(z) + h\underline{\phi_1}(z) + \underline{\phi_1}(z)g(\overline{\phi_2}(z)) - h(\psi * \underline{\phi_1})(z) & \\ & \leq -c\mu_1\epsilon_1 e^{\epsilon_1 z} + D_S\mu_1\epsilon_1^2 e^{\epsilon_1 z} + h(S_0 - \mu_1 e^{\epsilon_1 z}) + g'(0)(S_0 - \mu_1 e^{\epsilon_1 z})\overline{\phi_2}(z) - h(\psi * \underline{\phi_1})(z), \\ & \leq -c\mu_1\epsilon_1 e^{\epsilon_1 z} + D_S\mu_1\epsilon_1^2 e^{\epsilon_1 z} - h\mu_1 e^{\epsilon_1 z} + g'(0)e^{\lambda_1 z}(S_0 - \mu_1 e^{\epsilon_1 z}) + \mu_1 h \frac{(1-\eta)e^{D_p\tau\epsilon_1^2+(z-c\tau)\epsilon_1}}{1-\eta e^{D_p\tau\epsilon_1^2-c\tau\epsilon_1}}, \\ & \leq e^{\epsilon_1 z} \left[ -c\mu_1\epsilon_1 + D_S\mu_1\epsilon_1^2 - h\mu_1 + g'(0)S_0 e^{(\lambda_1-\epsilon_1)z} + \mu_1 h \frac{(1-\eta)e^{D_p\tau\epsilon_1^2-c\tau\epsilon_1}}{1-\eta e^{D_p\tau\epsilon_1^2-c\tau\epsilon_1}} \right], \\ & \leq e^{\epsilon_1 z} \left[ -c\mu_1\epsilon_1 + D_S\mu_1\epsilon_1^2 + g'(0)S_0 \left( \frac{S_0}{\mu_1} \right)^{(\lambda_1-\epsilon_1)/\epsilon_1} + \mu_1 h \frac{-1+e^{D_p\tau\epsilon_1^2-c\tau\epsilon_1}}{1-\eta e^{D_p\tau\epsilon_1^2-c\tau\epsilon_1}} \right]. \end{aligned}$$

The purpose now is to choose $\epsilon_1$ and $\mu_1 > S_0$ such that

$$-c\mu_1\epsilon_1 + D_S\mu_1\epsilon_1^2 + g'(0)S_0 \left(\frac{S_0}{\mu_1}\right)^{(\lambda_1-\epsilon_1)/\epsilon_1} + \mu_1 h \frac{-1 + e^{D_p\tau\epsilon_1^2-c\tau\epsilon_1}}{1 - \eta e^{D_p\tau\epsilon_1^2-c\tau\epsilon_1}} < 0.$$

In fact, for $\epsilon_1 = 1/\mu_1$ and $\mu_1$ sufficiently large, we have

$$-c + \epsilon_1 D_S + g'(0)S_0(\epsilon_1 S_0)^{(\lambda_1 - \epsilon_1)/\epsilon_1} + h \frac{-1 + e^{D_p \tau \epsilon_1^2 - c\tau \epsilon_1}}{\epsilon_1(1 - \eta e^{D_p \tau \epsilon_1^2 - c\tau \epsilon_1})} < 0.$$

For example, we can choose

$$\epsilon_1 < \min \left\{ \frac{c^*}{D_p}, \frac{1}{S_0} \right\}.$$

Therefore, for $z < z_0$, we arrive to

$$c\underline{\phi_1}'(z) \leq D_S\underline{\phi_1}''(z) - h\underline{\phi_1}(z) - \underline{\phi_1}(z)g(\overline{\phi_2}(z)) + h(\psi * \underline{\phi_1})(z).$$

Now we have to prove that, for $z \in \mathbb{R}$,

$$c\underline{\phi_2}'(z) \leq D_I\underline{\phi_2}''(z) + [G_E^1 \star \underline{\phi_1}g(\underline{\phi_2})](z) - r\underline{\phi_2}(z). \quad (18)$$
Recall that

$$ \underline{\phi}_2(z) = \begin{cases} 0, & z \ge z_1, \\ e^{\lambda_1 z} - \mu_2 e^{(\lambda_1 + \epsilon_2)z}, & z < z_1, \end{cases} $$

where $z_1 = -(1/\epsilon_2)\ln(\mu_2) < 0$. When $z \ge z_1$ it is obvious that (18) is verified. Indeed, $\underline{\phi}_2(z) = \underline{\phi}_2'(z) = \underline{\phi}_2''(z) = 0$ and $\underline{\phi}_1, \underline{\phi}_2$ are nonnegative. So, we have

$$ -[G_E^1 \star \underline{\phi}_1 g(\underline{\phi}_2)](z) = - \int_0^{+\infty} \int_{-\infty}^{+\infty} G_E^1(s, y) \underline{\phi}_1(z - y - cs) g(\underline{\phi}_2(z - y - cs)) dyds < 0. $$

Consequently, for $z \ge z_1$,

$$ c\underline{\phi}_2'(z) \le D_I\underline{\phi}_2''(z) + [G_E^1 \star \underline{\phi}_1 g(\underline{\phi}_2)](z) - r\underline{\phi}_2(z). $$

Now we suppose that $z < z_1$. Since

$$ g(\underline{\phi}_2(z)) \ge g'(0)\underline{\phi}_2(z) - ag'(0)\underline{\phi}_2^2(z), $$

and

$$ \underline{\phi}_2(z) \ge e^{\lambda_1 z} - \mu_2 e^{(\lambda_1 + \epsilon_2)z}, \quad \underline{\phi}_2^2(z) \le e^{2\lambda_1 z}, $$

we obtain the expression

$$ \begin{aligned} & c\underline{\phi}_2'(z) - D_I\underline{\phi}_2''(z) - [G_E^1 \star \underline{\phi}_1 g(\underline{\phi}_2)](z) + r\underline{\phi}_2(z) \\ & \le c\underline{\phi}_2'(z) - D_I\underline{\phi}_2''(z) + r\underline{\phi}_2(z) \\ & \quad - \int_0^{+\infty} \int_{-\infty}^{+\infty} G_E^1(s, y) (S_0 - \mu_1 e^{\epsilon_1(z-y-cs)}) g(\underline{\phi}_2(z-y-cs)) dyds, \\ & \le c\underline{\phi}_2'(z) - D_I\underline{\phi}_2''(z) + r\underline{\phi}_2(z) \\ & \quad + \int_0^{+\infty} \int_{-\infty}^{+\infty} G_E^1(s, y) (S_0 - \mu_1 e^{\epsilon_1(z-y-cs)}) (-g'(0)\underline{\phi}_2(z-y-cs) + ag'(0)e^{2\lambda_1(z-y-cs)}) dyds, \\ & = c\lambda_1 e^{\lambda_1 z} - \mu_2 c(\lambda_1 + \epsilon_2)e^{(\lambda_1 + \epsilon_2)z} - D_I\lambda_1^2 e^{\lambda_1 z} + D_I\mu_2(\lambda_1 + \epsilon_2)^2 e^{(\lambda_1 + \epsilon_2)z} + re^{\lambda_1 z} - \mu_2 re^{(\lambda_1 + \epsilon_2)z} \\ & \quad + \int_0^{+\infty} \int_{-\infty}^{+\infty} G_E^1(s, y) (S_0 - \mu_1 e^{\epsilon_1(z-y-cs)}) \\ & \qquad \times (-g'(0)e^{\lambda_1(z-y-cs)} + \mu_2 g'(0)e^{(\lambda_1 + \epsilon_2)(z-y-cs)} + ag'(0)e^{2\lambda_1(z-y-cs)}) dyds. \end{aligned} $$
This leads to

$$
\begin{aligned}
& c\underline{\phi_2}'(z) - D_I\underline{\phi_2}''(z) - [G_E^1 \star \underline{\phi_1}g(\underline{\phi_2})](z) + r\underline{\phi_2}(z) \\
& \le e^{\lambda_1 z} \Delta(\lambda_1) - \mu_2 c(\lambda_1 + \epsilon_2) e^{(\lambda_1 + \epsilon_2)z} + D_I \mu_2 (\lambda_1 + \epsilon_2)^2 e^{(\lambda_1 + \epsilon_2)z} - \mu_2 r e^{(\lambda_1 + \epsilon_2)z} \\
& \quad + ag'(0) \int_0^{+\infty} \int_{-\infty}^{+\infty} G_E^1(s, y) (S_0 - \mu_1 e^{\epsilon_1(z-y-cs)}) e^{2\lambda_1(z-y-cs)} dyds \\
& \quad + \int_0^{+\infty} \int_{-\infty}^{+\infty} G_E^1(s, y) [S_0 \mu_2 g'(0) e^{(\lambda_1+\epsilon_2)(z-y-cs)} + \mu_1 g'(0) e^{(\lambda_1+\epsilon_1)(z-y-cs)}] dyds \\
& \quad + \int_0^{+\infty} \int_{-\infty}^{+\infty} G_E^1(s, y) [-\mu_1 \mu_2 g'(0) e^{(\lambda_1+\epsilon_1+\epsilon_2)(z-y-cs)}] dyds.
\end{aligned}
$$

Then, we have

$$
\begin{aligned}
& c\underline{\phi_2}'(z) - D_I\underline{\phi_2}''(z) - [G_E^1 \star \underline{\phi_1}g(\underline{\phi_2})](z) + r\underline{\phi_2}(z) \\
& \le -\mu_2 e^{(\lambda_1+\epsilon_2)z} \Delta(\lambda_1 + \epsilon_2) + \int_0^{+\infty} \int_{-\infty}^{+\infty} G_E^1(s, y) (ag'(0)S_0 e^{2\lambda_1(z-y-cs)}) dyds \\
& \quad + \int_0^{+\infty} \int_{-\infty}^{+\infty} G_E^1(s, y) [\mu_1 g'(0) e^{(\lambda_1+\epsilon_1)(z-y-cs)} - \mu_1 \mu_2 g'(0) e^{(\lambda_1+\epsilon_1+\epsilon_2)(z-y-cs)}] dyds, \\
& = e^{(\lambda_1+\epsilon_2)z} \left[ -\mu_2 \Delta(\lambda_1 + \epsilon_2) + \int_0^{+\infty} \int_{-\infty}^{+\infty} G_E^1(s, y) ag'(0) S_0 e^{(\lambda_1-\epsilon_2)z + 2\lambda_1(-y-cs)} dyds \right] \\
& \quad + e^{(\lambda_1+\epsilon_2)z} \left[ \int_0^{+\infty} \int_{-\infty}^{+\infty} G_E^1(s, y) (\mu_1 g'(0) e^{(\lambda_1+\epsilon_1)(-y-cs)} + (\epsilon_1-\epsilon_2)z) dyds \right] \\
& \quad + e^{(\lambda_1+\epsilon_2)z} \left[ \int_0^{+\infty} \int_{-\infty}^{+\infty} G_E^1(s, y) (-\mu_1 \mu_2 g'(0) e^{\epsilon_1(z-y-cs)} + (\lambda_1+\epsilon_2)(-y-cs)) dyds \right].
\end{aligned}
$$

The double integral

$$
\int_0^{+\infty} \int_{-\infty}^{+\infty} G_E^1(s, y) e^{2\lambda_1(-y-cs)} dyds = \frac{\mu}{\mu + 2c\lambda_1 - 4D_E\lambda_1^2},
$$
is well defined (see Remark 2). Since $\epsilon_2$ is chosen such that $\epsilon_2 < \epsilon_1 < \lambda_1$, then we have

$$c\underline{\phi}_2'(z) - D_I\underline{\phi}_2''(z) - [G_E^1 \star \underline{\phi}_1 g(\underline{\phi}_2)](z) + r\underline{\phi}_2(z) \\ \le e^{(\lambda_1+\epsilon_2)z} \left[ -\mu_2\Delta(\lambda_1+\epsilon_2) + \int_0^{+\infty} \int_{-\infty}^{+\infty} G_E^1(s,y)ag'(0)S_0 e^{2\lambda_1(-y-cs)}dyds \right. \\ \left. + \int_0^{+\infty} \int_{-\infty}^{+\infty} G_E^1(s,y)(\mu_1 g'(0)e^{(\lambda_1+\epsilon_1)(-y-cs)})dyds \right].$$

Moreover, we are considering $\epsilon_2 < \lambda_2 - \lambda_1$. Hence, $\Delta(\lambda_1 + \epsilon_2) > 0$. Accordingly, if we choose

$$\mu_2 > \frac{1}{\Delta(\lambda_1 + \epsilon_2)} \left( \int_0^{+\infty} \int_{-\infty}^{+\infty} G_E^1(s,y)ag'(0)S_0 e^{2\lambda_1(-y-cs)}dyds \right. \\ \left. + \int_0^{+\infty} \int_{-\infty}^{+\infty} G_E^1(s,y)(\mu_1 g'(0)e^{(\lambda_1+\epsilon_1)(-y-cs)})dyds \right),$$

we arrive to

$$-\mu_2\Delta(\lambda_1 + \epsilon_2) + \int_0^{+\infty} \int_{-\infty}^{+\infty} G_E^1(s,y)ag'(0)S_0 e^{2\lambda_1(-y-cs)}dyds \\ + \int_0^{+\infty} \int_{-\infty}^{+\infty} G_E^1(s,y)(\mu_1 g'(0)e^{(\lambda_1+\epsilon_1)(-y-cs)})dyds < 0.$$

Consequently, for $z < z_1$, we have proven that

$$c\underline{\phi}_2'(z) \le D_I\underline{\phi}_2''(z) + [G_E^1 \star \underline{\phi}_1 g(\underline{\phi}_2)](z) - r\underline{\phi}_2(z).$$

□

# 5 Existence and non-existence of traveling wave solutions

In this section, we focus on the study of the existence and non-existence of solutions to the problem (13)-(14). The classical approach is to formulate the system as a fixed-point problem, in order to obtain appropriate solutions. This guarantees the existence of traveling wave solutions. First, we define the profile set of traveling wave solutions,

$$X = \{(\phi_1, \phi_2) \in C(\mathbb{R}, \mathbb{R}^2) / \forall z \in \mathbb{R}: \underline{\phi}_1(z) \le \phi_1(z) \le \overline{\phi}_1(z), \underline{\phi}_2(z) \le \phi_2(z) \le \overline{\phi}_2(z)\},$$

where $(\underline{\phi}_1, \underline{\phi}_2)$ and $(\overline{\phi}_1, \overline{\phi}_2)$ are the upper and lower solutions obtained in Lemma 3. Note that $X$ is non-empty since $(\underline{\phi}_1, \underline{\phi}_2)$, $(\overline{\phi}_1, \overline{\phi}_2) \in X$. In the sequel, we consider the operators $H_1, H_2$
defined by

$$H_1: \quad \mathbb{X} \to C(\mathbb{R}, \mathbb{R})$$

$$(\phi_1, \phi_2) \mapsto H_1(\phi_1, \phi_2)(z) = A_1\phi_1(z) - \frac{1}{g'(0)}\phi_1(z)g(\phi_2(z)) + \frac{h}{g'(0)}(\psi * \phi_1)(z),$$

and

$$H_2: \quad \mathbb{X} \to C(\mathbb{R}, \mathbb{R})$$

$$(\phi_1, \phi_2) \mapsto H_2(\phi_1, \phi_2)(z) = \frac{1}{g'(0)}[G_E^1 \star \phi_1 g(\phi_2)](z) + \frac{A_2 - r}{g'(0)}\phi_2(z).$$

The constants $A_1$ and $A_2$ are chosen large enough to satisfy

$$A_1 > \max \left\{ S_0, \sup_{y \ge 0} \left[ \frac{g(y)}{g'(0)} \right] \right\} \quad \text{and} \quad A_2 > r. \qquad (19)$$

Now, we define the operator

$$(F_1, F_2) : \mathbb{X} \to C(\mathbb{R}, \mathbb{R}^2),$$

by

$$F_1(\phi_1, \phi_2)(z) = \frac{g'(0)}{\delta_1} \int_{-\infty}^{z} e^{\gamma_{11}(z-x)} H_1(\phi_1, \phi_2)(x)dx + \frac{g'(0)}{\delta_1} \int_{z}^{+\infty} e^{\gamma_{12}(z-x)} H_1(\phi_1, \phi_2)(x)dx,$$

and

$$F_2(\phi_1, \phi_2)(z) = \frac{g'(0)}{\delta_2} \int_{-\infty}^{z} e^{\gamma_{21}(z-x)} H_2(\phi_1, \phi_2)(x)dx + \frac{g'(0)}{\delta_2} \int_{z}^{+\infty} e^{\gamma_{22}(z-x)} H_2(\phi_1, \phi_2)(x)dx,$$

where

$$\gamma_{11} := \frac{c - \sqrt{c^2 + 4D_S(g'(0)A_1 + h)}}{2D_S} < 0, \quad \gamma_{12} := \frac{c + \sqrt{c^2 + 4D_S(g'(0)A_1 + h)}}{2D_S} > 0,$$

$$\gamma_{21} := \frac{c - \sqrt{c^2 + 4D_I A_2}}{2D_I} < 0, \quad \gamma_{22} := \frac{c + \sqrt{c^2 + 4D_I A_2}}{2D_I} > 0,$$

and

$$\delta_1 = D_S(\gamma_{12} - \gamma_{11}) > 0 \quad \text{and} \quad \delta_2 = D_I(\gamma_{22} - \gamma_{21}) > 0.$$

We need to have $\lambda_1 < -\gamma_{11}$ and $\lambda_1 < -\gamma_{21}$. This is guaranteed by the choice of $A_1$ and $A_2$ large enough. Note that, the operator $(F_1, F_2)$ is well defined. We can check that a fixed point of this operator $\phi_1 = F_1(\phi_1, \phi_2)$ and $\phi_2 = F_2(\phi_1, \phi_2)$ is a solution of the wave system (13). In the next results, we give some properties about the wave solutions of (13).

**Lemma 4.** Assume that $(\phi_1, \phi_2)$ is a solution of System (13)-(14). Then, we have

$$\int_{-\infty}^{+\infty} \phi_1(s)g(\phi_2(s))ds < +\infty, \quad \int_{-\infty}^{+\infty} [G_E^1 \star \phi_1 g(\phi_2)](z)dz < +\infty,$$

and

$$\int_{-\infty}^{+\infty} \phi_2(z)dz < +\infty.$$
Proof. Suppose that a traveling wave solution exists. So, the problem (13)-(14) holds for some $(\phi_1, \phi_2)$. We arrive to the following expression, for $z \in \mathbb{R}$,

$$ \phi_2(z) = \frac{1}{\tilde{\delta}_2} \int_{-\infty}^{z} e^{\tilde{\gamma}_{21}(z-x)} [G_E^1 \star \phi_1 g(\phi_2)](x) dx + \frac{1}{\tilde{\delta}_2} \int_{z}^{+\infty} e^{\tilde{\gamma}_{22}(z-x)} [G_E^1 \star \phi_1 g(\phi_2)](x) dx, $$

where

$$ \tilde{\delta}_2 = D_I(\tilde{\gamma}_{22} - \tilde{\gamma}_{21}), \quad \tilde{\gamma}_{21} = \frac{c}{2D_I} \left( 1 - \sqrt{1 + \frac{4D_I r}{c^2}} \right) \quad \text{and} \quad \tilde{\gamma}_{22} = \frac{c}{2D_I} \left( 1 + \sqrt{1 + \frac{4D_I r}{c^2}} \right). $$

By doing a change of variables, we get

$$ \phi_2(z) = \frac{1}{\tilde{\delta}_2} \int_{0}^{+\infty} e^{\tilde{\gamma}_{21}x} [G_E^1 \star \phi_1 g(\phi_2)](z-x) dx + \frac{1}{\tilde{\delta}_2} \int_{-\infty}^{0} e^{\tilde{\gamma}_{22}x} [G_E^1 \star \phi_1 g(\phi_2)](z-x) dx. \quad (20) $$

It is not difficult to see that we can write $\phi_1$ as

$$ \phi_1(z) = \frac{g'(0)}{\delta_1} \int_{-\infty}^{z} e^{\gamma_{11}(z-x)} H_1(\phi_1, \phi_2)(x) dx + \frac{g'(0)}{\delta_1} \int_{z}^{+\infty} e^{\gamma_{12}(z-x)} H_1(\phi_1, \phi_2)(x) dx. $$

In addition, we have

$$ \phi_1'(z) = \gamma_{11} \frac{g'(0)}{\delta_1} \int_{-\infty}^{z} e^{\gamma_{11}(z-x)} H_1(\phi_1, \phi_2)(x) dx + \gamma_{12} \frac{g'(0)}{\delta_1} \int_{z}^{+\infty} e^{\gamma_{12}(z-x)} H_1(\phi_1, \phi_2)(x) dx. $$

By doing a change of variables, we obtain

$$ \phi_1'(z) = \gamma_{11} \frac{g'(0)}{\delta_1} \int_{0}^{+\infty} e^{\gamma_{11}\omega} H_1(\phi_1, \phi_2)(z-\omega) d\omega + \gamma_{12} \frac{g'(0)}{\delta_1} \int_{-\infty}^{0} e^{\gamma_{12}\omega} H_1(\phi_1, \phi_2)(z-\omega) d\omega. \quad (21) $$

Since $\lim_{z \to -\infty} \phi_1(z) = S_0$ and $\lim_{z \to -\infty} \phi_2(z) = 0$, we get, for any $\omega > 0$,

$$ \begin{aligned} & \lim_{z \to -\infty} H_1(\phi_1, \phi_2)(z - \omega) \\ &= \lim_{z \to -\infty} \left[ A_1 \phi_1(z - \omega) - \frac{1}{g'(0)} \phi_1(z - \omega) g(\phi_2(z - \omega)) + \frac{h}{g'(0)} (\psi * \phi_1)(z - \omega) \right], \\ &= \left( A_1 + \frac{h}{g'(0)} \right) S_0. \end{aligned} $$
Recall that $\gamma_{11} < 0$ and $\gamma_{12} > 0$. Therefore, (21) implies that

$$
\begin{align*}
\lim_{z \to -\infty} \phi_1'(z) &= \gamma_{11} \frac{g'(0)}{\delta_1} \int_0^{+\infty} e^{\gamma_{11}\omega} \lim_{z \to -\infty} H_1(\phi_1, \phi_2)(z-\omega) d\omega \\
&\quad + \gamma_{12} \frac{g'(0)}{\delta_1} \int_{-\infty}^0 e^{\gamma_{12}\omega} \lim_{z \to -\infty} H_1(\phi_1, \phi_2)(z-\omega) d\omega, \\
&= \left( A_1 + \frac{h}{g'(0)} \right) S_0 \frac{g'(0)}{\delta_1} [e^{\gamma_{11}\omega}|_0^{+\infty} + e^{\gamma_{12}\omega}|_{-\infty}^0], \\
&= 0.
\end{align*}
$$

We integrate the first equation of (13) from $-\infty$ to $z$, to obtain

$$
c[\phi_1(z) - S_0] = D_S \phi_1'(z) - \int_{-\infty}^{z} \phi_1(s)g(\phi_2(s))ds + h \int_{-\infty}^{z} (\psi * \phi_1)(y) - \phi_1(y)dy. \quad (22)
$$

Then, we used the fact that $\lim_{z \to -\infty} \phi_1'(z) = 0$. We set

$$
J(z) = h \int_{-\infty}^{z} (\psi * \phi_{1})(y) - \phi_{1}(y) dy.
$$

Despite that $\phi_1$ is not integrable on $(-\infty, r]$, the function $J$ is well defined. In fact, we set $t < z$
to obtain

$$
\int_t^z [(\psi * \phi_1)(y) - \phi_1(y)] dy = \int_t^z \left[ \int_{-\infty}^{+\infty} \psi(s)\phi_1(y-s)ds - \phi_1(y) \right] dy.
$$

Since

$$
\int_{-\infty}^{+\infty} \psi(s) ds = 1.
$$

Then, we get

$$
\int_t^z \left[ \int_{-\infty}^{+\infty} \psi(s) \phi_1(y-s) ds - \phi_1(y) \right] dy = \int_t^z \int_{-\infty}^{+\infty} \psi(s) [\phi_1(y-s) - \phi_1(y)] ds dy.
$$

On another hand, we have

$$
\phi_1(y-s) - \phi_1(y) = - \int_0^1 s \phi'_1(y - \eta s) d\eta.
$$

So, we obtain

$$
\int_t^z [(\psi * \phi_1)(y) - \phi_1(y)] dy = - \int_t^z \int_{-\infty}^{+\infty} s \psi(s) \int_0^1 \phi'_1(y - \eta s) d\eta ds dy.
$$
According to Fubini's theorem, we get

$$ \int_t^z [(\psi * \phi_1)(y) - \phi_1(y)] dy = - \int_{-\infty}^{+\infty} s\psi(s) \int_0^1 \int_t^z \phi'_1(y - \eta s) dy d\eta ds. $$

Therefore, the dominated convergence theorem implies that

$$ \int_{-\infty}^{z} [(\psi * \phi_{1})(y) - \phi_{1}(y)] dy = - \int_{-\infty}^{+\infty} s\psi(s) \int_{0}^{1} \left[ \lim_{t \to -\infty} \int_{t}^{z} \phi'_{1}(y - \eta s) dy \right] d\eta ds. $$

Recall that $\lim_{x \to -\infty} \phi_1(x) = S_0$. Thus, we have

$$ \lim_{t \to -\infty} \int_t^z \phi'_1(y - \eta s) dy = \phi_1(z - \eta s) - S_0. $$

As a result, we obtain

$$ \int_{-\infty}^{z} [(\psi * \phi_{1})(y) - \phi_{1}(y)] dy = - \int_{-\infty}^{+\infty} y\psi(y) \int_{0}^{1} \phi_{1}(z - \eta y) - S_{0}d\eta dy. $$

We know that the function $\phi_1$ is bounded on $\mathbb{R}$. Then, the function

$$ (s, y) \in \mathbb{R} \times \mathbb{R} \mapsto \int_{0}^{1} \phi_{1}(y - \eta s) - S_{0}d\eta \text{ is bounded.} $$

Moreover, we have

$$ \int_{-\infty}^{+\infty} |y\psi(y)| dy < +\infty. $$

Consequently, we get

$$ \int_{-\infty}^{z} |(\psi * \phi_{1})(y) - \phi_{1}(y)| dy < +\infty, \quad z \in \mathbb{R}. $$

If $\int_{-\infty}^{+\infty} \phi_1(s)g(\phi_2(s))ds = +\infty$, then from (22), we obtain

$$ c[\phi_1(z) - S_0] = D_S \phi'_1(z) - \int_{-\infty}^{z} \phi_1(s)g(\phi_2(s))ds + J(z), \quad z \in \mathbb{R}. $$

This implies that

$$ \lim_{z \to +\infty} \phi'_{1}(z) = +\infty. $$

We obtain a contradiction, since for all $z \in \mathbb{R}$ we have $\phi_1(z) \le S_0$. We finally conclude that

$$ \int_{-\infty}^{+\infty} \phi_1(s)g(\phi_2(s))ds < +\infty. \qquad (23) $$
On the other hand, from (23) and since $\int_0^{+\infty} \int_{-\infty}^{+\infty} G_E^1(s, y) ds dy = 1$, we conclude that

$$ \int_{-\infty}^{+\infty} [G_E^1 \star \phi_1 g(\phi_2)](x) dx = \int_0^{+\infty} \int_{-\infty}^{+\infty} G_E^1(s, y) \left[ \int_{-\infty}^{+\infty} (\phi_1 g(\phi_2))(x - y - cs) dx \right] dy ds < +\infty. $$

Moreover, we can show from (20) that $\lim_{z \to -\infty} \phi_2'(z) = 0$ and $\lim_{z \to +\infty} \phi_2'(z) = 0$. Therefore, by considering the second equation of (13) integrated from $-\infty$ to $+\infty$, we get

$$ \int_{-\infty}^{+\infty} \phi_2(z) dz = \frac{1}{r} \int_{-\infty}^{+\infty} [G_E^1 \star \phi_1 g(\phi_2)](z) dz < +\infty. $$

□

In the following, we show that the condition (17) is necessary to ensure the existence of traveling wave solutions for the system (6) (or (13)).

**Theorem 5.** Suppose that the condition (17) is not satisfied. This means,

$$ \frac{S_0 g'(0)}{r} < 1, \quad \left( \text{or} \quad \frac{(1-\eta)g'(0)u_0}{rh} < 1 \right). $$

Then, for any $c > 0$, there is no traveling wave solution $(\phi_1, \phi_2)$ for the system (6)-(8), such that $0 \le \phi_1(z) \le S_0$ and $0 \le \phi_3(z) \le u_0$ for $z \in \mathbb{R}$.

*Proof.* We will proceed by contradiction. Assume that there exists a pair $(\phi_1, \phi_2)$ satisfying (13)-(14). From the wave system (13), we have the expression, for $z \in \mathbb{R}$,

$$ \phi_2(z) = \frac{1}{\tilde{\delta}_2} \int_0^{+\infty} e^{\tilde{\gamma}_{21} x} [G_E^1 \star \phi_1 g(\phi_2)](z-x) dx + \frac{1}{\tilde{\delta}_2} \int_{-\infty}^{0} e^{\tilde{\gamma}_{22} x} [G_E^1 \star \phi_1 g(\phi_2)](z-x)(x) dx. $$

Then, with the help of Lemma 4, we obtain

$$ \begin{aligned} \int_{-\infty}^{+\infty} \phi_2(z) dz &= \frac{1}{\tilde{\delta}_2} \int_0^{+\infty} e^{\tilde{\gamma}_{21} x} \int_{-\infty}^{+\infty} [G_E^1 \star \phi_1 g(\phi_2)](z-x) dz dx \\ &\quad + \frac{1}{\tilde{\delta}_2} \int_{-\infty}^{0} e^{\tilde{\gamma}_{22} x} \int_{-\infty}^{+\infty} [G_E^1 \star \phi_1 g(\phi_2)](z-x) dz dx, \\ &= \frac{1}{\tilde{\delta}_2} \left[ \int_0^{+\infty} e^{\tilde{\gamma}_{21} x} dx + \int_{-\infty}^{0} e^{\tilde{\gamma}_{22} x} dx \right] \int_{-\infty}^{+\infty} [G_E^1 \star \phi_1 g(\phi_2)](z) dz, \\ &\le \frac{S_0}{\tilde{\delta}_2} \left[ \frac{1}{\tilde{\gamma}_{22}} - \frac{1}{\tilde{\gamma}_{21}} \right] \int_{-\infty}^{+\infty} [G_E^1 \star g(\phi_2)](z) dz. \end{aligned} $$
Since $\int_{0}^{+\infty} \int_{-\infty}^{+\infty} G_{E}^{1}(s, y) dy ds = 1$, then

$$ 
\begin{aligned}
\int_{-\infty}^{+\infty} \phi_2(z) dz &\le \frac{g'(0)S_0}{\tilde{\delta}_2} \left[ \frac{1}{\tilde{\gamma}_{22}} - \frac{1}{\tilde{\gamma}_{21}} \right] \int_{-\infty}^{+\infty} \phi_2(z) dz, \\
&= \frac{g'(0)S_0}{r} \int_{-\infty}^{+\infty} \phi_2(z) dz < \int_{-\infty}^{+\infty} \phi_2(z) dz.
\end{aligned}
$$

By the above strict inequality, we get a contradiction. This completes the proof. □

Now, we turn to show that our system admits non-trivial traveling wave solutions with speed $c \in (c^*, +\infty)$. We assume that

$$ \frac{S_0 g'(0)}{r} > 1. $$

The approach used below is to work with a suitable invariant set for the operator $(F_1, F_2)$ and then apply Schauder's fixed point theorem. Consider the constant $\alpha \in (\lambda_1, \min\{-\gamma_{11}, -\gamma_{21}\})$. Remember that we choose the constants $A_1$ and $A_2$ large enough. Consider the Banach space

$$ \mathbb{Q}_{\alpha}(\mathbb{R}, \mathbb{R}^2) = \left\{ (\phi_1, \phi_2) \in C(\mathbb{R}, \mathbb{R}^2) : \sup_{z \in \mathbb{R}} |\phi_1(z)|e^{-\alpha|z|} < +\infty, \quad \sup_{z \in \mathbb{R}} |\phi_2(z)|e^{-\alpha|z|} < +\infty \right\}, $$

for the exponential decay norm $|\cdot|_{\alpha}$ defined by

$$ |(\phi_1, \phi_2)|_{\alpha} = \max \left\{ \sup_{z \in \mathbb{R}} |\phi_1(z)| e^{-\alpha|z|}, \quad \sup_{z \in \mathbb{R}} |\phi_2(z)| e^{-\alpha|z|} \right\}. $$

In the following lemmas, we establish several results ensuring the existence of a fixed point for the operator $(F_1, F_2)$ in the subset $\mathbb{X}$ of $\mathbb{Q}_{\alpha}(\mathbb{R}, \mathbb{R}^2)$.

**Lemma 6.** 1. The set $\mathbb{X}$ is nonempty, bounded, closed and convex in $(\mathbb{Q}_{\alpha}(\mathbb{R}, \mathbb{R}^2), |\cdot|_{\alpha})$.

2. The operator $(F_1, F_2)$ maps $\mathbb{X}$ into $\mathbb{X}$.

*Proof.* It is not difficult to establish the statements of the point 1. We just focus on the proof of point 2. We consider $(\phi_1, \phi_2) \in \mathbb{X}$ and we verify the following inequalities

$$ \underline{\phi}_1(z) \le F_1(\phi_1, \phi_2)(z) \le \overline{\phi}_1(z) \quad \text{and} \quad \underline{\phi}_2(z) \le F_2(\phi_1, \phi_2)(z) \le \overline{\phi}_2(z). $$

First, we show that $\underline{\phi}_1(z) \le F_1(\phi_1, \phi_2)(z) \le \overline{\phi}_1(z)$, for all $z \in \mathbb{R}$. Recall that $F_1$ is defined by

$$ F_1(\phi_1, \phi_2)(z) = \frac{g'(0)}{\delta_1} \int_{-\infty}^{z} e^{\gamma_{11}(z-x)} H_1(\phi_1, \phi_2)(x) dx + \frac{g'(0)}{\delta_1} \int_{z}^{+\infty} e^{\gamma_{12}(z-x)} H_1(\phi_1, \phi_2)(x) dx. $$

Since $(\phi_1, \phi_2) \in \mathbb{X}$, then $\phi_1(z) \le \overline{\phi}_1(z) = S_0$, for all $z \in \mathbb{R}$. Hence, we get

$$ H_1(\phi_1, \phi_2)(z) = A_1\phi_1(z) - \frac{1}{g'(0)}\phi_1(z)g(\phi_2(z)) + \frac{h}{g'(0)}(\psi * \phi_1)(z) \le \left(A_1 + \frac{h}{g'(0)}\right) S_0. $$
This leads to

$$
\begin{align*}
F_1(\phi_1, \phi_2)(z) &\le \frac{(g'(0)A_1 + h)S_0}{\delta_1} \left[ \int_{-\infty}^{z} e^{\gamma_{11}(z-s)} ds + \int_{z}^{+\infty} e^{\gamma_{12}(z-s)} ds \right], \\
&= \frac{(g'(0)A_1 + h)S_0}{\delta_1} \left( \frac{1}{\gamma_{12}} - \frac{1}{\gamma_{11}} \right) = S_0 = \bar{\phi}_1(z).
\end{align*}
$$

Moreover, by using the fact that $\underline{\phi}_1(z) \le \phi_1(z)$ and $\phi_2(z) \le \bar{\phi}_2(z)$ for all $z \in \mathbb{R}$, and since $A_1 > 0$ satisfies (19), we have for $z \neq z_1$,

$$
\begin{align*}
c\underline{\phi}_1'(z) &\le D_S\underline{\phi}_1''(z) - h\underline{\phi}_1(z) - \underline{\phi}_1(z)g(\bar{\phi}_2(z)) + h(\psi * \underline{\phi}_1)(z), \\
&\le D_S\underline{\phi}_1''(z) - h\underline{\phi}_1(z) - \underline{\phi}_1(z)g(\phi_2(z)) + h(\psi * \phi_1)(z), \\
&= D_S\underline{\phi}_1''(z) - h\underline{\phi}_1(z) - A_1g'(0)\underline{\phi}_1(z) + g'(0)\underline{\phi}_1(z) \left(A_1 - \frac{1}{g'(0)}g(\phi_2(z))\right) + h(\psi * \underline{\phi}_1)(z), \\
&\le D_S\underline{\phi}_1''(z) - h\underline{\phi}_1(z) - A_1g'(0)\underline{\phi}_1(z) + g'(0)\underline{\phi}_1(z) \left(A_1 - \frac{1}{g'(0)}g(\phi_2(z))\right) + h(\psi * \underline{\phi}_1)(z).
\end{align*}
$$

We can rewrite this inequality as

$$
\begin{align*}
g'(0)H_1(\phi_1, \phi_2)(z) &= g'(0)\phi_1(z) \left( A_1 - \frac{1}{g'(0)}g(\phi_2(z)) \right) + h(\psi * \phi_1)(z), \\
&\geq c\phi_1'(z) - D_S\phi_1''(z) + h\phi_1(z) + A_1g'(0)\phi_1(z).
\end{align*}
$$

So, we obtain that

$$
F_1(\phi_1, \phi_2)(z) \geq \frac{1}{\delta_1} \int_{-\infty}^{z} e^{\gamma_{11}(z-s)} (c\underline{\phi}_1'(s) - D_S\underline{\phi}_1''(s) + (g'(0)A_1 + h)\underline{\phi}_1(s))ds \\
+ \frac{1}{\delta_1} \int_{z}^{+\infty} e^{\gamma_{12}(z-s)} (c\underline{\phi}_1'(s) - D_S\underline{\phi}_1''(s) + (g'(0)A_1 + h)\underline{\phi}_1(s))ds.
$$

If $z \ge z_0$, we have $\underline{\phi}_1(z) = 0$. Then, we obtain

$$
\begin{align*}
F_1(\phi_1, \phi_2)(z) &\ge \frac{1}{\delta_1} \int_{-\infty}^{z_0} e^{\gamma_{11}(z-s)} [c\underline{\phi}_1'(s) - D_S\underline{\phi}_1''(s) + (g'(0)A_1 + h)\underline{\phi}_1(s)] ds \\
&\quad + \frac{1}{\delta_1} \int_{z_0}^{z} e^{\gamma_{11}(z-s)} [c\underline{\phi}_1'(s) - D_S\underline{\phi}_1''(s)] ds + \frac{1}{\delta_1} \int_{z}^{+\infty} e^{\gamma_{12}(z-s)} [c\underline{\phi}_1'(s) - D_S\underline{\phi}_1''(s)] ds, \\
&= \frac{1}{\delta_1} e^{\gamma_{11}(z-z_0)} [c\underline{\phi}_1(z_0) - D_S\underline{\phi}_1'(z_0 - 0) - c\underline{\phi}_1(z_0) + D_S\underline{\phi}_1'(z_0 + 0)] + \frac{1}{\delta_1} [-D_S\underline{\phi}_1'(z) + D_S\underline{\phi}_1'(z)] \\
&\quad + \frac{1}{\delta_1} \int_{-\infty}^{z_0} e^{\gamma_{11}(z-s)} [\gamma_{11}c\underline{\phi}_1(s) - \gamma_{11}D_S\underline{\phi}_1'(s) + (g'(0)A_1 + h)\underline{\phi}_1(s)] ds \\
&\quad - \frac{1}{\delta_1} \int_{z_0}^{z} e^{\gamma_{11}(z-s)} \gamma_{11}D_S\underline{\phi}_1'(s) ds - \frac{1}{\delta_1} \int_{z}^{+\infty} e^{\gamma_{12}(z-s)} \gamma_{12}D_S\underline{\phi}_1'(s) ds.
\end{align*}
$$
Since $\underline{\phi_1}'(z_0 + 0) = 0 > \underline{\phi_1}'(z_0 - 0) = -\mu_1 \epsilon_1 e^{\epsilon_1 z_0}$, we obtain

$$F_1(\phi_1, \phi_2)(z) \geq \frac{D_S}{\delta_1} e^{\gamma_{11}(z-z_1)} [\underline{\phi_1}'(z_0 + 0) - \underline{\phi_1}'(z_0 - 0)] \\ + \frac{1}{\delta_1} \int_{-\infty}^{z_0} e^{\gamma_{11}(z-s)} [c\gamma_{11} - D_S\gamma_{11}^2 + g'(0)A_1 + h]\underline{\phi_1}(s)ds \geq 0 = \underline{\phi_1}(z).$$

If $z < z_0$, by the same procedure, we get

$$F_1(\phi_1, \phi_2)(z) \geq \frac{D_S}{\delta_1} e^{\gamma_{12}(z-z_1)} [\underline{\phi_1}'(z_0 + 0) - \underline{\phi_1}'(z_0 - 0)] \\ - \frac{D_S}{\delta_1} \gamma_{11} \underline{\phi_1}(z) + \frac{1}{\delta_1} \int_{-\infty}^{z} e^{\gamma_{11}(z-s)} [c\gamma_{11} - D_S\gamma_{11}^2 + g'(0)A_1 + h]\underline{\phi_1}(s)ds \\ + \frac{D_S}{\delta_1} \gamma_{12} \underline{\phi_1}(z) + \frac{1}{\delta_1} \int_{z}^{z_0} e^{\gamma_{12}(z-s)} [c\gamma_{12} - D_S\gamma_{12}^2 + g'(0)A_1 + h]\underline{\phi_1}(s)ds, \\ \geq \frac{D_S(\gamma_{12} - \gamma_{11})}{\delta_1} \underline{\phi_1}(z) = \underline{\phi_1}(z).$$

Next, we show that, for all $z \in \mathbb{R}$,

$$\underline{\phi_2}(z) \leq F_2(\phi_1, \phi_2)(z) \leq \overline{\phi_2}(z).$$

As the function $g$ is monotone increasing and, for all $z \in \mathbb{R}$, we have $\phi_1(z) \geq \underline{\phi_1}(z)$ and $\phi_2(z) \geq \underline{\phi_2}(z)$, then

$$c\underline{\phi_2}'(z) \leq D_I\underline{\phi_2}''(z) + \int_{0}^{+\infty} \int_{-\infty}^{+\infty} G_E^1(s, y)\underline{\phi_1}(z - y - cs)g(\underline{\phi_2}(z - y - cs))dyds - r\underline{\phi_2}(z), \\ \leq D_I\underline{\phi_2}''(z) - A_2\underline{\phi_2}(z) \\ + \int_{0}^{+\infty} \int_{-\infty}^{+\infty} G_E^1(s, y)\phi_1(z - y - cs)g(\phi_2(z - y - cs))dyds + (A_2 - r)\phi_2(z).$$
Consequently, we obtain

$$
\begin{align*}
F_2(\phi_1, \phi_2)(z) &= \frac{g'(0)}{\delta_2} \int_{-\infty}^{z} e^{\gamma_{21}(z-x)} \left[ \frac{1}{g'(0)} [G_E^1 \star \phi_1 g(\phi_2)](x) + \frac{A_2 - r}{g'(0)} \phi_2(x) \right] dx \\
&\quad + \frac{g'(0)}{\delta_2} \int_{z}^{+\infty} e^{\gamma_{22}(z-x)} \left[ \frac{1}{g'(0)} [G_E^1 \star \phi_1 g(\phi_2)](x) + \frac{A_2 - r}{g'(0)} \phi_2(x) \right] dx, \\
&= \frac{1}{\delta_2} \int_{-\infty}^{z} e^{\gamma_{21}(z-x)} \left[ (G_E^1 \star \phi_1 g(\phi_2))(x) + (A_2 - r)\phi_2(x) \right] dx \\
&\quad + \frac{1}{\delta_2} \int_{z}^{+\infty} e^{\gamma_{22}(z-x)} \left[ (G_E^1 \star \phi_1 g(\phi_2))(x) + (A_2 - r)\phi_2(x) \right] dx, \\
&\geq \frac{1}{\delta_2} \int_{-\infty}^{z} e^{\gamma_{21}(z-x)} \left[ c\phi_2'(x) - D_I\phi_2''(x) + A_2\phi_2(x) \right] dx \\
&\quad + \frac{1}{\delta_2} \int_{z}^{+\infty} e^{\gamma_{22}(z-x)} \left[ c\phi_2'(x) - D_I\phi_2''(x) + A_2\phi_2(x) \right] dx = \phi_2(z).
\end{align*}
$$

On another hand, for all $z \in \mathbb{R}$, we have

$$
\phi_1(z) \leq \bar{\phi}_1(z) = S_0 \quad \text{and} \quad \phi_2(z) \leq \bar{\phi}_2(z).
$$

So, $\overline{\phi_2}$ satisfies

$$
\begin{align*}
c\overline{\phi_2'}(z) &\geq D_I\overline{\phi_2}''(z) + S_0[G_E^1 \star g(\overline{\phi_2})](z) - r\overline{\phi_2}(z), \\
&\geq D_I\overline{\phi_2}''(z) - A_2\overline{\phi_2}(z) \\
&\quad + \int_0^{+\infty} \int_{-\infty}^{+\infty} G_E^1(s,y)\phi_1(z-y-cs)g(\phi_2(z-y-cs))dyds + (A_2-r)\phi_2(z).
\end{align*}
$$

Therefore, we obtain

$$
\begin{align*}
F_2(\phi_1, \phi_2)(z) &= \frac{1}{\delta_2} \int_{-\infty}^{z} e^{\gamma_{21}(z-x)} \left[ (G_E^1 \star \phi_1 g(\phi_2))(x) + (A_2 - r)\phi_2(x) \right] dx \\
&\quad + \frac{1}{\delta_2} \int_{z}^{+\infty} e^{\gamma_{22}(z-x)} \left[ (G_E^1 \star \phi_1 g(\phi_2))(x) + (A_2 - r)\phi_2(x) \right] dx, \\
&\leq \frac{1}{\delta_2} \int_{-\infty}^{z} e^{\gamma_{21}(z-x)} \left[ c\overline{\phi_2'}(x) - D_I\overline{\phi_2}''(x) + A_2\overline{\phi_2}(x) \right] dx \\
&\quad + \frac{1}{\delta_2} \int_{z}^{+\infty} e^{\gamma_{22}(z-x)} \left[ c\overline{\phi_2'}(x) - D_I\overline{\phi_2}''(x) + A_2\overline{\phi_2}(x) \right] dx = \overline{\phi_2}(z).
\end{align*}
$$

This completes the proof.

<figure><img src="image_3.png" /></figure>
We can also state the following lemma.

**Lemma 7.** • The operator $(F_1, F_2)$ is continuous with respect to the norm $|\cdot|_{\alpha}$ in $Q_{\alpha}(\mathbb{R}, \mathbb{R}^2)$.

• The operator $(F_1, F_2) : \mathbb{X} \to \mathbb{X}$ is compact with respect to the norm $|\cdot|_{\alpha}$ in $Q_{\alpha}(\mathbb{R}, \mathbb{R}^2)$.

*Proof.* We set for all $(\phi_1, \phi_2) \in \mathbb{X}$,

$$ \varphi(\phi_1, \phi_2) = \phi_1 g(\phi_2). $$

By the assumptions on $g$ and according to the mean value theorem, for all $\Phi = (\phi_1, \phi_2) \in \mathbb{X}$ and $\Psi = (\tilde{\phi}_1, \tilde{\phi}_2) \in \mathbb{X}$ and for some $M > 0$, we obtain

$$ |\varphi(\Phi) - \varphi(\Psi)| \leq M(|\phi_1 - \tilde{\phi}_1| + |\phi_2 - \tilde{\phi}_2|). \quad (24) $$

From the expression of the operator $H_1$, we have for some $M' > 0$,

$$ \left| \left[ A_1 \phi_1 - \frac{1}{g'(0)} \varphi(\phi_1, \phi_2) + \frac{h}{g'(0)} (\psi * \phi_1) \right] - \left[ A_1 \tilde{\phi}_1 - \frac{1}{g'(0)} \varphi(\tilde{\phi}_1, \tilde{\phi}_2) + \frac{h}{g'(0)} (\psi * \tilde{\phi}_1) \right] \right| \\ \leq M' (|\phi_1 - \tilde{\phi}_1| + |\phi_2 - \tilde{\phi}_2|) + \frac{h}{g'(0)} |(\psi * \phi_1) - (\psi * \tilde{\phi}_1)|. $$

Hence, from the expression of the operator $F_1$, we get ($M_1 > 0$)

$$ \begin{aligned} & |F_1(\phi_1, \phi_2)(z) - F_1(\tilde{\phi}_1, \tilde{\phi}_2)(z)| e^{-\alpha|z|} \\ & \leq \frac{M_1 g'(0)}{\delta_1} \left[ \int_{-\infty}^{z} e^{\gamma_{11}(z-x)+\alpha|x|-\alpha|z|} dx + \int_{z}^{+\infty} e^{\gamma_{12}(z-x)+\alpha|x|-\alpha|z|} dx \right] \\ & \times \left( |\phi_1 - \tilde{\phi}_1|_{\alpha} + |\phi_2 - \tilde{\phi}_2|_{\alpha} + |(\psi * \phi_1) - (\psi * \tilde{\phi}_1)|_{\alpha} \right). \end{aligned} $$

We put

$$ \mathcal{P}(z) := \int_{-\infty}^{z} e^{\gamma_{11}(z-x)+\alpha|x|-\alpha|z|} dx + \int_{z}^{+\infty} e^{\gamma_{12}(z-x)+\alpha|x|-\alpha|z|} dx > 0. $$

Recall that $0 < \alpha < -\gamma_{11} < \gamma_{12}$ and $0 < \alpha < -\gamma_{21} < \gamma_{22}$. So, we have

$$ \mathcal{P}(-\infty) = \frac{1}{\alpha + \gamma_{12}} - \frac{1}{\alpha + \gamma_{11}} > 0 \quad \text{and} \quad \mathcal{P}(+\infty) = \frac{1}{\alpha - \gamma_{11}} + \frac{1}{\gamma_{12} - \alpha}. $$

This means that $\mathcal{P}$ is uniformly bounded on $\mathbb{R}$. It is sufficient to prove that $|\phi_1 - \tilde{\phi}_1|_{\alpha} \to 0$ implies $|(\psi * \phi_1) - (\psi * \tilde{\phi}_1)|_{\alpha} \to 0$. In fact, since $\psi \in L^1(\mathbb{R}, \mathbb{R}^+)$, there exists $N > 0$ such that

$$ \int_{\mathbb{R}/[-N,N]} \psi(y) dy < \epsilon. $$

We note that $|\phi_1(z) - \tilde{\phi}_1(z)| \leq S_0$, for all $z \in \mathbb{R}$. Moreover, we have

$$ e^{-\alpha|z|} \int_{-N}^{N} |\phi_1(z-y) - \tilde{\phi}_1(z-y)| \psi(y) dy \leq |\phi_1 - \tilde{\phi}_1|_{\alpha} e^{-\alpha|z|} e^{\alpha(|z|+N)} \int_{-N}^{N} \psi(y) dy \leq \epsilon, $$
when $|\phi_1 - \tilde{\phi}_1|_\alpha \le \epsilon e^{-\alpha N}$. Consequently, $|(\psi * \phi_1) - (\psi * \tilde{\phi}_1)|_\alpha \to 0$ as $|\phi_1 - \tilde{\phi}_1|_\alpha \to 0$. This leads to the continuity of the component $F_1$ with respect to the norm $|\cdot|_\alpha$. Next, we focus on $F_2$. From (24) and the expression of $H_2$, we get for some $M_2 > 0$, the inequality

$$ \left| \left[ \frac{1}{g'(0)} [G_E^1 \star \varphi(\phi_1, \phi_2)] + \frac{A_2 - r}{g'(0)} \phi_2 \right] - \left[ \frac{1}{g'(0)} [G_E^1 \star \varphi(\tilde{\phi}_1, \tilde{\phi}_2)] + \frac{A_2 - r}{g'(0)} \tilde{\phi}_2 \right] \right| \\ \le M_2 (|\phi_2 - \tilde{\phi}_2| + G_E^1 \star |\phi_1 - \tilde{\phi}_1| + G_E^1 \star |\phi_2 - \tilde{\phi}_2|). $$

Consequently, we obtain that

$$ |F_2(\phi_1, \phi_2)(z) - F_2(\tilde{\phi}_1, \tilde{\phi}_2)(z)|e^{-\alpha|z|} \\ \le \frac{M_2 g'(0)}{\delta_2} \left[ \int_{-\infty}^{z} e^{\gamma_{21}(z-x)} e^{-\alpha|z|} (|\phi_2(x) - \tilde{\phi}_2(x)| + (G_E^1 \star |\phi_1 - \tilde{\phi}_1|)(x) + (G_E^1 \star |\phi_2 - \tilde{\phi}_2|)(x)) dx \right] \\ + \frac{M_2 g'(0)}{\delta_2} \left[ \int_{z}^{+\infty} e^{\gamma_{22}(z-x)} e^{-\alpha|z|} (|\phi_2(x) - \tilde{\phi}_2(x)| + (G_E^1 \star |\phi_1 - \tilde{\phi}_1|)(x) + (G_E^1 \star |\phi_2 - \tilde{\phi}_2|)(x)) dx \right]. $$

Therefore,

$$ |F_2(\phi_1, \phi_2)(z) - F_2(\tilde{\phi}_1, \tilde{\phi}_2)(z)|e^{-\alpha|z|} \le \frac{M_2 g'(0)}{\delta_2} (|\phi_2 - \tilde{\phi}_2|_\alpha + |G_E^1 \star |\phi_1 - \tilde{\phi}_1||_\alpha + |G_E^1 \star |\phi_2 - \tilde{\phi}_2||_\alpha) \mathcal{K}(z), $$

with

$$ \mathcal{K}(z) := \int_{-\infty}^{z} e^{\gamma_{21}(z-x)+\alpha|x|-\alpha|z|} dx + \int_{z}^{+\infty} e^{\gamma_{22}(z-x)+\alpha|x|-\alpha|z|} dx. $$

Clearly, we have

$$ \mathcal{K}(-\infty) = \frac{1}{\alpha + \gamma_{22}} - \frac{1}{\alpha + \gamma_{21}} > 0 \quad \text{and} \quad \mathcal{K}(+\infty) = \frac{1}{\gamma_{22} - \alpha} + \frac{1}{\alpha - \gamma_{21}} > 0. $$

Then, $\mathcal{K}$ is uniformly bounded on $\mathbb{R}$. On another hand, we note that the operator $F_2$ is almost the same as the corresponding one in [31]. So, we can show as in the proof of Lemma 2.8 of [31] and as we did above, that $|G_E^1 \star |\phi_1 - \tilde{\phi}_1||_\alpha \to 0$ and $|G_E^1 \star |\phi_2 - \tilde{\phi}_2||_\alpha \to 0$ as $|\phi_1 - \tilde{\phi}_1|_\alpha \to 0$ and $|\phi_2 - \tilde{\phi}_2|_\alpha \to 0$, respectively. Thus, the component $F_2$ is continuous with respect to the norm $|\cdot|_\alpha$. The conclusion is that the operator $(F_1, F_2)$ is continuous with respect to the norm $|\cdot|_\alpha$ in $\mathbb{Q}_\alpha(\mathbb{R}, \mathbb{R}^2)$.

For the second point of the lemma, the proof of compactness can be easily adapted from the corresponding result in [31], with only a few modifications. We therefore omit the details of this proof. The proof is completed. □

Taking into account all the previous lemmas, we obtain the following theorem, which states the existence of a traveling wave solution for the system (6)-(8).

**Theorem 8.** Suppose that

$$ \frac{S_0 g'(0)}{r} > 1, $$
and let $c^*$ (minimum wave speed) be given by Lemma 1. For every $c > c^*$, the system (6) admits a traveling wave solution,

$$ (S(t, x), I(t, x), u(t, x)) = (\phi_1(x + ct), \phi_2(x + ct), \phi_3(x + ct)) $$

satisfying the asymptotic boundary conditions (8).

*Proof.* Suppose that $c > c^*$. Then, according to Schauder's fixed point theorem, Lemma 6 and Lemma 7, we ensure the existence of a pair $(\phi_1, \phi_2) \in \mathcal{X}$, which is a fixed point of the operator $(F_1, F_2)$. We have

$$ \phi_1(z) = \frac{g'(0)}{\delta_1} \int_{-\infty}^{z} e^{\gamma_{11}(z-x)} H_1(\phi_1, \phi_2)(x)dx + \frac{g'(0)}{\delta_1} \int_{z}^{+\infty} e^{\gamma_{12}(z-x)} H_1(\phi_1, \phi_2)(x)dx, $$

and

$$ \phi_2(z) = \frac{g'(0)}{\delta_2} \int_{-\infty}^{z} e^{\gamma_{21}(z-x)} H_2(\phi_1, \phi_2)(x)dx + \frac{g'(0)}{\delta_2} \int_{z}^{+\infty} e^{\gamma_{22}(z-x)} H_2(\phi_1, \phi_2)(x)dx. $$

By a direct verification, we can check that $(\phi_1, \phi_2)$ satisfies the wave problem (13). As shown previously, this is sufficient to get the existence of a non-negative and nontrivial solution for the wave system (6). Next, we show the asymptotic behavior (8). In fact, we know that for all $z \in \mathbb{R}$,

$$ \underline{\phi_1}(z) \le \phi_1(z) \le \overline{\phi_1}(z) \quad \text{and} \quad \underline{\phi_2}(z) \le \phi_2(z) \le \overline{\phi_2}(z). $$

Then, we get

$$ \lim_{z \to -\infty} \phi_1(z) = S_0 \quad \text{and} \quad \lim_{z \to -\infty} \phi_2(z) = 0. $$

Furthermore, $0 \le \phi_1(z) \le S_0$ for $z \in \mathbb{R}$. Now, it remains to prove that

$$ \lim_{z \to +\infty} \phi_2(z) = 0 \quad \text{and} \quad \lim_{z \to +\infty} \phi_1(z) = S_{\infty} < S_0. $$

From the wave system (13), we have, for $z \in \mathbb{R}$,

$$ \phi_2(z) = \frac{1}{\tilde{\delta}_2} \int_{0}^{+\infty} e^{\tilde{\gamma}_{21}x} [G_E^1 \star \phi_1 g(\phi_2)](z-x)dx + \frac{1}{\tilde{\delta}_2} \int_{-\infty}^{0} e^{\tilde{\gamma}_{22}x} [G_E^1 \star \phi_1 g(\phi_2)](z-x)(x)dx, $$

where

$$ \tilde{\gamma}_{21} := \frac{c}{2D_I} \left( 1 - \sqrt{1 + \frac{4D_I r}{c^2}} \right), \quad \tilde{\gamma}_{22} := \frac{c}{2D_I} \left( 1 - \sqrt{1 + \frac{4D_I r}{c^2}} \right) \quad \text{and} \quad \tilde{\delta}_2 := D_I(\tilde{\gamma}_{22} - \tilde{\gamma}_{21}). $$

From Lemma 4, we have

$$ \int_{-\infty}^{+\infty} \phi_1(s)g(\phi_2(s))ds = \tilde{M} < +\infty. $$
Then, by applying Fubini's theorem, we obtain

$$ \int_{-\infty}^{+\infty} \phi_2(z) dz = \frac{1}{\tilde{\delta}_2} \int_{0}^{+\infty} \int_{-\infty}^{+\infty} G_E^1(s, y) dy ds \left[ \int_{0}^{+\infty} e^{\tilde{\gamma}_{21}s} ds + \int_{-\infty}^{0} e^{\tilde{\gamma}_{22}s} ds \right] \tilde{M} < +\infty. $$

Also, we can show that $\phi_2'(z)$ is uniformly bounded and with the combination with the fact that $\phi_2(z) > 0$ is integrable on $\mathbb{R}$, we get

$$ \lim_{z \to +\infty} \phi_2(z) = 0. $$

Next, we show that

$$ \lim_{z \to +\infty} \phi_1(z) = S_{\infty} < S_0. $$

Recall from the proof of Theorem 5 that

$$ |J(z)| = h| \int_{-\infty}^{z} (\psi * \phi_1)(s) - \phi_1(s) ds|. $$

We can show easily that $\lim_{z \to +\infty} J(z) = 0$. By integrating the equation of $\phi_1$ from $-\infty$ to $z$, we get the following limit

$$ \lim_{z \to +\infty} [c\phi_1(z) - D_S\phi_1'(z)] = cS_0 - \int_{-\infty}^{+\infty} \phi_1(s)g(\phi_2(s))ds := \mathcal{L}. \quad (25) $$

Hence, for any $\epsilon > 0$, there exists $z^* \in \mathbb{R}$ such that

$$ c\phi_1(z) - D_S\phi_1'(z) \geq \mathcal{L} - \epsilon, \quad z \geq z^*. $$

This means that, for any $\xi \geq z^*$, we have

$$ \begin{aligned} \phi_1(z) &\leq \phi_1(\xi)e^{\frac{c}{D_S}(z-\xi)} - \frac{1}{D_S}\int_{\xi}^{z} e^{\frac{c}{D_S}(z-s)}[\mathcal{L}-\epsilon]ds, \\ &= \left[\phi_1(\xi) - \frac{\mathcal{L}-\epsilon}{c}\right]e^{\frac{c}{D_S}(z-\xi)} + \frac{\mathcal{L}-\epsilon}{c}. \end{aligned} $$

The above inequality implies that

$$ \phi_1(\xi) \geq \frac{\mathcal{L} - \epsilon}{c}. $$

Because if we suppose that $\phi_1(\xi) < (\mathcal{L} - \epsilon)/c$, then $\lim_{z \to +\infty} \phi_1(z) = -\infty$, and we find a contradiction. As $\xi \geq z^*$ is chosen arbitrarily, we have

$$ \phi_1(z) \geq \frac{\mathcal{L} - \epsilon}{c}, \quad z \geq z^*. $$

Moreover, (25) implies that there exists $\tilde{z} \geq z^*$ such that, for all $z \geq \tilde{z}$,

$$ c\phi_1(z) - D_S\phi_1'(z) \leq \mathcal{L} + \epsilon. $$
On the other hand, we have, for any $\xi \ge z^*$,

$$ \phi_1(z) \ge \left[ \phi_1(\xi) - \frac{\mathcal{L} + \epsilon}{c} \right] e^{\frac{c}{D_s}(z-\xi)} + \frac{\mathcal{L} + \epsilon}{c}. $$

Then, we get

$$ \phi_1(\xi) \le \frac{\mathcal{L} + \epsilon}{c}. $$

Otherwise, we obtain $\lim_{z \to +\infty} \phi_1(z) = +\infty$, which is again a contradiction. Since $\xi \ge \tilde{z}$ is arbitrary, we have, $\phi_1(\xi) \le (\mathcal{L} + \epsilon)/c$ for all $z \ge \tilde{z}$. As $\tilde{z} \ge z^*$, then we have $z \ge \tilde{z}$, $(\mathcal{L} - \epsilon)/c \le \phi_1(z) \le (\mathcal{L} + \epsilon)/c$. Recall that $\epsilon > 0$ is also arbitrary, then, the previous inequalities give $\lim_{z \to +\infty} \phi_1(z) = \mathcal{L}/c =: S_\infty$. As we know that $\phi_1(z) \ge \underline{\phi_1}(z) \ge 0$ for all $z \in \mathbb{R}$, then we have $S_\infty \ge 0$. Finally, by

$$ S_\infty = \frac{\mathcal{L}}{c} = S_0 - \frac{1}{c} \int_{-\infty}^{+\infty} \phi_1(s)g(\phi_2(s))ds, $$

we get $S_\infty < S_0$. The proof is completed. $\square$

# 6 Discussion and numerical simulations

This work is devoted to the study of the existence of nontrivial traveling wave solutions for a system of reaction-diffusion equations with nonlocal terms in time and space. The presence or absence of nontrivial traveling wave solutions is important in epidemiology, as it indicates whether an epidemic has the potential to spread through the population, and also gives the speed of the propagation. The model considered is new because, in addition to having nonlocal terms, it is written in a hybrid form as a coupling between a reaction-diffusion system and a difference equation. We have already investigated this type of hybrid problem in cell population dynamics [3] and epidemiological problem [4]. Nevertheless, the mathematical difficulties associated with this new model are more significant. This is mainly due to the fact of having a double convolution due to the nonlocal terms that take into account the movement of individuals during both the protection phase and the exposure phase.

We start by reducing the main three-equation traveling wave model, (7), to a system with only two equations but two nonlocal terms, (13). For this reduced system (13), we construct an upper and a lower solution and use them to establish a framework for the application of Schauder's fixed point theorem (see Theorem 8 and Lemma 3). Then, under the following classical threshold

$$ \frac{S_0 g'(0)}{r} > 1, $$

we prove the existence of $c^*$ (minimum wave speed) such that for $c \in (c^*, +\infty)$ we have the existence of traveling wave solutions connecting the disease-free equilibrium to an endemic equilibrium. We performed numerical simulations, Figures 2 and 3, that confirm the theoretical results on wave propagation in our model. A traveling wave solution appears under some initial conditions. We note that both $\tau$ (duration of the protection period) and $\eta$ (rate of re-entry into the protection compartment) affect the monotonicity of the S-wave and the u-wave. In many other works, see for instance [31], the authors proved that the traveling wave solutions
<figure><img src="image_4.png" /><figcaption>Figure 2: The solution of System (3) (or equivalently (4)) is plotted. The saturated incidence rate used is $g(I) = \beta I/(1+aI)$ and the parameters are arbitrarily fixed as follows: $h = 0.1$, $\beta = 2$, $a = 0.5$, $\eta = 0.1$, $D_S = D_p = D_I = D_E = 0.05$, $r = 0.8$, $\mu = 1$, and $S_0 = 1$. In this case, $S_0 g'(0)/r = 2.5 > 1$. For this figure, we choose a sufficiently small delay: $\tau = 5$. We can see that a small delay induces the monotonicity of the S-wave and the u-wave.</figcaption></figure>

are monotone. This is not the case in our model due to the difference equation, which is a consequence of the protection phase (time-delay). For other examples where the authors obtained monotone S-waves, we can mention [17, 31, 36, 40]. For non-monotone S-waves, we can refer to the papers [4, 12, 13]. For $c \in (0, c^*)$, we have the nonexistence of traveling wave solutions. The proof of this result is classical, based on the Laplace transform, and follows the same approach as [31] (we have not reproduced it here), see also [8, 9, 30]. The case of $c = c^*$ is more difficult to treat and is left open for the moment.

For the continuation of our research, we plan to study hybrid systems including other types of incidence functions (bilinear, for example). The method used in this article is no longer appropriate. We plan to develop other techniques inspired by other epidemic models that can be adapted to the hybrid model to prove the existence of traveling wave solutions.

## Acknowledgements

A.C thanks grant, PRFU: C00L03UN29012022002, from DGRSDT of Algeria. A.C. thanks also the IMU-CDC and Simons Foundation for partially supporting this work with a grant.

## Conflict of interest statement

This work does not have any conflicts of interest.
<figure><img src="image_5.png" /><figcaption>Figure 3: The solution of System (3) (or equivalently (4)) is plotted. The saturated incidence rate used is $g(I) = \beta I/(1 + aI)$ and the parameters are arbitrarily set as follows: $h = 0.1$, $\beta = 2$, $a = 0.5$, $\eta = 0.1$, $D_S = D_p = D_I = D_E = 0.05$, $r = 0.8$, $\mu = 1$, and $S_0 = 1$. Using these values, we obtain $S_0 g'(0)/r = 2.5 > 1$. For this figure, we choose a sufficiently large delay: $\tau = 20$. We can see that the S-wave and the u-wave are not always monotone due to the protection phase (with a large delay).</figcaption></figure>

References

[1] M. Adimy, A. Chekroun, C. Dugourd-Camus, and H. Meghelli. Global asymptotic stability of a hybrid differential-difference system describing SIR and SIS epidemic models with a protection phase and a nonlinear force of infection. Qual. Theory Dyn. Syst., 23(34):1–31, 2024.

[2] M. Adimy, A. Chekroun, and C. P. Ferreira. Global dynamics of a differential-difference system: a case of Kermack-McKendrick SIR model with age-structured protection phase. Mathematical Biosciences and Engineering, 17(2):1329, 2020.

[3] M. Adimy, A. Chekroun, and B. Kazmierczak. Traveling waves in a coupled reaction-diffusion and difference model of hematopoiesis. Journal of Differential Equations, 262(7):4085–4128, 2017.

[4] M. Adimy, A. Chekroun, and T. Kuniya. Traveling waves of a differential-difference diffusive Kermack-McKendrick epidemic model with age-structured protection phase. Journal of Mathematical Analysis and Applications, 505(1):125464, 2022.

[5] M. Adimy, A. Chekroun, and T.-M. Touaoula. Global asymptotic stability for an age-structured model of hematopoietic stem cell dynamics. Applicable Analysis, 96(3):429–440, 2016.
[6] F. Brauer and C. Castillo-Chavez. *Mathematical models in population biology and epidemiology*. Texts in Applied Mathematics, Springer, New York, 2nd edition, 2012.

[7] V. Capasso and G. Serio. A generalization of the Kermack-McKendrick deterministic epidemic model. *Mathematical Biosciences*, 42(1):43–61, 1978.

[8] J. Carr and A. Chmaj. Uniqueness of travelling waves for nonlocal monostable equations. *Proceedings of the American Mathematical Society*, 132(8):2433–2439, 2004.

[9] A. Chekroun. Existence and asymptotics of traveling wave fronts for a coupled nonlocal diffusion and difference system with delay. *Electron. J. Qual. Theory Differ. Equ.*, 85:1–23, 2019.

[10] O. Diekmann and J. Heesterbeek. *Mathematical Epidemiology of Infective Diseases: Model Building, Analysis and Interpretation*. Wiley, New York, 2000.

[11] O. Diekmann, J. A. P. Heesterbeek, and J. A. J. Metz. On the definition and the computation of the basic reproduction ratio $R_0$ in models for infectious diseases in heterogeneous populations. *Journal of Mathematical Biology*, 28:365–382, 1990.

[12] A. Ducrot and P. Magal. Travelling wave solutions for an infection-age structured model with diffusion. *Proceedings of the Royal Society of Edinburgh Section A: Mathematics*, 139(3):459–482, 2009.

[13] H. El Mamouni, K. Hattaf, and N. Yousfi. Existence of traveling waves by means of fixed point theory for an epidemic model with hattaf-yousfi incidence rate and temporary immunity acquired by vaccination. *Adv. Fixed Point Theory*, 13(3):1–20, 2023.

[14] S. Gao, L. Chen, J. J. Nieto, and A. Torres. Analysis of a delayed epidemic model with pulse vaccination and saturation incidence. *Vaccine*, 24(35):6037–6045, 2006.

[15] S. Gao, L. Chen, and Z. Teng. Pulse vaccination of an SEIR epidemic model with time delay. *Nonlinear Analysis: Real World Applications*, 9(2):599–607, 2008.

[16] S.-A. Gourley and J. Wu. Delayed non-local diffusive systems in biological invasion and disease spread. *Nonlinear Dynamics and Evolution Equations, Amer. Math. Soc*, 48:137–200, 2006.

[17] Y. Hosono and B. Ilyas. Traveling waves for a simple diffusive epidemic model. *Mathematical Models and Methods in Applied Sciences*, 05(07):935–966, 1995.

[18] H.-F. Huo, S.-K. Hu, and H. Xiang. Traveling wave solution for a diffusion SEIR epidemic model with self-protection and treatment. *Electronic Research Archive*, 29(3):2325–2358, 2021.

[19] W. O. Kermack and A. G. McKendrick. A contribution to the mathematical theory of epidemics. *Proc. R. Soc. Lond. A*, 115(772):700–721, 1927.

[20] M. Y. Li. *An Introduction to Mathematical Modeling of Infectious Diseases*. Mathematics of Planet Earth, Mathematics of Planet Earth Collection. Springer International Publishing, 2018.
[21] W.-T. Li, G. Lin, C. Ma, and F.-Y. Yang. Traveling wave solutions of a nonlocal delayed SIR model without outbreak threshold. *Discrete & Continuous Dynamical Systems B*, 19:467–484, 2014.

[22] M. Martcheva. *An Introduction to Mathematical Epidemiology*. Springer Publishing Company, Incorporated, 1st edition, 2015.

[23] J. D. Murray. *Mathematical Biology: I. An Introduction*. Springer, New York, 3rd edition edition, 2002.

[24] C. V. Pao. *Nonlinear Parabolic and Elliptic Equations*. Springer New York, NY, 1992.

[25] S. Ruan. *Spatial-Temporal Dynamics in Nonlocal Epidemiological Models*, pages 97–122. Springer Berlin Heidelberg, Berlin, Heidelberg, 2007.

[26] S. Ruan and J. Wu. Modeling spatial spread of communicable diseases involving animal hosts, in: S. Cantrell, C. Cosner, S. Ruan (Eds.), *Spatial Ecology, Chapman and Hall/CRC, New York*. pages 293–316, 2009.

[27] H. Shu, X. Pan, B. Wade, and X.-S. Wang. Traveling waves of nonlocal delayed disease models: critical wave speed and propagation speed. *Applicable Analysis*, 102(2):385–405, 2023.

[28] B. Tian and R. Yuan. Traveling waves for a diffusive SEIR epidemic model with non-local reaction. *Applied Mathematical Modelling*, 50:432–449, 2017.

[29] A. I. Volpert, V. A. Volpert, and V. A. Volpert. *Traveling Wave Solutions of Parabolic Systems*. American Mathematical Soc., 1994.

[30] Z.-C. Wang, W.-T. Li, and S. Ruan. Traveling fronts in monostable equations with nonlocal delayed effects. *J Dyn Diff Equat*, 20(3):573–607, 2008.

[31] Z.-C. Wang and J. Wu. Travelling waves of a diffusive Kermack-McKendrick epidemic model with non-local delayed transmission. *Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences*, 466(2113):237–261, 2010.

[32] P. Weng and X.-Q. Zhao. Spreading speed and traveling waves for a multi-type SIS epidemic model. *Journal of Differential Equations*, 229(1):270–296, 2006.

[33] J. Wu. *Theory and Applications of Partial Functional Differential Equations*. Springer, 1996.

[34] X. Wu and Z. Ma. Wave propagation in a diffusive SEIR epidemic model with nonlocal transmission and a general nonlinear incidence rate. *Bound Value Probl*, 87:1–32, 2021.

[35] X. Wu, B. Tian, and R. Yuan. Traveling waves in an SEIR epidemic model with a general nonlinear incidence rate. *Applicable Analysis*, 99(1):133–157, 2020.

[36] Z. Xu. Traveling waves in a Kermack-McKendrick epidemic model with diffusion and latent period. *Nonlinear Analysis: Theory, Methods & Applications*, 111:66–81, 2014.
[37] Z. Xu. Traveling waves for a diffusive SEIR epidemic model. *Communications on Pure & Applied Analysis*, 15(3):871, 2016.

[38] X.-Q. Zhao and W. Wang. Fisher waves in an epidemic model. *Discrete and Continuous Dynamical Systems. Series B*, 4(4):1117–1128, 2004.

[39] X.-Q. Zhao and D. Xiao. The asymptotic speed of spread and traveling waves for a vector disease model. *J Dyn Diff Equat*, 18(4):1001–1019, 2006.

[40] Z. Zhen, J. Wei, L. Tian, J. Zhou, and W. Chen. Wave propagation in a diffusive SIR epidemic model with spatiotemporal delay. *Mathematical Methods in the Applied Sciences*, 41(16):7074–7098, 2018.