On a model for epidemic spread with interpopulation contact and repellent taxis

TL;DR

This paper introduces a PDE model for epidemic spread considering interpopulation contact and repellent taxis, analyzing the existence of solutions and illustrating dynamics through numerical simulations.

Contribution

It develops a PDE framework incorporating diffusion and repellent taxis for susceptibles and nonlinear diffusion for infecteds, with existence proofs and simulations.

Findings

Existence of weak-strong solutions in 1D

Existence of supersolutions in higher dimensions

Numerical simulations illustrating space-time dynamics

Abstract

We study a PDE model for dynamics of susceptible-infected interactions. The dispersal of susceptibles is via diffusion and repellent taxis as they move away from the increasing density of infected. The diffusion of infected is a nonlinear, possibly degenerating term in nondivergence form. We prove the existence of so-called weak-strong solutions in 1D for a positive susceptible initial population. For dimension $N\geq 2$ and nonnegative susceptible initial density we show the existence of supersolutions. Numerical simulations are performed for different scenarios and illustrate the space-time behaviour of solutions.

Tables

Table 1. Table 1: Parameters

$C_1$	$C_2$	$C_3$	$ϵ$	${\lambda }_S$	${\lambda }_I$	${\mu }_I$	${\mu }_S$
$0.1$	$0.2$	$0.3$	$0.25$	$0.5$	$0.5$	$0.05$	$0.01$

Equations

$${\partial }_{t}S=\nabla \cdot \left(\nabla S+\chi (S)S\nabla I\right)+f(S,I)$$

$${\partial }_{t}I=S\Delta I+g(S,I)$$

$${\partial }_{\nu }S={\partial }_{\nu }I=0$$

$$S(0,\cdot )=S_0,I(0,\cdot )=I_0$$

$$\chi (S)=K(1-S),$$

$$f(S,I)=-{\lambda }_S\frac{SI}{S+I}+{\mu }_{S}S(1-S),g(S,I)={\lambda }_I\frac{SI}{S+I}-{\mu }_{I}I,$$

$${\partial }_{t}R={\mu }_{I}I.$$

$${\partial }_t{S}_{\epsilon }=\nabla \cdot \left(\nabla S_{\epsilon }+\chi (S_{\epsilon })S_{\epsilon }\nabla I_{\epsilon }\right)+f(S_{\epsilon },I_{\epsilon })$$

$${\partial }_t{I}_{\epsilon }=(\epsilon +S_{\epsilon })\Delta I_{\epsilon }+g(S_{\epsilon },I_{\epsilon })$$

$${\partial }_{\nu }{S}_{\epsilon }={\partial }_{\nu }{I}_{\epsilon }=0$$

$$S_{\epsilon }(0,\cdot )=S_0,I_{\epsilon }(0,\cdot )=I_0$$

$$S_0\in L^{\infty }(\Omega ),$$

$$I_0\in W^{1,\infty }(\Omega ),$$

$$N=1\text{and}\underset{\Omega }{\inf }{S}_0\in (0,1]$$

$$\Vert \nabla I_{\epsilon }{\Vert }_{L^{\infty }({\mathbb{R}}^{+};L^2(\Omega ))}\le \text{C},$$

$$\Vert I_{\epsilon }{\Vert }_{L^2(0,T;L^2(\Omega ))}\le \text{C}(T),$$

$$\sqrt{\epsilon }{\Vert \Delta I_{\epsilon }\Vert }_{L^2(0,T;L^2(\Omega ))}+{\Vert \sqrt{S_{\epsilon }}\Delta I_{\epsilon }\Vert }_{L^2(0,T;L^2(\Omega ))}\le \text{C}(T)$$

$$\Vert {\partial }_t{I}_{\epsilon }{\Vert }_{L^2(0,T;L^2(\Omega ))}\le \text{C}(T),$$

$$\Vert \nabla S_{\epsilon }{\Vert }_{L^2(0,T;L^2(\Omega ))}\le \text{C}(T),$$

$$\Vert {\partial }_t{S}_{\epsilon }{\Vert }_{L^2(0,T;(H^1(\Omega ){)}^{\prime })}\le \text{C}(T)$$

$$\frac{1}{2}\frac{d}{dt}\Vert \nabla I_{\epsilon }{\Vert }_{L^2(\Omega )}^2=$$

$$\le$$

$$\frac{1}{2}\frac{d}{dt}\Vert I_{\epsilon }{\Vert }_{L^2(\Omega )}^2$$

$$=$$

$$\le$$

$$\frac{1}{2}\frac{d}{dt}\Vert S_{\epsilon }{\Vert }_{L^2(\Omega )}^2=$$

$$\le$$

$$\nabla I_{{\epsilon }_n}\underset{n\to \infty }{\stackrel{\ast }{\rightharpoonup }}\nabla I$$

$${\epsilon }_n\Delta I_{{\epsilon }_n}\underset{n\to \infty }{\to }0$$

$${\partial }_t{I}_{{\epsilon }_n}\underset{n\to \infty }{\rightharpoonup }{\partial }_{t}I$$

$$I_{{\epsilon }_n}\underset{n\to \infty }{\to }I$$

$$\nabla S_{{\epsilon }_n}\underset{n\to \infty }{\rightharpoonup }\nabla S$$

$${\partial }_t{S}_{{\epsilon }_n}\underset{n\to \infty }{\rightharpoonup }{\partial }_{t}S$$

$$S_{{\epsilon }_n}\underset{n\to \infty }{\to }S$$

$$⟨{\partial }_t{S}_{\epsilon },\phi ⟩=-{\int }_{\Omega }(\nabla S_{\epsilon }+\chi (S_{\epsilon })S_{\epsilon }\nabla I_{\epsilon })\cdot \nabla \phi dx+{\int }_{\Omega }f(S_{\epsilon },I_{\epsilon })\phi dx\text{a.e. in }{\mathbb{R}}_0^{+},$$

$${\partial }_t{I}_{\epsilon }=S_{\epsilon }\Delta I_{\epsilon }+g(S_{\epsilon },I_{\epsilon })\text{a.e. in }{\mathbb{R}}_0^{+}\times \Omega ,$$

$${\partial }_{\nu }{I}_{\epsilon }=0\text{a.e. in }{\mathbb{R}}_0^{+}\times \partial \Omega ,$$

$$S_{\epsilon }(0,\cdot )=S_0,I_{\epsilon }(0,\cdot )=I_0\text{a.e. in }\Omega .$$

$${\partial }_t{I}_{\epsilon }=(\epsilon +{\bar{S}}_{\epsilon })\Delta I_{\epsilon }+g({\bar{S}}_{\epsilon },I_{\epsilon })$$

$${\partial }_{\nu }{I}_{\epsilon }=0$$

$$I_{\epsilon }(0,\cdot )=I_0$$

$${\partial }_t{S}_{\epsilon }=\nabla \cdot \left(\nabla S_{\epsilon }+\chi (S_{\epsilon })S_{\epsilon }\nabla I_{\epsilon }\right)+f(S_{\epsilon },I_{\epsilon })$$

$${\partial }_{\nu }{S}_{\epsilon }=0$$

Full text

Keywords: epidemic model, nonlinear diffusion with interpopulation contact, repellent taxis, non-divergence form, weak solution, weak supersolution

MSC 2010: 35Q92, 35K55, 92D30

1 Introduction

There exists by now a vast literature about models of epidemic spread, most of which take the form of ODE systems. The heterogeneity of space, however, can play an essential role in the dynamics of an infectious disease, as the environmental conditions can differ from one site to the other and moreover the individuals in the populations can move in space. Therefore various settings accounting for both space and time variability have been proposed, among the first being the contact model with diffusion in [14], followed by [9, 8, 10, 13] and many others, mainly performing traveling wave analysis. 333In this context we are interested only in reaction-diffusion models; population balance models accounting for further structures and featuring integro-differential PDEs (see e.g. [18] and the references therein) are not addressed here. During the last two decades the investigation of well-posedness and qualitative properties of solutions to reaction-diffusion PDE systems describing epidemics has attracted increasing interest, see e.g. [1, 6, 11, 21, 23]; we also refer to [7, 22] for earlier works. Some of the more recent models [5, 4, 12, 15] were extended to account for at least one of the interacting populations having a motility bias in a certain direction (e.g., due to environmental influences like fluid or air flow), which leads to a drift term supplementing the diffusion and the source/decay terms.

A model with linear cross-diffusion of susceptibles has been considered in [20], while [3] proposed a numerical scheme to handle a nonlinear one. Such terms are included into the epidemic model in order to account for the response of an active population of susceptibles toward the other population’s degree of infectiveness. The former would typically try to avoid the latter by biasing its movement in the direction opposite to the gradient of infected population density. How effective this avoidance is depends, of course, on the amount of susceptibles and infected in the system and their respective ratio: A small amount of infected in an overwhelming susceptible population will pass -at least for a while- unnoticed (carelessness). In the other extreme, a large density of infected in a comparatively rather small population of susceptibles will lead to fatalism (very low possibility of avoidance). At moderate densities and ratios the avoidance mechanism can be quite effective, leading to patterns and (local) phase separation, see the simulations in [3, 20]. The corresponding advection term occurring in the PDE for the density of susceptibles can be interpreted as a repellent taxis term, in analogy with models of (chemo)repellence involving (chemo)taxis terms with a sign opposite to that usual for Keller-Segel type models.

In this work we propose and investigate a reaction-diffusion model for infection spread with contact and with repellent taxis, which combines cross-diffusion (in the sense mentioned above) with the influence of contacts between susceptibles and infected on the self-diffusion of the latter. To our knowledge, previous epidemic models have either one or the other of these features, see [14, 19, 20]; the model in [3] involves self-diffusion of infected with self-contact. The mentioned studies in [19, 20, 3] mostly focus upon computational aspects. They do not include any proofs of existence of solutions to such settings. Here we aim at addressing this analytical issue for the announced model extension.

The paper is structured as follows: Section 2 specifies the model to be investigated, along with the mathematical challenges arising from its structure, and introduces a family of regularized problems approximating the actual PDE system of interest. Section 3 is dedicated to the analysis of this regularized problem, followed in Section 4 by the existence proof for solutions to the full model in 1D and with positive initial density of susceptibles. Section 5 contains the existence proof for supersolutions of our model in space dimension $N\geq 2$ and for nonnegative initial density of susceptibles. Finally, Section 6 provides some numerical simulations to illustrate the behavior of the model for a few different scenarios of the epidemics. Some concluding remarks are provided as well.

2 Model setup

We consider the following model for the evolution of susceptibles ($S$) and infected ($I$):

[TABLE]

Here $\Omega$ is a bounded domain in $\mathbb{R}^{N}$, $N\in\mathbb{N}$, with the corresponding outer normal unit vector $\nu$ on $\partial\Omega$.

The equations describe the interactions of the two populations, whereby $S$ performs linear diffusion and repellent taxis in the sense mentioned in Section 1: Susceptibles tend to avoid the infected. The efficiency of avoidance is characterized by the function $\chi$ which in analogy to chemotaxis models will be called in the following tactic sensitivity. It may in fact depend both on $S$ and $I$, however for our analysis we choose it in the form

[TABLE]

which accounts for the crowding effect: amidst a large mass of susceptibles their awareness of infectives is reduced (’drowned’). In particular, the threshold value for $S$ corresponding to a tight packing state is assumed to be normalised to $S_{max}=1$.

The first term on the right hand side of (2.1b) describes self-diffusion with interpopulation contact, similarly to the model in [14]. Further, we modify the interaction terms having the usual form $SI$ to account for an infectiveness threshold with a limitation given by the total population. Thus, following e.g., [20] or [1] we choose

[TABLE]

where the second terms in these expressions describe as usual logistic growth of susceptibles and linear removal of infected, respectively. 444The function $f_{1}(S,I)=\frac{SI}{S+I}$ is Lipschitz with respect to $S$ and $I$ in the open first quadrant, hence its definition can be extended to the closure of that set by letting it be zero when either $S=0$ or $I=0$. The removed (including dead and recovered) population can be described by

[TABLE]

This equation is decoupled from (2.1), thus not contributing to the dynamics of $(S,I)$ and will therefore be ignored in the following.

Analytical challenges.

System (2.1) combines several effects which jointly make the analysis challenging:

1. (i)

(2.1a) is in divergence form, while (2.1b) is not, and they are strongly coupled; 2. (ii)

equation (2.1a) for $S$ includes a potentially destabilising chemotaxis transport term, in this case in the direction opposite to $\nabla I$; 3. (iii)

equation (2.1b) for $I$ features a non-standard degeneracy occurring on the zero level set of the variable $S$.

System (2.1) can be seen as a formal limit as $\varepsilon\rightarrow 0$ of the following family of regularised problems:

[TABLE]

Thereby a small number $\varepsilon\in(0,1]$ added to the diffusion coefficient of $I$ eliminates the degeneracy issue. Standard tools can be used (see Section 3 below) in order to prove the existence of solutions for (2.2). This observation naturally leads to an attempt to apply the standard compactness method which is based on establishing uniform w.r.t. $\varepsilon$ a priori estimates for $S_{\varepsilon}$, $I_{\varepsilon}$, and their derivatives in suitable Bochner spaces and utilising some known compact embeddings and other necessary results in order to prove the existence of a sequence which converges to some weak solution of the non-perturbed system. In general, however, owing to the sort of degeneracy present in the original system, it seems difficult, if not impossible, to get a priori bounds which would allow to pass rigorously to the limit in the term $S_{\varepsilon}\Delta I_{\varepsilon}$ in order to regain $S\Delta I$.

In order to obtain our existence results, we assume that

[TABLE]

In this work we consider first the special situation when

[TABLE]

and prove in Section 4 that under these assumptions a solution does exist. In Section 5 we then turn to the general case of an arbitrary space dimension and without a positive lower bound for $S_{0}$. In this case we are able to establish the existence of a weak supersolution (see Definition 5.1 below).

Remark 2.1** (Notation).**

Throughout the paper we make the following useful conventions:

For any index $i$, a quantity $C_{i}$ denotes a positive constant or function; 2. 2.

Dependence upon such parameters as: the space dimension $N$, domain $\Omega$, constants $K,\lambda_{I},\lambda_{S},\mu_{I},\mu_{S}$, and the norms of the initial data $S_{0}$ and $I_{0}$ is mostly not indicated in an explicit way; 3. 3.

We assume the reader to be familiar with the conventional Lebesgue, Sobolev, and Bochner spaces and standard results concerning them. We denote $\left<\cdot,\cdot\right>$ the duality paring between $H^{1}(\Omega)$ and its dual $(H^{1}(\Omega))^{\prime}$.

3 Analysis of the regularized system (2.2)

3.1 A priori estimates and compactness

To begin with, we establish several necessary a priori estimates for (2.2).

Lemma 3.1**.**

Let a pair of measurable functions $(S_{\varepsilon},I_{\varepsilon}):\mathbb{R}_{0}^{+}\times\overline{\Omega}\rightarrow[0,1]\times\mathbb{R}_{0}^{+}$ be a sufficiently regular solution to system (2.2). Then it satisfies the following estimates:

[TABLE]

Proof.

Testing (2.2b) with $-\Delta I_{\varepsilon}$ and using Young’s inequality we obtain that

[TABLE]

Using the Gronwall lemma and integrating with respect to time when necessary we conclude from (3.7) that estimates (3.1) and (3.3) hold. Testing (2.2b) with $I_{\varepsilon}$ and integrating over $\Omega$ we obtain with the Cauchy-Schwarz and Young inequalities that

[TABLE]

Combining the Gronwall lemma with (3.3) and (3.8) we obtain (3.2). Altogether, estimates (3.1)-(3.3) allow to estimate the right-hand side of (2.2b) yielding (3.4).

Next, we test (2.2a) with $S_{\varepsilon}$ and integrate over $\Omega$, by parts where necessary. We thus obtain by using the Cauchy-Schwarz and Young inequalities and the estimate (3.1) that

[TABLE]

Integrating (3.9) over $(0,T)$ we obtain (3.5). Finally, (3.1) and (3.5) allow to estimate the right-hand side of (2.2a) yielding (3.6). ∎

Lemma 3.2**.**

Let for each $\varepsilon\in(0,1]$ the pair $(S_{\varepsilon},I_{\varepsilon}):\mathbb{R}_{0}^{+}\times\overline{\Omega}\rightarrow[0,1]\times\mathbb{R}_{0}^{+}$ be a sufficiently regular solution to system (2.2). Then there exists a sequence $\varepsilon_{n}\rightarrow 0$ and a pair $(S,I)$ such that for any $T>0$

[TABLE]

Proof.

This is a direct consequence of uniform estimates (3.1)-(3.6) combined with the Lions-Aubin lemma and the Banach-Alaoglu theorem. ∎

3.2 Existence of solutions

Since system (2.2) couples two parabolic equations, one of which is in divergence form, while the other is not, the standard theory, e.g. from [16] or [2], seems not to be directly applicable. Still, existence of solutions can be obtained by using the well-established procedure based on the Schauder fixed point theorem. For the convenience of the reader, we state the corresponding existence result and sketch its proof. We choose the following notion of a solution:

Definition 3.3** (Weak-strong solution).**

We call a pair of measurable functions $(S_{\varepsilon},I_{\varepsilon}):\mathbb{R}_{0}^{+}\times\overline{\Omega}\rightarrow[0,1]\times\mathbb{R}_{0}^{+}$ a weak-strong solution to system (2.2) if:

$\nabla S_{\varepsilon}\in L^{2}_{loc}(\mathbb{R}_{0}^{+};L^{2}(\Omega))$, $\partial_{t}S_{\varepsilon}\in L^{2}_{loc}(\mathbb{R}_{0}^{+};(H^{1}(\Omega))^{\prime})$; 2. 2.

$I_{\varepsilon}\in L^{2}_{loc}(\mathbb{R}_{0}^{+};H^{2}(\Omega))$, $\nabla I_{\varepsilon}\in L^{\infty}(\mathbb{R}_{0}^{+};L^{2}(\Omega))$, $\partial_{t}I_{\varepsilon}\in L^{2}_{loc}(\mathbb{R}_{0}^{+};L^{2}(\Omega))$; 3. 3.

the pair $(S_{\varepsilon},I_{\varepsilon})$ is a weak solution to (2.1a) and a strong solution to (2.1b), i.e., for all $\varphi\in H^{1}(\Omega)$ it holds that

[TABLE]

Theorem 3.4** (Existence of a weak-strong solution).**

Let (2.3) hold. Then, system (2.2) possesses a weak-strong solution in terms of Definition 3.3.

Proof.

(Sketch) To begin with, we decouple the equations for the two components:

[TABLE]

and

[TABLE]

For smooth $\bar{S}_{\varepsilon}$ and $I_{0}$ standard parabolic theory [16] insures the existence of a unique classical solution $I_{\varepsilon}$ to (3.18). Similarly, such $I_{\varepsilon}$ and smooth $S_{0}$ lead to a unique classical solution $S_{\varepsilon}$ to (3.19). Moreover, it is clear from the proofs that results of Lemmas 3.1-3.2 continue to hold. They allow to obtain solutions with the regularity as stated in Definition 3.3 under assumption (2.3) by means of an approximation procedure.

Uniqueness can be established in both cases in the usual way by considering equations for differences of two solutions, testing with suitable test functions, performing estimates, and, finally, applying the Gronwall lemma. Here we only check uniqueness for (3.19): Let $S_{\varepsilon}^{(1)}$ and $S_{\varepsilon}^{(2)}$ be two solutions corresponding to some $I_{\varepsilon}$ and $S_{0}$. Set

[TABLE]

In this notation we have for $U$ the linear equation

[TABLE]

We need to verify that $U=0$ a.e. Observe that

[TABLE]

We are going to test (3.20) with

[TABLE]

and then pass to the limit as $\delta\rightarrow 0$. Since $\operatorname{sign}_{\delta}\in W^{1,\infty}(\mathbb{R})$, it holds that

[TABLE]

i.e., it is a valid test function. Using the weak chain and product rules and (3.21)-(3.23) where necessary, we thus compute that

[TABLE]

[TABLE]

[TABLE]

Combining (3.24)-(LABEL:est5) we obtain that

Finally, applying the Gronwall lemma to (LABEL:est10), we conclude that $\|U\|_{L^{1}(\Omega)}\equiv 0$.

Altogether, we have a well-defined operator $\Phi$ in the following setting:

[TABLE]

Thanks to Lemma 3.2 the image $\Phi({\cal M})$ is precompact in $X$. Moreover, this Lemma together with fact that both equations are uniquely solvable imply the continuity of $\Phi$. Therefore, the Schauder fixed point theorem implies the existence of a fixed point and, as a result, of a solution to (2.2).

∎

4 Existence of solutions to (2.1) for $N=1$ and $S_{0}>0$

Throughout this section we assume that

[TABLE]

and

[TABLE]

4.1 An a priori lower bound

To begin with, we return to the regularised problem (2.2) and establish an a priori uniform positive lower bound for $S_{\varepsilon}$.

Lemma 4.1**.**

Solutions to (2.2) satisfy

[TABLE]

Proof.

We use the standard method of propagation of $L^{p}$ bounds in order to derive a finite uniform upper bound for $S_{\varepsilon}^{-1}$. Its inverse gives a uniform positive lower bound for $S_{\varepsilon}$. Let $p\geq 1$. Multiplying (2.2a) by $-pS_{\varepsilon}^{-p-1}$ and integrating by parts over $\Omega$ we obtain using the Hölder, Young, and Gagliardo-Nirenberg inequalities as well as estimate (3.1) where necessary that

[TABLE]

Using the Gronwall lemma we conclude from (4.5) that

[TABLE]

Consequently,

[TABLE]

Further, integrating (4.4) over $(0,t)$ we obtain that

[TABLE]

Consequently,

[TABLE]

For $n\in\mathbb{N}$ we introduce

[TABLE]

Due to estimate (4.7) we have that

[TABLE]

whereas (4.6) implies that

[TABLE]

Using the standard recursive result [16, Chapter 2, Lemma 5.6] we conclude with (4.8) that for all $n\in\mathbb{N}$

[TABLE]

Combining (4.9)-(4.10) we conclude that

[TABLE]

Since $\|\cdot\|_{2^{n+1}}\rightarrow\|\cdot\|_{\infty}$ as $n\rightarrow\infty$, (4.11) implies (4.3) . ∎

4.2 Existence of solutions

Having obtained the estimate (4.3) we can now prove the existence of weak-strong solutions to (2.1). The definition is as follows:

Definition 4.2** (Weak-strong solution).**

We call a pair of measurable functions $(S,I):\mathbb{R}_{0}^{+}\times\overline{\Omega}\rightarrow[0,1]\times\mathbb{R}_{0}^{+}$ a weak-strong solution to system (2.1) if:

$\nabla S\in L^{2}_{loc}(\mathbb{R}_{0}^{+};L^{2}(\Omega))$, $\partial_{t}S\in L^{2}_{loc}(\mathbb{R}_{0}^{+};(H^{1}(\Omega))^{\prime})$; 2. 2.

$I\in L^{2}_{loc}(\mathbb{R}_{0}^{+};H^{2}(\Omega))$, $\nabla I\in L^{\infty}(\mathbb{R}_{0}^{+};L^{2}(\Omega))$, $\partial_{t}I\in L^{2}_{loc}(\mathbb{R}_{0}^{+};L^{2}(\Omega))$; 3. 3.

the pair $(S,I)$ is a weak solution to (2.1a) and a strong solution to (2.1b), i.e., for all $\varphi\in H^{1}(\Omega)$ it holds that

[TABLE]

Theorem 4.3** (Existence of a weak-strong solution).**

Under assumptions (2.3) and (4.1)-(4.2) system (2.1) possesses a weak-strong solution in terms of Definition 4.2.

Proof.

Our starting point is the weak formulation from Definition 3.3 and Lemma 3.2 on convergence. Thanks to the uniform estimates (3.3) and (4.3) and the Banach-Alaoglu theorem we may assume that the sequence from that Lemma is chosen in such a way that

[TABLE]

holds as well. Using (3.1)-(3.6) and (4.13) together with the dominated convergence theorem and compensated compactness, we can pass to the limit in (3.17) along the sequence $\varepsilon_{n}\rightarrow 0$ and thus obtain that $(S,I)$ satisfies all conditions from Definition 4.2. ∎

5 Existence of supersolutions to (2.1) for $N\geq 2$ and $S_{0}\geq 0$

In this Section we consider the general case of an arbitrary space dimension, assume $0\leq S_{0}\leq 1$, and prove the existence of a weak supersolution to (2.1). The definition is as follows:

Definition 5.1** (Weak supersolution).**

We call a pair of measurable functions $(S,I):\mathbb{R}_{0}^{+}\times\overline{\Omega}\rightarrow[0,1]\times\mathbb{R}_{0}^{+}$ a weak supersolution to system (2.1) if:

$\nabla S\in L^{2}_{loc}(\mathbb{R}_{0}^{+};L^{2}(\Omega))$, $\partial_{t}S\in L^{2}_{loc}(\mathbb{R}_{0}^{+};(H^{1}(\Omega))^{\prime})$; 2. 2.

$I\in L^{2}_{loc}(\mathbb{R}_{0}^{+};L^{2}(\Omega))$, $\nabla I\in L^{\infty}(\mathbb{R}_{0}^{+};L^{2}(\Omega))$, $\partial_{t}I\in L^{2}_{loc}(\mathbb{R}_{0}^{+};L^{2}(\Omega))$; 3. 3.

the pair $(S,I)$ is a weak solution to (2.1a) and a weak supersolution to (2.1b), i.e., for all $\varphi\in H^{1}(\Omega)$ and $0\leq\psi\in W^{1,\infty}(\Omega)$ it holds that

[TABLE]

Remark 5.2**.**

Recently weak (generalised) supersolutions in the form of a variational inequality have been used in order to provide a solution concept for models with positive chemotaxis, see e.g. [17, 25, 24]. In those cases, however, both equations are in divergent form, while equation (2.1b) is not. The latter precludes the possibility to close (5.1) by imposing a suitable mass control from above. As a result, even a smooth supersolution $(S,I)$ is not automatically a subsolution to (2.1).

Theorem 5.3** (Existence of a weak supersolution).**

Let (2.3) hold. Then, system (2.1) possesses a weak supersolution in terms of Definition 5.1.

Proof.

Once again, our starting point is the weak formulation from Definition 3.3 and Lemma 3.2 on convergence. To begin with, we construct a suitable reformulation of (2.2b). Since this equation is not in divergence form, the standard approach based on testing and integration by parts is not a good foundation for a limit procedure. It turns out useful to construct instead a variational identity which combines equations for both solution components.

Let $0\leq\eta\in L^{\infty}(0,T;W^{1,\infty}(\Omega))$, so that $\eta I_{\varepsilon}\in L^{2}(0,T;H^{1}(\Omega))$. Testing (3.17a) and (2.2b) with $-\eta I_{\varepsilon}$ and $\eta$, respectively, adding the results, integrating over $(0,T)$, and using the chain rule where necessary, we compute

[TABLE]

Using (3.10)-(3.16) together with the dominated convergence theorem and the compensated compactness, we can pass to the limit in (3.17a) and (3.17d) along the sequence $\varepsilon_{n}\rightarrow 0$ which yields (5.1a) and (5.1c), respectively. Owing to the presence of the quadratic term $|\nabla I_{\varepsilon}|^{2}$ in one of the integrals, we cannot justify the equality while passing to the limit in (5.2). Instead, we take limit inferior on both sides, use the above mentioned convergences and theorems, as well as the weak lower semicontinuity of a norm, and, finally, the chain rule where necessary and thus arrive at

[TABLE]

Plugging ${\varphi=}\eta I\in L^{2}(0,T;H^{1}(\Omega))$ into (3.17a), integrating over $(0,T)$, passing to the limit for $\varepsilon_{n}\to 0$, and then adding the result to (5.3) finally yields (5.1b). ∎

6 Numerical simulations and discussion

In order to illustrate the solution behavior we present in this section some 2D numerical simulation results for the system (2.1). A finite difference scheme with a first order upwind discretization of the repellent taxis term was used to produce them. Here we show contour plots for the populations of susceptibles/infected at various time points in a square domain $\Omega=[0,10]^{2}$. We consider complementary initial densities

[TABLE]

and use parameters given in Table 1.

Figure 1 shows simulations of the model in the cases with ($\chi=15$) and without ($\chi=0$) repellent taxis, respectively. The former exhibits a spread of the infection comparable with the latter case, but with a reduced suppression of the susceptibles and an overall slightly higher infected population density. Due to the infectives’ diffusion with contact, however, the avoidance efficiency is diminished, which for even larger values of $\chi$ leads to a more effective spread of the invasion. Further numerical simulations (not shown in this paper) suggest the existence of a critical value of $\chi$ above which the repellent taxis actually triggers the opposite effect. Moreover, here the sensitivity of susceptibles towards infected was taken to be linearly decreasing with $S$ which also contributed to the mentioned infection enhancement, but it is reasonable to assume its dependence also on $I$, more precisely on the interactions between the two populations. Analytically determining the critical $\chi$ range and its effect on the epidemic spread would be an interesting problem in the framework of travelling wave analysis.