\documentclass[reqno]{amsart} \usepackage{hyperref} \usepackage{graphicx} \AtBeginDocument{{\noindent\small \emph{Electronic Journal of Differential Equations}, Vol. 2016 (2016), No. 284, pp. 1--14.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu} \thanks{\copyright 2016 Texas State University.} \vspace{8mm}} \begin{document} \title[\hfilneg EJDE-2016/284\hfil Global stability of SVIR models] {Global asymptotic stability of a diffusive SVIR epidemic model with immigration of individuals} \author[S. Abdelmalek, S. Bendoukha \hfil EJDE-2016/284\hfilneg] {Salem Abdelmalek, Samir Bendoukha} \address{Salem Abdelmalek \newline Department of mathematics, University of Tebessa 12002, Algeria} \email{sallllm@gmail.com} \address{Samir Bendoukha \newline Electrical Engineering Department, College of Engineering at Yanbu, Taibah University, Saudi Arabia} \email{sbendoukha@taibahu.edu.sa} \thanks{Submitted May 10, 2016. Published October 24, 2016.} \subjclass[2010]{35K45, 35K57} \keywords{Reaction diffusion systems; SVIR; immigration; Lyapunov functional; \hfill\break\indent large time behavior} \begin{abstract} In this article, we consider a spatially SVIR model of infectious disease epidemics which allows for continuous immigration of all classes of individuals. We show that the proposed model has a unique steady state that is asymptotically stable. Using an appropriately constructed Lyapunov functional, we establish its global asymptotic stability. Numerical results obtained through Matlab simulations are presented to confirm the results. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{lemma}[theorem]{Lemma} \newtheorem{proposition}[theorem]{Proposition} \allowdisplaybreaks \section{Introduction} In this article, we are concerned with reaction-diffusion models of disease epidemics. Of the many models available in the literature, see \cite{Busenberg1993}, we will deal with one of the suceptible-vaccinated-infectious-recovered (SVIR) type, which as the name suggests takes into consideration four classes of individuals according to their relation to the disease. Numerous recent publications can be found in the literature regarding the subject. In the following is a brief description of the most relevant of these studies. Liu et al.\ \cite{Liu2008} presented two different models to represent the two vaccination strategies: continuous and pulse and showed that the dynamics of both models depend on the basic reproduction number. The study of Kuniya \cite{Kuniya2013} considered a multi-group SVIR model that allows for the heterogeneity of the population and the effect of immunity induced by the vaccination. Results showed that the long time behaviour of the model depends on the basic reproductive number. In \cite{Duan2014}, Duan et al.\ examined an ODE SVIR model which allows for the vaccinated individuals to become suceptible again after a certain period of time as the vaccine loses its cover. They studied the dynamics of the model based on LaSalle's invariance principle and appropriately constructed Lyapunov functionals and showed that the global stability of the equilibriums depend only upon the basic reproductive number. In \cite{Henshaw2015}, Henshaw and McCluskey studied the local and global asymptotic stability of an ODE SVIR model with immigration of individuals. The model they proposed is the basis of the work that will be presented in this paper. Our aim is to show that the inclusion of spatial spreading in the model does not affect the asymptotic stability of the equilibrium. The work carried out here is analogous to that of Abdelmalek et al.\ \cite{Abdelmalek2015}, where they studied the asymptotic stability of an SEI model including immigration of all classes of individuals. The remainder of this paper consists of three sections. Section \ref{SecModel} will present the proposed system model and identify its main characteristics and the conditions on the parameters. Section \ref{SecSteadyState} will examine the main properties of the steady state solutions. Section \ref{SecGlobal} will prove that the unique steady state of the model is globally asymptotically stable using an appropriate Lyapunov functional. \section{System model} \label{SecModel} In this article, we study the SVIR epidemic model with immigration of individuals, \begin{equation} \begin{gathered} \partial_{t}u-d_1\Delta u=\Lambda_1-uf(w) -(\mu+\alpha) u:=f_1(u,v,w) \quad \text{in }\mathbb{R}^{+}\times \Omega,\\ \partial_{t}v-d_2\Delta v=\Lambda_2+\alpha u-vg(w) -( \mu+\beta) v:=f_2(u,v,w) \quad \text{in }\mathbb{R}^{+}\times \Omega,\\ \partial_{t}w-d_3\Delta w=\Lambda_3+uf(w) +vg(w) -(\mu+\gamma+\delta) w :=f_3(u,v,w) \quad \text{in }\mathbb{R}^{+}\times \Omega,\\ \partial_{t}R-d_{4}\Delta R=\Lambda_{4}+\beta v+\delta w-\mu R :=f_{4}(v,w,R) \quad \text{in }\mathbb{R}^{+}\times \Omega, \end{gathered}\label{2.1} \end{equation} where $\Omega$ is an open bounded subset of $\mathbb{R}^{n}$ with piecewise smooth boundary $\partial\Omega$. We assume the initial conditions \begin{equation} u_0(x) =u(x,0),\quad v_0(x) =v(x,0) ,\quad w_0(x) =w(x,0),\quad R_0(x) =R(x,0) ,\quad \text{in } \Omega,\label{2.2} \end{equation} where $u_0(x), v_0(x), w_0(x),R_0(x) \in C^2(\Omega) \cap C^{0}(\overline{\Omega}) $, and homogoneous Neumann boundary conditions \begin{equation} \frac{\partial u}{\partial\nu} =\frac{\partial v}{\partial\nu} =\frac{\partial w}{\partial\nu}=\frac{\partial R}{\partial\nu} =0\quad \text{on }\mathbb{R}^{+}\times\partial\Omega,\label{2.3} \end{equation} with $\nu$ being the unit outer normal to $\partial\Omega$. We will also assume that the initial conditions $u_0(x) ,v_0(x) ,w_0(x) ,R_0(x) \in \mathbb{R}_{\geq0}$. Note that this model is similar to that proposed in \cite{Henshaw2015} but with the inclusion of spatial diffusion. In the proposed model, the positive functions $u(x,t) ,v(x,t) ,w(x,t) ,R(x,t) \geq0$ represent the population distributions of four classes of people: suceptible, vaccinated, infectious, and recovered, respectively. However, since the recovered class $R$ does not have an impact on the remaining classes, it will be omitted in the sequel. The parameters $\Lambda_i>0$ denote the growth of the different classes of individuals whether through birth or immigration and migration. The parameter $\alpha$ denotes the rate at which the suceptible population is vaccinated. In this model, death can either be attributed to the infectious disease or to other reasons. The per capita death rate for the former is denoted by $\gamma$, whereas the latter is denoted by $\mu>0$. Since in reality, it takes a while for the vaccinated individual to develop full immunity, the parameter $\beta$ has been introduced here indicating an average duration $\frac{1}{\beta}$. The parameter $\delta$ is introduced to allow for some of the infected individuals to recover on their own after a duration $\frac{1}{\delta}$. We will assume that $\alpha,\beta,\gamma,\delta\geq0$. The transfer diagram shown in Figure \ref{Sys_Model} presents a summary of the proposed model. The model \eqref{2.1}--\eqref{2.3} includes the spatial spreading of the individuals. The parameters $d_i\geq0$ represent the diffusivity constants modelling the movement of a certain class as a result of its distribution. \begin{figure}[htb] \label{Sys_Model} \begin{center} \includegraphics[width = 3.5in]{fig1} % Paper9_Model.pdf \end{center} \caption{A transfer diagram of the proposed system model.} \end{figure} The functions $f(w) $ and $g(w)$ are known as the incidence functions allowing for a nonlinear relation between the first three classes of individuals. We will assume that the incidence functions satisfy the following conditions for all $w\geq$ $0$: \begin{itemize} \item[(H1)] $f(w) ,g(w) \geq 0$ with equality if and only if $w=0$, \item[(H2)] $f'(w),g'(w) \geq0$, \item[(H3)] $f''(w) ,g''(w) \leq0$, \item[(H4)] $g(w) \leq f(w) $. \end{itemize} In addition, note that for $(u,v,w) \in\mathbb{R}_{\geq0}^{3}$, we have \begin{gather*} f_1(0,v,w) =\Lambda_1\geq0,\\ f_2(u,0,w) =\Lambda_2+uf(w) \geq0,\\ f_3(u,v,0) =\Lambda_3+\beta v\geq0. \end{gather*} Hence, the function $(f_1,f_2,f_2) ^{T}$ is \emph{essentially nonnegative}. Then, the non-negative octant $\mathbb{R}_{\geq0}^{3}$ is an invariant set (see \cite[Proposition 2.1]{Haddad2010} and \cite[page 288]{Quitner2007}). \section{Steady states and stability\label{SecSteadyState}} \subsection{ODE Case} Before we determine the steady state solutions to our proposed model \eqref{2.1}--\eqref{2.3} and their asymptotic stability, let us recall the results obtained by Henshaw and McCluskey in \cite{Henshaw2015}. We mentioned previously that the fourth equation of the system \eqref{2.1}--\eqref{2.3} will be omitted as it has no impact on the remaining three. In the absence of diffusion, the proposed system reduces to \begin{equation} \begin{gathered} \partial_{t}u=\Lambda_1-uf(w) -(\mu+\alpha) u,\\ \partial_{t}v=\Lambda_2+\alpha u-vg(w) -(\mu +\beta) v,\\ \partial_{t}w=\Lambda_3+uf(w) +vg(w) -(\mu+\gamma+\delta) w. \end{gathered} \label{3.1} \end{equation} First, let us define \[ \Lambda=\Lambda_1+\Lambda_2+\Lambda_3, \] and for any $\epsilon\geq0$, \begin{equation} D_{\epsilon}=\big\{ (u,v,w) :u,v,w>\epsilon\text{ and } u+v+w\leq\frac{\Lambda}{\mu}\big\} . \end{equation} System \eqref{3.1} was shown in \cite{Henshaw2015} to have the positively invariant non-negative octant $\mathbb{R}_{\geq0}^{3}$ and that there exists a number $\epsilon>0$ such that $D_{\epsilon}$ is non-empty, attracting and positively invariant. Henshaw and McCluskey also showed that the system has a unique equilibrium in the attraction region $(u^{\ast},v^{\ast},w^{\ast}) \in D_{\epsilon}$. This equilibrium is the solution of the system \begin{equation} \begin{gathered} \Lambda_1=u^{\ast}f(w^{\ast}) +(\mu+\alpha) u^{\ast}\\ \Lambda_2=-\alpha u^{\ast}+v^{\ast}g(w^{\ast}) +(\mu+\beta) v^{\ast}\\ (\mu+\gamma+\delta) =\frac{\Lambda_3+u^{\ast}f(w^{\ast}) +vg(w^{\ast}) }{w^{\ast}}. \end{gathered} \label{3.3} \end{equation} To determine the local stability of this unique equilibrium, we need to examine the eigenvalues of the Jacobian. The Jacobian and its second additive compound (see \eqref{App2}) are \[ J=\begin{pmatrix} -H_0-\mu-\alpha & 0 & -F_2\\ \alpha & -H_1-\mu-\beta & -G_2\\ H_0 & H_1 & -H_2 \end{pmatrix} \] and \[ J^{[2] }=\begin{pmatrix} -H_0-H_1-\alpha-\beta-2\mu & -G_2 & F_2\\ H_1 & -H_0-H_2-\mu-\alpha & 0\\ -H_0 & \alpha & -H_1-H_2-\mu-\beta \end{pmatrix}, \] respectively, where \begin{equation} \begin{gathered} F_2=u^{\ast}f'(w^{\ast}) \geq0,\quad G_2=v^{\ast}g'(w^{\ast}) \geq0,\\ H_0=f(w^{\ast}) \geq0,\quad H_1=g(w^{\ast}) \geq0,\\ H_2=(\mu+\gamma+\delta) -u^{\ast}f'(w^{\ast }) -v^{\ast}g'(w^{\ast}) \geq0. \end{gathered} \label{3.4} \end{equation} The positivity of the terms $F_2,G_2,H_0,H_1$ is trivial. The term $H_2$, however, requires a careful attention. Note that by applying Proposition \ref{Proposition} in the Appendix to $f$ and $g$ and using the third equation of \eqref{3.3},\ we have \begin{align*} H_2 & \geq(\mu+\gamma+\delta) -u^{\ast}\frac{f( w^{\ast}) }{w^{\ast}}-v^{\ast}\frac{g(w^{\ast}) } {w^{\ast}}\\ & =\frac{(\mu+\gamma+\delta) w^{\ast}-u^{\ast}f(w^{\ast }) -v^{\ast}g(w^{\ast}) }{w^{\ast}}\\ & =\frac{\Lambda_3}{w^{\ast}}\geq0. \end{align*} For information about the meaning and properties of additive compounds, we refer to \cite{Muldowney1990}. The local stability of the equilibrium can be examined by looking at the determinant of the Jacobian $\det(J)$, its trace $\operatorname{tr}(J) $, and the determinant of its second additive compound $\det(J^{[2]}) $ and ensuring that they are all negative (see Proposition \ref{PropositionM}). We have \begin{gather} \begin{aligned} \det J&=-\alpha F_2H_1-(H_0+\mu+\alpha) [(H_1+\mu+\beta) H_2+G_2H_1] \\ &\quad -F_2(H_1+\mu +\beta) H_0, \end{aligned} \label{3.5} \\ \operatorname{tr}J=-(H_0+H_1+H_2+\alpha+\beta+2\mu) ,\label{3.6} \\ \begin{aligned} \det J^{[2] } & =F_2[\alpha H_1-H_0( H_0+H_2+\mu+\alpha) ] -G_2H_1(H_1+H_2 +\mu+\beta) \\ &\quad -(H_0+H_1+\alpha+\beta+2\mu) (H_0+H_2 +\mu+\alpha) (H_1+H_2+\mu+\beta) . \end{aligned} \label{3.7} \end{gather} It is evident that $\det J<0$ and $\operatorname{tr}J<0$. However, for $\det J^{[2] }$, the term $\alpha H_1-H_0(H_0+H_2 +\mu+\alpha) $ needs to be examined. Using condition (H4), we have $H_1\leq H_0$, leading to \begin{align*} \alpha H_1-H_0(H_0+H_2+\mu+\alpha) & \leq\alpha H_0-H_0(H_0+H_2+\mu+\alpha) \\ & =-H_0(H_0+H_2+\mu) \leq0. \end{align*} Therefore, we see that $\det J^{[2] }<0$. Hence, as shown in \cite{Henshaw2015}, the unique equilibrium $(u^{\ast},v^{\ast},w^{\ast}) $ is in fact locally asymptotically stable. \subsection{Properties of the steady states} In this subsection, we shall discuss the basic properties of the non-homogeneous steady states of the proposed epidemic model \eqref{2.1}--\eqref{2.3}. In the presence of diffusion, the steady state solution satisfies \begin{equation} \begin{gathered} d_1\Delta u+\Lambda_1-u^{\ast}f(w^{\ast}) -( \mu+\alpha) u^{\ast}=0,\\ d_2\Delta v+\Lambda_2+\alpha u^{\ast}-v^{\ast}g(w^{\ast}) -(\mu+\beta) v^{\ast}=0,\\ d_3\Delta w+\Lambda_3+u^{\ast}f(w^{\ast}) +v^{\ast}g( w^{\ast}) -(\mu+\gamma+\delta) w^{\ast}=0. \end{gathered} \label{4.1} \end{equation} subject to the homogeneous Neumann boundary condition $\frac{\partial u}{\partial\nu}=\frac{\partial v}{\partial\nu} =\frac{\partial w}{\partial\nu}=0$ for all\ $x\in\partial\Omega$. Let $0=\lambda_0<\lambda_1\leq\lambda_2\leq\dots$. be the sequence of eigenvalues for the elliptic operator ($-\Delta$) subject to the homogeneous Neumann boundary condition on $\Omega$, where each $\lambda_i $ has multiplicity $m_i\geq1$. Also let $\Phi_{ij},1\leq j\leq m_i $, (recall that $\Phi_0=$ const and $\lambda_i\to\infty$ at $i\to \infty$) be the normalized eigenfunctions corresponding to $\lambda_i$. That is, $\Phi_{ij}$ and $\lambda_i$ satisfy $-\Delta\Phi_{ij}=\lambda_i \Phi_{ij}$ in $\Omega$, with $\frac{\partial\Phi_{ij}}{\partial\nu}=0$ in $\partial\Omega$, and $\int_{\Omega}\Phi_{ij}^2(x) dx=1$. \begin{theorem}\label{PropLocal} The constant steady state $(u^{\ast},v^{\ast},w^{\ast}) $ is asymptotically stable. \end{theorem} \begin{proof} Let us define the linearizing operator \[ \mathcal{L}=\begin{pmatrix} -d_1\Delta-(H_0+\mu+\alpha) & 0 & -F_2\\ \alpha & -d_2\Delta-(H_1+\mu+\beta) & -G_2\\ H_0 & H_1 & -d_3\Delta-H_2 \end{pmatrix}. \] Similar to the ODE case, the asymptotic stability of the steady state solution $(u^{\ast},v^{\ast},w^{\ast}) $ can be determined by examining the eigenvalues of the operator $\mathcal{L}$. That is the solution is asymptotically stable if all the eigenvalues of $\mathcal{L}$ have negative real parts. In order to achieve that, suppose $(\phi(x),\psi(x) ,\Upsilon(x) ) $ is an eigenfunction of $\mathcal{L}$ corresponding to an eigenvalue $\xi$. By definition, we have \[ \mathcal{L}(\phi(x) ,\psi(x) ,\Upsilon (x) ) ^{t} =\xi(\phi(x) ,\psi(x) ,\Upsilon(x) ) ^{t}, \] leading to \[ (\mathcal{L-}\xi I) \begin{pmatrix} \phi\\ \psi\\ \Upsilon \end{pmatrix} =\begin{pmatrix} 0\\ 0\\ 0 \end{pmatrix}. \] This can be rearranged to the form \[ \sum_{0\leq i\leq\infty,1\leq j\leq m_i}(A_i\mathcal{-}\xi I) \begin{pmatrix} a_{ij}\\ b_{ij}\\ c_{ij} \end{pmatrix} \Phi_{ij} =\begin{pmatrix} 0\\ 0\\ 0 \end{pmatrix}, \] where \[ \phi=\sum_{0\leq i\leq\infty,1\leq j\leq m_i}a_{ij}\Phi_{ij}, \quad \psi=\sum_{0\leq i\leq\infty,1\leq j\leq m_i}b_{ij}\Phi_{ij}, \quad \Upsilon=\sum_{0\leq i\leq\infty,1\leq j\leq m_i}c_{ij}\Phi_{ij}, \] and \[ A_i=\begin{pmatrix} -d_1\lambda_i-(H_0+\mu+\alpha) & 0 & -F_2\\ \alpha & -d_2\lambda_i-(H_1+\mu+\beta) & -G_2\\ H_0 & H_1 & -d_3\lambda_i-H_2 \end{pmatrix}. \] The stability of the steady state now reduces to examining the eigenvalues of the matrices $A_i$. The negativity of the real parts of every eigenvalue is ensured if the trace and determinant of $A_i$ and the determinant of its second additive compound $A_i^{[2] }$ are all negative. The trace of $A_i$ is given by \[ \operatorname{tr}A_i=-(d_1+d_2+d_3) \lambda_i+\operatorname{tr}J, \] which is clearly negative for all $i\geq0$\ since $\operatorname{tr}J<0$ (see \eqref{3.6}). The determinant of $A_i$ can be shown to be \begin{equation} \det A_i=-d_1d_2d_3\lambda_i^{3}-B_{A}\lambda_i^2-C_{A} \lambda_i+\det J,\label{4.5} \end{equation} where \begin{gather*} B_{A}=H_2d_1d_2+(H_0+\mu+\alpha) d_2d_3+(H_1+\mu+\beta) d_1d_3>0, \\ \begin{aligned} C_{A} & =(G_2H_1+(H_1+\mu+\beta) H_2) d_1 +((H_0+\alpha+\mu) H_2+F_2H_0) d_2\\ &\quad +(\beta+\mu+H_1) (\alpha+\mu+H_0) d_3>0. \end{aligned} \end{gather*} Clearly, $\det A_i$ for all $i\geq0$ since $\det J<0$. The last thing is to examine $\det A_i^{[2] }<0$. The matrix $A_i^{[2] }$ is the second additive compound of $A_i$ given by \begin{equation} A_i^{[2] }=\begin{pmatrix} -(d_1+d_2) \lambda_i-A & -G_2 & F_2\\ H_1 & -(d_1+d_3) \lambda_i-B & 0\\ -H_0 & \alpha & -(d_2+d_3) \lambda_i-C \end{pmatrix},\label{4.6} \end{equation} where \begin{gather*} A =H_0+H_1+2\mu+\alpha+\beta>0\\ B =H_0+H_2+\mu+\alpha>0\\ C =H_1+H_2+\mu+\beta>0. \end{gather*} Therefore, \begin{equation} \det A_i^{[2] }=-(d_2+d_3) ( d_1+d_3) (d_1+d_2) \lambda_i^{3}-B_{A^{[ 2] }}\lambda_i^2-C_{A^{[2] }}\lambda_i+\det J^{[2] },\label{4.7} \end{equation} with \begin{gather*} B_{A^{[2] }}=(B+C) d_1d_2+( A+B+C) d_2d_3+(A+B+C) d_1d_3+Ad_3^2 +Bd_2^2+Cd_1^2, \\ \begin{aligned} C_{A^{[2] }} & =(AC+BC+F_2H_0) d_1+(AB+BC+G_2H_1) d_2 \\\ &\quad + (AB+F_2H_0+AC+G_2H_1) d_3. \end{aligned} \end{gather*} We can see that $B_{A^{[2] }},C_{A^{[2] }}>0$, and since $\det J^{[2] }<0$, it follows that $\det A_i^{[2] }<0$ for all $i\geq0$. Hence, the steady state solution is locally asymptotically stable. This concludes the proof of the Proposition. \end{proof} \section{Global asymptotic stability\label{SecGlobal}} In this section, we study the global asymptotic stability of the steady state solutions for the proposed system \eqref{2.1}--\eqref{2.3}. In the ODE case, Henshaw and McCluskey \cite{Henshaw2015} established the global asymptotic stability of the unique equilibrium using an appropriate Lyapunov functional. The aim here is to show that in the presence of diffusion, every solution of the system \eqref{2.1}--\eqref{2.3} with a positive initial value that is different from the equilibrium point will converge to the equilibrium. First, let \begin{equation} L(x) =x-1-\ln(x) \label{4.10} \end{equation} for $x>0$. \begin{theorem} Let \[ V(t) =\int_{\Omega}[u^{\ast}L(\frac{u}{u^{\ast} }) +u_2^{\ast}L(\frac{v}{v^{\ast}}) +u_3^{\ast }L(\frac{w}{w^{\ast}}) ] dx. \] Then, $V(t) $ is non-negative and is strictly minimized at the unique equilibrium $(u^{\ast},v^{\ast},w^{\ast}) $, i.e. it is a valid Lyapunov functional. Hence, $(u^{\ast},v^{\ast},w^{\ast}) $ is globally asymptotically stable. \end{theorem} \begin{proof} To prove that the steady state solution $(u^{\ast},v^{\ast },w^{\ast}) $ is globally asymptotically stable, we need to establish that $V(t) $ is a Lyapunov functional. First, we differentiate $V(t) $ with respect to time \[ \frac{dV}{dt}=\int_{\Omega}\big[(1-\frac{u^{\ast}}{u}) \frac{du}{dt}+(1-\frac{v^{\ast}}{v}) \frac{dv}{dt}+( 1-\frac{w^{\ast}}{w}) \frac{dw}{dt}\big] dx. \] Substituting the time derivatives with their values from \eqref{2.1} yields \begin{align*} \frac{dV}{dt} & =\int_{\Omega}(1-\frac{u^{\ast}}{u}) [ d_1\Delta u+\Lambda_1-uf(w) -(\mu+\alpha)u] dx\\ &\quad +\int_{\Omega}(1-\frac{v^{\ast}}{v}) [d_2\Delta v+\Lambda_2+\alpha u-vg(w) -(\mu+\beta)v] dx\\ &\quad +\int_{\Omega}(1-\frac{w^{\ast}}{w}) [d_3\Delta w+\Lambda_3+uf(w) +vg(w) -(\mu +\gamma+\delta) w] dx\\ & =I+J. \end{align*} The first part is \begin{equation} I=I_1+I_2+I_3,\label{4.12} \end{equation} where \begin{gather*} I_1=\int_{\Omega}d_1(1-\frac{u^{\ast}}{u}) \Delta u\,dx, \\ I_2=\int_{\Omega}d_2(1-\frac{v^{\ast}}{v}) \Delta v\,dx, \\ I_3=\int_{\Omega}d_3(1-\frac{w^{\ast}}{w}) \Delta w\,dx. \end{gather*} The second part of the derivative is \begin{equation} \label{4.11} \begin{aligned} J & =\int_{\Omega}(1-\frac{u^{\ast}}{u}) [\Lambda _1-uf(w) -(\mu+\alpha) u] dx \\ &\quad +\int_{\Omega}(1-\frac{v^{\ast}}{v}) [\Lambda _2+\alpha u-vg(w) -(\mu+\beta) v] dx \\ &\quad +\int_{\Omega}(1-\frac{w^{\ast}}{w}) [\Lambda _3+uf(w) +vg(w) -(\mu+\gamma +\delta) w] dx, \end{aligned} \end{equation} We start by looking at $I$. Using Green's formula and assuming the Neumann boundary conditions in \eqref{2.3}, we obtain \begin{align*} I_1 & =\int_{\Omega}d_1(1-\frac{u^{\ast}}{u}) \Delta udx\\ & =-d_1\int_{\Omega}\nabla(1-\frac{u^{\ast}}{u}) \nabla udx\\ & =-d_1\int_{\Omega}\frac{u^{\ast}}{u^2}| \nabla u| ^2dx, \end{align*} \[ I_2 =\int_{\Omega}d_2(1-\frac{v^{\ast}}{v}) \Delta vdx\\ =-d_2\int_{\Omega}\frac{v^{\ast}}{v^2}| \nabla v| ^2dx, \] and \[ I_3 =\int_{\Omega}d_3(1-\frac{w^{\ast}}{w}) \Delta wdx\\ =-d_3\int_{\Omega}\frac{w^{\ast}}{w^2}| \nabla w| ^2dx. \] Therefore, by \eqref{4.12}, we have \[ I=-\int_{\Omega}\big[d_1\frac{u^{\ast}}{u^2}| \nabla u| ^2+d_2\frac{v^{\ast}}{v^2}| \nabla v| ^2+d_3\frac{w^{\ast}}{w^2}| \nabla w| ^2\big] dx<0. \] The second part of the derivative $J$ can be simplified by replacing $\Lambda_1$, $\Lambda_2$, and $(\mu+\gamma+\delta) $ with their values from \eqref{3.3} and rearranging to the form \begin{align*} J & =\int_{\Omega}(1-\frac{u^{\ast}}{u}) [u^{\ast }f(w^{\ast}) +(\mu+\alpha) u^{\ast}-uf( w) -(\mu+\alpha) u] dx\\ &\quad +\int_{\Omega}(1-\frac{v^{\ast}}{v}) [v^{\ast}g( w^{\ast}) +(\mu+\beta) v^{\ast}-\alpha u^{\ast}+\alpha u-vg(w) -(\mu+\beta) v] dx\\ &\quad +\int_{\Omega}(1-\frac{w^{\ast}}{w}) \Big[\Lambda _3+uf(w) +vg(w) -\frac{\Lambda_3+u^{\ast }f(w^{\ast}) +v^{\ast}g(w^{\ast}) }{w^{\ast}}w\Big] dx\\ & =\int_{\Omega}(1-\frac{u^{\ast}}{u}) \Big[u^{\ast}f( w^{\ast}) \Big(1-\frac{uf(w) }{u^{\ast}f( w^{\ast}) }\Big) +(\mu+\alpha) u^{\ast}( 1-\frac{u}{u^{\ast}}) \Big] dx\\ &\quad +\int_{\Omega}(1-\frac{v^{\ast}}{v}) \Big[v^{\ast}g( w^{\ast}) \Big(1-\frac{vg(w) }{v^{\ast}g( w^{\ast}) }\Big) +(\mu+\beta) v^{\ast}( 1-\frac{v}{v^{\ast}}) +\alpha u^{\ast}(\frac{u}{u^{\ast} }-1) \Big] dx\\ &\quad +\int_{\Omega}(1-\frac{w^{\ast}}{w}) \Big[\Lambda _3(1-\frac{w}{w^{\ast}}) +u^{\ast}f(w^{\ast}) \Big(\frac{uf(w) }{u^{\ast}f(w^{\ast}) } -\frac{w}{w^{\ast}}\Big) \\ &\quad +v^{\ast}g(w^{\ast}) \Big(\frac{vg(w) }{v^{\ast}g(w^{\ast}) } -\frac{w}{w^{\ast}}\Big) \Big] dx. \end{align*} Further simplification yields \begin{align} J & =\int_{\Omega}(\mu+\alpha) u^{\ast}(1-\frac {u}{u^{\ast}}) (1-\frac{u^{\ast}}{u}) +\Lambda_3( 1-\frac{w}{w^{\ast}}) (1-\frac{w^{\ast}}{w}) \\ &\quad +u^{\ast}f(w^{\ast}) \Big[\Big(\frac{uf(w) }{u^{\ast}f(w^{\ast}) }-\frac{w}{w^{\ast}}\Big) ( 1-\frac{w^{\ast}}{w}) +(1-\frac{u^{\ast}}{u}) \Big( 1-\frac{uf(w) }{u^{\ast}f(w^{\ast}) }\Big)\Big] \\ &\quad +v^{\ast}g(w^{\ast}) \Big[\Big(1-\frac{vg( w) }{v^{\ast}g(w^{\ast}) }\Big) (1-\frac {v^{\ast}}{v}) +\Big(\frac{vg(w) }{v^{\ast}g( w^{\ast}) }-\frac{w}{w^{\ast}}\Big) (1-\frac{w^{\ast}} {w}) \Big] \\ &\quad +(\mu+\beta) v^{\ast}(1-\frac{v}{v^{\ast}}) (1-\frac{v^{\ast}}{v}) +\alpha u^{\ast}(\frac{u}{u^{\ast }}-1) (1-\frac{v^{\ast}}{v}) dx.\label{4.13} \end{align} Now, to show that $J$ is negative, we observe the following equalities \begin{gather*} L(\frac{u}{u^{\ast}}) +L(\frac{u^{\ast}}{u}) =-(1-\frac{u}{u^{\ast}}) (1-\frac{u^{\ast}}{u}) ,\\ \begin{aligned} &-L(\frac{u^{\ast}}{u}) +L\big(\frac{f(w) }{f(w^{\ast}) }\big) -L(\frac{w}{w^{\ast}}) -L\big(\frac{uf(w) w^{\ast}}{u^{\ast}f(w^{\ast}) w}\big) \\ &= \Big(\frac{uf(w) }{u^{\ast}f(w^{\ast}) }-\frac{w}{w^{\ast}}\Big) (1-\frac{w^{\ast}}{w}) +(1-\frac{u^{\ast}}{u}) \big(1-\frac{uf(w) }{u^{\ast}f(w^{\ast})}\big) , \end{aligned} \\ \begin{aligned} &-L\big(\frac{w}{w^{\ast}}\big) -L(\frac{vg(w) w^{\ast}}{v^{\ast}g(w^{\ast}) w}) -L(\frac {v^{\ast}}{v}) +L\big(\frac{g(w) }{g(w^{\ast }) }\big) \\ &= \Big(1-\frac {vg(w) }{v^{\ast}g(w^{\ast}) }\Big) (1-\frac{v^{\ast}}{v}) +\Big(\frac{vg(w) }{v^{\ast }g(w^{\ast}) }-\frac{w}{w^{\ast}}\Big) (1-\frac {w^{\ast}}{w}) , \end{aligned} \\ L(\frac{u}{u^{\ast}}) -L(\frac{uv^{\ast}}{u^{\ast} v}) +L(\frac{v^{\ast}}{v}) =(\frac{u}{u^{\ast} }-1) (1-\frac{v^{\ast}}{v}) . \end{gather*} Substituting these in \eqref{4.13} leads to \begin{align*} J & =-\int_{\Omega}(\mu+\alpha) u^{\ast}\Big[L( \frac{u}{u^{\ast}}) +L(\frac{u^{\ast}}{u}) \Big] dx-\int_{\Omega}\Lambda_3\frac{(w-w^{\ast}) ^2}{ww^{\ast} }dx\\ &\quad -\int_{\Omega}u^{\ast}f(w^{\ast}) \Big[L( \frac{u^{\ast}}{u}) -L(\frac{f(w) }{f( w^{\ast}) }) +L(\frac{w}{w^{\ast}}) +L\Big( \frac{uf(w) w^{\ast}}{u^{\ast}f(w^{\ast}) w}\Big) \Big] dx\\ &\quad -\int_{\Omega}v^{\ast}g(w^{\ast}) \Big[L(\frac {w}{w^{\ast}}) +L\Big(\frac{vg(w) w^{\ast}}{v^{\ast }g(w^{\ast}) w}\Big) +L(\frac{v^{\ast}}{v}) -L\Big(\frac{g(w) }{g(w^{\ast}) }\Big) \Big] dx\\ &\quad -(\mu+\beta) \int_{\Omega}v^{\ast}\big[L(\frac {v}{v^{\ast}}) +L(\frac{v^{\ast}}{v}) \big] dx\\ &\quad +\alpha\int_{\Omega}u^{\ast}\big[L(\frac{u}{u^{\ast}}) -L(\frac{uv^{\ast}}{u^{\ast}v}) +L(\frac{v^{\ast}} {v}) \big] dx \end{align*} Now, using Proposition \ref{HenshawProp} and simplification similar to \cite{Henshaw2015} yields the inequality \begin{align*} J & \leq-\int_{\Omega}(\mu+\alpha) u^{\ast}\Big[L( \frac{u}{u^{\ast}}) +L(\frac{u^{\ast}}{u}) \Big]dx -\int_{\Omega}\Lambda_3\frac{(w-w^{\ast}) ^2}{ww^{\ast}}dx\\ &\quad -\int_{\Omega}u^{\ast}f(w^{\ast}) \Big[L( \frac{u^{\ast}}{u}) +L\Big(\frac{uf(w) w^{\ast} }{u^{\ast}f(w^{\ast}) w}\Big) \Big] dx\\ &\quad -\int_{\Omega}v^{\ast}g(w^{\ast}) \Big[L\Big( \frac{vg(w) w^{\ast}}{v^{\ast}g(w^{\ast}) w}\Big) \Big] dx\\ &\quad -(\mu+\beta) \int_{\Omega}v^{\ast}L(\frac{v}{v^{\ast} }) dx-\alpha\int_{\Omega}u^{\ast}L(\frac{uv^{\ast}}{u^{\ast} v}) dx. \end{align*} It is clear that $J\leq0$, which leads to $\frac{dV}{dt}\leq0$; $\frac{dV}{dt}=0$ only at the steady state $(u^{\ast},v^{\ast},w^{\ast})$. Therefore, by Lyapunov's direct method, the steady state solution $(u^{\ast},v^{\ast},w^{\ast}) $ is globally asymptotically stable. \end{proof} \section{Numerical examples} In this section, we present two numerical examples that illustrate and confirm the findings of this study. The parameters utilized for the examples are stated in Table \ref{Tab1}. \begin{table}[htb] \caption{Simulation parameters for the stated numerical examples.} \label{Tab1} \begin{center} \begin{tabular}{|c|c|c|}\hline Parameters & Example 1 & Example 2\\ \hline $f(x) $ & $\frac{x}{x+1}$ & $\frac{x}{x+1}$\\ \hline $g(x) $ & $\frac{x}{2x+2}$ & $\frac{x}{2x+2}$\\ \hline $\Lambda_1$ & $1$ & $1.5$\\ \hline $\Lambda_2$ & $1.2$ & $1.2$\\ \hline $\Lambda_3$ & $0.95$ & $0.95$\\ \hline $\Lambda_3$ & $0.1$ & $0.1$\\ \hline $\mu$ & $0.2$ & $0.2$\\ \hline $\alpha$ & $0.3$ & $0.3$\\ \hline $\beta$ & $0.45$ & $0.45$\\ \hline $\gamma$ & $0.02$ & $0.005$\\ \hline $\delta$ & $0.1$ & $0.1$\\ \hline $d_1$ & $1$ & $100$\\ \hline $d_2$ & $1.5$ & $600$\\ \hline $d_3$ & $1.3$ & $1000$\\ \hline $d_{4}$ & $1$ & $500$\\ \hline $u_0$ & $50\operatorname{sinc}[0.2(x^2+y^2) ] $ & $50\operatorname{sinc}[0.2(x^2+y^2) ] $\\ \hline $v_0$ & $15\operatorname{sinc}[0.8(x^2+y^2) ] $ & $15\operatorname{sinc}[0.8(x^2+y^2) ] $\\ \hline $w_0$ & $10\operatorname{sinc}[0.8(x^2+y^2) ] $ & $10\operatorname{sinc}[0.8(x^2+y^2) ] $\\ \hline $r_0$ & $0.1\operatorname{sinc}[0.8(x^2+y^2) ] $ & $0.1\operatorname{sinc}[0.8(x^2+y^2) ] $\\ \hline \end{tabular} \end{center} \end{table} \subsection{First Example} We use the parameters from the first column of Table \ref{Tab1}. Solving the system of equations \eqref{3.3} numerically yields the equilibrium solution $(u^{\ast},v^{\ast},w^{\ast}) =(0.7296,1.3073,6.7323,6.8076) $. In the ODE case, the initial data in the ODE case is simply $(50,15,10) $ and the equilibrium can be clearly seen to be asymptotically stable as seen in Figure \ref{Ex1_0D}. Figure \ref{Ex1_2D} shows the solutions in the two-dimensional diffusion case and the steady state solution is again asymptotically stable. We see that over-time, the solutions tend to the steady state $(u^{\ast},v^{\ast},w^{\ast},r^{\ast}) $ and become close to uniformly distributed in space. \begin{figure}[htb] \begin{center} \includegraphics[width = 4in]{fig2} % Pap9_Ex1_0D.pdf \end{center} \caption{Solutions of the proposed system \eqref{2.1} in the ODE case using parameters from the first column of Table \ref{Tab1}.} \label{Ex1_0D} \end{figure} \begin{figure}[ptb] \begin{center} \includegraphics[width = \textwidth]{fig3} % Pap9_Ex1_2D.pdf \end{center} \caption{Solutions of the proposed system \eqref{2.1} in the two-dimensional PDE diffusion case using parameters from the first column of Table \ref{Tab1}. The snapshots from top to bottom are taken at times $t=0$, $t=1$, and $t=10$, respectively.} \label{Ex1_2D} \end{figure} \subsection{Second Example} The aim of this example is to show that high diffusivity constants do not affect the asymptotic stability of the solutions. The system parameters are shown in the second column of Table \ref{Tab1}. Figures \ref{Ex2_0D} and \ref{Ex2_2D} show the solutions in the ODE and two-dimensional cases, respectively. Observe that due to the high diffusivities, the solutions reach the equilibrium $(u^{\ast},v^{\ast},w^{\ast},r^{\ast}) =(1.0772,1.3895,8.2997,7.7761) $ in a shorter time and that the solutions remain stable in both scenarios. \begin{figure}[htb] \begin{center} \includegraphics[width = 4in]{fig4} % Pap9_Ex2_0D.pdf \end{center} \caption{Solutions of the proposed system \eqref{2.1} in the ODE case using parameters from the second column of Table \ref{Tab1}.} \label{Ex2_0D} \end{figure} \begin{figure}[htb] \begin{center} \includegraphics[width = \textwidth]{fig5} % Pap9_Ex2_2D.pdf \end{center} \caption{Solutions of the proposed system \eqref{2.1} in the two-dimensional PDE diffusion case using parameters from the second column of Table \ref{Tab1}. The snapshots from top to bottom are taken at times $t=0$, $t=1$, and $t=10$, respectively.} \label{Ex2_2D} \end{figure} \section{Appendix} \begin{lemma}[\cite{McCluskey2004}] \label{PropositionM} Let $M$ be a $3\times3$ real matrix. If $\operatorname{tr}(M) $, $\det(M) $, and $\det(M^{[2] }) $ are all negative, then all of the eigenvalues of $M$ have negative real parts, where (see \cite{Muldowney1990}) \begin{equation} M=\begin{pmatrix} a_{11} & a_{12} & a_{13}\\ a_{21} & a_{22} & a_{23}\\ a_{31} & a_{32} & a_{33} \end{pmatrix}, \quad M^{[2] }=\begin{pmatrix} a_{11}+a_{22} & a_{23} & -a_{13}\\ a_{32} & a_{11}+a_{33} & -a_{12}\\ -a_{31} & a_{21} & a_{22}+a_{33} \end{pmatrix}.\label{App2} \end{equation} \end{lemma} \begin{proof} Let $\lambda_{j}$, $j=1,2,3$ be the eigenvalues of $M$ with $\Re(\lambda_1) \leq\Re(\lambda_2) \leq\Re(\lambda_3) $. It follows from $\det(M) <0$ that $\lambda_1\lambda_2\lambda_3<0$. Thus, either $\Re(\lambda_{j}) <0$ for $j=1,2,3$ (which would prove the lemma) or $\Re(\lambda_1) <0\leq\Re(\lambda_2) \leq\Re(\lambda_3) $. Suppose that the second set of inequalities holds. Since $\operatorname{tr}(M) <0$, it follows that $\lambda_1 +\lambda_2+\lambda_3<0$, which implies that $\Re(\lambda _1+\lambda_2) <0$ and $\Re(\lambda_1+\lambda_3)<0$. The eigenvalues of $M^{[2]}$ are $\lambda_i+\lambda_{j}$, $1\leq i0$. \end{proposition} \begin{proof} Let $w>0$. Since $f(w) $ is continuous on $[0,w]$ and differentiable on $(0,w) $, the mean value theorem implies that there exists $c\in(0,w) $ such that $f'(c) \leq\frac{f(w) -f(0) }{w-0}$. By (H1), we have $f'(c) =\frac{f(w) }{w}$. From \textbf{(H3)}, $f'$ is monotone decreasing. Thus, $f'(w) \leq f'(c) =\frac{f(w) }{w}$. The same can be said about $g$. \end{proof} \begin{proposition}[\cite{Henshaw2015}] \label{HenshawProp} Suppose the incidence functions $f$ and $g$ satisfy the criteria in (H1)--(H4). It follows that if $w>0$, then \begin{gather} L\Big(\frac{f(w) }{f(w^{\ast}) }\Big) \leq L\big(\frac{w}{w^{\ast}}\big), \label{App3} \\ L\Big(\frac{g(w) }{g(w^{\ast}) }\Big) \leq L\big(\frac{w}{w^{\ast}}\Big). \label{App4} \end{gather} \end{proposition} \begin{proof} In this proof, we will only establish the property \eqref{App3}. However, the same can be said about \eqref{App4}. Let $w\geq w^{\ast}$ and $m(w) =\frac{f(w) }{w}$. It follows that \[ m'(w) =\frac{f'(w)w-f(w) }{w^2}\leq\frac{f(w) -f(w)}{w^2}=0. \] Therefore, we conclude that $m$ is decreasing, which leads to $m(w) \leq m(w^{\ast})$, i.e., \[ \frac{f(w) }{w}\leq\frac{f(w^{\ast}) }{w^{\ast}}, \] and so \[ \frac{f(w) }{f(w^{\ast}) }\leq\frac{w}{w^{\ast}}. \] Since $f$ is increasing, we have \[ 1\leq\frac{f(w) }{f(w^{\ast}) }\leq\frac{w}{w^{\ast}}. \] Note that $L(x) =1-\frac{1}{x}$. Hence, $L$ is increasing for $x>1$, and \[ L\Big(\frac{f(w) }{f(w^{\ast}) }\Big) \leq L(\frac{w}{w^{\ast}}) . \] \end{proof} \subsection*{Acknowledgements} The authors would like to thank Prof. M. Kirane and the anonymous referee for their most insightful suggestions that improved the quality of this article. \begin{thebibliography}{99} \bibitem{Abdelmalek2013} S. Abdelmalek, M. Kirane, A. Youkana; \emph{A Lyapunov functional for a triangular reaction--diffusion system with nonlinearities of exponential growth}, Math. Methods in the Appl. Sci., Vol. 36, 80--85, (2013). \bibitem{Abdelmalek2015} S. Abdelmalek, S. Bendoukha; \emph{Global Asymptotic Stability for an SEI Reaction-Diffusion Model of Infectious Diseases with Immigration}, submitted. \bibitem{Busenberg1993} S. Busenberg, K. Cooke; \emph{Vertically Transmitted Diseases, Methods and Dynamics}. Biomathematics, Vol. 23. Springer Verlag: New York, (1993). \bibitem{Casten1977} R. G. Casten, C. J. Holland; \emph{Stability properties of solutions to systems of reaction-diffusion equations}, SIAM J. Appl. Math. 33, 353--364 (1977). \bibitem{Duan2014} X. Duan, S. Yuan, X. Li; \emph{Global stability of an SVIR model with age of vaccination}, Applied Math. and Comp., Vol. 226, 528--540, (2014). \bibitem{Haddad2010} W. M. Haddad, V. Chellaboina, Q. Hui; \emph{Nonnegative and Compartmental Dynamical Systems}, Princeton University Press, Princeton, New Jersey, (2010). \bibitem{Henshaw2015} S. Henshaw, C. C. McCluskey; \emph{Global stability of a vaccination model with immigration}, Electron. J. of Diff. Eqs., Vol. 2015, No. 92, pp. 1--10, (2015). \bibitem{Kuniya2013} T. Kuniya; \emph{Global stability of a multi-group SVIR epidemic model}, Nonlinear Anal.: Real World Apps., Vol. 14, 1135--1143, (2013). \bibitem{Liu2008} X. Liu, Y. Takeuchib, S. Iwami; \emph{SVIR epidemic models with vaccination strategies}, J. Theor. Bio., Vol. 253, 1--11, (2008). \bibitem{McCluskey2004} C. C. McCluskey, P. van den Driessche; \emph{Global analysis of two tuberculosis models}. J. Dynam. Diff. Eqs., Vol. 16(1): 139-166, (2004). \bibitem{Muldowney1990} J. S. Muldowney; \emph{Compound matrices and ordinary differential equations}, Rocky Mountain J. Math., Vol. 20: 857--872, (1990). \bibitem{Quitner2007} P. Quittner, Ph. Souplet; \emph{Superlinear Parabolic Problems. Blow-up, Global Existence and Steady States}, Birkhauser Verlag AG, Basel. Boston., Berlin, (2007). \bibitem{Sigdel2014} Ram P. Sigdel, C. Connell McCluskey; \emph{Global stability for an SEI model of infectious disease with immigration}, Applied Math. and Comp., Vol. 243: 684--689, (2014). \end{thebibliography} \end{document}