Geographical network model for COVID-19 spread among dynamic epidemic regions

Roman Zúñiga Macías; Humberto Gutiérrez-Pulido; Edgar Alejandro Guerrero Arroyo; Abel Palafox González; Roman Zúñiga Macías; Humberto Gutiérrez-Pulido; Edgar Alejandro Guerrero Arroyo; Abel Palafox González

doi:10.3934/mbe.2022196

Mathematical Biosciences and Engineering

2022, Volume 19, Issue 4: 4237-4259. doi: 10.3934/mbe.2022196

Previous Article Next Article

Research article Special Issues

Geographical network model for COVID-19 spread among dynamic epidemic regions

Universidad de Guadalajara, CUCEI, Blvd. Marcelino García Barragán 1421, 44430, Guadalajara, Jal., México

Received: 11 December 2021 Revised: 25 January 2022 Accepted: 08 February 2022 Published: 23 February 2022

Pandemic due to SARS-CoV-2 (COVID-19) has affected to world in several aspects: high number of confirmed cases, high number of deaths, low economic growth, among others. Understanding of spatio-temporal dynamics of the virus is helpful and necessary for decision making, for instance to decide where, whether and how, non-pharmaceutical intervention policies are to be applied. This point has not been properly addressed in literature since typical strategies do not consider marked differences on the epidemic spread across country or large territory. Those strategies assume similarities and apply similar interventions instead. This work is focused on posing a methodology where spatio-temporal epidemic dynamics is captured by means of dividing a territory in time-varying epidemic regions, according to geographical closeness and infection level. In addition, a novel Lagrangian-SEIR-based model is posed for describing the dynamic within and between those regions. The capabilities of this methodology for identifying local outbreaks and reproducing the epidemic curve are discussed for the case of COVID-19 epidemic in Jalisco state (Mexico). The contagions from July 31, 2020 to March 31, 2021 are analyzed, with monthly adjustments, and the estimates obtained at the level of the epidemic regions present satisfactory results since Relative Root Mean Squared Error RRMSE is below 15% in most of regions, and at the level of the whole state outstanding with RRMSE below 5%.

Keywords:

Citation: Roman Zúñiga Macías, Humberto Gutiérrez-Pulido, Edgar Alejandro Guerrero Arroyo, Abel Palafox González. Geographical network model for COVID-19 spread among dynamic epidemic regions[J]. Mathematical Biosciences and Engineering, 2022, 19(4): 4237-4259. doi: 10.3934/mbe.2022196

Related Papers:

[1]	Haoyang Ji, Wenxiu Ma . Phase transition for piecewise linear fibonacci bimodal map. AIMS Mathematics, 2023, 8(4): 8403-8430. doi: 10.3934/math.2023424
[2]	Skander Hachicha, Najmeddine Attia . Probabilistic approaches to exploring Binet's type formula for the Tribonacci sequence. AIMS Mathematics, 2025, 10(5): 11957-11975. doi: 10.3934/math.2025540
[3]	Waleed Mohamed Abd-Elhameed, Omar Mazen Alqubori . New expressions for certain polynomials combining Fibonacci and Lucas polynomials. AIMS Mathematics, 2025, 10(2): 2930-2957. doi: 10.3934/math.2025136
[4]	Ümit Tokeşer, Tuğba Mert, Yakup Dündar . Some properties and Vajda theorems of split dual Fibonacci and split dual Lucas octonions. AIMS Mathematics, 2022, 7(5): 8645-8653. doi: 10.3934/math.2022483
[5]	Moussa Benoumhani . Restricted partitions and convex topologies. AIMS Mathematics, 2025, 10(4): 10187-10203. doi: 10.3934/math.2025464
[6]	Taja Yaying, S. A. Mohiuddine, Jabr Aljedani . Exploring the $q$ -analogue of Fibonacci sequence spaces associated with $c$ and $c_0$ . AIMS Mathematics, 2025, 10(1): 634-653. doi: 10.3934/math.2025028
[7]	Man Chen, Huaifeng Chen . On ideal matrices whose entries are the generalized $k-$ Horadam numbers. AIMS Mathematics, 2025, 10(2): 1981-1997. doi: 10.3934/math.2025093
[8]	Rajiniganth Pandurangan, Sabri T. M. Thabet, Imed Kedim, Miguel Vivas-Cortez . On the Generalized $\overline{\theta({\tt{t}})}$ -Fibonacci sequences and its bifurcation analysis. AIMS Mathematics, 2025, 10(1): 972-987. doi: 10.3934/math.2025046
[9]	Kritkhajohn Onphaeng, Prapanpong Pongsriiam . Exact divisibility by powers of the integers in the Lucas sequence of the first kind. AIMS Mathematics, 2020, 5(6): 6739-6748. doi: 10.3934/math.2020433
[10]	Faik Babadağ, Ali Atasoy . On hyper-dual vectors and angles with Pell, Pell-Lucas numbers. AIMS Mathematics, 2024, 9(11): 30655-30666. doi: 10.3934/math.20241480

Abstract

1. Introduction

Competition is the fundamental relation of ecological systems in which species of one kind have to fight for limited resources in some particular habitat. These sorts of competitions become severe when the potential biological needs exceed the resources. The predator-prey populations are reasonably the building blocks for dynamic ecosystems. In ecology, among different kinds of relationships between species, the most important relationship is the predator-prey relationship. Mathematical models governed by differential equations are more appropriate for modeling the interaction in which populations are overlapped. Lotka and Volterra first put forward the fundamental predator-prey model in the early twentieth century. This basic model comprises two first-order nonlinear ordinary differential equations. Since then, researchers have presented many models covering several issues regarding the complex natural relationship. Consider

$\frac{dU}{dt} = {\alpha }_{1}U-{\alpha }_{2}UV$

(1)

$\frac{dV}{dt} = {\alpha }_{3}UV-{\alpha }_{4}V$

(2)

where $U$ and $V$ respectively indicate the populations of prey and predator; ${\alpha }_{1}, {\alpha }_{2}, {\alpha }_{3}$ and ${\alpha }_{4}$ are the positive real parameters. Holling introduced the density-dependent response ^[1], after which several researchers came forward to contribute their valuable research using density-dependent responses. The main concern in the population dynamical models is the stability of fixed points. Many mathematical models have been presented and fixed-point stability has been studied ^[2]. The stability of the predator-prey model has been brought into the discussion by many researchers ^[3,4,5]. The other targeted area is the positivity of the solution of the prey-predator relation ^[6]. During the last few decades, many species have faced extinction due to limited resources, over-exploitation, pollution and predation exercise. Mostly, the extinction of some species occurs due to transgression of the environment or ecological structure. To impede species from extinction, external supporting trends are the refuge and restriction of the population to a specific area ^[7,8,9]. Considering the issue of extinction, many researchers have formulated models to deal with the problem and discussed the factors in detail. A mathematical model of two prey and one predator was considered by Takeuchi and Adachi ^[10]. Equilibrium points were investigated, and a mathematical analysis of equilibrium points was performed. Another important issue is bifurcation, which is the critical behavior of fixed points. Bifurcation analysis for interacting predator-prey models has been provided by many investigators ^{[10,11,12,13]}. In ^[14], the authors investigated the influence of nonlocal competition among the species. The study was further enhanced to investigate the effects of such competition on the stability of equilibrium points. Due to the overuse of chemicals for the higher production of crops and other factors, infectious diseases are spreading among the human population and animals. Researchers have published several articles on the problem of population density dynamics. Whenever a disease spreads in an environment, it seriously affects the population dynamics and disturbs the whole ecosystem. A model was presented by Haque ^[15] to discuss the impact of the diseased predator on the population dynamics of predator and prey. He showed that the absence of prey has a moderately strong impact on the predator population when there exists an infection in the predator. The mathematical model for infected predator and prey, as presented by the author, is as follows:

$\frac{dU}{dt} = aU(1-\frac{U}{M})-c{V}_{1}U-cq{V}_{2}$

(3)

$\frac{dF}{dt} = rF(1-\frac{F}{K})+ec{V}_{1}U+ecq{V}_{2}U$

(4)

$\frac{d{V}_{1}}{dt} = (b-\frac{\theta rF}{K})F-(d+\frac{(1-\theta )rF}{k}{)V}_{1}-\frac{\beta {V}_{1}{V}_{2}}{F}+\gamma {V}_{2}+ec{V}_{1}U$

(5)

$\frac{d{V}_{2}}{dt} = F{V}_{2}(1-\frac{{V}_{2}}{K})+\frac{\beta {V}_{1}{V}_{2}}{F}-\gamma {V}_{2}-(d+\frac{\left(1-\theta \right)rF}{k}){V}_{2}+ecq{V}_{2}U$

(6)

Here, $P$ and $F = {V}_{1}+{V}_{2}$ respectively represent the prey and predator populations (susceptible and infected). Fang and Wang discussed the interaction between two species when the predator and prey use common food resources. Moreover, the predator is also consuming the prey species. The model suggested by them is as follows ^[16]:

$\frac{dU}{dt} = U(\frac{af}{aU+bV}+{\alpha }_{1}V-m)$

(7)

$\frac{dV}{dt} = V(\frac{bf}{aU+bV}-{\alpha }_{1}U-n)$

(8)

where $U$ and $V$ respectively represent the population of predator and prey, $f$ is the common resource, $a$ and $b$ are the consumption rates of resources by predator and prey, respectively, and $m$ and $n$ are the natural death rates of predator and prey, respectively. Fang and his co-author discussed the global stability of the above dynamical system with the additional effects of diffusion. A detailed discussion was provided on the effects of diffusion by the authors. They discussed the instability of equilibrium points due to the involvement of diffusion. Moreover, to carry out the stability of the system the Lyapunov method was employed. Inseba et al. studied the nonlinear interaction of multi-species, taking into account diffusion and prey taxis. The authors also investigated the linearized stability of the system ^[17]. Chen and Yu discussed the multi-species predator-prey diffusive system. They showed the pattern formation of the bifurcating system in which the conversion rate was taken as a bifurcation parameter ^[18]. Djilali discussed the effect of the prey's social behavior on the diffusive predator-prey system. He also discussed the bifurcation and stability analysis ^[19]. In ^[20], Ferreira et al. studied the stability of the cross-diffusive system for the three-species interacting model. The equilibrium points were examined to perform a local stability analysis. Kant and Kumar studied the interacting model of species in which predator and prey are diseased. They showed that, for the basic reproductive number ${R}_{0} > 1$ , the disease spread ^[21]. Owolabi formulated a numerical scheme to deal with the fractional diffusion system of interacting species. He applied the Caputo fractional operator rather than the first-order time derivative ^[22]. Song et al. presented a qualitative approach to study the diffusive predator-prey system. They also discussed the system's instability when cross-diffusion was introduced ^[23]. Zhang et al. discussed the influence of diffusion on predator-prey predators' hyperbolic mortality. The turning instability region was also found by linear stability ^[24]. Ghosh et al. investigated the stability and performed a bifurcation analysis of Bozykin's prey-predator model. The authors also discussed some properties of the solution, like uniqueness and boundedness ^[25]. In ^[26], Owolabi formulated a numerical scheme for the computation of fractional (time) derivatives by using the finite difference and Fourier algorithm. Many researchers have discussed the effects of diseased species on the evolution of population dynamics when prey, predator or both have some infection. The authors have investigated the models by using core characteristics like the existence and positivity of equilibria and local and global stability. Some authors have also discussed the bifurcation parameters in the model. Such studies have been supported by simulating the models via numerical schemes ^{[27,28,29,30,31]}.

Recent years have seen a rise in the application of fractional differential calculus to solve critical and practical problems. In ^[32], a prey-predator model with three species is examined within the framework of a fractional operator. Two species, both of which grow logistically, are represented in the model. A competitor is considered to fall behind the third group due to their Holling type II functional response. To represent the interplay between tumor growth and the immune system, the authors of ^[33] provided some new approximate solutions to a computational formulation by using numerous fractional and fractal operators.

In ^[34], the researcher looked at a computer model to investigate the spread of a viral infectious disease that is more common in children, i.e., hand, foot, and mouth syndrome. A contribution of ^[35] is the development of a new, time-saving method for obtaining exact fractional solutions to local fractional equations. Equation of Gardner on Cantor sets defined by efficient numerical methods are studied.

In ^[36], the authors study the stability of a disease model with susceptible prey and infected predators around an internal steady state. The Mittag-Leffler kernels from the Liouville-Caputo idea, which are used to calculate fractional derivatives, have been considered for this purpose. In ^[37], the researchers looked at one realistic application from the current state of edge detection research. To achieve this goal, they first suggested two overarching structures that can be used to create brand-new fractional masks. Then, the Atangana-Baleanu operator, a fractional integral, was used to assess the roles of the various parts in these two architectures.

In ^[38], the authors suggested a predator-prey model with a Michaelis–Menten functional response and split prey into susceptible and infected subpopulations. The researchers presented an eco-epidemiological model of an infected predator-prey system ^[39]. Incorporating prey refuge ensures that a portion of the diseased prey is available for ingestion by the predator. The authors of ^[40] suggested a predator-prey model in which the prey population is affected by a disease. Here, healthy prey species exhibited defense mechanisms in response to a predator attack. A Leslie-Gower predator-prey model incorporating disease in the predator has been developed ^[41]. The genetic repercussions of the Lesli-Gower model allow it to be considered an evolutionary version of the Lotka-Volterra model. In ^[42], the authors created an eco-epidemic model with two prey populations and one predator population, where only the first prey population is affected by an infectious disease. More literature on predator-prey models can be seen in ^{[43,44,45,46,47]}.

In the present research work, we formulate a multi-species eco-epidemiological model. Two major aspects of ecology are discussed, i.e., food resources and infection in the species. These factors bring drastic changes in the population density of the species under consideration. It is shown that this kind of interaction leads to system stability even under self-diffusion. It is also proved that the system becomes unstable when cross-diffusion is considered. A two-stage numerical scheme is also applied in the present work, and numerical simulations support the theoretical findings. The stability of the proposed scheme is presented by using the von Neumann stability criterion. It is proved that the scheme is conditionally stable.

2. Model formulation

The ecosystem refers to a natural complex structure in which interacting species affect the population density of one another on a large scale. At the same time, the interacting species are in danger of infectious diseases in the real environment. Numerous mathematical models present the interaction among species with infection in some species. The infections are so critical that they disturb the whole ecosystem. The population density of all species is affected badly, as they depend on one another directly or indirectly. This problem has motivated several investigators to develop their model to study population dynamics effectively. Here, we propose a model to study the predator-prey system in which both species use the same food resources and there is a fatal infectious disease in the prey population. Suppose that $U, V$ respectively represent the predator and prey populations at any time t. Here, we assume that $V = {V}_{S}+{V}_{I}$ , where ${V}_{S}\wedge {V}_{I}$ respectively denote the susceptible and infected prey.

$\frac{dU}{dt} = U(\frac{af}{aU+b{V}_{S}+c{V}_{I}}+{\alpha }_{1}{V}_{S}+{\alpha }_{2}{V}_{I}-m)$

(9)

$\frac{d{V}_{S}}{dt} = {V}_{S}(\frac{bf}{aU+b{V}_{S}+c{V}_{I}}-\gamma {V}_{I}-{\alpha }_{1}U-{n}_{1}$

(10)

$\frac{d{V}_{I}}{dt} = {V}_{I}(\frac{cf}{aU+b{V}_{S}+c{V}_{I}}+\gamma {V}_{S}-{\alpha }_{2}U-{n}_{2})$

(11)

where

${\phi }_{0} = \left\{U\right(0) = {U}_{0}\ge 0, {V}_{S}(0) = {V}_{S, 0}\ge 0, {V}_{I}(0) = {V}_{I, 0}\ge 0\}$

(12)

is the set of initial conditions. Here, m, ${n}_{1}$ and ${n}_{2}$ are the death rates of the predator and preys, and ${n}_{2} = {n}_{1}+\delta$ displays the sum for the related natural and infectious disease.

3. Stability for the predator-prey system of equations

The present section deals with the computation and stability of the equilibria of the system described by (9)–(11). By solving the system when the rate of change of population becomes zero, it is easy to obtain the equilibria of the system, which are seven in number for the present case. The equilibria of the system are described as follows:

${E}_{1} = (0, \frac{f}{{n}_{1}}, \ 0), {E}_{2} = (0, \ \ \ \ 0, \frac{f}{{n}_{2}}), {E}_{3} = (\frac{f}{m}, \ 0, \ 0)$

${E}_{4} = ({U}_{4}, {V}_{S, 4}, {V}_{I, 4}), {E}_{5} = ({U}_{5}, {V}_{S, 5}, {V}_{I, 5}), {E}_{6} = ({U}_{6}, {V}_{S, 6}, {V}_{I, 6}), {E}_{7} = ({U}^{\mathrm{*}}, {V}_{S}^{\mathrm{*}}, {V}_{I}^{\mathrm{*}})$

(13)

where

${U}_{4} = 0, {V}_{S, 4} = -\frac{-b{{n}_{2}}^{2}+cf\gamma +c{n}_{1}{n}_{2}}{\gamma \left(b{n}_{2}-c{n}_{1}\right)}, {V}_{I, 4} = \frac{bf\gamma -b{n}_{1}{n}_{2}+c{{n}_{1}}^{2}}{\gamma \left(b{n}_{2}-c{n}_{1}\right)}$

${U}_{5} = -\frac{a{{n}_{1}}^{2}+bf{\alpha }_{1}-bm{n}_{1}}{{\alpha }_{1}\left(a{n}_{1}-bm\right)}, {V}_{S, 5} = \frac{af{\alpha }_{1}+am{n}_{1}-b{m}^{2}}{{\alpha }_{1}\left(a{n}_{1}-bm\right)}, {V}_{I, 5} = 0$

${U}_{6} = -\frac{a{{n}_{2}}^{2}+cf{\alpha }_{2}-cm{n}_{2}}{{\alpha }_{2}\left(a{n}_{2}-cm\right)}, {V}_{S, 6} = 0, {V}_{I, 6} = \frac{af{\alpha }_{2}+am{n}_{2}-c{m}^{2}}{{\alpha }_{2}\left(a{n}_{2}-cm\right)}$

${U}^{\mathrm{*}} = \frac{c{n}_{1}-b{n}_{2}}{a\gamma -c{\alpha }_{1}+b{\alpha }_{2}}+\frac{f\gamma }{m\gamma -{n}_{2}{\alpha }_{1}+{n}_{1}{\alpha }_{2}}$

(14)

${V}_{S}^{\mathrm{*}} = \frac{-cm+a{n}_{2}}{a\gamma -c{\alpha }_{1}+b{\alpha }_{2}}+\frac{f{\alpha }_{2}}{m\gamma -{n}_{2}{\alpha }_{1}+{n}_{1}{\alpha }_{2}}$

(15)

${V}_{I}^{\mathrm{*}} = \frac{bm-a{n}_{1}}{a\gamma -c{\alpha }_{1}+b{\alpha }_{2}}-\frac{f{\alpha }_{1}}{m\gamma -{n}_{2}{\alpha }_{1}+{n}_{1}{\alpha }_{2}}$

(16)

In the present section, we study the stability of the fixed point, which deals with the coexistence of all species. The existence condition for the coexistence equilibrium point ${E}_{7}$ is

$f 𝛾( {\chi }_{1} + a\gamma ) > {\delta }_{1}(\gamma m+{\xi }_{1}) , f{\alpha }_{2} ( {\chi }_{2} + b{\alpha }_{2}) > {-\delta }_{2}({n}_{1}{\alpha }_{2}+{\xi }_{2}) ,\\ f{\alpha }_{1}\left({\chi }_{3}+c{\alpha }_{1}\right) > -{\delta }_{3}\left({n}_{2}{\alpha }_{1}-{\xi }_{3}\right) \ \ and \ \ \alpha \beta > 0$

(17)

where

${\delta }_{1} = b{n}_{2}-c{n}_{1} , {\delta }_{2} = a{n}_{2}-c\mathrm{m} , {\delta }_{3} = a{n}_{1}-bm, {\xi }_{1} = {n}_{1}{\alpha }_{2}-{n}_{2}{\alpha }_{1} ,$

${\xi }_{2} = \gamma m-{n}_{2}{\alpha }_{1} , {\xi }_{3} = \gamma m+{n}_{1}{\alpha }_{2} , \alpha = \gamma m+{\xi }_{1} , \beta = a\gamma +{\chi }_{1} ,$

${\chi }_{1} = b{\alpha }_{2}-c{\alpha }_{1}, {\chi }_{2} = a\gamma -c{\alpha }_{1} , {\chi }_{3} = -b{\alpha }_{2}-a\gamma$

Theorem 1: Let $a, b, c, f, {\alpha }_{1}, {\alpha }_{2}, m, {n}_{1}$ and ${n}_{2}$ be positive; then, $({U}^{\mathrm{*}}, {V}_{S}^{\mathrm{*}}, {V}_{I}^{\mathrm{*}})$ is globally stable in $\omega = \{\left(U, {V}_{S}, {V}_{I}\right):U > 0, {V}_{S} > 0, {V}_{I} > 0\} .$

Proof: We construct the Lyapunov function to determine the stability of the coexistence equilibrium point. Consider the following

$W(U, {V}_{S}, {V}_{I}) = {\int }_{{U}^{\mathrm{*}}}^{U}\text{‍}\frac{\eta -{U}^{\mathrm{*}}}{\eta }d\eta +{\int }_{{V}_{S}^{\mathrm{*}}}^{{V}_{S}}\text{‍}\frac{\kappa -{V}_{S}^{\mathrm{*}}}{\kappa }d\kappa +{\int }_{{V}_{I}^{\mathrm{*}}}^{{V}_{I}}\text{‍}\frac{\xi -{z}^{\mathrm{*}}}{\xi }d\xi$

(18)

By taking derivative with respect to t, we have

$\frac{dW}{dt} = \frac{U-{U}^{\mathrm{*}}}{U}\frac{dU}{dt}+\frac{{V}_{S}-{V}_{S}^{\mathrm{*}}}{{V}_{S}}\frac{d{V}_{S}}{dt}+\frac{{V}_{I}-{V}_{I}^{\mathrm{*}}}{{V}_{I}}\frac{d{V}_{I}}{dt}$

(19)

$\frac{dW}{dt} = \left(U-{U}^{\mathrm{*}}\right)\left(\frac{af}{aU+b{V}_{S}+c{V}_{I}}+{\alpha }_{1}{V}_{S}+{\alpha }_{2}{V}_{I}-m\right)$

$+\left({V}_{S}-{V}_{{S}^{\mathrm{*}}}\right)\left(\frac{bf}{aU+b{V}_{S}+c{V}_{I}}-{\alpha }_{1}U-\gamma {V}_{I}-{n}_{1}\right)$

$+({V}_{I}-{V}_{{I}^{\mathrm{*}}})(\frac{cf}{aU+b{V}_{S}+c{V}_{I}}-{\alpha }_{2}U+\gamma {V}_{S}-{n}_{2})$

(20)

After some simplification, we get

$\frac{dW}{dt} = \left(U-{U}^{\mathrm{*}}\right)af\frac{a{U}^{\mathrm{*}}+b{V}_{S}^{\mathrm{*}}+c{V}_{I}^{\mathrm{*}}-aU-b{V}_{S}-c{V}_{I}}{\left(aU+b{V}_{S}+c{V}_{I}\right)\left(a{U}^{\mathrm{*}}+b{V}_{S}^{\mathrm{*}}+c{V}_{I}^{\mathrm{*}}\right)}$

$+\left({V}_{S}-{V}_{S}^{\mathrm{*}}\right)bf\frac{a{U}^{\mathrm{*}}+b{V}_{S}^{\mathrm{*}}+c{V}_{I}^{\mathrm{*}}-aU-b{V}_{S}-c{V}_{I}}{\left(aU+b{V}_{S}+c{V}_{I}\right)\left(a{U}^{\mathrm{*}}+b{V}_{S}^{\mathrm{*}}+c{V}_{I}^{\mathrm{*}}\right)}$

$+({V}_{I}-{V}_{I}^{\mathrm{*}})cf\frac{a{U}^{\mathrm{*}}+b{V}_{S}^{\mathrm{*}}+c{V}_{I}^{\mathrm{*}}-aU-b{V}_{S}-c{V}_{I}}{(aU+b{V}_{S}+c{V}_{I})(a{U}^{\mathrm{*}}+b{V}_{S}^{\mathrm{*}}+c{V}_{I}^{\mathrm{*}})}$

$\frac{dW}{dt} = -f\frac{(aU+b{V}_{S}+c{V}_{I}-(a{U}^{\mathrm{*}}+b{V}_{S}^{\mathrm{*}}+c{V}_{I}^{\mathrm{*}}){)}^{2}}{(aU+b{V}_{S}+c{V}_{I})(a{U}^{\mathrm{*}}+b{V}_{S}^{\mathrm{*}}+c{V}_{I}^{\mathrm{*}})}\le 0$

(21)

$\frac{dW}{dt} = 0 \text{ iff } aU+b{V}_{S}+c{V}_{I} = a{U}^{\mathrm{*}}+b{V}_{S}^{\mathrm{*}}+c{V}_{I}^{\mathrm{*}}$

Using the values from (16)–(18) and applying further simplification of (9) leads to the following result

$\frac{dU}{dt} = U(\frac{{\alpha }_{1}a}{b}\left({U}^{\mathrm{*}}-\mathrm{U}\right)+(\frac{{\alpha }_{1}c}{b}-{\alpha }_{2}\left)\right({V}_{I}^{\mathrm{*}}-{V}_{I}) )$

(22)

Then $U$ → ${U}^{\mathrm{*}}$ and ${V}_{I}$ → ${V}_{I}^{\mathrm{*}}$ when t→∞ provided that $U\left(0\right)\ge 0.$ Similarly ${V}_{S}$ → ${V}_{S}^{\mathrm{*}};$ hence, the largest invariant set { $U$ , ${V}_{I}, {V}_{S}ϵ\omega :$ $\frac{dW}{dt} = 0$ } has the unique positive equilibrium $({U}^{\mathrm{*}}, {V}_{S}^{\mathrm{*}}, {V}_{I}^{\mathrm{*}})$ . According to the LaSalle theorem $({U}^{\mathrm{*}}, {V}_{S}^{\mathrm{*}}, {V}_{I}^{\mathrm{*}})$ is globally stable in $\omega$ .

Corollary 1. For the equilibrium points ${E}_{i}$ $\left(i = \mathrm{1, 2}, 3, \dots , 6\right)$ , we have the following results.

(a). ${E}_{1}$ , ${E}_{2}$ and ${E}_{3}$ always exist and are globally stable.

(b). ${E}_{4}$ exists if ${{n}_{2}\delta }_{1} > cf\gamma$ and ${{n}_{2}\delta }_{1} > cf\gamma$ , and it is globally stable.

(c). ${E}_{5}$ and ${E}_{6}$ do not exist.

Proof:

(a). The results are straightforward from Eqs (13) and (22).

(b). Consider ${E}_{4} = \left({U}_{4}, {V}_{S, 4}, {V}_{I, 4}\right)$ , where

${V}_{S, 4} = -\frac{-b{{n}_{2}}^{2}+cf\gamma +c{n}_{1}{n}_{2}}{\gamma \left(b{n}_{2}-c{n}_{1}\right)}$

$= \frac{b{{n}_{2}}^{2}-cf\gamma -c{n}_{1}{n}_{2}}{\gamma \left(b{n}_{2}-c{n}_{1}\right)}$

We can rewrite as

${V}_{S, 4} = \frac{{{n}_{2}\delta }_{1}-cf\gamma }{\gamma {\delta }_{1}} \Rightarrow {{n}_{2}\delta }_{1} > cf\gamma$

Similarly,

${V}_{I, 4} = \frac{bf\gamma -b{n}_{1}{n}_{2}+c{{n}_{1}}^{2}}{\gamma \left(b{n}_{2}-c{n}_{1}\right)} \Rightarrow {{bf\gamma > n}_{1}\delta }_{1}$

The global stability is obvious from Eq (22).

(c). Consider ${E}_{5} = ({U}_{5}, {V}_{S, 5}, {V}_{I, 5})$ , where

${U}_{5} = -\frac{a{{n}_{1}}^{2}+bf{\alpha }_{1}-bm{n}_{1}}{{\alpha }_{1}\left(a{n}_{1}-bm\right)} \Rightarrow {U}_{5} = {{-bf{\alpha }_{1}-n}_{1}\delta }_{3} < 0$

Hence, the result, similarly, ${E}_{6}$ does not exist.

4. Global stability for a diffusive system

In the dynamical population model, random walking plays a central role in the structure of the habitat. The movement of species can vary from region to region depending on the distribution of food resources in the habitat. In addition, the species can vary with respect to density in the area where they live. The following section deals with the influence of diffusion on the coexistence of fixed points.

4.1. Global stability for the self-diffusive system

The self-diffusive system for (9)–(11) can be written as

$\frac{dU}{dt} = {D}_{11}△U+U(\frac{af}{aU+b{V}_{S}+c{V}_{I}}+{\alpha }_{1}{V}_{S}+{\alpha }_{2}{V}_{I}-m) , \ \ \ \ \ \ \ \ \ \ \text{ x ϵ Ω, t > 0 }$

(23)

$\frac{d{V}_{S}}{dt} = {D}_{22}△{V}_{S}+{V}_{S}(\frac{bf}{aU+b{V}_{S}+c{V}_{I}}-\gamma {V}_{I}-{\alpha }_{1}U-{n}_{1}) , \ \ \ \ \ \ \ \ \ \ \text{ x ϵ Ω, t > 0 }$

(24)

$\frac{d{V}_{I}}{dt} = {D}_{33}△{V}_{I}+{V}_{I}(\frac{cf}{aU+b{V}_{S}+c{V}_{I}}+\gamma {V}_{S}-{\alpha }_{2}U-{n}_{2}) , \ \ \ \ \ \ \ \ \ \ \text{ x ϵ Ω, t > 0}$

(25)

$\frac{\partial U}{\partial \vartheta } = \frac{\partial {V}_{S}}{\partial \vartheta } = \frac{\partial {V}_{I}}{\partial \vartheta } = 0, \mathrm{x}ϵ\partial \mathrm{\Omega }, \ \mathrm{t} > 0$

(26)

with the following initial condition

$\phi \left(0\right) = \left\{U\right(0) = {U}_{0}\ge 0, {V}_{S}(0) = {V}_{S, 0}\ge 0, {V}_{I}(0) = {V}_{I, 0}\ge 0\}$

(27)

where the self-diffusion coefficients ${D}_{11}, {D}_{22}, {D}_{33}$ are non-negative. According to Hollis' results, the solution to the above system has global existence and boundedness ^[47]. The self-diffusive system has an equilibrium point ${E}_{7} = ({U}^{\mathrm{*}}, {V}_{S}^{\mathrm{*}}, {V}_{I}^{\mathrm{*}})$ under the assumption of (17)

Theorem 2: Let $a, b, c, {D}_{11}, {D}_{22}, {D}_{33}, f, g, m, {n}_{1}, {n}_{2}, {\alpha }_{1}, {\alpha }_{2}$ and $\gamma$ be positive; then, $({U}^{\mathrm{*}}, {V}_{S}^{\mathrm{*}}, {V}_{I}^{\mathrm{*}})$ is globally stable in $ꙍ = \{\left(U, {V}_{S}, {V}_{I}\right):U > 0, {V}_{S}, {V}_{I} > 0\}$ .

Proof. It is obvious that $\{\left(U, {V}_{S}, {V}_{I}\right):U = 0\}$ , $\{\left(U, {V}_{S}, {V}_{I}\right):{V}_{S} = 0\}$ and $\{\left(U, {V}_{S}, {V}_{I}\right):{V}_{I} = 0\}$ are the invariant manifolds. We establish the Lyapunov function to determine the global stability of the self-diffusive system as follows

$W(U, {V}_{S}, {V}_{I}) = {\int }_{\mathrm{\Omega }}\text{‍}{\int }_{{U}^{\mathrm{*}}}^{U}\text{‍}\frac{\eta -{U}^{\mathrm{*}}}{\eta }d\eta dX+{\int }_{\mathrm{\Omega }}\text{‍}{\int }_{{N}_{S}^{\mathrm{*}}}^{{N}_{S}}\text{‍}\frac{\kappa -{V}_{S}^{\mathrm{*}}}{\kappa }d\kappa dX+{\int }_{\mathrm{\Omega }}\text{‍}{\int }_{{N}_{I}^{\mathrm{*}}}^{{N}_{I}}\text{‍}\frac{\xi -{V}_{I}^{\mathrm{*}}}{\xi }d\xi dX$

(28)

Now, taking derivative with respect to t on both sides, we get

$\frac{dW}{dt} = {\int }_{\mathrm{\Omega }}\text{‍}\frac{U-{U}^{\mathrm{*}}}{U}\frac{\partial U}{dt}dX+{\int }_{\mathrm{\Omega }}\text{‍}\frac{{V}_{S}-{V}_{S}^{\mathrm{*}}}{{V}_{S}}\frac{\partial {V}_{S}}{dt}dX+{\int }_{\mathrm{\Omega }}\text{‍}\frac{{V}_{I}-{V}_{I}^{\mathrm{*}}}{{V}_{I}}\frac{\partial {V}_{I}}{dt}dX$

$= {\int }_{\mathrm{\Omega }}\text{‍}\frac{U-{U}^{\mathrm{*}}}{U}\left({D}_{11}△U+\mathrm{U}\left(\frac{af}{aU+b{V}_{S}+c{V}_{I}}+{\alpha }_{1}{V}_{S}+{\alpha }_{2}{V}_{I}-m\right)\right)$

$+{\int }_{\mathrm{\Omega }}\text{‍}\frac{{V}_{S}-{V}_{S}^{\mathrm{*}}}{{V}_{S}}\left({D}_{22}△{V}_{S}+{V}_{S}\left(\frac{bf}{aU+b{V}_{S}+c{V}_{I}}-\gamma {V}_{I}-{\alpha }_{1}U-{n}_{1}\right)\right)$

$+{\int }_{\mathrm{\Omega }}\text{‍}\frac{{V}_{I}-{V}_{I}^{\mathrm{*}}}{{V}_{I}}{(D}_{33}△{V}_{I}+{V}_{I}(\frac{cf}{aU+b{V}_{S}+c{V}_{I}}+\gamma {V}_{S}-{\alpha }_{2}U-{n}_{2}))$

$\frac{dW}{dt}≕{M}_{1}+{M}_{2}+{M}_{3}+{M}_{4}$

(29)

where

${M}_{1} = -{D}_{11}{U}^{\mathrm{*}}{\int }_{\mathrm{\Omega }}\text{‍}\frac{1}{{U}^{2}}|\nabla {U|}^{2}dX\le 0$

(30)

${M}_{2} = -{D}_{22}{V}_{S}^{\mathrm{*}}{\int }_{\mathrm{\Omega }}\text{‍}\frac{1}{{V}_{S}^{2}}|\nabla {V}_{S}^{2}|dX\le 0$

(31)

${M}_{3} = -{D}_{33}{V}_{I}^{\mathrm{*}}{\int }_{\mathrm{\Omega }}\text{‍}\frac{1}{{V}_{I}^{2}}|\nabla {V}_{I}^{2}|dX\le 0$

(32)

${M}_{4} = {\int }_{\mathrm{\Omega }}\text{‍}\left(U-{U}^{\mathrm{*}}\right)\left[\left(\frac{af}{aU+b{V}_{S}+c{V}_{I}}+{\alpha }_{1}{V}_{S}+{\alpha }_{2}{V}_{I}-m\right)\right]dX$

$+{\int }_{\mathrm{\Omega }}\text{‍}\left({V}_{S}-{V}_{S}^{\mathrm{*}}\right)\left[\left(\frac{bf}{aU+b{V}_{S}+c{V}_{I}}-\gamma {V}_{I}-{\alpha }_{1}U-{n}_{1}\right)\right]dX$

$+{\int }_{\mathrm{\Omega }}\text{‍}({V}_{I}-{V}_{I}^{\mathrm{*}})\left[\right(\frac{cf}{aU+b{V}_{S}+c{V}_{I}}+\gamma {V}_{S}-{\alpha }_{2}U-{n}_{2}\left)\right]dX\le 0$

(33)

${M}_{4} = -{\int }_{\mathrm{\Omega }}\text{‍}f\frac{(aU+b{V}_{S}+c{V}_{I}-(a{U}^{\mathrm{*}}+b{V}_{S}^{\mathrm{*}}+c{V}_{I}^{\mathrm{*}}){)}^{2}}{(aU+b{V}_{S}+c{V}_{I})(a{U}^{\mathrm{*}}+b{V}_{S}^{\mathrm{*}}+c{V}_{I}^{\mathrm{*}})}\le 0$

(34)

It is obvious that

$\frac{dW}{dt}\le 0$

(35)

$\frac{dW}{dt} = 0 \text{ iff } aU+b{V}_{S}+c{V}_{I} = a{U}^{\mathrm{*}}+b{V}_{S}^{\mathrm{*}}+c{V}_{I}^{\mathrm{*}}$

Further simplification of (23) by using the values of the coexistence equilibrium point $\left({U}^{\mathrm{*}}, {V}_{S}^{\mathrm{*}}, {V}_{I}^{\mathrm{*}}\right)$ leads to the following

$\frac{dU}{dt} = {D}_{11}+U(\frac{{\alpha }_{1}a}{b}\left({U}^{\mathrm{*}}-\mathrm{U}\right)+\left(\frac{{\alpha }_{1}c}{b}-{\alpha }_{2}\right)\left({V}_{I}^{\mathrm{*}}-{V}_{I}\right))$

It is well known that the above equation represents a gradient system, where every orbit converges to a unique steady state $\mathrm{U} = {U}^{\mathrm{*}}$ and ${V}_{I} = {V}_{I}^{\mathrm{*}}$ ^[47]. Thus, ${V}_{S}$ → ${V}_{S}^{\mathrm{*}}$ since $\left(U, {V}_{S}, {V}_{I}\right)$ is on $\Gamma ≔\left\{\right(U, {V}_{S}, {V}_{I}):aU+b{V}_{S}+c{V}_{I} = a{U}^{\mathrm{*}}+b{V}_{S}^{\mathrm{*}}+c{V}_{I}^{\mathrm{*}}, U\ge 0, {V}_{S}\ge 0, {V}_{I}\ge 0\}$ . Therefore, positive solution $\left({U}^{\mathrm{*}}, {V}_{S}^{\mathrm{*}}, {V}_{I}^{\mathrm{*}}\right)$ is globally asymptotically stable.

4.2. Induced stability for the cross-diffusive system

The predator is attracted to its prey in search of food resources. This attraction causes prey taxis, which is the predator's movement as controlled by the density of prey, and on the contrary, the prey is repelled by the predator. This results in the cross-diffusion of species in the habitat. Consider the following system of equations:

$\frac{dU}{dt} = {D}_{11}△U+{D}_{12}△{V}_{S}+{D}_{13}△{V}_{I}+U\left(\frac{af}{aU+b{V}_{S}+c{V}_{I}}+{\alpha }_{1}{V}_{S}+{\alpha }_{2}{V}_{I}-m\right) ,$

$xϵ\varOmega , t > 0$

(36)

$\frac{d{V}_{S}}{dt} = {D}_{21}△U+{D}_{22}△{V}_{S}+{V}_{S}\left(\frac{bf}{aU+b{V}_{S}+c{V}_{I}}-\gamma {V}_{I}-{\alpha }_{1}U-{n}_{1}\right) ,$

$xϵ\varOmega , t > 0$

(37)

$\frac{d{V}_{I}}{dt} = {D}_{31}△U+{D}_{33}△{V}_{I}+{V}_{I}\left(\frac{cf}{aU+b{V}_{S}+c{V}_{I}}+\gamma {V}_{S}-{\alpha }_{2}U-{n}_{2}\right),$

$xϵ\varOmega , t > 0$

(38)

$\frac{\partial U}{\partial \vartheta } = \frac{\partial {V}_{S}}{\partial \vartheta } = \frac{\partial {V}_{I}}{\partial \vartheta } = 0 , xϵ\partial \varOmega , t > 0$

$\phi \left(0\right) = \left\{U\right(0) = {U}_{0}\ge 0, {V}_{S}(0) = {V}_{S, 0}\ge 0, {V}_{I}(0) = {V}_{I, 0}\ge 0\} ,\text{ x ϵ Ω, t > 0 }$

(39)

One can see that the coexisting equilibrium point also holds for the cross-diffusive system. The self-diffusive system does not induce instability, but the Turing instability can be seen in the following section. The above system, as described by Eqs (36)–(38), can be written in the following way:

$R\mathrm{\Phi } = D△\mathrm{\Phi }+J\mathrm{\Phi }$

(40)

$\mathrm{\Phi }\left(t, x\right) = \left[\begin{array}{c}\phi \left(t, x\right)\\ \varphi \left(t, x\right)\\ \psi \left(t, x\right)\end{array}\right], D = \left[\begin{array}{ccc}{D}_{11}& {D}_{12}& {D}_{13}\\ {D}_{21}& {D}_{22}& 0\\ {D}_{31}& 0& {D}_{33}\end{array}\right] \text{ and } J = \left[\begin{array}{ccc}{A}_{11}& {A}_{12}& {A}_{13}\\ {A}_{21}& {A}_{22}& {A}_{23}\\ {A}_{31}& {A}_{32}& {A}_{33}\end{array}\right]$

where

${A}_{11} = -\frac{{a}^{2}f{U}^{\mathrm{*}}}{{\left(a{U}^{\mathrm{*}}+c{V}_{I}^{\mathrm{*}}+b{V}_{S}^{\mathrm{*}}\right)}^{2}},$

${A}_{12} = -\frac{abf{U}^{\mathrm{*}}}{{\left(a{U}^{\mathrm{*}}+c{V}_{I}^{\mathrm{*}}+b{V}_{S}^{\mathrm{*}}\right)}^{2}}+{\alpha }_{1}{U}^{\mathrm{*}},$

${A}_{13} = -\frac{acf{U}^{\mathrm{*}}}{{\left(a{U}^{\mathrm{*}}+c{V}_{I}^{\mathrm{*}}+b{V}_{S}^{\mathrm{*}}\right)}^{2}}+{\alpha }_{2}{U}^{\mathrm{*}},$

${A}_{21} = -\frac{abf{V}_{S}^{\mathrm{*}}}{{\left(a{U}^{\mathrm{*}}+c{V}_{I}^{\mathrm{*}}+b{V}_{S}^{\mathrm{*}}\right)}^{2}}-{\alpha }_{1}{V}_{S}^{\mathrm{*}},$

${A}_{22} = -\frac{{b}^{2}f{V}_{S}^{\mathrm{*}}}{{\left(a{U}^{\mathrm{*}}+c{V}_{I}^{\mathrm{*}}+b{V}_{S}^{\mathrm{*}}\right)}^{2}},$

${A}_{23} = -\frac{bcf{V}_{S}^{\mathrm{*}}}{{\left(a{U}^{\mathrm{*}}+c{V}_{I}^{\mathrm{*}}+b{V}_{S}^{\mathrm{*}}\right)}^{2}}-\gamma {V}_{S}^{\mathrm{*}},$

${A}_{31} = -\frac{acf{V}_{I}^{\mathrm{*}}}{{\left(aU+c{V}_{I}+b{V}_{S}\right)}^{2}}-{\alpha }_{2}{V}_{I}^{\mathrm{*}},$

${A}_{32} = -\frac{bcf{V}_{I}^{\mathrm{*}}}{{\left(a{U}^{\mathrm{*}}+c{V}_{I}^{\mathrm{*}}+b{V}_{S}^{\mathrm{*}}\right)}^{2}}+\gamma {V}_{I}^{\mathrm{*}},$

${A}_{33} = -\frac{{c}^{2}f{V}_{I}^{\mathrm{*}}}{{\left(a{U}^{\mathrm{*}}+c{V}_{I}^{\mathrm{*}}+b{V}_{S}^{\mathrm{*}}\right)}^{2}}$

(41)

If all of the eigenvalues of the operator $R$ are negative, then $({U}^{\mathrm{*}}, {V}_{S}^{\mathrm{*}}, {V}_{I}^{\mathrm{*}})$ is asymptotically stable ^[47]. Next, it is clear that $-△$ under the Neumann boundary condition has eigenvalues

$0 = {\mu }_{0} < {\mu }_{1}\le {\mu }_{2}\le \dots$

and ${\mu }_{i}$ → ∞ when t → ∞. It is also known that $\nu$ is the eigenvalue of the operator $R$ if and only if it is the eigenvalue of ${J}_{i} = -{\mu }_{i}D+J$ for $i$ ≥ 0. Moreover, we obtain

$\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{e}\left({J}_{i}\right) = -{\mu }_{i}({D}_{11}+{D}_{22}+{D}_{33})+\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{e}(J )$

It is obvious from (41) that

$\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{e}(J ) = {A}_{11}+{A}_{22}+{A}_{33} < 0.$

Next, we prove Turing instability for the cross-diffusive system described by (36)–(38).

Theorem 3: Let $a, b, c, {D}_{11}, {D}_{22}, {D}_{33}, f, g, m, {n}_{1}, {n}_{2}, {\alpha }_{1}, {\alpha }_{2}$ and $\gamma$ be positive and

$det\left(D\right) = \left({D}_{11}{D}_{22}{D}_{33}-{D}_{12}{D}_{21}{D}_{33}+{D}_{13}{D}_{22}{D}_{31}\right) > 0$

then, there exists an unbounded region where ${s}_{1} < {s}_{2}+{s}_{3}+{s}_{4}$

such that $({U}^{\mathrm{*}}, {V}_{S}^{\mathrm{*}}, {V}_{I}^{\mathrm{*}})$ is unstable in the region where the above condition holds.

Additionally,

${s}_{1} = -[\frac{1}{3}\frac{F\left(J, D\right)-\sqrt{3 \ \mathrm{d}\mathrm{e}\mathrm{t}\left(D\right)f\left(J, D\right)+{\left(F\left(J, D\right)\right)}^{2}}}{\mathrm{d}\mathrm{e}\mathrm{t}\left(D\right)}{]}^{3}\mathrm{d}\mathrm{e}\mathrm{t}\left(D\right)$

${s}_{2} = [\frac{1}{3}\frac{F\left(J, D\right)-\sqrt{3 \ \mathrm{d}\mathrm{e}\mathrm{t}\left(D\right)f\left(J, D\right)+{\left(F\left(J, D\right)\right)}^{2}}}{\mathrm{d}\mathrm{e}\mathrm{t}\left(D\right)}{]}^{2}F\left(J, D\right)$

${s}_{3} = [1/3\frac{F\left(J, D\right)-\sqrt{3 \ \mathrm{d}\mathrm{e}\mathrm{t}\left(D\right)f\left(J, D\right)+{\left(F\left(J, D\right)\right)}^{2}}}{\mathrm{d}\mathrm{e}\mathrm{t}\left(D\right)}]f\left(J, D\right)$

${s}_{4} = {A}_{11}{A}_{22}{A}_{33}-{A}_{11}{A}_{23}{A}_{32}-{A}_{12}{A}_{21}{A}_{33}$

$+{A}_{12}{A}_{23}{A}_{31}+{A}_{13}{A}_{21}{A}_{32}-{A}_{13}{A}_{22}{A}_{31}$

Proof. It is clear from the value of ${J}_{i}$ that $trace\left({J}_{i}\right)$ is negative. To show the instability of $({U}^{\mathrm{*}}, {V}_{S}^{\mathrm{*}}, {V}_{I}^{\mathrm{*}})$ , we shall prove that $det\left({J}_{i}\right) > 0$ for some $i$ , where $i\ge 1$ . As

$\mathrm{det}\left({J}_{0}\right) = -{A}_{11}{A}_{21}{A}_{33}+{A}_{11}{A}_{22}{A}_{33}$

$+{A}_{11}{A}_{23}{A}_{31}-{A}_{11}{A}_{23}{A}_{32}$

$+{A}_{13}{A}_{21}{A}_{32}-{A}_{13}{A}_{22}{A}_{31}$

(42)

$det\left({J}_{i}\right) = \left[\begin{array}{lll}-{\mu }_{i}{D}_{11}+{A}_{11}& -{\mu }_{i}{D}_{12}+{A}_{12}& -{\mu }_{i}{D}_{13}+{A}_{13}\\ -{\mu }_{i}{D}_{21}+{A}_{21}& -{\mu }_{i}{D}_{22}+{A}_{22}& {A}_{23}\\ -{\mu }_{i}{D}_{31}+{A}_{31}& {A}_{32}& -{\mu }_{i}{D}_{33}+{A}_{33}\end{array}\right]$

(43)

$det\left({J}_{i}\right) = -{D}_{11}{D}_{22}{D}_{33}{{\mu }_{i}}^{3}+{A}_{11}{D}_{22}{D}_{33}{{\mu }_{i}}^{2}+{A}_{22}{D}_{11}{D}_{33}{{\mu }_{i}}^{2}+{A}_{33}{D}_{11}{D}_{22}{{\mu }_{i}}^{2}\\ -{A}_{11}{A}_{22}{D}_{33}{\mu }_{i}-{A}_{11}{A}_{33}{D}_{22}{\mu }_{i}-{A}_{22}{A}_{33}{D}_{11}{\mu }_{i}+{A}_{23}{A}_{32}{D}_{11}{\mu }_{i}\\ +{D}_{12}{D}_{21}{D}_{33}{{\mu }_{i}}^{3}-{A}_{12}{D}_{21}{D}_{33}{{\mu }_{i}}^{2}-{A}_{21}{D}_{12}{D}_{33}{{\mu }_{i}}^{2}+{A}_{32}{D}_{13}{D}_{21}{{\mu }_{i}}^{2}\\ -{A}_{33}{D}_{12}{D}_{21}{{\mu }_{i}}^{2}+{A}_{12}{A}_{21}{D}_{33}{\mu }_{i}+{A}_{12}{A}_{33}{D}_{21}{\mu }_{i}-{A}_{13}{A}_{32}{D}_{21}{\mu }_{i}\\ -{A}_{21}{A}_{32}{D}_{13}{\mu }_{i}+{A}_{21}{A}_{33}{D}_{12}{\mu }_{i}+{D}_{13}{D}_{22}{D}_{31}{{\mu }_{i}}^{3}-{A}_{13}{D}_{22}{D}_{31}{{\mu }_{i}}^{2}\\ -{A}_{22}{D}_{13}{D}_{31}{{\mu }_{i}}^{2}+{A}_{23}{D}_{12}{D}_{31}{{\mu }_{i}}^{2}-{A}_{31}{D}_{13}{D}_{22}{{\mu }_{i}}^{2}-{A}_{12}{A}_{23}{D}_{31}{\mu }_{i}\\ +{A}_{13}{A}_{22}{D}_{31}{\mu }_{i}+{A}_{13}{A}_{31}{D}_{22}{\mu }_{i}+{A}_{22}{A}_{31}{D}_{13}{\mu }_{i}-{A}_{23}{A}_{31}{D}_{12}{\mu }_{i}\\ -{A}_{12}{A}_{21}{A}_{33}+{A}_{13}{A}_{21}{A}_{32}+{A}_{12}{A}_{23}{A}_{31}-{A}_{13}{A}_{22}{A}_{31}+{A}_{11}{A}_{22}{A}_{33}\\ -{A}_{11}{A}_{23}{A}_{32}$

By collecting the terms containing ${\mu }_{i}$ , we obtain

$det\left({J}_{i}\right) = \left(-{D}_{11}{D}_{22}{D}_{33}+{D}_{12}{D}_{21}{D}_{33}+{D}_{13}{D}_{22}{D}_{31}\right){{\mu }_{i}}^{3}\\ +\left({A}_{11}{D}_{22}{D}_{33}+{A}_{22}{D}_{11}{D}_{33}+{A}_{33}{D}_{11}{D}_{22}-{A}_{12}{D}_{21}{D}_{33}-{A}_{21}{D}_{12}{D}_{33}\\ +{A}_{32}{D}_{13}{D}_{21}-{A}_{33}{D}_{12}{D}_{21}-{A}_{13}{D}_{22}{D}_{31}-{A}_{22}{D}_{13}{D}_{31}+{A}_{23}{D}_{12}{D}_{31}\\ -{A}_{31}{D}_{13}{D}_{22}\right){{\mu }_{i}}^{2}\\ +\left(-{A}_{11}{A}_{22}{D}_{33}-{A}_{11}{A}_{33}{D}_{22}-{A}_{22}{A}_{33}{D}_{11}+{A}_{23}{A}_{32}{D}_{11}+{A}_{12}{A}_{21}{D}_{33}\\ +{A}_{12}{A}_{33}{D}_{21}-{A}_{13}{A}_{32}{D}_{21}-{A}_{21}{A}_{32}{D}_{13}+{A}_{21}{A}_{33}{D}_{12}-{A}_{12}{A}_{23}{D}_{31}\\ +{A}_{13}{A}_{22}{D}_{31}+{A}_{13}{A}_{31}{D}_{22}+{A}_{22}{A}_{31}{D}_{13}-{A}_{23}{A}_{31}{D}_{12}\right){\mu }_{i}+{A}_{11}{A}_{22}{A}_{33}\\ -{A}_{11}{A}_{23}{A}_{32}-{A}_{12}{A}_{21}{A}_{33}+{A}_{12}{A}_{23}{A}_{31}+{A}_{13}{A}_{21}{A}_{32}-{A}_{13}{A}_{22}{A}_{31}$

$\mathrm{d}\mathrm{e}\mathrm{t}\left({J}_{i}\right) = -{{\mu }_{i}}^{3}\mathrm{d}\mathrm{e}\mathrm{t}\left(D\right)+{{\mu }_{i}}^{2}F\left(J, D\right)$

$+{\mu }_{i}f\left(J, D\right)+{A}_{11}{A}_{22}{A}_{33}-{A}_{11}{A}_{23}{A}_{32}-{A}_{12}{A}_{21}{A}_{33}$

$+{A}_{12}{A}_{23}{A}_{31}+{A}_{13}{A}_{21}{A}_{32}-{A}_{13}{A}_{22}{A}_{31}$

(44)

where

$det\left(D\right) = \left({D}_{11}{D}_{22}{D}_{33}-{D}_{12}{D}_{21}{D}_{33}+{D}_{13}{D}_{22}{D}_{31}\right) > 0$

$F(J, D) = \left(\begin{array}{c}{A}_{11}{D}_{22}{D}_{33}+{A}_{22}{D}_{11}{D}_{33}\\ +{A}_{33}{D}_{11}{D}_{22}-{A}_{12}{D}_{21}{D}_{33}-{A}_{21}{D}_{12}{D}_{33}\\ +{A}_{32}{D}_{13}{D}_{21}-{A}_{33}{D}_{12}{D}_{21}-{A}_{13}{D}_{22}{D}_{31}-{A}_{22}{D}_{13}{D}_{31}\\ +{A}_{23}{D}_{12}{D}_{31}-{A}_{31}{D}_{13}{D}_{22}\end{array}\right)$

(45)

$f\left(J, D\right) = -{A}_{11}{A}_{22}{D}_{33}-{A}_{11}{A}_{33}{D}_{22}-{A}_{22}{A}_{33}{D}_{11}$

$+{A}_{23}{A}_{32}{D}_{11}+{A}_{12}{A}_{21}{D}_{33}$

$+{A}_{12}{A}_{33}{D}_{21}-{A}_{13}{A}_{32}{D}_{21}-{A}_{21}{A}_{32}{D}_{13}$

$+{A}_{21}{A}_{33}{D}_{12}-{A}_{12}{A}_{23}{D}_{31}+{A}_{13}{A}_{22}{D}_{31}$

$+{A}_{13}{A}_{31}{D}_{22}+{A}_{22}{A}_{31}{D}_{13}-{A}_{23}{A}_{31}{D}_{12}$

(46)

The above equality of (44) leads to its minimum value ${min}_{{\mu }_{i}\mathrm{ϵ}{\mathrm{R}}^{+}}$ $\mathrm{d}\mathrm{e}\mathrm{t}\left({J}_{i}\right)$ at

$\mathrm{\mu } = \ 1/3\frac{F\left(J, D\right)-\sqrt{3 \ \mathrm{d}\mathrm{e}\mathrm{t}\left(D\right)f\left(J, D\right)+{\left(F\left(J, D\right)\right)}^{2}}}{\mathrm{d}\mathrm{e}\mathrm{t}\left(D\right)}$

After substituting the above value in (44), we get

$\mathrm{D}\mathrm{e}\mathrm{t}\left({J}_{i}\right) = \ -[1/3\frac{F\left(J, D\right)-\sqrt{3 \ \mathrm{d}\mathrm{e}\mathrm{t}\left(D\right)f\left(J, D\right)+{\left(F\left(J, D\right)\right)}^{2}}}{\mathrm{d}\mathrm{e}\mathrm{t}\left(D\right)}{]}^{3}\mathrm{d}\mathrm{e}\mathrm{t}\left(D\right)$

$+[1/3\frac{F\left(J, D\right)-\sqrt{3 \ \mathrm{d}\mathrm{e}\mathrm{t}\left(D\right)f\left(J, D\right)+{\left(F\left(J, D\right)\right)}^{2}}}{\mathrm{d}\mathrm{e}\mathrm{t}\left(D\right)}{]}^{2}F\left(J, D\right)$

$+\left[\frac{1}{3}\frac{F\left(J, D\right)-\sqrt{3 \ \mathrm{d}\mathrm{e}\mathrm{t}\left(D\right)f\left(J, D\right)+{\left(F\left(J, D\right)\right)}^{2}}}{\mathrm{d}\mathrm{e}\mathrm{t}\left(D\right)}\right]f\left(J, D\right)$

$+{A}_{11}{A}_{22}{A}_{33}-{A}_{11}{A}_{23}{A}_{32}-{A}_{12}{A}_{21}{A}_{33}$

$+{A}_{12}{A}_{23}{A}_{31}+{A}_{13}{A}_{21}{A}_{32}-{A}_{13}{A}_{22}{A}_{31}$

$\mathrm{d}\mathrm{e}\mathrm{t}\left({J}_{i}\right) = -{s}_{1}+{s}_{2}+{s}_{3}+{s}_{4}$

(47)

For the instability of the system, we must have $det\left({J}_{i}\right) > 0.$ Therefore above the equilibrium point $({U}^{\mathrm{*}}, {V}_{S}^{\mathrm{*}}, {V}_{I}^{\mathrm{*}})$ is unstable in the region where ${s}_{1} < {s}_{2}+{s}_{3}+{s}_{4}$ .

5. Construction of numerical scheme

A numerical scheme is constructed for solving Eqs (23)–(27). The scheme's construction discretizes the time variable, and suitable difference formula can be adopted to discretize the space variable. The scheme is the two-stage explicit scheme. The scheme is constructed on three-time levels. Consider the following equation:

$\frac{\partial V}{\partial t} = G\left({V}_{1}\frac{{\partial }^{2}V}{\partial {x}^{2}}\right)$

(48)

The first stage of the scheme is written as

${\stackrel{-}{V}}_{i}^{n+1} = {V}_{i}^{n}+\Delta t{\left(\frac{\partial V}{\partial t}\right)}_{i}^{n}$

(49)

The second stage of the scheme contains three unknowns. Their values will be found later. The second stage is written as

${V}_{i}^{n+1} = \frac{1}{2}({\stackrel{-}{V}}_{i}^{n+1}+{V}_{i}^{n-1})+\mathrm{\Delta }t\left\{{a}_{1}{\left(\frac{\partial V}{\partial t}\right)}_{i}^{n+1}+{b}_{1}{\left(\frac{\partial V}{\partial t}\right)}_{i}^{n}+{c}_{1}{\left(\frac{\partial V}{\partial t}\right)}_{i}^{n-1}\right\}$

(50)

To find the values of unknowns ${a}_{1}, {b}_{1}$ and ${c}_{1}$ , the Taylor series expansions for ${V}_{i}^{n+1}, {V}_{i}^{n-1}, {\left(\frac{\partial V}{\partial t}\right)}_{i}^{n+1}$ and ${\left(\frac{\partial V}{\partial t}\right)}_{i}^{n-1}$ are given as

${V}_{i}^{n+1} = {V}_{i}^{n}+\mathrm{\Delta }t{\left(\frac{\partial V}{\partial t}\right)}_{i}^{n}+\frac{{\left(\mathrm{\Delta }t\right)}^{2}}{2}{\left(\frac{{\partial }^{2}V}{\partial {t}^{2}}\right)}_{i}^{n}+\frac{{\left(\mathrm{\Delta }t\right)}^{3}}{6}{\left(\frac{{\partial }^{3}v}{\partial {t}^{3}}\right)}_{i}^{n}+O\left({\left(\mathrm{\Delta }t\right)}^{4}\right)$

(51)

${V}_{i}^{n-1} = {V}_{i}^{n}-\mathrm{\Delta }t{\left(\frac{\partial V}{\partial t}\right)}_{i}^{n}+\frac{{\left(\mathrm{\Delta }t\right)}^{2}}{2}{\left(\frac{{\partial }^{2}V}{\partial {t}^{2}}\right)}_{i}^{n}-\frac{{\left(\mathrm{\Delta }t\right)}^{3}}{6}{\left(\frac{{\partial }^{3}V}{\partial {t}^{3}}\right)}_{i}^{n}+O\left({\left(\mathrm{\Delta }t\right)}^{4}\right)$

(52)

${\left(\frac{\partial V}{\partial t}\right)}_{i}^{n+1} = {\left(\frac{\partial V}{\partial t}\right)}_{i}^{n}+\mathrm{\Delta }t{\left(\frac{{\partial }^{2}V}{\partial {t}^{2}}\right)}_{i}^{n}+\frac{{\left(\mathrm{\Delta }t\right)}^{2}}{2}{\left(\frac{{\partial }^{3}V}{\partial {t}^{3}}\right)}_{i}^{n}+O\left({\left(\mathrm{\Delta }t\right)}^{3}\right)$

(53)

${\left(\frac{\partial V}{\partial t}\right)}_{i}^{n-1} = {\left(\frac{\partial V}{\partial t}\right)}_{i}^{n}-\mathrm{\Delta }t{\left(\frac{{\partial }^{2}V}{\partial {t}^{2}}\right)}_{i}^{n}+\frac{{\left(\mathrm{\Delta }t\right)}^{2}}{2}{\left(\frac{{\partial }^{3}V}{\partial {t}^{3}}\right)}_{i}^{n}+O\left({\left(\mathrm{\Delta }t\right)}^{3}\right)$

(54)

Substituting Taylor series expansions (51)–(54) into Eq (50), we obtain the following

Comparison of coefficients of ${u}_{i}^{n}, \mathrm{\Delta }t{\left(\frac{\partial u}{\partial t}\right)}_{i}^{n}, {\left(\mathrm{\Delta }t\right)}^{2}{\left(\frac{{\partial }^{2}u}{\partial {t}^{2}}\right)}_{i}^{n}$ and ${\left(\mathrm{\Delta }t\right)}^{3}{\left(\frac{{\partial }^{3}u}{\partial {t}^{3}}\right)}_{i}^{n}$ leads to the following

$\left\{\begin{array}{c}1 = {a}_{1}+{b}_{1}+{c}_{1}\\ \frac{1}{2} = \frac{1}{4}+{a}_{1}-{c}_{1}\\ \frac{1}{6} = -\frac{1}{12}+\frac{{c}_{1}}{2}\end{array}\right.$

(55)

Solving the above system of equations gives the values of the unknown parameters ${a}_{1}, {b}_{1}$ and ${c}_{1}$ as

${a}_{1} = \frac{3}{4}, {b}_{1} = -\frac{1}{4}, {c}_{1} = \frac{1}{2}$

(56)

5.1. Stability analysis

The stability condition of the proposed scheme for the parabolic Eq (48) von Neumann stability criterion is applied. For applying the stability analysis, some transformations are substituted into difference equations. For the considered problem, the transformations are given as follows.

Using a second-order central difference formula for the diffusion term, the first stage of the proposed scheme for the linearized equation is given as

${\stackrel{-}{V}}_{i}^{n+1} = {V}_{i}^{n}+\frac{1}{4}{u}_{i}^{n-1}+d({V}_{i+1}^{n}-2{V}_{i}^{n}-{V}_{i-1}^{n})$

(57)

where

$d = \frac{\Delta \mathrm{t}}{{(\Delta \mathrm{x})}^{2}}$

According to von Neumann's stability criterion, consider the following transformations

${V}_{i}^{n+1} = {\stackrel{-}{E}}^{n+1}{e}^{iI\mathrm{\theta }}, {V}_{i\pm 1}^{n} = {E}^{n}{e}^{(i\pm 1)I\mathrm{\theta }}, {u}_{i\pm 1}^{n-1} = {E}^{n+1}{e}^{(i\pm 1)I\mathrm{\theta }}$

$\text{ where } I = \sqrt{-1}$

(58)

Substituting the transformations of (50) into Eq (49), the following is obtained:

${\stackrel{-}{E}}^{n+1} = {E}^{n}+d({e}^{I\theta }-2+{e}^{-I\theta }){E}^{n}$

${\stackrel{-}{E}}^{n+1} = (1+2d\left(cos\theta -1\right)){E}^{n}$

(59)

By using the second-order central difference formula for the diffusion term, the second stage of the proposed scheme for the linearized equation is given as

${V}_{i}^{n+1} = \frac{1}{2}\left({\stackrel{-}{V}}_{i}^{n+1}+{V}_{i}^{n-1}\right)+d\left\{\begin{array}{c}{a}_{1}\left({\stackrel{-}{V}}_{i+1}^{n+1}-2{\stackrel{-}{V}}_{i}^{n+1}{+\stackrel{-}{V}}_{i-1}^{n+1}\right)\\ +{b}_{1}\left({V}_{i+1}^{n}-2{V}_{i}^{n}+{V}_{i-1}^{n}\right)\\ +{c}_{1}({V}_{i+1}^{n-1}-2\end{array}\right\}$

${V}_{i}^{n-1}+{V}_{i-1}^{n-1}) \}$

(60)

Substituting the transformations of (58) into (60) and dividing the resulting equation by ${e}^{iI\mathrm{\theta }}$ , we get

${E}^{n+1} = \frac{1}{2}\left({\stackrel{-}{E}}_{i}^{n+1}+{E}_{i}^{n-1}\right)+d\left\{\begin{array}{c}{a}_{1}({e}^{I\theta }-2+{e}^{-I\theta }){\stackrel{-}{E}}^{n+1}+\\ {b}_{1}({e}^{I\theta }-2+{e}^{-I\theta }){E}^{n}+\\ {c}_{1}({e}^{I\theta }-2+{e}^{-I\theta }){E}^{n-1}\end{array}\right\}$

Collecting the coefficients of ${E}^{n+1}$ on the left-hand side of Eq (61) provides

${E}^{n+1} = \left(\frac{1}{2}+2{a}_{1}d\left(cos\theta -1\right)\right){\stackrel{-}{E}}^{n+1}+2bd\left(cos\theta -1\right){E}^{n}+\left(\frac{1}{2}+2{c}_{1}d\left(cos\theta -1\right)\right){E}^{n-1}$

${E}^{n+1} = {A}_{1}{E}^{n}+{B}_{1}{E}^{n-1}$

(61)

${A}_{1} = \left(\frac{1}{2}+2{a}_{1}d\left(cos\theta -1\right)\right)(1+2d\left(cos\theta -1\right)) + 2{b}_{1}d\left(cos\theta -1\right)$

${B}_{1} = \frac{1}{2}+2{c}_{1}d\left(cos\theta -1\right)$

${E}^{n} = {E}^{n}+0{E}^{n-1}$

(62)

Equations (61) and (62) can be written in matrix form as

$⌈\begin{array}{c}{E}^{n+1}\\ {E}^{n}\end{array}⌉ = \left[\begin{array}{cc}{A}_{1}& {B}_{1}\\ 1& 0\end{array}\right]\left[\begin{array}{c}{E}^{n}\\ {E}^{n-1}\end{array}\right]$

(63)

The stability conditions can be imposed on the eigenvalue of the coefficient matrix, and these conditions are expressed as

$\left|\frac{{A}_{1}+\sqrt{{{A}_{1}}^{2}+4{B}_{1}}}{2}\right|\le 1 \text{ and } \left|\frac{{A}_{1}-\sqrt{{{A}_{1}}^{2}+4{B}_{1}}}{2}\right|\le 1$

(64)

The proposed scheme will be stable if it satisfies inequalities of (64).

6. Numerical simulations

The spatial terms in (23)–(25) are discretized by using the classical central difference formula to apply the proposed scheme for the considered diffusive model. The discretized equations for the model are given as

${\stackrel{-}{U}}_{i}^{n+1} = {U}_{i}^{n}+\Delta t\left\{\begin{array}{c}{D}_{11}\left({U}_{i+1}^{n}-2{U}_{i}^{n}+{U}_{i-1}^{n}\right)/{\left(\Delta x\right)}^{2}\\ +{U}_{i}^{n}\left(\frac{af}{a{U}_{i}^{n}+b{V}_{S, i}^{n}+c{V}_{I, i}^{n}}+{\alpha }_{1}{V}_{S, i}^{n}+{\alpha }_{2}{V}_{I, i}^{n}-m\right)\end{array}\right\}$

(65)

${\stackrel{-}{V}}_{S, i}^{n+1} = {V}_{S, i}^{n}+\Delta t\left\{\begin{array}{c}{D}_{22}\left({V}_{S, i+1}^{n}-2{V}_{S, i}^{n}+{V}_{S, i-1}^{n}\right)/{\left(\Delta x\right)}^{2}\\ +{U}_{S, i}^{n}\left(\frac{bf}{a{U}_{i}^{n}+b{V}_{S, i}^{n}+c{V}_{I, i}^{n}}-{\alpha }_{1}{U}_{i}^{n}-\mathrm{\gamma }{V}_{I, i}^{n}-{n}_{1}\right)\end{array}\right\}$

(66)

${\stackrel{-}{V}}_{S, i}^{n+1} = {V}_{I, i}^{n}+\Delta t\left\{\begin{array}{c}{D}_{33}\left({V}_{I, i+1}^{n}-2{V}_{I, i}^{n}+{V}_{I, i-1}^{n}\right)/{\left(\Delta x\right)}^{2}\\ +{U}_{I, i}^{n}\left(\frac{cf}{a{U}_{i}^{n}+b{V}_{S, i}^{n}+c{V}_{I, i}^{n}}-{\alpha }_{2}{U}_{i}^{n}-\mathrm{\gamma }{V}_{S, i}^{n}-{n}_{2}\right)\end{array}\right\}$

(67)

${U}_{i}^{n+1} = {U}_{i}^{n}+\Delta t\left[\begin{array}{c}{a}_{1}\left\{\begin{array}{c}{D}_{11}\left({\stackrel{-}{U}}_{i+1}^{n+1}-2{\stackrel{-}{U}}_{i}^{n+1}+{\stackrel{-}{U}}_{i-1}^{n+1}\right)/{\left(\Delta x\right)}^{2}+\\ {\stackrel{-}{U}}_{i}^{n+1}\left(\frac{af}{a{\stackrel{-}{U}}_{i}^{n+1}+b{\stackrel{-}{V}}_{S, i}^{n+1}+c{\stackrel{-}{V}}_{I, i}^{n+1}}+{\alpha }_{1}{\stackrel{-}{V}}_{S, i}^{n+1}+{\alpha }_{2}{\stackrel{-}{V}}_{I, i}^{n+1}-m\right)\end{array}\right\}+\\ {b}_{1}\left\{\begin{array}{c}{D}_{11}\left({U}_{i+1}^{n}-2{U}_{i}^{n}+{U}_{i-1}^{n}\right)/{\left(\Delta x\right)}^{2}+\\ {U}_{i}^{n}\left(\frac{af}{a{U}_{i}^{n}+b{V}_{S, i}^{n}+c{V}_{I, i}^{n}}+{\alpha }_{1}{V}_{S, i}^{n}+{\mathrm{\alpha }}_{2}{V}_{I, i}^{n}-m\right)\end{array}\right\}+\\ {c}_{1}\left\{\begin{array}{c}{D}_{11}\left({U}_{i+1}^{n-1}-2{U}_{i}^{n-1}+{U}_{i-1}^{n-1}\right)/{\left(\Delta x\right)}^{2}+\\ {U}_{i}^{n-1}\left(\frac{af}{a{U}_{i}^{n-1}+b{U}_{i}^{n-1}+c{V}_{I, i}^{n-1}}+{\alpha }_{1}{V}_{S, i}^{n-1}+{\mathrm{\alpha }}_{2}{V}_{I, i}^{n-1}-m\right)\end{array}\right\}\end{array}\right]$

${V}_{S, i}^{n+1} = {V}_{S, i}^{n+1}+\Delta t\left[\begin{array}{c}{a}_{1}\left\{\begin{array}{c}{D}_{22}\left({\stackrel{-}{V}}_{S, i+1}^{n+1}-2{\stackrel{-}{V}}_{S, i}^{n+1}+{\stackrel{-}{V}}_{S, i-1}^{n+1}\right)/{\left(\Delta x\right)}^{2}+\\ {\stackrel{-}{U}}_{i}^{n+1}\left(\frac{bf}{a{\stackrel{-}{U}}_{i}^{n+1}+b{\stackrel{-}{V}}_{S, i}^{n+1}+c{\stackrel{-}{V}}_{I, i}^{n+1}}+{\alpha }_{1}{\stackrel{-}{V}}_{S, i}^{n+1}+{\alpha }_{2}{\stackrel{-}{V}}_{I, i}^{n+1}-{n}_{1}\right)\end{array}\right\}+\\ {b}_{1}\left\{\begin{array}{c}{D}_{22}\left({V}_{S, i+1}^{n}-2{V}_{S, i}^{n}+{V}_{S, i-1}^{n}\right)/{\left(\Delta x\right)}^{2}+\\ {U}_{i}^{n}\left(\frac{bf}{a{U}_{i}^{n}+b{V}_{S, i}^{n}+c{V}_{I, i}^{n}}+{\alpha }_{1}{V}_{S, i}^{n}+{\mathrm{\alpha }}_{2}{V}_{I, i}^{n}-{n}_{2}\right)\end{array}\right\}+\\ {c}_{1}\left\{\begin{array}{c}{D}_{22}\left({V}_{S, i+1}^{n-1}-2{V}_{S, i}^{n-1}+{V}_{S, i-1}^{n-1}\right)/{\left(\Delta x\right)}^{2}+\\ {U}_{i}^{n-1}\left(\frac{bf}{a{U}_{i}^{n-1}+b{U}_{i}^{n-1}+c{V}_{I, i}^{n-1}}+{\alpha }_{1}{V}_{S, i}^{n-1}+{\mathrm{\alpha }}_{2}{V}_{I, i}^{n-1}-{n}_{2}\right)\end{array}\right\}\end{array}\right]$

${V}_{I, i}^{n+1} = {V}_{I, i}^{n}+\Delta t\left[\begin{array}{c}{a}_{1}\left\{\begin{array}{c}{D}_{33}\left({\stackrel{-}{V}}_{I, i+1}^{n+1}-2{\stackrel{-}{V}}_{I, i}^{n+1}+{\stackrel{-}{V}}_{I, i-1}^{n+1}\right)/{\left(\Delta x\right)}^{2}+\\ {\stackrel{-}{V}}_{I, i}^{n+1}\left(\frac{cf}{a{\stackrel{-}{U}}_{i}^{n+1}+b{\stackrel{-}{V}}_{S, i}^{n+1}+c{\stackrel{-}{V}}_{I, i}^{n+1}}+\mathrm{\gamma }{\stackrel{-}{V}}_{S, i}^{n+1}-{\alpha }_{2}{U}_{i}^{n+1}-{n}_{2}\right)\end{array}\right\}+\\ {b}_{1}\left\{\begin{array}{c}{D}_{33}\left({V}_{I, i+1}^{n}-2{V}_{i, i}^{n}+{V}_{I, i-1}^{n}\right)/{\left(\Delta x\right)}^{2}+\\ {V}_{I, i}^{n}\left(\frac{cf}{a{U}_{i}^{n}+b{V}_{S, i}^{n}+c{V}_{I, i}^{n}}+\mathrm{\gamma }{V}_{S, i}^{n}-{\mathrm{\alpha }}_{2}{U}_{i}^{n}-{n}_{2}\right)\end{array}\right\}+\\ {c}_{1}\left\{\begin{array}{c}{D}_{33}\left({V}_{I, i+1}^{n-1}-2{V}_{I, i}^{n-1}+{V}_{I, i-1}^{n-1}\right)/{\left(\Delta x\right)}^{2}+\\ {V}_{I, i}^{n-1}\left(\frac{cf}{a{U}_{i}^{n-1}+b{V}_{S, i}^{n-1}+c{V}_{I, i}^{n-1}}+{\alpha }_{1}{V}_{S, i}^{n-1}+{\mathrm{\alpha }}_{2}{V}_{I, i}^{n-1}-{n}_{2}\right)\end{array}\right\}\end{array}\right]$

7. Results and discussion

T he present research yielded an eco-epidemiological model for the predator-prey system. We considered two classes of prey species, i.e., susceptible and infected prey. It is supposed in this study that the transmission of infection is non-vertical. Predation activity influences the biological environment on a large scale. Here, we considered that there is no disease transmission due to predation. The populations of two classes of prey species face disturbance due to infection. and show the stability of the coexisting equilibrium point for the original and self-diffusive systems. The parameters used for the simulation in are $\bf \boldsymbol{a} = 0.45, \boldsymbol{b} = 0.64,$ $\bf \boldsymbol{c} = 0.53,$ $\bf \boldsymbol{f} = 0.53, {\boldsymbol{\alpha }}_{1} = 0.14,$ $\bf {\boldsymbol{\alpha }}_{2} = 0.27, \boldsymbol{m} = 0.56,$ $\bf {\boldsymbol{n}}_{1} = 0.55,$ $\bf {\boldsymbol{n}}_{2} = 0.49,$ $\bf \boldsymbol{\gamma } = 0.49,$ whereas, in , the values of all of the parameters are the same as in , with the diffusion coefficients $\bf {\boldsymbol{D}}_{11} = 0.10,$ $\bf {\boldsymbol{D}}_{22} = 0.20\wedge {\boldsymbol{D}}_{33} = 0.10$ . shows the solution of the system for $\bf \boldsymbol{\gamma } = 0.20, 0.30, 0.35, \wedge 0.45$ . Other parameters have the same values as in . – show time-series solutions for the self-diffusive system with $\bf {\boldsymbol{D}}_{11} = 0.01,$ $\bf {\boldsymbol{D}}_{22} = 0.06\wedge {\boldsymbol{D}}_{33} = 0.07$ for different resource values ( $\bf \boldsymbol{f} = 0.40, 0.45, \wedge 0.50$ ); the values of the other parameters were the same as in Figure 1. For Figure 3, the initial condition was (0.1, 0.1, 0.1), whereas, for all Figures 4–6, the initial condition was (0.6096, 0.2060, 0.1533). The Figures 4– results indicate that the oscillation produced in the system dies out, resulting in the stability of the species' coexistence. – present contour plots for the predator, susceptible prey and infected prey, respectively. The values of the parameters were $\bf \boldsymbol{a} = 0.91,$ $\bf \boldsymbol{b} = 0.50,$ $\bf \boldsymbol{c} = 0.55, \boldsymbol{f} = 0.91,$ $\bf {\boldsymbol{\alpha }}_{1} = 0.77,$ $\bf {\boldsymbol{\alpha }}_{2} = 0.95, \boldsymbol{m} = 0.06, {\boldsymbol{n}}_{1} = 0.85, {\boldsymbol{n}}_{2} = 0.055, \boldsymbol{\gamma } = 0.60$ . The values of the diffusion parameters were $\bf {\boldsymbol{D}}_{11} = 0.02, {\boldsymbol{D}}_{12} = 0.84, {\boldsymbol{D}}_{13} = 0.5$ , $\bf {\boldsymbol{D}}_{21} = 0.4,$ $\bf {\boldsymbol{D}}_{22} = 0.07,$ $\bf {\boldsymbol{D}}_{23} = 0$ , $\bf {\boldsymbol{D}}_{31} = 0.2,$ $\bf {\boldsymbol{D}}_{32} = 0\wedge {\boldsymbol{D}}_{33} = 0.01$ . In Figures 1 and 2, it is evident that the systems described by (9)–(11) and (23)–(25) turn to the coexistence of the fixed point, which is biologically important. Figure 3 shows the direct impact of 𝛾 on the infected prey population. When the value of 𝛾 is raised from 0.20 to 0.45, there is a significant increase in the said population. In the time-series plots, it is interesting to note that increased values of the food parameter f lead to a remarkable increase in the predator and susceptible prey populations. Still, a reverse effect is obvious on the infected population of prey species. – show a comparison among the systems presented in (9)–(11), (23)–(25) and (36)–(38). The population density of all species is remarkably affected by the inclusion of self- and cross-diffusion. The values of the parameters taken in are $\bf \boldsymbol{a} = 0.1, \boldsymbol{b} = 0.50,$ $\bf \boldsymbol{c} = 0.05, \boldsymbol{f} = 0.01,$ $\bf {\boldsymbol{\alpha }}_{1} = 0.0177,$ $\bf {\boldsymbol{\alpha }}_{2} = 0.95,$ $\bf \boldsymbol{m} = 0.06,$ $\bf {\boldsymbol{n}}_{1} = 0.185,$ $\bf {\boldsymbol{n}}_{2} = 0.055,$ $\bf \wedge \boldsymbol{\gamma } = 0.60$ . In , all of the parameters have the same value, except $\bf {\boldsymbol{\alpha }}_{1} = 0.0177$ , and the diffusion coefficients are $\bf {\boldsymbol{D}}_{11} = 0.46,$ $\bf {\boldsymbol{D}}_{22} = 0.7\wedge {\boldsymbol{D}}_{33} = 0.71$ . All of the parameters have the same values as those applied in , except $\bf {\boldsymbol{D}}_{12} = 0.31,$ $\bf {\boldsymbol{D}}_{13} = 0.75\wedge {\boldsymbol{D}}_{21} = 0.14,$ $\bf {\boldsymbol{D}}_{23} = 0$ , $\bf {\boldsymbol{D}}_{31} = 0.65\wedge {\boldsymbol{D}}_{32} = 0$ .

Figure 1. Phase portrait for the systems (9–11) at

$U\left(0\right) = 0.60, {V}_{S}\left(0\right) = 0.06, {V}_{I}\left(0\right) = 0.73$ .

DownLoad: Full-Size Img PowerPoint

Figure 2. Phase portrait for the self-diffusive system at

$U\left(0\right) = 0.53, {V}_{S}\left(0\right) = 0.0101, {V}_{I}\left(0\right) = 0.70$ .

DownLoad: Full-Size Img PowerPoint

Figure 3. Solution for

$t = 5$ and

$\gamma = 0.20, \ 0.30, \ 0.35, \ 0.45$ .

DownLoad: Full-Size Img PowerPoint

Figure 4. Time series solution for t = 3500 and

$f = 0.40$ .

DownLoad: Full-Size Img PowerPoint

Figure 5. Time series solution for

$t = 3500$ and

$f = 0.40$ .

DownLoad: Full-Size Img PowerPoint

Figure 6. Time series solution for t = 3500 and

$f = 0.40$ .

DownLoad: Full-Size Img PowerPoint

Figure 7. Contour plot for the predator.

DownLoad: Full-Size Img PowerPoint

Figure 8. Contour plot for the susceptible prey.

DownLoad: Full-Size Img PowerPoint

Figure 9. Contour plot for the infected prey.

DownLoad: Full-Size Img PowerPoint

Figure 10. Population density for non-diffusive system.

DownLoad: Full-Size Img PowerPoint

Figure 11. Population density for self-diffusive system.

DownLoad: Full-Size Img PowerPoint

Figure 12. Population density for self-diffusive system.

DownLoad: Full-Size Img PowerPoint

8. Conclusions

In the present article, we have formulated a dynamical nonlinear eco-epidemiological population model considering two species, i.e., the predator and prey, taking into account the species' dependence on the same food resources. Moreover, it is assumed that a fatal infectious, transmissible disease is present in the prey species. The equilibria of the proposed model have been found. We analyzed the global stability of the coexistence equilibrium point for the non-diffusive and diffusive cases by establishing the Lyapunov function. We have proved in Theorems 1 and 2 that the coexisting equilibrium point $({U}^{\mathrm{*}}, {V}_{S}^{\mathrm{*}}, {V}_{I}^{\mathrm{*}})$ is globally stable for the original and self-diffusive models, respectively. Also, it is proved in Theorem 3 that Turing instability of the system for the cross-diffusive case exists conditionally. Phase portraits were drawn to support the theoretical results and show the equilibrium point's stability for the non-diffusive and self-diffusive cases. Moreover, a two-stage explicit numerical scheme was constructed, and the stability of the proposed scheme was evaluated by using the von Neumann stability criterion. The impacts of parameters like the disease transmission parameter 𝛾 and food resource f have been analyzed through the use of plots. It is obvious from the plots that increasing the disease transmission parameter 𝛾 enhances the infected population density, and vice versa. The solution of the self-diffusive system for different values of the resource parameter f has been obtained, and it is shown that the increment in the food resource affects the population size.

Acknowledgments

The authors wish to express their gratitude to Prince Sultan University for facilitating the publication of this article through the Theoretical and Applied Sciences Lab. Funding Statement: The authors would like to acknowledge the support of Prince Sultan University for paying the Article Processing Charges for this publication.

Conflict of Interest

The authors declare that they have no conflict of interest to report regarding the present study

Appendix A

For the proof of (22):

$\frac{dU}{dt} = U(\frac{{\alpha }_{1}a}{b}\left({U}^{\mathrm{*}}-\mathrm{U}\right)+(\frac{{\alpha }_{1}c}{b}-{\alpha }_{2}\left)\right({V}_{I}^{\mathrm{*}}-{V}_{I}) )$

We consider

$aU+b{V}_{S}+c{V}_{I} = a{U}^{\mathrm{*}}+b{V}_{S}^{\mathrm{*}}+c{V}_{I}^{\mathrm{*}}$

$aU+b{V}_{S}+c{V}_{I} = a\left(\frac{af{\gamma }^{2}+bf\gamma {\alpha }_{2}-b\gamma m{n}_{2}-b{n}_{1}{n}_{2}{\alpha }_{2}+b{{n}_{2}}^{2}{\alpha }_{1}-cf\gamma {\alpha }_{1}+c\gamma m{n}_{1}+c{{n}_{1}}^{2}{\alpha }_{2}-c{n}_{1}{n}_{2}{\alpha }_{1}}{\left(a\gamma +b{\alpha }_{2}-c{\alpha }_{1}\right)\left(\gamma m+{n}_{1}{\alpha }_{2}-{n}_{2}{\alpha }_{1}\right)}\right)$

$+b\left(\frac{af\gamma {\alpha }_{2}+a\gamma m{n}_{2}+a{n}_{1}{n}_{2}{\alpha }_{2}-a{{n}_{2}}^{2}{\alpha }_{1}+bf{{\alpha }_{2}}^{2}-cf{\alpha }_{1}{\alpha }_{2}-c\gamma {m}^{2}-cm{n}_{1}{\alpha }_{2}+cm{n}_{2}{\alpha }_{1}}{\left(a\gamma +b{\alpha }_{2}-c{\alpha }_{1}\right)\left(\gamma m+{n}_{1}{\alpha }_{2}-{n}_{2}{\alpha }_{1}\right)}\right)$

$+c(-\frac{af\gamma {\alpha }_{1}+a\gamma m{n}_{1}+a{{n}_{1}}^{2}{\alpha }_{2}-a{n}_{1}{n}_{2}{\alpha }_{1}+bf{\alpha }_{1}{\alpha }_{2}-b\gamma {m}^{2}-bm{n}_{1}{\alpha }_{2}+bm{n}_{2}{\alpha }_{1}-cf{{\alpha }_{1}}^{2}}{\left(a\gamma +b{\alpha }_{2}-c{\alpha }_{1}\right)\left(\gamma m+{n}_{1}{\alpha }_{2}-{n}_{2}{\alpha }_{1}\right)})$

$= \frac{f\left(a\gamma +b{\alpha }_{2}-c{\alpha }_{1}\right)}{(\gamma m+{n}_{1}{\alpha }_{2}-{n}_{2}{\alpha }_{1)}}$

$b{V}_{S} = \frac{f\left(a\gamma +b{\alpha }_{2}-c{\alpha }_{1}\right)}{(\gamma m+{n}_{1}{\alpha }_{2}-{n}_{2}{\alpha }_{1)}}-aU-c{V}_{I}$

${\alpha }_{1}{V}_{S} = \frac{{\alpha }_{1}f\left(a\gamma +b{\alpha }_{2}-c{\alpha }_{1}\right)}{(\gamma m+{n}_{1}{\alpha }_{2}-{n}_{2}{\alpha }_{1)}}-\frac{{\alpha }_{1}aU}{b}-\frac{{\alpha }_{1}c{V}_{I}}{b}$

Using above value in (9) and applying simplification leads us to the following result:

$\frac{dU}{dt} = U(\frac{{\alpha }_{1}a}{b}\left({U}^{\mathrm{*}}-\mathrm{U}\right)+(\frac{{\alpha }_{1}c}{b}-{\alpha }_{2}\left)\right({V}_{I}^{\mathrm{*}}-{V}_{I}) )$

Appendix B

To prove the value of ${M}_{1} = -{D}_{11}{U}^{\mathrm{*}}{\int }_{\mathrm{\Omega }}\text{‍}\frac{1}{{U}^{2}}|\nabla {U|}^{2}dX\le 0$ in (28), we proceed as follows:

$\text{ As } \ \ \ {M}_{1} = {\int }_{\mathrm{\Omega }}\text{‍}\frac{U-{U}^{\mathrm{*}}}{U}{D}_{11}△UdX$

By using Green's first identity, we have

${\int }_{\mathrm{\Omega }}\text{‍}\frac{U-{U}^{\mathrm{*}}}{U}{D}_{11}△UdX$

$= {D}_{11} ( {\int }_{\mathrm{\Omega }}\text{‍}\frac{U-{U}^{\mathrm{*}}}{U}\frac{\partial U}{\partial \vartheta }dX - {\int }_{\mathrm{\Omega }}\text{‍}{\nabla }_{\mathrm{x}}\frac{U-{U}^{\mathrm{*}}}{U}.{\nabla }_{\mathrm{x}}UdX )$

(-)

$= - {D}_{11}{\int }_{\mathrm{\Omega }}\text{‍}{\nabla }_{\mathrm{x}}\frac{U-{U}^{\mathrm{*}}}{U}.{\nabla }_{\mathrm{x}}UdX$

$= - {D}_{11}{\int }_{\mathrm{\Omega }}\text{‍}{\nabla }_{\mathrm{x}}\left(1-\frac{{U}^{\mathrm{*}}}{U}\right){\nabla }_{\mathrm{x}}UdX$

Further simplification leads to the following:

$= -{D}_{11}{\int }_{\mathrm{\Omega }}\text{‍}\frac{{U}^{\mathrm{*}}}{{U}^{2}}|\nabla {U|}^{2}dX$

So, we have

${M}_{1} = -{D}_{11}{U}^{\mathrm{*}}{\int }_{\mathrm{\Omega }}\text{‍}\frac{1}{{U}^{2}}|\nabla {U|}^{2}dX\le 0$

Similarly, (31) and (32) can be proved.

References

[1]	E. Dong, H. Du, L. Gardner, An interactive web-based dashboard to track COVID-19 in real time, Lancet Infect. Dis., 20 (2020), 533–534.
[2]	The World Bank, World Bank open data, 2021. Available from: https://databank.worldbank.org/home.aspx.
[3]	SEGOB, Acuerdo por el que se establecen acciones extraordinarias para atender la emergencia sanitaria generada por el virus SARS-CoV2, 2020. Available from: https://dof.gob.mx/nota_detalle.php?codigo=5590914&fecha=31%2F03%2F2020.
[4]	INEGI, Censo de población y vivienda 2020, características de las localidades, tabulados predefinidos, 2021. Available from: https://www.inegi.org.mx/programas/ccpv/2020/.
[5]	I. Ibarra-Nava, J. A. Cardenas-De La Garza, R. E. Ruiz-Lozano, R. G. Salazar-Montalvo, Mexico and the COVID-19 response, Disaster Med. Public Health Prep., 14 (2020), e17–e18. https://doi.org/10.1017/dmp.2020.260 doi: 10.1017/dmp.2020.260
[6]	F. Knaul, H. Arreola-Ornelas, T. Porteny, M. Touchton, M. Sánchez-Talanquer, Ó. Méndez, et al., Not far enough: Public health policies to combat COVID-19 in Mexico's states, PLoS One, 16 (2021), e0251722. https://doi.org/10.1371/journal.pone.0251722 doi: 10.1371/journal.pone.0251722
[7]	M. Sánchez-Talanquer, E. González-Pier, J. Sepúlveda, L. Abascal-Miguel, J. Fieldhouse, C. Del Rio, et al., Mexico's Response to COVID-19: A Case Study, Institute for Global Health Sciences, University of California, San Francisco, 2021. Available from: https://globalhealthsciences.ucsf.edu/sites/globalhealthsciences.ucsf.edu/files/mexico-covid-19-case-study-english.pdf.
[8]	I. Franch-Pardo, B. M. Napoletano, F. Rosete-Verges, L. Billa, Spatial analysis and GIS in the study of COVID-19. A review, Sci. Total Environ., 739 (2020), 140033. https://doi.org/10.1016/j.scitotenv.2020.140033 doi: 10.1016/j.scitotenv.2020.140033
[9]	C. Scarpone, S. T. Brinkmann, T. Große, D. Sonnenwald, M.Fuchs, B. B. Walker, A multimethod approach for county-scale geospatial analysis of emerging infectious diseases: a cross-sectional case study of COVID-19 incidence in Germany, Int. J. Health Geographics, 19 (2020), 1–17. https://doi.org/10.1186/s12942-020-00225-1 doi: 10.1186/s12942-020-00225-1
[10]	N. Arinaminpathy, J. Das, T. H. McCormick, P. Mukhopadhyay, N. Sircar, Quantifying heterogeneity in SARS-CoV-2 transmission during the lockdown in India, Epidemics, 36 (2021), 100477. https://doi.org/10.1101/2020.09.10.20190017 doi: 10.1101/2020.09.10.20190017
[11]	Y. Ge, Z. Chen, A. Handel, L. Martinez, Q. Xiao, C. Li, et al., The impact of social distancing, contact tracing, and case isolation interventions to suppress the COVID-19 epidemic: A modeling study, Epidemics, 36 (2021), 100483. https://doi.org/10.1016/j.epidem.2021.100483 doi: 10.1016/j.epidem.2021.100483
[12]	S. Triambak, D. P. Mahapatra, N. Mallick, R. Sahoo, A new logistic growth model applied to COVID-19 fatality data, Epidemics, (2021), 100515. https://doi.org/10.1016/j.epidem.2021.100515 doi: 10.1016/j.epidem.2021.100515
[13]	P. Pincheira-Brown, A. Bentancor, Forecasting COVID-19 infections with the semi-unrestricted generalized growth model, Epidemics, 37 (2021), 100486. https://doi.org/10.1016/j.epidem.2021.100486 doi: 10.1016/j.epidem.2021.100486
[14]	M. Fatima, K. J. O'Keefe, W. Wei, S. Arshad, O. Gruebner, Geospatial analysis of COVID-19: A scoping review, Int. J. Environ. Res. Public Health, 18 (2021), 2336. https://doi.org/10.3390/ijerph18052336 doi: 10.3390/ijerph18052336
[15]	S. He, S. Tang, L. Rong, A discrete stochastic model of the COVID-19 outbreak: Forecast and control, Math. Biosci. Eng., 17 (2020), 2792–2804. https://doi.org/10.3934/mbe.2020153 doi: 10.3934/mbe.2020153
[16]	A. K. Srivastav, P. K. Tiwari, P. K. Srivastava, M. Ghosh, Y. Kang, A mathematical model for the impacts of face mask, hospitalization and quarantine on the dynamics of COVID-19 in India: deterministic vs. stochastic, Math. Biosci. Eng., 18 (2021), 182–213. https://doi.org/10.3934/mbe.2021010 doi: 10.3934/mbe.2021010
[17]	D. Olabode, J. Culp, A. Fisher, A. Tower, D. Hull-Nye, X. Wang, Deterministic and stochastic models for the epidemic dynamics of COVID-19 in Wuhan, China, Math. Biosci. Eng., 18 (2021), 950–967. https://doi.org/10.3934/mbe.2021050 doi: 10.3934/mbe.2021050
[18]	T. S. Brown, K. Engø-Monsen, M. V. Kiang, A. S. Mahmud, R. J. Maude, C. O. Buckee, The impact of mobility network properties on predicted epidemic dynamics in Dhaka and Bangkok, Epidemics, 35 (2021), 100441. https://doi.org/10.1016/j.epidem.2021.100441 doi: 10.1016/j.epidem.2021.100441
[19]	W. Chiang, X. Liu, G. Mohler, Hawkes process modeling of COVID-19 with mobility leading indicators and spatial covariates, Int. J. Forecast., 38 (2021), 505–520. https://doi.org/10.1016/j.ijforecast.2021.07.001 doi: 10.1016/j.ijforecast.2021.07.001
[20]	R. Engbert, M. M. Rabe, R. Kliegl, S. Reich, Sequential data assimilation of the stochastic SEIR epidemic model for regional COVID-19 dynamics, Bull. Math. Biol., 83 (2021), 1–16. https://doi.org/10.1007/s11538-020-00834-8 doi: 10.1007/s11538-020-00834-8
[21]	V. Gopalakrishnan, S. Pethe, S. Kefayati, S. Raman, P. Hake, A. Deshpande, et al., Globally local: Hyper-local modeling for accurate forecast of COVID-19, Epidemics, 37 (2021), 100510. https://doi.org/10.1016/j.epidem.2021.100510 doi: 10.1016/j.epidem.2021.100510
[22]	R. Huang, M. Liu, Y. Ding, Spatial-temporal distribution of COVID-19 in China and its prediction: A data-driven modeling analysis, J. Infect. Dev. Countries, 14 (2020), 246–253. https://doi.org/10.3855/jidc.12585 doi: 10.3855/jidc.12585
[23]	Y. Mammeri, A reaction-diffusion system to better comprehend the unlockdown: Application of SEIR-type model with diffusion to the spatial spread of COVID-19 in France, Comput. Math. Biophys., 8 (2020), 102–113. https://doi.org/10.1515/cmb-2020-0104 doi: 10.1515/cmb-2020-0104
[24]	X. Hou, S. Gao, Q. Li, Y. Kang, N. Chen, K. Chen, et al., Intracounty modeling of COVID-19 infection with human mobility: Assessing spatial heterogeneity with business traffic, age, and race, Proc. Natl. Acad. Sci., 118 (2021). https://doi.org/10.1073/pnas.2020524118 doi: 10.1073/pnas.2020524118
[25]	T. W. Alleman, J. Vergeynst, L. De Visscher, M. Rollier, E. Torfs, I. Nopens, et al., Assessing the effects of non-pharmaceutical interventions on SARS-CoV-2 transmission in Belgium by means of an extended SEIQRD model and public mobility data, Epidemics, 37 (2021), 100505. https://doi.org/10.1016/j.epidem.2021.100505 doi: 10.1016/j.epidem.2021.100505
[26]	R. Pastor-Satorras, C. Castellano, P. Van Mieghem, A. Vespignani, Epidemic processes in complex networks, Rev. Mod. Phys., 87 (2015), 925. https://doi.org/10.1103/RevModPhys.87.925 doi: 10.1103/RevModPhys.87.925
[27]	Y. Wang, J. Cao, Final size of network epidemic models: properties and connections, Sci. China Inf. Sci, 64 (2021), 179201. https://doi.org/10.1007/s11432-019-2656-2 doi: 10.1007/s11432-019-2656-2
[28]	D. O'Sullivan, M. Gahegan, D. J. Exeter, B. Adams, Spatially explicit models for exploring COVID-19 lockdown strategies, Trans. GIS, 24 (2020), 967–1000. https://doi.org/10.1111/tgis.12660 doi: 10.1111/tgis.12660
[29]	D. Laroze, E. Neumayer, T. Plümper, COVID-19 does not stop at open borders: Spatial contagion among local authority districts during England's first wave, Social Sci. Med., 270 (2021). https://doi.org/10.1016/j.socscimed.2020.113655 doi: 10.1016/j.socscimed.2020.113655
[30]	I. Franch-Pardo, M. R. Desjardins, I. Barea-Navarro, A. Cerdà, A review of GIS methodologies to analyze the dynamics of COVID-19 in the second half of 2020, Trans. GIS, 25 (2021), 2191–2239. https://doi.org/10.1111/tgis.12792 doi: 10.1111/tgis.12792
[31]	A. Gandolfi, Percolation methods for SEIR epidemics on graphs, in Dynamic Models of Infectious Diseases, Springer, New York, NY, (2013), 31–58. https://doi.org/10.1007/978-1-4614-9224-5_2
[32]	R. Wallace, L. Bergmann, L. Hogerwerf, R. Kock, R. G. Wallace, Ebola in the hog sector: Modeling pandemic emergence in commodity livestock, in Neoliberal Ebola, Springer, Cham, (2016), 13–53. https://doi.org/10.1007/978-3-319-40940-5_2
[33]	R. Wallace, D. Wallace, H. Andrews, R. Fullilove, M. T. Fullilove, The spatiotemporal dynamics of AIDS and TB in the New York metropolitan region from a sociogeographic perspective: Understanding the linkages of central city and suburbs, Environ. Plann. A, 27 (1995), 1085–1108. https://doi.org/10.1068/a271085 doi: 10.1068/a271085
[34]	U. N. Raghavan, R. Albert, S. Kumara, Near linear time algorithm to detect community structures in large-scale networks, Phys. Rev. E, 76 (2007), 036106. https://doi.org/10.1103/PhysRevE.76.036106 doi: 10.1103/PhysRevE.76.036106
[35]	M. Martcheva, An introduction to Mathematical Epidemiology, Springer, 61 (2015). https://doi.org/10.1007/978-1-4899-7612-3
[36]	Center for Disease Control and Prevention, Scientific brief: SARS-CoV-2 transmission, 2022. Available from: https://www.cdc.gov/coronavirus/2019-ncov/science/science-briefs/sars-cov-2-transmission.html.
[37]	J. M. Brauner, S. Mindermann, M. Sharma, D. Johnston, J. Salvatier, T. Gavenčiak, et al., Inferring the effectiveness of government interventions against COVID-19, Science, 371 (2021). https://doi.org/10.1126/science.abd9338 doi: 10.1126/science.abd9338
[38]	N. Imai, K. A. M. Gaythorpe, S. Abbott, S. Bhatia, S. van Elsland, K. Prem, et al., Adoption and impact of non-pharmaceutical interventions for COVID-19, Wellcome Open Res., 5 (2020). https://doi.org/10.12688/wellcomeopenres.15808.1 doi: 10.12688/wellcomeopenres.15808.1
[39]	N. L. Komarova, A. Azizi, D. Wodarz, Network models and the interpretation of prolonged infection plateaus in the COVID19 pandemic, Epidemics, 35 (2021), 100463. https://doi.org/10.1016/j.epidem.2021.100463 doi: 10.1016/j.epidem.2021.100463
[40]	D. Pradhan, P. Biswasroy, P. Kumar Naik, G. Ghosh, G. Rath, A review of current interventions for COVID-19 prevention, Arch. Med. Res., 51 (2020), 363–374. https://doi.org/10.1016/j.arcmed.2020.04.020 doi: 10.1016/j.arcmed.2020.04.020
[41]	K. Soltesz, F. Gustafsson, T. Timpka, J. Jaldén, C. Jidling, A. Heimerson, et al., The effect of interventions on COVID-19, Nature, 588 (2020), E26–E28. https://doi.org/10.1038/s41586-020-3025-y doi: 10.1038/s41586-020-3025-y
[42]	P. Bacchetti, N. P. Jewell, Nonparametric estimation of the incubation period of AIDS based on a prevalent cohort with unknown infection times, Biometrics, 47 (1991), 947–960. https://doi.org/10.2307/2532651 doi: 10.2307/2532651
[43]	B. Cazelles, S. Hales, Infectious diseases, climate influences, and nonstationarity, PLoS Med., 3 (2006), e328. https://doi.org/10.1371/journal.pmed.0030328 doi: 10.1371/journal.pmed.0030328
[44]	J. Lu, A. Lin, C. Jiang, A. Zhang, Z. Yang, Influence of transportation network on transmission heterogeneity of COVID-19 in China, Transp. Res. Part C Emerging Technol., 129 (2021), 103231. https://doi.org/10.1016/j.trc.2021.103231 doi: 10.1016/j.trc.2021.103231
[45]	F. Brauer, P. V. den Driessche, J. Wu, Mathematical Epidemiology, Berlin, Germany, Springer, 1945 (2008). https://doi.org/10.1007/978-3-540-78911-6
[46]	J. Chen, S. Levin, S. Eubank, H. Mortveit, S. Venkatramanan, A. Vullikanti, et al., Networked epidemiology for COVID-19, SIAM News, 53 (2020). Available from: https://par.nsf.gov/servlets/purl/10213734.
[47]	M. Marathe, A. K. S. Vullikanti, Computational epidemiology, Commun. ACM, 56 (2013), 88–96. https://doi.org/10.1145/2483852.2483871 doi: 10.1145/2483852.2483871
[48]	W. O. Kermack, A. G. McKendrick, A contribution to the mathematical theory of epidemics, Proc. R. Soc. London, Ser. A, 115 (1927), 700–721. https://doi.org/10.1098/rspa.1927.0118 doi: 10.1098/rspa.1927.0118
[49]	F. Arrigoni, A. Pugliese, Limits of a multi-patch SIS epidemic model, J. Math. Biol., 45 (2002), 419–440. https://doi.org/10.1007/s002850200155 doi: 10.1007/s002850200155
[50]	C. Wolf, M. Langlais, D. Ponter, A multi-patch epidemic model with periodic demography, direct and indirect transmission and variable maturation rate, Math. Popul. Stud., 13 (2006), 153–177. https://doi.org/10.1080/08898480600788584 doi: 10.1080/08898480600788584
[51]	J. Arino, J. M. Miller, P. Van den Driessche, A multi-species epidemic model with spatial dynamics, Math. Med. Biol., 22 (2005), 129–142. https://doi.org/10.1093/imammb/dqi003 doi: 10.1093/imammb/dqi003
[52]	R. K. McCormack, L. J. S. Allen, Multi-patch deterministic and stochastic models for wildlife diseases, J. Biol. Dyn., 1 (2007), 63–85. https://doi.org/10.1080/17513750601032711 doi: 10.1080/17513750601032711
[53]	D. Bichara, A. Iggidr, Multi-patch and multi-group epidemic models: a new framework, J. Math. Biol, 77 (2018), 107–134. https://doi.org/10.1007/s00285-017-1191-9 doi: 10.1007/s00285-017-1191-9
[54]	M. Girvan, M. E. J. Newman, Community structure in social and biological networks, Proc. Natl. Acad. Sci., 99 (2002), 7821–7826. https://doi.org/10.1073/pnas.122653799 doi: 10.1073/pnas.122653799
[55]	F. Parés, D. G. Gasulla, A. Vilalta, J. Moreno, E. Ayguadé, J. Labarta, et al., Fluid communities: A competitive, scalable and diverse community detection algorithm, in International Conference on Complex Networks and Their Applications, Springer, Cham, (2017), 229–240. https://doi.org/10.1007/978-3-319-72150-7_19
[56]	S. Fortunato, Community detection in graphs, Phys. Rep., 486 (2010), 75–174. https://doi.org/10.1016/j.physrep.2009.11.002 doi: 10.1016/j.physrep.2009.11.002
[57]	A. Clauset, M. E. J. Newman, C. Moore, Finding community structure in very large networks, Phys. Rev. E, 70 (2004), 066111. https://doi.org/10.1103/PhysRevE.70.066111 doi: 10.1103/PhysRevE.70.066111
[58]	M. Bastian, S. Heymann, M. Jacomy, Gephi: An open source software for exploring and manipulating networks, in Proceedings of the Third International Conference on Weblogs and Social Media, 2009. https://doi.org/10.13140/2.1.1341.1520
[59]	R. Wallace, Y. S. Huang, P. Gould, D. Wallace, The hierarchical diffusion of AIDS and violent crime among US metropolitan regions: inner-city decay, stochastic resonance and reversal of the mortality transition, Social Sci. Med., 44 (1997), 935–947. https://doi.org/10.1016/s0277-9536(96)00197-9 doi: 10.1016/s0277-9536(96)00197-9
[60]	Secretaría de Salud Jalisco, Datos de Radar Jalisco, 2021. Available from: https://coronavirus.jalisco.gob.mx/datos-abiertos.
[61]	H. Gutiérrez-Pulido, D. Gómez-Alvarez, Estudio técnico para la adecuación de las regiones de Jalisco, Technical report, Gobierno de Jalisco, 2016. Available from: https://seplan.app.jalisco.gob.mx/biblioteca/archivo/verDocumento/939.
[62]	K. Ejima, K. S. Kim, C. Ludema, A. I. Bento, S. Iwanami, Y. Fujita, et al., Estimation of the incubation period of COVID-19 using viral load data, Epidemics, 35 (2021), 100454. https://doi.org/10.1016/j.epidem.2021.100454 doi: 10.1016/j.epidem.2021.100454
[63]	H. Najafimehr, K. M. Ali, S. Safari, M. Yousefifard, M. Hosseini, Estimation of basic reproduction number for COVID-19 and the reasons for its differences, Int. J. Clin. Pract., 74 (2020), e13518. https://doi.org/10.1111/ijcp.13518 doi: 10.1111/ijcp.13518
[64]	W. He, G. Y. Yi, Y. Zhu, Estimation of the basic reproduction number, average incubation time, asymptomatic infection rate, and case fatality rate for COVID-19: Meta-analysis and sensitivity analysis, J. Med. Virol., 92 (2020), 2543–2550. https://doi.org/10.1002/jmv.26041 doi: 10.1002/jmv.26041
[65]	S. Zhao, B. Tang, S. S. Musa, S. Ma, J. Zhang, M. Zeng, et al., Estimating the generation interval and inferring the latent period of COVID-19 from the contact tracing data, Epidemics, 36 (2021), 100482. https://doi.org/10.1016/j.epidem.2021.100482 doi: 10.1016/j.epidem.2021.100482
[66]	J. A. Backer, D. Klinkenberg, J. Wallinga, Incubation period of 2019 novel coronavirus (2019-NCOV) infections among travellers from Wuhan, China, 20–28 January 2020, Eurosurveillance, 25 (2020), 2000062. https://doi.org/10.2807/1560-7917.ES.2020.25.5.2000062 doi: 10.2807/1560-7917.ES.2020.25.5.2000062
[67]	R. J. Hyndman, A. B. Koehler, Another look at measures of forecast accuracy, Int. J. Forecast., 22 (2006), 679–688. https://doi.org/10.1016/j.ijforecast.2006.03.001 doi: 10.1016/j.ijforecast.2006.03.001

mbe-19-04-196-supplementary.pdf

This article has been cited by:

Skander Hachicha, Najmeddine Attia, Probabilistic approaches to exploring Binet's type formula for the Tribonacci sequence, 2025, 10, 2473-6988, 11957, 10.3934/math.2025540

Reader Comments

Your name:*

Email:*
© 2022 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)

通讯作者: 陈斌, bchen63@163.com

1.
沈阳化工大学材料科学与工程学院沈阳 110142

Mathematical Biosciences and Engineering

4.4

Metrics

Article views(3870) PDF downloads(130) Cited by(1)

Preview PDF

Download XML

Export Citation

Article outline

Show full outline

Figures and Tables

Figures(4)

Mathematical Biosciences and Engineering

Geographical network model for COVID-19 spread among dynamic epidemic regions

Related Papers:

Abstract

1. Introduction

2. Model formulation

3. Stability for the predator-prey system of equations

4. Global stability for a diffusive system

4.1. Global stability for the self-diffusive system

4.2. Induced stability for the cross-diffusive system

5. Construction of numerical scheme

5.1. Stability analysis

6. Numerical simulations

7. Results and discussion

8. Conclusions

Acknowledgments

Conflict of Interest

Appendix A

Appendix B

References

This article has been cited by:

Reader Comments

通讯作者: 陈斌, bchen63@163.com

Metrics

Figures and Tables

Other Articles By Authors

Catalog

Mathematical Biosciences and Engineering

Geographical network model for COVID-19 spread among dynamic epidemic regions

Related Papers:

Abstract

1. Introduction

2. Model formulation

3. Stability for the predator-prey system of equations

4. Global stability for a diffusive system

4.1. Global stability for the self-diffusive system

4.2. Induced stability for the cross-diffusive system

5. Construction of numerical scheme

5.1. Stability analysis

6. Numerical simulations

7. Results and discussion

8. Conclusions

Acknowledgments

Conflict of Interest

Appendix A

Appendix B

References

This article has been cited by:

Reader Comments

通讯作者: 陈斌, bchen63@163.com

Metrics

Figures and Tables

Other Articles By Authors

Related pages

Tools

Export File

Citation

Format

Content

Catalog