Loading [MathJax]/jax/output/SVG/jax.js

Mathematical modelling of a mushy region formation during sulphation of calcium carbonate

  • The subject of the present paper is the derivation and asymptotic analysis of a mathematical model for the formation of a mushy region during sulphation of calcium carbonate. The model is derived by averaging, with the use of the multiple scales method, applied on microscopic moving - boundary problems. The latter problems describe the transformation of calcium carbonate into gypsum on the microscopic scale. The derived macroscopic model is solved numerically with the use of a finite element method. The results of some simulations and a relevant discussion are also presented.

    Citation: Christos V. Nikolopoulos. Mathematical modelling of a mushy region formation during sulphation of calcium carbonate[J]. Networks and Heterogeneous Media, 2014, 9(4): 635-654. doi: 10.3934/nhm.2014.9.635

    Related Papers:

    [1] Chang-Yuan Cheng, Shyan-Shiou Chen, Xingfu Zou . On an age structured population model with density-dependent dispersals between two patches. Mathematical Biosciences and Engineering, 2019, 16(5): 4976-4998. doi: 10.3934/mbe.2019251
    [2] Yun Kang, Sourav Kumar Sasmal, Komi Messan . A two-patch prey-predator model with predator dispersal driven by the predation strength. Mathematical Biosciences and Engineering, 2017, 14(4): 843-880. doi: 10.3934/mbe.2017046
    [3] Tingting Zheng, Huaiping Zhu, Zhidong Teng, Linfei Nie, Yantao Luo . Patch model for border reopening and control to prevent new outbreaks of COVID-19. Mathematical Biosciences and Engineering, 2023, 20(4): 7171-7192. doi: 10.3934/mbe.2023310
    [4] Fen-fen Zhang, Zhen Jin . Effect of travel restrictions, contact tracing and vaccination on control of emerging infectious diseases: transmission of COVID-19 as a case study. Mathematical Biosciences and Engineering, 2022, 19(3): 3177-3201. doi: 10.3934/mbe.2022147
    [5] Nancy Azer, P. van den Driessche . Competition and Dispersal Delays in Patchy Environments. Mathematical Biosciences and Engineering, 2006, 3(2): 283-296. doi: 10.3934/mbe.2006.3.283
    [6] Cheng-Cheng Zhu, Jiang Zhu, Xiao-Lan Liu . Influence of spatial heterogeneous environment on long-term dynamics of a reaction-diffusion SVIR epidemic model with relaps. Mathematical Biosciences and Engineering, 2019, 16(5): 5897-5922. doi: 10.3934/mbe.2019295
    [7] Qianqian Cui, Zhipeng Qiu, Ling Ding . An SIR epidemic model with vaccination in a patchy environment. Mathematical Biosciences and Engineering, 2017, 14(5&6): 1141-1157. doi: 10.3934/mbe.2017059
    [8] Maria Guadalupe Vazquez-Peña, Cruz Vargas-De-León, Jorge Velázquez-Castro . Global stability for a mosquito-borne disease model with continuous-time age structure in the susceptible and relapsed host classes. Mathematical Biosciences and Engineering, 2024, 21(11): 7582-7600. doi: 10.3934/mbe.2024333
    [9] Ali Mai, Guowei Sun, Lin Wang . The impacts of dispersal on the competition outcome of multi-patch competition models. Mathematical Biosciences and Engineering, 2019, 16(4): 2697-2716. doi: 10.3934/mbe.2019134
    [10] Minjuan Gao, Lijuan Chen, Fengde Chen . Dynamical analysis of a discrete two-patch model with the Allee effect and nonlinear dispersal. Mathematical Biosciences and Engineering, 2024, 21(4): 5499-5520. doi: 10.3934/mbe.2024242
  • The subject of the present paper is the derivation and asymptotic analysis of a mathematical model for the formation of a mushy region during sulphation of calcium carbonate. The model is derived by averaging, with the use of the multiple scales method, applied on microscopic moving - boundary problems. The latter problems describe the transformation of calcium carbonate into gypsum on the microscopic scale. The derived macroscopic model is solved numerically with the use of a finite element method. The results of some simulations and a relevant discussion are also presented.


    For some infectious diseases, recovered individuals may relapse after some time in recovery class, reverting them back into the infectious class. Actually, such recurrence of disease is an important feature of some animal and human diseases, for example, tuberculosis, including human and bovine [2,1], and herpes [2,3,4]. In general, a recovered individual may or may not relapse; and in the former case, the relapse time varies from individuals to individuals, following certain type of distributions. In order to describe the above mentioned individual variance, van den Driessche and Zou [4] proposed an approach in the form of integro-differential equations involving a probability function to track the recovered individuals and their possible relapses. To briefly review this approach, we let $ I(t) $ denote the population of infectious class and $ \gamma $ denote the recovery rate. Considering that the relapse times for recovered individuals may differ from individual to individual, a function $ P(t) $ is introduced in [4] to denote the probability that a recovered individual still remains in the recovered class $ t $ time units after recovery. By the meaning of $ P(t) $, it is then assumed to satisfy the following property:

    (A) $ P: [0, \infty) \to [0, \infty) $ is differentiable except at possibly finite many points where it may have jump discontinuities, non-increasing and satisfies $ P(0) = 1, \, \, \, P(\infty) = 0 $ and $ \int_0^\infty P(t)\, dt $ positive and finite.

    Then, the population of the recovered class at time $ t $ is given by

    $ R(t)=t0γI(ξ)ed(tξ)P(tξ)dξ $ (1.1)

    where $ d $ is the death rate of the recovered class.

    In order to fit into a general differential equation model for a disease that relapse, we differentiate (1.1) with respect to $ t $ to obtain

    $ R(t)=γI(t)dR(t)+t0γI(ξ)ed(tξ)dtP(tξ)dξ. $ (1.2)

    Here the first term represents the new entry into the recovered class from infectious class and the second term explains the deaths of the recovered individuals, while the third term is nothing but the rate at which recovered individuals revert into infectious class. Thus, fitting (1.2) into a model with susceptible population $ S(t) $ and infectious population $ I(t) $ leads to the following model system with general distribution for the relapse time reflected by the probability function $ P(t) $:

    $ {S(t)=KdS(t)λS(t)I(t),I(t)=λS(t)I(t)dI(t)γI(t)t0γI(ξ)ed(tξ)dtP(tξ)dξ,R(t)=γI(t)dR(t)+t0γI(ξ)ed(tξ)dtP(tξ)dξ. $ (1.3)

    Here, a simple demographic dynamics $ S'(t) = K-dS(t) $ and a mass action infection mechanism $ \lambda S(t) I(t) $ are adopted.

    We would particularly mention two special forms of $ P(t) $:

    (I) exponential decay function, i.e., $ P(t) = e^{-rt} $ for $ t\ge 0 $ where $ r > 0 $ is a constant;

    (II) step function, i.e., $ P(t) = 1 $ for $ t \in [0, \tau) $ and $ P(t) = 0 $ for $ t\ge \tau $, where $ \tau > 0 $ is a constant.

    We remark that the choice (II) is a reasonable choice for those diseases for which recovered individuals have a relatively concentrated relapse time which is approximated by $ \tau > 0 $.

    With choice (I), the integral term in (1.2) becomes

    $ \int_0^t \gamma I(\xi) e^{-d(t-\xi)} d_t P(t-\xi) d\xi = -r R(t), $

    and accordingly, (1.3) reduces to

    $ {S(t)=KdS(t)λS(t)I(t),I(t)=λS(t)I(t)dI(t)γI(t)+rR(t),R(t)=γI(t)(d+r)R(t). $ (1.4)

    With choice (II), $ H(t): = 1-P(t) $ is the Heaviside function at $ \tau $ whose derivative is the Dirac delta function at $ \tau $, i.e., $ H'(t) = \delta(t-\tau) $. Hence $ d_t P(t-\xi) = \, - \delta (t-\xi -\tau) $ and therefore, the integral in (1.2) becomes

    $ \int_0^t \gamma I(\xi) e^{-d(t-\xi)} d_t P(t-\xi) d\xi = -\gamma I(t-\tau) e^{-d\tau}, $

    and accordingly, (1.3) splits to

    $ {S(t)=KdS(t)λS(t)I(t),I(t)=λS(t)I(t)dI(t)γI(t),R(t)=γI(t)dR(t).fort[0,τ], $ (1.5)

    and

    $ {S(t)=KdS(t)λS(t)I(t),I(t)=λS(t)I(t)dI(t)γI(t)+edτγI(tτ),R(t)=γI(t)dR(t)edτγI(tτ).fort>τ. $ (1.6)

    On the other hand, the world is highly connected nowadays, and travels between different regions/cities are more and more frequent and common. In order to model the transmission dynamics of infectious diseases, patch models are typically used. There have been plenty of patch models for transmission dynamics of diseases of various types in the literature, including SI, SIS, SIR, SEIR types and even vector-borne diseases. See, e.g., [5, 6, 7, 8, 9] and the references there in; particularly the more recent works [10,11] which contains more recent references on patch models for diseases dynamics. For our concerns in this paper, because of the travels of human beings, an individual recovered from an infectious disease in one region/city may be in another region when he reverts back into the infectious class. Thus, in order to describe the transmission dynamics of a disease over $ n $ patches (e.g., regions or cities) that may relapse, one cannot simply add dispersion terms in the set of $ n $ subsystems of the forms of (1.3) (or (1.4) or (1.6)) indexed by $ i $ with $ i = 1, 2, \cdots, n $, as in the aforementioned references. Instead, one needs to carefully track the dispersals of the recovered individuals among all matches to accurately evaluate the reverting rate at each patch, and this would lead to a phenomenon of "non-locality". For this purpose, the approach reviewed above faces a big challenge, if not impossible. Hence, it seems that an alternative approach needs to be sought to achieve the aforementioned goal.

    This work is motivated by [9,12,13,14] for the notion of "non-locality". In [14], based on the basic McKendrick-von Foerster equation with the structure variable being the natural age and assuming the immature individuals may disperse between patches, a non-local population dynamics model is derived. In [9,12,13], adopting the infection age as the structure variable in the McKendrick-von Foerster equation, some patch models for transmission dynamics of infections diseases with a fixed latency are derived and explored. In this paper, we will follow the framework in [9,12,13,14] but using the recovery age in the McKendrick-von Foerster equation to derive a patch model with non-local reverting in each patch, meaning that the reverting rate in each patch is actually a result of dispersals of the individuals recovered in all patches. We point out that the notion of "recovery age" is also used in [15] to derive a non-patch model for influenza disease; while in [16], another structure variable, immunity level, which is similar to the recovery age, is introduced to track the rate of recovered individuals returning to the susceptible class to derive a non-patch model.

    The rest of the paper is organized as follows. In the next section, we present the model formulation for a two-patch environment. Section 3 is devoted to confirming the well-posedness of the model obtained in Section 2. In Sections 4 and 5, we deal with the situation when all dispersal rate matrices are irreducible, in which the disease extinction/persistence, existence and stability of the disease-free equilibrium and endemic equilibrium are analysed. In Section 6, we are concerned with the situation when the irreducibility of the travel matrices for recovered class and infectious class does not hold. We only consider three special cases, which enable us to obtain more detailed results on the joint impact of relapse time and the mobility of the individuals. In Section 7, we summarize the main results of the paper, and discuss some implications of the mathematical results.

    In order to avoid the main ideas to be hidden from the complicated notations, we only consider a two patch environment. Consider a population that lives in the two patches (e.g., cities). Let $ S_i (t) $, $ I_i (t) $, $ R_i (t) $ be the sub-populations of the susceptible, infectious and recovered classes on patch $ i $, $ i = 1, \; 2 $ at time $ t $, respectively. These two patches are connected in the sense that individuals can disperse between these two patches. To track the dispersals of the recovered individuals during the recovered period, we denote recovery age by $ a $, which is the time elapsed since recovery.

    Let $ r_i (t, a) $ be the density (with respect to the recovery age $ a $) of recovered individuals at time $ t $ in patch $ i \; (i = 1, \; 2) $. We assume that all recovered individual relapse to infected class at the recover age $ a = \tau $. This assumption is in the line of choice (II) for the probability function $ P(t) $. We admit that, strictly speaking, this assumption is not that realistic, however, it makes our main idea mathematically trackable and the resulting model workable so that we are able to explore, to some extent, the impact of travels on the disease dynamics. With this assumption the total number of recovered individuals in patch $ i $ is then given by

    $ Ri(t)=τ0ri(t,a)da. $ (2.1)

    Similar to the equation governing the evolution of a population with natural age structure (see [24]), the densities $ r_i (t, a) $, $ i = 1, \; 2 $ are described by the following system of first-order partial differential equations:

    $ {r1(t,a)t+r1(t,a)a=d1r1(t,a)+DR12r2(t,a)DR21r1(t,a),r2(t,a)t+r2(t,a)a=d2r2(t,a)+DR21r1(t,a)DR12r2(t,a),t>0,a(0,τ]. $ (2.2)

    Here $ D^R_{i j}\, r_j (t, a) $ corresponds to the dispersal of the recovered individuals at the recovery age $ a $ from patch $ j $ to patch $ i $; constant $ d_i > 0 $ denotes the natural death rate in patch $ i $ which is independent of the recovery age and the disease status. In addition, we assume that there is no delay in the dispersal between patches and there is no loss during migration from patch $ j $ to patch $ i $, that is, all of those who leave patch $ j $ for patch $ i $ arrive at patch $ i $ safely. Without loss of generality, we set $ t = 0 $ to be the time when the disease epidemics starts, and hence, initially there is no recovered individuals, meaning that $ r_i(0, a) = 0 $ for $ a \in (0, \tau] $.

    Obviously, $ r_i (t, 0) $ corresponds to the new recovery individuals in patch $ i $ who come from the infectious individuals. Assuming a constant recovery rate $ \gamma_i > 0 $ in patch $ i $, we then have

    $ ri(t,0)=γiIi(t). $ (2.3)

    Now, integrating (2.2) with respect to $ a $ from 0 to $ \tau $ and making use of (2.1) and (2.3) leads to

    $ {dR1(t)dt=d1R1(t)+DR12R2(t)DR21R1(t)+γ1I1(t)r1(t,τ),dR2(t)dt=d2R2(t)+DR21R1(t)DR12R2(t)+γ2I2(t)r2(t,τ). $ (2.4)

    As in (1.3) as well as in [13], we adopt the simplest demographic structure of the population under consideration, in which we assume that there is a constant recruitment of susceptible individuals denoted by $ K_i > 0 $ in patch $ i $, $ i = 1, \; 2 $, and a constant natural death rate for each class denoted still by $ d_i > 0 $ and assume that the disease does not transmit vertically. With these assumptions, the disease dynamics can be described by the following equations:

    $ {dS1(t)dt=K1d1S1(t)+DS12S2(t)DS21S1(t)λ1S1(t)I1(t),dS2(t)dt=K2d2S2(t)+DS21S1(t)DS12S2(t)λ2S2(t)I2(t),dI1(t)dt=d1I1(t)+DI12I2(t)DI21I1(t)+λ1S1(t)I1(t)γ1I1(t)+r1(t,τ),dI2(t)dt=d2I2(t)+DI21I1(t)DI12I2(t)+λ2S2(t)I2(t)γ2I2(t)+r2(t,τ),dR1(t)dt=d1R1(t)+DR12R2(t)DR21R1(t)+γ1I1(t)r1(t,τ),dR2(t)dt=d2R2(t)+DR21R1(t)DR12R2(t)+γ2I2(t)r2(t,τ), $ (2.5)

    where $ D^S_{i j} $ is the rate at which susceptible individuals migrate from patch $ j $ to patch $ i $, $ i \neq j $, and $ D^I_{i j} $ is the rate at which infectious individuals migrate from patch $ j $ to patch $ i $, $ i \neq j $. Here $ \lambda_i > 0 $, $ i = 1, \; 2 $ are the transmission rate in patch $ i $. Note that the equations for the recovered class $ R_i $, $ i = 1, \; 2 $ are decoupled from the equations for $ S_i $ and $ I_i $, $ i = 1, \; 2 $. Thus we only need to consider the 4 equations for $ S_i $ and $ I_i $, $ i = 1, \; 2 $ in (2.5).

    Obviously, $ r_i (t, \tau) $ is the rate at which patch $ i $ gains relapse individuals, which can be determined below in terms of $ I_j (t) $ for $ j = 1, \; 2 $.

    For fixed $ \xi\geq0 $, let

    $ V^\xi_i (t) = r_i (t, t-\xi), \; \; \; for \; \; \xi\leq t \leq \xi+\tau\; \; \; and\; \; i = 1, \; 2. $

    Then for $ 1 \leq i \neq j \leq 2 $,

    $ ddtVξi(t)=tri(t,a)a=tξ+ari(t,a)a=tξ=diri(t,tξ)αi(tξ)ri(t,tξ)+DRijrj(t,tξ)DRjiri(t,tξ)=diri(t,tξ)+DRijrj(t,tξ)DRjiri(t,tξ)=diVξi(t)+DRijVξj(t)DRjiVξi(t). $ (2.6)

    Denote $ \mathbf{V}^\xi(t) = (V^\xi_1(t), V^\xi_2(t))^\top $, where $ \top $ represents the transpose of a vector. Then $ \mathbf{V}^\xi(t) $ satisfies

    $ ddtVξ(t)=BVξ(t), $ (2.7)

    where

    $ \mathbf{B} = \left( d1DR21DR12DR21d2DR12 \right)\\ . $

    Integrating (2.7) with respect to $ t $ from $ \xi $ to $ t $, we have

    $ Vξ(t)=exp(B(tξ))(Vξ1(ξ),Vξ2(ξ)),ξtξ+τ. $ (2.8)

    By using the definition of $ V_i^\xi(t) $ and (2.3),

    $ Vξ(t)=exp(B(tξ))(r1(ξ,0),r2(ξ,0)),ξtξ+τ=exp(B(tξ))(γ1I1(ξ),γ2I2(ξ)). $ (2.9)

    For $ t \geq\tau $ (hence $ t-\tau\geq 0 $), letting $ \mathbf{r}(t, \tau) = (r_1(t, \tau), r_2(t, \tau))^\top $, we obtain

    $ r(t,τ)=Vtτ(t)=exp(Bτ)(γ1I1(tτ),γ2I2(tτ)). $ (2.10)

    Denoting $ [b_{i j} (\tau)]_{2\times 2} : = \mathrm{exp}(\mathbf{B}\tau) $, it follows that

    $ ri(t,τ)=2j=1bij(τ)γjIj(tτ). $ (2.11)

    Substituting $ r_i(t, \tau) $ back into the $ I_i $ equations in (2.5) and taking out the first 4 equations for $ S_i $, and $ I_i $, $ i = 1, \; 2 $, results in the following new model:

    $ {dS1(t)dt=K1d1S1(t)+DS12S2(t)DS21S1(t)λ1S1(t)I1(t),dS2(t)dt=K2d2S2(t)+DS21S1(t)DS12S2(t)λ2S2(t)I2(t),t>τ,dI1(t)dt=d1I1(t)+DI12I2(t)DI21I1(t)+λ1S1(t)I1(t)γ1I1(t)+2j=1b1j(τ)γjIj(tτ),dI2(t)dt=d2I1(t)+DI21I1(t)DI12I2(t)+λ2S2(t)I2(t)γ2I2(t)+2j=1b2j(τ)γjIj(tτ). $ (2.12)

    For $ 0 < t \le \tau $, there is no relapsed individual reverting to the infectious class yet, and thus, the dynamics of $ S $ and $ I $ classes are governed by the following system of ordinary differential equations:

    $ {dS1(t)dt=K1d1S1(t)+DS12S2(t)DS21S1(t)λ1S1(t)I1(t),dS2(t)dt=K2d2S2(t)+DS21S1(t)DS12S2(t)λ2S2(t)I2(t),dI1(t)dt=d1I1(t)+DI12I2(t)DI21I1(t)+λ1S1(t)I1(t)γ1I1(t),dI2(t)dt=d2I2(t)+DI21I1(t)DI12I2(t)+λ2S2(t)I2(t)γ2I2(t).tτ. $ (2.13)

    The last term on the right side of the $ I_i $ equation in (2.12) accounts for non-local reverting force, reflecting how the individuals recovered $ \tau $ time units ago in all patches contribute to the growth of the infectious population in patch $ i $ through reverting to the infectious class. As is clear from the structure of the matrix $ \mathbf{B} $ and the expression (2.1), such a non-locality is caused by the mobility of the individuals during the recovered period.

    If we further assume $ d_1 = d_2 = d $, then

    $ \begin{array}{rcl} \mathbf{B} = \left( \begin{array}{cc} -d-D^R_{21} & D^R_{12}\\ D^R_{21} & -d-D^R_{12} \end{array} \right) = \left( d00d \right) +\left( DR21DR12DR21DR12 \right), \\ \end{array} $

    and we can obtain $ [b_{i j} (\tau)] = \mathrm{exp}(\mathbf{B}\tau) $ explicitly as

    $ b11(τ)=edτ(1δ1(τ)),b12(τ)=edτδ2(τ),b22(τ)=edτ(1δ2(τ)),b21(τ)=edτδ1(τ), $ (2.14)

    where

    $ δi(τ)=DRjiDRji+DRij(1e(DRji+DRij)τ),for1ij2. $ (2.15)

    Hence the model becomes

    $ {dS1(t)dt=K1d1S1(t)+DS12S2(t)DS21S1(t)λ1S1(t)I1(t),dS2(t)dt=K2d2S2(t)+DS21S1(t)DS12S2(t)λ2S2(t)I2(t),dI1(t)dt=d1I1(t)+DI12I2(t)DI21I1(t)+λ1S1(t)I1(t)γ1I1(t)+edτ(1δ1(τ))γ1I1(tτ)+edτδ2(τ)γ2I2(tτ),dI2(t)dt=d2I2(t)+DI21I1(t)DI12I2(t)+λ2S2(t)I2(t)γ2I2(t)+edτ(1δ2(τ))γ2I2(tτ)+edτδ1(τ)γ1I1(tτ). $ (2.16)

    From this model, it is seen that the dispersion of the individuals in the recovered period plays a different role from that of the susceptible and infectious individuals. The explanation for those instantaneous terms in (2.16) is quite straightforward, and we now explain those delayed terms in the model. The probability that an individual recovered in patch 1 can survive the relapse period is $ e^{-d\tau} $. Due to the mobility during the recovered period between the two patches, $ \tau $ time units later, a survived recovered individual may relapse in patch 1 with probability $ (1-\delta_1(\tau)) $ or in patch 2 with probability $ \delta_1(\tau) $. This explains the term $ e^{-d\tau}[1-\delta_1(\tau)]\gamma_1I_1(t-\tau) $ in the $ I_1 $ equation and the term $ e^{-d\tau}\delta_1(\tau)\gamma_1I_1(t-\tau) $ in the $ I_2 $ equation. The terms $ e^{-d\tau}(1-\delta_2(\tau))\gamma_2I_2(t-\tau) $ in $ I_2 $ equation and the term $ e^{-d\tau}\delta_2(\tau)\gamma_2I_2(t-\tau) $ in $ I_1 $ equation are explained similarly. Alternatively, we may explain these terms in light of fractions as below: among the individuals recovered in the first patch $ \tau $ time units ago, a fraction $ e^{-d\tau} $ can survive the relapse period, and a fraction $ (1-\delta_1(\tau)) $ of these survived individuals is now still in patch 1 while the remaining fraction $ \delta_1(\tau) $ of them has now moved to patch 2.

    For $ b_{ij}(\tau) $, one should expect the relation $ 0 \leq b_{ij}(\tau) \leq 1 $, and this relation will be used later in Section 5 to prove the persistence of the disease. Now we prove the above expectation by a comparison argument and properties of nonnegative matrices.

    Lemma 2.1. Let

    $ \underline{d} = \mathrm{min}\{d_1, d_2\}, \; \; \; and\; \; \; \bar{d} = \mathrm{max}\{d_1, d_2\}. $

    Then

    $ eˉdτ2i=1bij(τ)ed_τ,forj=1,2. $ (2.17)

    Proof. Choose a constant $ H > 0 $ sufficiently large such that

    $ H \gt \mathrm{max}\{d_1 \tau+D^R_{21}\tau, d_2\tau+D^R_{12}\tau\}. $

    Write $ \mathbf{B}\tau $ as $ \mathbf{B}\tau = -H\mathbf{E} + H\mathbf{E} + \mathbf{D}_0 + \mathbf{D}_1 $, where $ \mathbf{E} $ is the $ 2\times 2 $ identity matrix and

    $ \mathbf{D}_0 = \left( d1τ00d2τ \right) , \; \; \; \; \mathbf{D}_1 = \left( DR21τDR12τDR21τDR12τ \right). $

    Let $ \underline{\mathbf{D}} = -\underline{d} \tau \mathbf{E} $. Then both $ H \mathbf{E}+\mathbf{D}_0 +\mathbf{D}_1 $ and $ H\mathbf{E}+ \underline{\mathbf{D}} +\mathbf{D}_1 $ are nonnegative matrices and

    $ H\mathbf{E}+\mathbf{D}_0 +\mathbf{D}_1\leq H\mathbf{E}+ \underline{\mathbf{D}} +\mathbf{D}_1. $

    Thus,

    $ exp(Bτ)=exp(HE+HE+D0+D1)=exp(HE)exp(HE+D0+D1)exp(HE)exp(HE+D_+D1)=exp(D1)exp(D_). $ (2.18)

    Let $ V = (1, 1) $. It is easy to verify that $ V\mathbf{D}_1 = \mathbf{0} $, and hence $ V \mathrm{exp}(\mathbf{D}_1) = V\mathbf{E} $. Therefore,

    $ Vexp(Bτ)Vexp(D1)exp(D_)=Vexp(D_), $ (2.19)

    leading to the right side inequalities in (2.17). The left side inequalities in (2.17) can be similarly proved, and the proof of the lemma is completed.

    Our new model consists of two parts: a system of ODEs (2.13) for $ t\in [0, \tau] $ and a system of DDEs (2.12) for $ t\geq\tau $. For biological reasons, the following non-negative initial value conditions should be posed for the model:

    $ Si(0)0andIi(0)0,i=1,2. $ (3.1)

    In order for the model to be biologically well-posed, we need to make sure that the model (2.13)–(2.12) with (3.1) has a unique solution which remains non-negative and bounded. The following theorem confirms these properties.

    Theorem 3.1. The initial value problem (2.13)-(2.12)-(3.1) has a unique solution which exists globally (i.e., for all $ t \geq 0 $), remains non-negative and is bounded.

    Proof. The standard theory of ODEs ensures that the initial value problem (2.13)–(3.1) has a unique solution $ (S^0_1(t), S^0_2(t), I^0_1(t), I^0_2(t)) $, which exists globally, remains non-negative and is bounded. Consider the restriction of this solution on $ [0, \tau] $ and denote its components by

    $ \phi_i(\theta) = S^0_i(\theta), \; \; \; and\; \; \; \psi_i(\theta) = I^0_i(\theta), \; \; \; for\; \; \; i = 1, \; 2, \; \; \; and\; \; \; \theta\in[0, \tau]. $

    Then, $ \phi_i(\theta) $ and $ \psi_i(\theta) $ are continuous and non-negative functions on $ [0, \tau] $. By the fundamental theory of delay differential equations (see, e.g., [18]), we know that the DDE system (2.12) with the initial conditions

    $ S_i(\theta) = \phi_i(\theta)\; \; \; and\; \; \; I_i(\theta) = \psi_i(\theta), \; \; \; for\; \; \; i = 1, \; 2, $

    has a unique solution $ (S(t, \phi, \psi), I(t, \phi, \psi)) $, which is well-defined on its maximal interval of existence $ [\tau, t_{max}(\phi, \psi)) $, where

    $ (S(t, \phi, \psi), I(t, \phi, \psi)): = (S_1(t, \phi, \psi), S_2(t, \phi, \psi), I_1(t, \phi, \psi), I_2(t, \phi, \psi)), $
    $ (\phi, \psi): = (\phi_1(\theta), \phi_2(\theta), \psi_1(\theta), \psi_2(\theta)). $

    Firstly, we show the non-negativity of the solution for $ t \in[\tau, t_{max}(\phi, \psi)) $. For this purpose, let us rewrite the system (2.12) as follows:

    $ ddtS(t)=K+D(t)S(t), $ (3.2)
    $ ddtI(t)=C(t)I(t)+AI(tτ),tτ, $ (3.3)

    where $ \mathbf{S}(t) = (S_1(t), S_2(t))^\top $, $ \mathbf{I}(t) = (I_1(t), I_2(t))^\top $ and $ \mathbf{K} = (K_1, K_2)^\top $, and

    $ \begin{array}{rcl} \mathbf{D}(t) = \left( \begin{array}{cc} -d_1-D^S_{21}-\lambda_1I_1(t) & D^S_{12}\\ D^S_{21} & -d_2-D^S_{12}-\lambda_2I_2(t) \end{array} \right), \end{array} \begin{array}{rcl} \mathbf{A} = \left( \begin{array}{cc} b_{11}(\tau)\gamma_1 & b_{12}(\tau)\gamma_2 \\ b_{21}(\tau)\gamma_1 & b_{22}(\tau)\gamma_2 \end{array} \right), \end{array} $
    $ \begin{array}{rcl} \mathbf{C}(t) = \left( \begin{array}{cc} -d_1-D^I_{21}+\lambda_1S_1(t)-\gamma_1 & D^I_{12}\\ D^I_{21} & -d_2-D^I_{12}+\lambda_2S_2(t)-\gamma_2 \end{array} \right). \end{array} \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; $

    Noting that the off-diagonal elements of matrix $ \mathbf{D}(t) $ are non-negative, we conclude that the entries of the matrix $ e^{\int_{\tau}^t \mathbf{D}(\xi)d\xi} $ are all nonnegative. Indeed, let

    $ G(t) = \mathrm{max}\left\{d_1+D^S_{21}+\lambda_1I_1(t)+1, d_2+D^S_{12}+\lambda_2I_2(t)+1\right\} $

    and rewrite $ \mathbf{D}(t) $ as

    $ D(t)=(G(t)00G(t))+(G(t)d1DS21λ1I1(t)DS12DS21G(t)d2DS12λ2I2(t))G(t)E+ˉD(t). $

    Then all entries of $ \mathbf{\bar{D}}(t) $ are nonnegative, and hence, so are the entries of $ e^{\int_{\tau}^t \mathbf{\bar{D}}(\xi)d\xi} $. It is obvious that

    $ \begin{array}{rcl} e^{\int_{\tau}^t \left(-G(\xi)\mathbf{E}\right)d\xi} = \left( \begin{array}{cc} e^{\int_{\tau}^t \left(-G(\xi)\right)d\xi} & 0\\ 0 & e^{\int_{\tau}^t \left(-G(\xi)\right)d\xi} \end{array} \right). \end{array} $

    Noting that the scalar matrix $ -G(t)\mathbf{E} $ commutes with any $ 2\times 2 $ matrix (hence with $ \mathbf{\bar{D}}(t) $), we have

    $ e^{\int_{\tau}^t \mathbf{D}(\xi)d\xi} = e^{\int_{\tau}^t \left(-G(\xi)\mathbf{E}\right)d\xi}\cdot e^{\int_{\tau}^t \mathbf{\bar{D}}(\xi)d\xi}, $

    implying that all entries of $ e^{\int_{\tau}^t \mathbf{D}(\xi)d\xi} $ are nonnegative. Now from (3.2), we have

    $ S(t)=etτD(ξ)dξS(τ)+tτKetsτD(ξ)dξds0,fort[τ,tmax(ϕ,ψ)). $ (3.4)

    Similarly, for any $ t\geq\tau $, all entries of $ e^{\int_{\tau}^t \mathbf{C}(\xi)d\xi} $ are nonnegative. Moreover, it is obvious that all entries of $ \mathbf{A} $ are all non-negative. Now, (3.3) leads to

    $ I(t)=etτC(ξ)dξI(τ)+tτetsτC(ξ)dξAI(sτ)ds,tτ, $ (3.5)

    implying $ \mathbf{I}(t) \geq 0 $ for $ t\in [\tau, 2\tau] $ from the initial condition $ I_i (\theta) \geq 0 $ for $ \theta \in [0, \tau] $ and $ i = 1, \; 2 $. This and (3.5) ensure $ \mathbf{I}(t) \geq 0 $ for $ t\in [2\tau, 3\tau] $. By induction, we then conclude that $ \mathbf{I}(t) \geq 0 $ for $ t\in [\tau, t_{max}(\phi, \psi)) $.

    Now, we show that $ S_i (t) $, $ I_i (t) $ and $ R_i (t) $ are bounded for $ t\in [\tau, t_{max}(\phi, \psi)) $ and $ i = 1, \; 2 $. Noting that, by using the method of characteristic lines for the model (2.2), we can derive that $ r_i(t, a)\geq0 $, as well as $ R_i(t)\geq0 $, $ i = 1, \; 2 $. Let $ N(t) = S_1(t)+S_2(t)+ I_1(t)+I_2(t)+R_1(t)+R_2(t) $. Then from System (2.5), we have

    $ ddtN(t)=K1+K2d1S1(t)d2S2(t)d1I1(t)d2I2(t)d1R1(t)d2R2(t)K1+K2min{d1,d2}N(t). $

    This implies that $ N(t) $ is bounded, and so are $ S_i (t) $, $ I_i (t) $ and $ R_i (t) $ for $ i = 1, \; 2 $ and $ t\in [\tau, t_{max}(\phi, \psi)) $. By the theory of continuation of solutions (see, e.g., [18]), we conclude that $ t_{max}(\phi, \psi) = \infty $, which means the solution $ (S(t, \phi, \psi), I (t, \phi, \psi)) $ exists globally. This together with the results on $ S^0_i (t) $ and $ I^0_i (t) $ on $ t\in [0, \tau) $ implies that all of the above results actually hold for all $ t \geq 0 $. This completes the proof of Theorem 3.1.

    Remark 3.1. From the proof of this theorem, we see that the $ S $-components of the solution to (2.12) with (3.1) actually remain positive. If we further assume $ \psi_i (0) > 0 $ for $ i = 1, \; 2 $, then the $ I $-components of the solution also remain positive.

    Remark 3.2. Although the new model consists of two parts, (2.13) only plays a role of generating the necessary initial functions on $ [0, \tau] $ for (2.12). The long term behavior of the solution to (2.12)-(2.13)-(3.1) is indeed determined by (2.12). Therefore we only consider (2.12) since we are only interested in the long term disease dynamics. Note that (2.12) is an autonomous system of delay differential equations, and hence, its long dynamics is independent of the initial time. Because of this and for convenience of notations in applying existing theories and results on long term dynamics delay differential equations, we move the initial time $ \tau $ to $ 0 $ and accordingly consider initial functions on the interval $ [-\tau, 0] $, leading to the following equivalent system for (2.12):

    $ {dS1(t)dt=K1d1S1(t)+DS12S2(t)DS21S1(t)λ1S1(t)I1(t),dS2(t)dt=K2d2S2(t)+DS21S1(t)DS12S2(t)λ2S2(t)I2(t),t>0,dI1(t)dt=d1I1(t)+DI12I2(t)DI21I1(t)+λ1S1(t)I1(t)γ1I1(t)+2j=1b1j(τ)γjIj(tτ),dI2(t)dt=d2I1(t)+DI21I1(t)DI12I2(t)+λ2S2(t)I2(t)γ2I2(t)+2j=1b2j(τ)γjIj(tτ). $ (3.6)

    with the initial conditions specified in $ [-\tau, 0] $ by

    $ Si(θ)=ϕi(θ)C([τ,0],R+)Ii(θ)=ψi(θ)C([τ,0],R+),fori=1,2. $ (3.7)

    In the rest of the paper, we will just explore the dynamics of (3.6)–(3.7).

    In this section, we assume that the travel rate matrices $ [D^S_{i j}] $, $ [D^I_{i j}] $ and $ [D^R_{i j}] $ are irreducible. As usual, we start by investigating disease-free equilibrium. A disease-free equilibrium (DFE) is a steady state solution of the system (3.6) with all infectious variables being zeros. A DFE for the model (3.6) is thus given by $ E_0 = (S_1^{(0)}, S_2^{(0)}, 0, 0) $ with $ \mathbf{S}^{(0)} = (S_1^{(0)}, S_2^{(0)})^\top $ satisfying the linear system $ \mathbf{M}\mathbf{S}^{(0)} = \mathbf{K} $, where

    $ \begin{array}{rcl} \mathbf{M} = \left( \begin{array}{cc} d_{1}+D^S_{21} & -D^S_{12}\\ -D^S_{21} & d_{2}+D^S_{12} \end{array} \right). \end{array} $

    Note that matrix $ \mathbf{M} $ is irreducible, has positive column sums and negative off-diagonal entries. Thus $ \mathbf{M} $ is a non-singular M-matrix (see [17], page 141) with $ \mathbf{M}^{-1} > 0 $, and therefore, the linear system has a unique solution given by $ \mathbf{S}^{(0)} = \mathbf{M}^{-1} \mathbf{K} > 0 $. Indeed, one can explicitly solve $ \mathbf{M}\mathbf{S}^{(0)} = \mathbf{K} $ to obtain

    $ S(0)1=DS12K1+DS12K2+d2K1d1d2+d1DS12+d2DS21,andS(0)2=DS21K1+DS21K2+d1K2d1d2+d1DS12+d2DS21. $ (4.1)

    Now we discuss the stability of $ E_0 $. To this end, we consider the linearization of (3.6) at $ E_0 $:

    $ {dS1(t)dt=d1S1(t)+DS12S2(t)DS21S1(t)λ1S(0)1I1(t),dS2(t)dt=d2S2(t)+DS21S1(t)DS12S2(t)λ2S(0)2(t)I2(t),dI1(t)dt=d1I1(t)+DI12I2(t)DI21I1(t)+λ1S(0)1I1(t)γ1I1(t)+b11(τ)γ1I1(tτ)+b12(τ)γ2I2(tτ),dI2(t)dt=d2I2(t)+DI21I1(t)DI12I2(t)+λ2S(0)2I2(t)γ2I2(t)+b21(τ)γ1I1(tτ)+b22(τ)γ2I2(tτ). $ (4.2)

    The characteristic equation of (4.2) is given by

    $ Δ1(z)Δ2(z,τ)=0, $ (4.3)

    where

    $ Δ1(z)=det[z+d1+DS21DS12DS21z+d2+DS12]=z2+(d1+d2+DS12+DS21)z+d1d2+d1DS12+d2DS21, $

    and

    $ Δ2(z,τ)=det[z+d1+DI21+γ1λ1S(0)1b11(τ)γ1ezτDI12b12(τ)γ2ezτDI21b21(τ)γ1ezτz+d2+DI12+γ2λ2S(0)2b22(τ)γ2ezτ]. $

    It is obvious that all roots of $ \Delta_1(z) $ have negative real parts. Thus, the stability of $ E_0 $ is fully determined by the roots of $ \Delta_2(z, \tau) $. Note that the $ I $ equations of the linearization (4.2) is decoupled from the $ S $ equations and $ \Delta_2(z, \tau) = 0 $ is nothing but precisely the characteristic equation of the decoupled $ I $ equations in (4.2). Write the $ I $-equations as

    $ I(t)=FI(t)+AI(tτ), $ (4.4)

    where $ \mathbf{A} $ is defined in Theorem 3.1, and

    $ \begin{array}{rcl} \mathbf{F} = \left( \begin{array}{cc} -d_1-D^I_{21}+\lambda_1S^0_1-\gamma_1 & D^I_{12}\\ D^I_{21} & -d_2-D^I_{12}+\lambda_2S^0_2-\gamma_2 \end{array} \right). \end{array} $

    Note that $ \mathbf{A} $ and $ \mathbf{F} $ are quasi-positive and irreducible matrices. Thus, a cooperative and irreducible system of ordinary differential equations can be associated with the system (4.4) by simply replacing $ I(t-\tau) $ with $ I(t) $ in (4.4). This leads to the system

    $ I(t)=(F+A)I(t)WI(t), $ (4.5)

    with

    $ \begin{array}{rcl} \mathbf{W} & = &\left( \begin{array}{cc} -d_1-D^I_{21}+\lambda_1S^{(0)}_1-\gamma_1 +b_{11}(\tau)\gamma_1 & D^I_{12}+b_{12}(\tau)\gamma_2\\ D^I_{21}+b_{21}(\tau)\gamma_1 & -d_2-D^I_{12}+\lambda_2S^{(0)}_2-\gamma_2+b_{22}(\tau)\gamma_2 \end{array} \right)\\ &\triangleq& \left( w11w12w21w22 \right). \end{array} $

    By using the stability criteria for the cooperative and irreducible systems (see Theorem 5.1 and Corollary 5.2 in [19]), we know that the linear stability of the trivial equilibrium for system (4.4) is equivalent to that for system (4.5). Therefore, we just need explore the roots for characteristic equation of (4.5). Noting that the off-diagonal elements of matrix $ \mathbf{W} $ are non-negative, we conclude that the characteristic equation of (4.5) has two real zeros $ z_2 < z_1 $:

    $ z_1 = \frac{w_{11}+w_{22}+\sqrt{(w_{11}-w_{22})^2+4w_{12}w_{21}}}{2}, \; \; \; z_2 = \frac{w_{11}+w_{22}-\sqrt{(w_{11}-w_{22})^2+4w_{12}w_{21}}}{2}. $

    Hence, if

    $ w11+w22+(w11w22)2+4w12w21<0, $ (4.6)

    the trivial solution of the system (4.4) is stable, and so is $ E_0 $ for (3.6); and when

    $ w11+w22+(w11w22)2+4w12w21>0, $ (4.7)

    the trivial solution of the system (4.4) and $ E_0 $ are both unstable and so is $ E_0 $ for (3.6).

    By estimating the trace and determinant of $ \mathbf{W} $, we find that the following three more explicit conditions (4.8) (4.9) and (4.10), directly in terms of the model parameters, imply that (4.6) hold.

    $ λ1S(0)1d1+DI21+γ1(1b11(τ))<1, $ (4.8)
    $ λ2S(0)2d2+DI12+γ2(1b22(τ))<1, $ (4.9)
    $ (d1DI21+λ1S(0)1γ1+b11(τ)γ1)(d2DI12+λ2S(0)2γ2+b22(τ)γ2)(DI12+b12(τ)γ2)(DI21+b21(τ)γ1)>1. $ (4.10)

    Based on the preceding discussion, we then have proved the following theorem.

    Theorem 4.1. If (4.6) holds, then the disease-free equilibrium $ E_0 = (S^{(0)}_1, S^{(0)}_2, 0, 0) $ of the system (3.6) is locally asymptotically stable.

    The next theorem shows that $ E_0 = (S^{(0)}_1, S^{(0)}_2, 0, 0, ) $ is actually globally asymptotically stable.

    Theorem 4.2. If (4.6) holds, then the disease-free equilibrium $ E_0 = (S^{(0)}_1, S^{(0)}_2, 0, 0) $ of the system (3.6) is globally asymptotically stable.

    Proof. By Theorem 4.1, we know that $ E_0 $ is locally asymptotically stable if (4.8)–(4.10) satisfied. It merely remains to prove that $ E_0 $ is globally attractive in the case (4.8)–(4.10) held, that is, for any non-negative solutions $ (S_1(t), S_2(t), I_1(t), I_2(t)) $ of (3.6), we will prove that $ \lim_{t\rightarrow+\infty}(S_1(t), S_1(t), I_1(t), I_2(t)) = (S^{(0)}_1, S^{(0)}_2, 0, 0) $. From the $ S $-equations in System (3.6) and the non-negativity of the solutions to the system (3.6) with (3.7), we have

    $ {dS1(t)dt=K1d1S1(t)+DS12S2(t)DS21S1(t)λ1S1(t)I1(t)K1d1S1(t)+DS12S2(t)DS21S1(t),dS2(t)dt=K2d2S2(t)+DS21S1(t)DS12S2(t)λ2S2(t)I2(t)K2d2S2(t)+DS21S1(t)DS12S2(t). $ (4.11)

    This suggests the following comparison system for the $ S $-equations of (3.6)

    $ {du1(t)dt=K1d1u1(t)+DS12u2(t)DS21u1(t),du2(t)dt=K2d2u2(t)+DS21u1(t)DS12u2(t). $ (4.12)

    We have seen that the system (4.12) admits a unique positive equilibrium $ (S^{(0)}_1, S^{(0)}_2) $. It is easy to see that the stability of $ (S^{(0)}_1, S^{(0)}_2) $ for (4.12) is precisely determined by $ \Delta_1(z) = 0 $ where $ \Delta_1(z) $ is as in (4.3). Since all roots have negative real parts, $ (S^{(0)}_1, S^{(0)}_2) $ is globally asymptotically stable (in a linear system, local stability is equivalent to global stability). By the comparison theorem (see, e.g., [19,20]), we then have

    $ Silim supt+Si(t)limt+ui(t)=S(0)i,i=1,2. $ (4.13)

    Thus, for any constant $ \epsilon > 0 $, there is a large enough $ T $ such that $ S_i(t)\leq S^{(0)}_i+\epsilon $, for all $ t\geq T $.

    Now, for $ t\geq T $, we construct the following comparison linear system for the $ I $ equations in (3.6):

    $ {dI1(t)dt=d1I1(t)+DI12I2(t)DI21I1(t)+λ1(S(0)1+ϵ)I1(t)γ1I1(t)+b11(τ))γ1I1(tτ)+b12(τ)γ2I2(tτ),dI2(t)dt=d2I2(t)+DI21I1(t)DI12I2(t)+λ2(S(0)2+ϵ)I2(t)γ2I2(t)+b21(τ))γ1I1(tτ)+b22(τ)γ2I2(tτ). $ (4.14)

    By the same argument as that for the stability of (4.4), we know that the trivial solution of this system is globally asymptotically stable, implying that all solutions of the linear system (4.14) tend to the trivial solution as $ t\rightarrow\infty $. Note that (4.14) is a cooperative system of delay differential equations. By the comparison theorem, we then conclude that all $ I $ components of the solution to (3.6) with (3.7) also tend to zeros as $ t\rightarrow\infty $. This in return implies that the $ S $ equation in (3.6) has (4.12) as its limiting system, which has the dynamics of global convergence to $ (S^{(0)}_1, S^{(0)}_2) $. Finally by the theory of asymptotically autonomous systems (see, e.g., [21,22]), we conclude that the $ S $ component of the solution to (3.6) with (3.7) also converges to $ (S^{(0)}_1, S^{(0)}_2) $. This confirms the global attractivity of $ E_0 $ for (3.6) under the condition (4.8)–(4.10) held, and hence completes the proof.

    In Section 4, under the assumption that the travel rate matrices $ [D^S_{i j}] $, $ [D^I_{i j}] $ and $ [D^R_{i j}] $ are irreducible, we have shown DFE $ E_0 $ is globally asymptotically stable if (4.6) is satisfied. One naturally wonders what happens when (4.7) holds instead. In this section, we still assume the irreducibility of all travel rate matrices, and we will prove that the disease is persistent in all patches when (4.7) is satisfied. This conclusion together with a well-known result for persistent systems actually implies the existence of an endemic equilibrium for (3.6).

    For the convenience of stating and proving the main results, we first introduce some notations. Let $ C: = C([-\tau, 0], \mathbb{R}^2) $ denote the set of all continuous functions from $ [-\tau, 0] $ to $ \mathbb{R}^2 $. As is customary, $ C_{+}: = C([-\tau, 0], \mathbb{R}_{+}^2) $ denotes the subset of $ C $ consisting of all non-negative functions. By Theorem 3.1 and Remark 3.1, for any $ (\phi, \psi)\in C_{+} \times C_{+} $, with $ \psi(0) > 0 $ there is a unique solution to (3.6), denoted by

    $ (S(t, \phi, \psi), I(t, \phi, \psi)): = (S_1(t, \phi, \psi), S_2(t, \phi, \psi), I_1(t, \phi, \psi), I_2(t, \phi, \psi)), $

    whose components are all positive and bounded for $ t > 0 $.

    Theorem 5.1. Assume that all three travel rate matrices $ [D^S_{i j}] $, $ [D^I_{i j}] $ and $ [D^R_{i j}] $ are irreducible. Suppose (4.7) hold, then there is an $ \varepsilon > 0 $ such that for every $ (\phi, \psi)\in C_{+} \times C_{+} $ with $ \psi(0) = (\psi_1(0), \psi_2(0)) > 0 $, meaning that $ \psi_i(0) \ge 0 $ for $ i = 1, 2 $ and $ \psi_1(0)+\psi_2(0) > 0 $, then the solution $ (S(t), I (t)) = (S(t, \phi, \psi), I (t, \phi, \psi)) $ of (3.6) satisfies

    $ \liminf\limits_{t\rightarrow\infty} I_i(t, \phi, \psi)\geq \varepsilon, \; \; \; \; i = 1, \; 2. $

    Moreover, the model (3.6) admits at least one (componentwise) positive equilibrium.

    Proof. Define $ X: = \{(\phi, \psi)\in C_{+} \times C_{+}\} $, $ X_0: = \{(\phi, \psi)\in X, \psi_i(0) > 0, \; i = 1, \; 2\} $ and $ \partial {X_0}: = X \setminus X_0 $. It then suffices to show that (3.6) is uniformly persistent with respect to $ (X_0, \partial {X_0}) $.

    Let $ \Phi(t) : X \rightarrow X $ be the solution semiflow of (3.6)-(3.7), that is, $ \Phi(t)(\phi, \psi) = (S_t(\phi, \psi), I_t(\phi, \psi)) $ holds for $ t \geq 0 $, with $ S(t) = \varphi(t) $, $ I(t) = \psi(t) $, $ t \in [-\tau, 0] $. By the fundamental theory for functional differential equations with bounded delays established in [18], the solution semin-flow $ \Phi(t) $ is actually compact for $ t \ge \tau $ (consequence of the Arzela-Ascoli Theorem.) By Theorem 3.1 and Remark 3.1, $ X $ and $ X_0 $ are positively invariant for $ \Phi(t) $. Clearly, $ \partial {X_0} = \{(\phi, \psi)\in X, \; \psi_i(0) = 0, \; for \; at \; least\; one \; i \in \{1, 2, \}\} $ and it is relatively closed in $ X $. Furthermore, system (3.6) is point dissipative in $ \mathbb{R}^2_+ $ since nonnegative solutions of (3.6) are ultimately bounded (see Theorem 3.1).

    Define $ \Omega_{\partial} = \{(\phi, \psi) \in X : (S_t(\phi, \psi), I_t (\phi, \psi)) \in \partial {X_0}, \; \forall t \geq 0\} $. We now show that

    $ Ω={(ϕ,ψ)X0:I(t,ϕ,ψ)=0,t0}. $ (5.1)

    Assume $ (\phi, \psi) \in \Omega_{\partial} $. It suffices to show that $ I(t, \phi, \psi) = 0 $, $ \forall t \geq 0 $. For the sake of contradiction, assume that there is a $ t_0 \geq 0 $ such that $ I_1 (t_0) > 0 $. Then by the definition of $ \Omega_{\partial} $, we can derive that $ I_2 (t_0, \phi, \psi) = 0 $. This leads to

    $ ddtI2(t,ϕ,ψ)t=t0=(d2+DI12+γ2)I2(t0,ϕ,ψ)+DI21I1(t0,ϕ,ψ)+λ2S2(t0,ϕ,ψ)I2(t0,ϕ,ψ)+b21(τ)γ1I1(t0τ,ϕ,ψ)+b22(τ)γ2I2(t0τ,ϕ,ψ)DI21I1(t0,ϕ,ψ)+b21(τ)γ1I1(t0τ,ϕ,ψ)+b22(τ)γ2I2(t0τ,ϕ,ψ)>0. $ (5.2)

    It follows that there is an $ \epsilon_0 > 0 $ such that $ I_2(t, \phi, \psi) > 0 $ and $ t_0 < t < t_0+\epsilon_0 $. Clearly, we can restrict $ \epsilon_0 > 0 $ small enough such that $ I_1 (t, \phi, \psi) > 0 $ for $ t_0 < t < t_0+\epsilon_0 $. This means that $ (S_t(\phi, \psi), I_t(\phi, \psi))\not\in \partial X_0 $ for $ t_0 < t < t_0+\epsilon_0 $, which contradicts the assumption that $ (\phi, \psi) \in \Omega_{\partial} $. This proves (5.1).

    Choose $ \xi > 0 $ small enough such that

    $ wξ11+wξ22+(wξ11wξ22)2+4w12w21>0, $ (5.3)

    where $ w_{12} $ and $ w_{21} $ are as in Section 4 and

    $ w^\xi_{11} = -d_1-D^I_{21}+\lambda_1(S^{(0)}_1-\xi)-\gamma_1 +b_{11}(\tau)\gamma_1, \; w^\xi_{22} = -d_2-D^I_{12}+\lambda_2(S^{(0)}_2-\xi)-\gamma_2+b_{22}(\tau)\gamma_2. $

    Let us consider the following linear system

    $ {dS1(t)dt=K1d1S1(t)+DS12S2(t)DS21S1(t)ελ1S1(t)=K1(d1+DS21+ελ1)S1(t)+DS12S2(t),dS2(t)dt=K2d2S2(t)+DS21S1(t)DS12S2(t)ελ2S2(t)=K2(d2+DS12+ελ2)S2(t)+DS21S1(t), $ (5.4)

    which is a perturbation of (4.12). Restrict $ \varepsilon > 0 $ small enough such that (5.4), just as (4.12), has a unique positive equilibrium $ (S^{(0)}_1(\varepsilon), S^{(0)}_2(\varepsilon)) $ which is globally asymptotically stable. By the implicit function theorem, it follows that $ (S^{(0)}_1(\varepsilon), S^{(0)}_2(\varepsilon)) $ is continuous in $ \varepsilon $. Thus, we can further restrict $ \varepsilon $ small enough such that $ (S^{(0)}_1(\varepsilon), S^{(0)}_2(\varepsilon)) > (S^{(0)}_1-\xi, S^{(0)}_2-\xi) $.

    Next for the solution $ (S(t, \phi, \psi), I (t, \phi, \psi)) $ of (3.6) through $ (\phi, \psi) $, we claim that

    $ lim suptmax{I1(t,ϕ,ψ),I2(t,ϕ,ψ)}>ε,forall(ϕ,ψ)X0. $ (5.5)

    Otherwise, there is a $ T_1 > 0 $ such that $ 0 < I_i(t, \phi, \psi) \leq \varepsilon $, $ i = 1, \; 2 $, for all $ t \geq T_1 $. Then for $ t \geq T_1 $, we have

    $ {dS1(t)dtK1d1S1(t)+DS12S2(t)DS21S1(t)λ1S1(t)ε=K1(d1+DS21+λ1ε)S1(t)+DS12S2(t),dS2(t)dtK2d2S2(t)+DS21S1(t)DS12S2(t)λ2S2(t)ε=K2(d2+DS12+λ2ε)S2(t)+DS21S1(t). $ (5.6)

    Since the equilibrium $ (S^{(0)}_1(\varepsilon), S^{(0)}_2(\varepsilon)) $ of (5.4) is globally asymptotically stable and $ (S^{(0)}_1(\varepsilon), S^{(0)}_2(\varepsilon)) > (S^{(0)}_1-\xi, S^{(0)}_2-\xi) $, there is a $ T_2 $ such that $ (S_1(t), S_2(t)) > (S^{(0)}_1-\xi, S^{(0)}_2-\xi) $ for $ t \geq T_1 + T_2 $. Consequently, for $ t \geq T_1 + T_2 $,

    $ {dI1(t)dtd1I1(t)+DI12I2(t)DI21I1(t)+λ1(S(0)1ξ)I1(t)γ1I1(t)+b11(τ)γ1I1(tτ)+b12(τ)γ2I2(tτ),dI2(t)dtd2I2(t)+DI21I1(t)DI12I2(t)+λ2(S(0)2ξ)I2(t)γ2I2(t)+b21(τ)γ1I1(tτ)+b22(τ)γ2I2(tτ). $ (5.7)

    By the same arguments as that for the stability and instability of the ODE (4.5) and the DDE (4.4) in Section. 4, we know that the condition (5.3) implies that the trivial solution of the linear system

    $ {dI1(t)dt=d1I1(t)+DI12I2(t)DI21I1(t)+λ1(S(0)1ξ)I1(t)γ1I1(t)+b11(τ)γ1I1(tτ)+b12(τ)γ2I2(tτ),dI2(t)dt=d2I2(t)+DI21I1(t)DI12I2(t)+λ2(S(0)2ξ)I2(t)γ2I2(t)+b21(τ)γ1I1(tτ)+b22(τ)γ2I2(tτ), $ (5.8)

    is unstable. This together with (5.7) and the comparison theorem implies that there is at least one $ i \in \{1, 2\} $ such that $ I_i (t) \rightarrow \infty $ as $ t\rightarrow\infty $, a contradiction to the boundedness of solutions. Therefore (5.5) holds.

    Note that $ (S^{(0)}_1, S^{(0)}_2) $ is globally asymptotically stable in $ \mathbb{R}^2_+\setminus \{0\} $ for system (4.12). By the aforementioned claim, it then follows that $ E_0 $ is an isolated invariant set in $ X $, and $ W^s (E_0) \cap X_0 = \emptyset $. Clearly, every orbit in $ \Omega_{\partial} $ converges to $ E_0 $, and $ E_0 $ is the only invariant set in $ \Omega_{\partial} $. By Theorem 4.6 in [25], we conclude that system (3.6) is indeed uniformly persistent with respect to $ (X_0, \partial X_0) $. Moreover, by Theorem 2.4 in [27], system (3.6) has an equilibrium $ (S^{\ast}_1, S^{\ast}_2, I^{\ast}_1, I^{\ast}_2) \in X_0 $, implying that $ (S^{\ast}_1, S^{\ast}_2) \in \mathbb{R}^2_+ $ and $ (I^{\ast}_1, I^{\ast}_2)\in \mathrm{int}(\mathbb{R}^2_+) $. We further claim that $ (S^{\ast}_1, S^{\ast}_2) \in \mathbb{R}^2_+\setminus \{0\} $. Suppose that $ (S^{\ast}_1, S^{\ast}_2) = 0 $. By the $ I $ -equations in (3.6), we then obtain

    $ 0=d1I1+DI12I2DI21I1γ1I1+b11(τ)γ1I1+b12(τ)γ2I2,0=d2I2+DI21I1DI12I2γ2I2+b21(τ)γ1I1+b22(τ)γ2I2. $

    and hence

    $ 0 = \left[-d_1+\left(b_{11}(\tau)+b_{21}(\tau)-1\right)\gamma_1\right]I_1^{\ast}+\left[-d_2+\left(b_{22}(\tau)+b_{12}(\tau)-1\right)\gamma_2\right]I_2^{\ast}. $

    By Lemma 2.1, we know that $ \sum_{i = 1}^{2}b_{ij}(\tau) \leq e^{-\underline{d}\tau} \leq1, \; for\; j = 1, \; 2 $, therefore, $ I_1^{\ast} = I_2^{\ast} = 0 $, a contradiction. By the $ S $-equation in (3.6) and the irreducibility of the cooperative matrix $ [D^S_{i j}] $, it follows that $ S^{\ast} = S(t, S^{\ast}, I^{\ast}) \in \mathrm{int}(\mathbb{R}^2_+) $ with $ S^{\ast}: = (S^{\ast}_1, S^{\ast}_2) $ and $ I^{\ast}: = (I^{\ast}_1, I^{\ast}_2) $, for $ \forall t > 0 $. Then $ (S^{\ast}, I^{\ast}) $ is a componentwise positive equilibrium of system (3.6), meaning that the system (3.6) has an epidemic equilibrium.

    The results in Sections 4 and 5 are obtained under the assumption that the travel rate matrices are all irreducible. In reality, these assumptions may not be satisfied. For example, when an infectious disease is reported in one or more cities, the health authorities in some or all cities may implement a ban against travel by the infected individuals. Such a measure may make some travel rate matrices reducible. In this section, we deal with cases allowing reducible rate matrices.

    For convenience of comparison later, we first consider the case when the two patches are fully disconnected by setting all dispersal rates to zero, implying that

    $ b12(τ)=b21(τ)=0,b11(τ)=ed1τ:=ϵ1,b22(τ)=ed2τ:=ϵ2. $ (6.1)

    Thus, (3.6) is decoupled to

    $ {dS1(t)dt=K1d1S1(t)λ1S1(t)I1(t),dI1(t)dt=d1I1(t)+λ1S1(t)I1(t)γ1I1(t)+ϵ1γ1I1(tτ), $

    for patch 1, and

    $ {dS2(t)dt=K2d2S2(t)λ2S2(t)I2(t),dI2(t)dt=d2I2(t)+λ2S2(t)I2(t)γ2I2(t)+ϵ2γ2I2(tτ), $

    for patch 2. By the results in [26], the disease dynamics in each patch in such a disconnected case is described by the corresponding basic reproduction number

    $ R_{i0}^{(0)}\triangleq \frac{K_i}{d_i}\frac{\lambda_i}{d_i+\gamma_i(1-\epsilon_i)}, \; \; \; i = 1, \; 2, $

    as summarized below.

    Theorem 6.1. If $ R_{i0}^{(0)} < 1 $, then the disease dies out in Patch $ i $ ($ i $ = 1, 2) in the sense that the disease-free equilibrium $ (\frac{K_i}{d_i}, 0) $ is globally asymptotically stable; if $ R_{i0}^{(0)} > 1 $, then the disease will persist in the population in the sense that the disease-free equilibrium is unstable and there is a unique endemic equilibrium

    $ (S_i^{\ast}, I_i^{\ast}) = \left(\frac{d_i+\gamma_i(1-\epsilon_i)}{\lambda_i}, \frac{K_i\lambda_i-d_i[d_i+\gamma_i(1-\epsilon_i)]}{\lambda_i[d_i+\gamma_i(1-\epsilon_i)]}\right), $

    which is asymptotically stable.

    In the rest of this section, we explore the impact of dispersals between the two patches on the disease dynamics of (3.6) in cases allowing reducible travel rate matrices. We only demonstrate three simpler scenarios that make the two patches connected: (i) only susceptible individuals disperse; (ii) the dispersals of recovered individuals are unidirectional; (iii) the dispersals of infected individuals are unidirectional.

    In this subsection, we assume that only susceptible individuals in the two patches travel. Such an assumption may account for the situation when all infectious and recovered individuals are prohibited (e.g., by health authorities) from traveling. This implies that $ D^S_{12} $ and $ D^S_{21} $are positive, but $ D^I_{12} = D^I_{21} = D^R_{12} = D^R_{21} = 0 $. Accordingly, one can compute to obtain the following:

    $ \mathbf{B} = \left( d100d2 \right) , \; \; \; and\; \; \; [b_{ij}(\tau)] = e^{\left(\mathbf{B}\tau\right)} = \left( ϵ100ϵ2 \right), $

    where $ \epsilon_i $, $ i = 1, \; 2 $ are defined in (6.1). In such a case, (3.6) reduces to

    $ {dS1(t)dt=K1d1S1(t)+DS12S2(t)DS21S1(t)λ1S1(t)I1(t),dS2(t)dt=K2d2S2(t)+DS21S1(t)DS12S2(t)λ2S2(t)I2(t),dI1(t)dt=d1I1(t)+λ1S1(t)I1(t)γ1I1(t)+ϵ1γ1I1(tτ),dI2(t)dt=d2I1(t)+λ2S2(t)I2(t)γ2I2(t)+ϵ2γ2I2(tτ). $ (6.2)

    We have seen that the DFE $ E_0 $ still exists and is given by (4.1), but its stability/instability can not be concluded from Theorems 4.1 and 4.2 as the irreducibility of $ [D^I_{i j}] $ and $ [D^R_{i j}] $ does not hold. Linearizing (6.2) at $ E_0 $ leads to

    $ {dS1(t)dt=d1S1(t)+DS12S2(t)DS21S1(t)λ1S(0)1I1(t),dS2(t)dt=d2S2(t)+DS21S1(t)DS12S2(t)λ2S(0)2(t)I2(t),dI1(t)dt=d1I1(t)+λ1S01I1(t)γ1I1(t)+ϵ1γ1I1(tτ),dI2(t)dt=d2I2(t)+λ2S02I2(t)γ2I2(t)+ϵ2γ2I2(tτ). $ (6.3)

    The characteristic equation of (6.3) is given by

    $ Δ1(z)Δ3(z,τ)Δ4(z,τ)=0, $

    where $ \Delta_1(z) $ is given by (4.3) and

    $ \Delta_3(z, \tau) = z+d_1+\gamma_1-\lambda_1S^{(0)}_1-\epsilon_1\gamma_1e^{-z\tau}, $
    $ \Delta_4(z, \tau) = z+d_2+\gamma_2-\lambda_2S^{(0)}_2-\epsilon_2\gamma_2e^{-z\tau}. $

    It is obvious that all roots of $ \Delta_1(z) $ have negative real parts. By the results on Hayes equation (see the Appendix in [18]), one knows that for $ i = 3, \; 4 $, all roots of $ \Delta_i(z, \tau) = 0 $ have negative real parts if and only if

    $ Ri0=λiS(0)idi+γi(1ϵi)<1. $

    Therefore, the DFE $ E_0 $ is asymptotically stable if $ \max\{R_{10}, R_{20}\} < 1 $ and it is unstable if $ \max\{R_{10}, R_{20}\} > \; 1 $.

    In the unstable case, we expect other equilibrium to appear. We start with looking for possible boundary equilibria, that is, equilibrium of the form $ E_1 = (S^{(1)}_1, S^{(1)}_2, I^{(1)}_1, 0) $ or $ E_2 = (S^{(2)}_1, S^{(2)}_2, 0, I^{(2)}_2) $ with $ I^{(1)}_1 > 0 $ for the former or $ I^{(2)}_2 > 0 $ for the latter. For $ E_1 $, we need to solve the algebraical equations

    $ {K1d1S(1)1+DS12S(1)2DS21S(1)1λ1S(1)1I(1)1=0,K2d2S(1)2+DS21S(1)1DS12S(1)2=0,d1I(1)1+λ1S(1)1I(1)1γ1I(1)1+ϵ1γ1I(1)1=0, $

    for positive $ S^{(1)}_1 $, $ S^{(1)}_2 $ and $ I^{(1)}_1 $ which are determined by

    $ S(1)1=d1+γ1(1ϵ1)λ1,S(1)2=1d2+DS12(K2+DS21d1+γ1(1ϵ1)λ1),I(1)1=1d1+γ1(1ϵ1)[K1(d1+DS21)(d1+γ1(1ϵ1)λ1+DS12d2+DS12(K2+DS21d1+γ1(1ϵ1)λ1)]=d1d2+d1DS12+d2DS21λ1(d2+DS12)(R101). $ (6.4)

    Thus, $ E_1 $ exists $ (I^{(1)}_1 > 0) $ if and only if

    $ R10=λ1S(0)1d1+γ1(1ϵ1)>1. $

    Similarly, for $ E_2 $ we have

    $ S(2)1=1d1+DS21(K1+DS12d2+γ2(1ϵ2)λ2),S(2)2=d2+γ2(1ϵ2)λ2,I(2)2=1d2+γ2(1ϵ2)[K2(d2+DS12)(d2+γ2(1ϵ2))λ2+DS21d1+DS21(K1+DS12d2+γ2(1ϵ2)λ2)]=d1d2+d1DS12+d2DS21λ2(d1+DS21)(R201). $

    Hence, $ E_2 $ exists $ (I^{(2)}_2 > 0) $ if and only if

    $ R20=λ2S(0)2d2+γ2(1ϵ2)>1. $

    Finally, an interior equilibrium is an equilibrium of the form $ E^* = (S^{\ast}_1, S^{\ast}_2, I^{\ast}_1, I^{\ast}_2) $ with all components positive, which can be determined from the following equations,

    $ {K1d1S1+DS12S2DS21S1λ1S1I1=0,K2d2S2+DS21S1DS12S2λ2S2I2=0,d1I1+λ1S1I1γ1I1+ϵ1γ1I1=0,d2I2+λ2S2I2γ2I2+ϵ2γ2I2=0. $

    Solving these equations for positive components leads to

    $ S1=d1+γ1(1ϵ1)λ1,S2=d2+γ2(1ϵ2)λ2,I1=1d1+γ1(1ϵ1)[K1(d1+DS21)(d1+γ1(1ϵ1))λ1+DS12d2+γ2(1ϵ2)λ2],I2=1d2+γ2(1ϵ2)[K2(d2+DS12)(d2+γ2(1ϵ2))λ2+DS21d1+γ1(1ϵ1)λ1]. $

    Define

    $ ˆR10=λ1S(2)1d1+γ1(1ϵ1),ˆR20=λ2S(1)2d2+γ2(1ϵ2). $

    By straightforward calculations we can further express $ I^{\ast}_1 $ and $ I^{\ast}_2 $ in terms of $ \hat{R}_{10} $ and $ \hat{R}_{20} $ as the following:

    $ I1=d1+DS21λ1(ˆR101), $
    $ I2=d2+DS12λ2(ˆR201). $

    Thus, the interior equilibrium $ E^* $ exists if and only if

    $ ˆR10>1andˆR20>1. $

    Remark 6.1. Direct computations show that

    $ \hat{R}_{10} \lt {R}_{10}\; \; \Leftrightarrow\; \; {R}_{20} \gt 1, \; \; \; \; \hat{R}_{20} \lt {R}_{20}\; \; \Leftrightarrow\; \; {R}_{10} \gt 1. $

    Moreover, $ \hat{R}_{10} < 1 < {R}_{10} $ and $ \hat{R}_{20} < 1 < {R}_{20} $ can not hold simultaneously.

    Summarizing the above analyses, we have obtained the following theorem for the model system (6.2), in terms of $ R_{i0} $ and $ \hat{R}_{i0} $ for $ i = 1, \; 2 $.

    Theorem 6.2. Consider the system (6.2)

    $\mathbf{(i)}$ If $ \max \{R_{10}, R_{20}\} < 1 $, then the disease-free equilibrium $ E_0 $ is locally asymptotically stable; if $ \max \{R_{10}, R_{20}\} > 1 $, then the disease-free equilibrium $ E_0 $ becomes unstable.

    $\mathbf{(ii)}$ If $ R_{10} > 1 $ and $ {R}_{20} < 1 $, then the boundary equilibrium $ E_1 $ exists and is asymptotically stable.

    $\mathbf{(iii)}$ If $ R_{20} > 1 $ and $ {R}_{10} < 1 $, then the boundary equilibrium $ E_2 $ exists and is asymptotically stable.

    $\mathbf{(iv)}$ Assume $ {R}_{10} > 1 $ and $ {R}_{20} > 1 $.

    (a) If $ {R}_{10} < {R}_{20} $ and $ \hat{R}_{10} < 1 $, then the boundary equilibrium $ E_2 $ is asymptotically stable;

    (b) If $ {R}_{10} < {R}_{20} $ and $ \hat{R}_{10} > 1 $, $ E_2 $ is unstable;

    (c) If $ {R}_{20} < {R}_{10} $ and $ \hat{R}_{20} < 1 $, then the boundary equilibrium $ E_1 $ is asymptotically stable;

    (d) If $ {R}_{20} < {R}_{10} $ and $ \hat{R}_{20} > 1 $, then $ E_1 $ is unstable.

    $\mathbf{(v)}$ If $ \hat{R}_{10} > 1 $ and $ \hat{R}_{20} > 1 $, then there is the interior equilibrium $ E_{\ast} $.

    In the above discussion, we have only shown the local asymptotical stability of the DFE $ E_0 $ when $ \max\{R_{10}, R_{20}\} < 1 $. By using the fluctuation lemma (see, e.g., [23]) and a comparison argument, we actually can prove that $ E_0 $ is indeed globally asymptotically stable for this case, as demonstrated below.

    Theorem 6.3. If $ \max\{R_{10}, R_{20}\} < 1 $, then the disease-free equilibrium $ E_0 $ is globally asymptotically stable for (6.2).

    Proof. We only need to show that every nonnegative solution of (6.2) converges to $ E_0 $. Following the convention, we use the following notations: for a continuous and bounded function $ f(t) $ defined on $ [0, \infty) $,

    $ f^{\infty}\triangleq \limsup\limits_{t\rightarrow\infty} f(t), \; \; \; and\; \; \; f_{\infty}\triangleq \liminf\limits_{t\rightarrow\infty} f(t). $

    Now, let $ (S_1(t), S_2(t), I_1(t), I_2(t)) $ be any non-negative solution of (6.2). Comparison theorem leads to (see (4.13) in Section 4)

    $ 0S1S1S(0)1,0S2S2S(0)2. $ (6.5)

    Also, by Theorem 3.1, we know that

    $ 0I1I1<,0I2I2<. $ (6.6)

    On the other hand, by the fluctuation lemma (see, e.g., [23]), there is a sequence $ t_n $ with $ t_n\rightarrow \infty $ as $ n \rightarrow \infty $ such that

    $ I_1(t_n)\rightarrow I_1^{\infty}\; \; \; and\; \; \; I'_1(t_n)\rightarrow 0, \; \; \; as \; n \rightarrow \infty. $

    Substituting the sequence $ t_n $ into the third equation of (6.2), letting $ n \rightarrow \infty $ and making use of (6.5), we obtain

    $ [d1+γ1(1ϵ1)]I1λ1I1S1<λ1I1S(0)1. $ (6.7)

    In a similar way, we can establish

    $ [d2+γ2(1ϵ2)]I2λ2I2S2<λ2I2S(0)2. $ (6.8)

    Under $ \max\{R_{10}, R_{20}\} < 1 $, (6.7)–(6.8) leads to $ I^{\infty}_i = 0 $, $ i = 1, \; 2 $. This together with (6.6) implies $ \lim _{t\rightarrow\infty} I_i (t) = I_{i\infty} = I^{\infty}_i = 0 $ for $ i = 1, \; 2 $. Finally, applying the theory of asymptotically autonomous systems (see, e.g., [21]) to the first and second equations of (6.2), we conclude that $ \lim _{t\rightarrow\infty} S_i (t) = S^{(0)}_i $, $ i = 1, \; 2 $. This completes the proof.

    In this subsection, we still assume positive $ D^S_{12} $ and $ D^S_{21} $. We consider a scenario that the travel of the recovered individuals is unidirectional. Without loss of generality, we assume that recovered individuals can travel from Patch $ 2 $ to Patch $ 1 $, but can not travel from Patch $ 1 $ to Patch $ 2 $. That is, we assume that $ D^I_{12} = D^I_{21} = D^R_{21} = 0 $, but $ D^R_{12} > 0 $. If the two patches are two cities, such a situation may occur when the two cities have different public health systems, or the health officials in the two cities disagree on the severity of an infectious disease, resulting in one city implementing a ban against the arrival of the recovered individuals from the other city but not vice-versa.

    In this case, the matrix $ \mathbf{B} $ is upper triangular, and so is $ [b_{ij}(\tau)] = e^{\left(\mathbf{B}\tau\right)} $, given by

    $ b11(τ)=ed1τ=ϵ1,b22(τ)=e(d2+DR12)τ,b12(τ)=ed1τ(1eDR12τ),b21(τ)=0. $

    Thus, (3.6) reduces to

    $ {dS1(t)dt=K1d1S1(t)+DS12S2(t)DS21S1(t)λ1S1(t)I1(t),dS2(t)dt=K2d2S2(t)+DS21S1(t)DS12S2(t)λ2S2(t)I2(t),dI1(t)dt=d1I1(t)+λ1S1(t)I1(t)γ1I1(t)+b11(τ)γ1I1(tτ)+b12(τ)γ2I2(tτ),dI2(t)dt=d2I2(t)+λ2S2(t)I2(t)γ2I2(t)+b22(τ)γ2I2(tτ). $ (6.9)

    The DFE $ E_0 $ is still given by (4.1). A possible boundary equilibrium of the form $ E_1 = (S^{(1)}_1, S^{(1)}_2, I^{(1)}_1, 0) $ is still given by (6.4). Hence, as is seen in Subsection 6.1, $ E_1 $ exists if and only if $ R_{10} > 1 $ where $ R_{10} $ is defined in Subsection 6.1. However, since $ b_{12}(\tau) > 0 $, a boundary equilibrium of the form $ E_2 = (S^{(2)}_1, S^{(2)}_2, 0, I^{(2)}_2) $ becomes impossible.

    For the convenience of discussing stability of the equilibria, we define

    $ R20=λ2S(0)2d2+γ2(1b22(τ)),ˆR20=λ2S(1)2d2+γ2(1b22(τ)). $

    Linearizing (6.9) at $ E_0 = (S^{(0)}_1, S^{(0)}_2, 0, 0) $ leads to

    $ {dS1(t)dt=(d1+DS21)S1(t)+DS12S2(t)λ1S(0)1I1(t),dS2(t)dt=(d2+DS12)S2(t)+DS21S1(t)λ2S(0)2(t)I2(t),dI1(t)dt=(d1+γ1)I1(t)+λ1S(0)1I1(t)+b11(τ)γ1I1(tτ)+b12(τ)γ2I2(tτ),dI2(t)dt=(d2+γ2)I2(t)+λ2S(0)2I2(t)+b22(τ)γ2I2(tτ). $ (6.10)

    The characteristic equation of (6.10)

    $ Δ1(z)Δ3(z,τ)ˆΔ4(z,τ)=0, $ (6.11)

    where $ \Delta_1(z) $ and $ \Delta_3(z, \tau) $ are as in Section 6.1, but $ \widehat{\Delta}_4(z, \tau) $ is a modification of $ \Delta_4(z, \tau) $ by the following formula:

    $ \widehat{\Delta}_4(z, \tau) = z+d_2+\gamma_2-\lambda_2S^{(0)}_2-b_{22}(\tau)\gamma_2e^{-z\tau}. $

    which is a result of replacing $ \epsilon_2 $ in Section 6.1 by $ b_{22}(\tau) $. Thus, by a similar argument to that for the stability/instability of $ E_0 $ in Section 6.1, we conclude that $ E_0 $ is locally asymptotically stable if $ \max\{R_{10}, R'_{20}\} < 1 $, and it becomes unstable if $ \max\{R_{10}, R'_{20}\} > 1 $. Actually, we can also further prove that $ E_0 $ is globally asymptotically stable if $ \max\{R_{10}, R'_{20}\} < 1 $ again by using the fluctuation lemma. In fact, for any nonnegative solution $ (S_1(t), S_2(t), I_1(t), I_2(t)) $ of (6.9), by argument similar to that in proof of Theorem 6.3, we have $ \lim _{t\rightarrow\infty} I_i (t) = 0 $ and $ \lim _{t\rightarrow\infty} S_i (t) = S^{(0)}_i $, $ i = 1, \; 2 $. This gives the globally asymptotically stability of $ E_0 $ for (6.9). Thus we have the following Theorem.

    Theorem 6.4. If $ \max\{R_{10}, R'_{20}\} < 1 $, then the disease-free equilibrium $ E_0 $ is globally asymptotically stable for (6.9); it is unstable if $ \max\{R_{10}, R'_{20}\} > 1 $.

    Next, we investigate what happens when $ \max\{R_{10}, R'_{20}\} > 1 $.

    $ Case \; 1 $: $ R_{10} > 1 $. We have seen above that in this case there is the boundary equilibrium $ E_1 $. To investigate the stability of $ E_1 $, we linearize (6.9) at $ E_1 $ to obtain

    $ {dS1(t)dt=(d1+DS21+λ1I(1)1)S1(t)+DS12S2(t)λ1S(1)1I1(t),dS2(t)dt=(d2+DS12)S2(t)+DS21S1(t)λ2S(1)2I2(t),dI1(t)dt=d1I1(t)+λ1S(1)1I1(t)+λ1I(1)1S1(t)γ1I1(t)+b11(τ)γ1I1(tτ)+b12(τ)γ2I2(tτ),dI2(t)dt=d2I2(t)+λ2S(1)2I2(t)γ2I2(t)+b22(τ)γ2I2(tτ). $ (6.12)

    Note that the $ I_1 $ and $ I_2 $ equations in (6.12) are decoupled and form a cooperative linear DDE system, and the stability of the trivial equilibrium of this subsystem is fully determined by the sign of $ -d_2+\lambda S_2^{(1)}- \gamma_2 +b_{22}(\tau) $, which is equivalently related to whether $ \widehat{R}'_{20} < 1 $ or $ \widehat{R}'_{20} > 1 $. Therefore, we actually have the following theorem.

    Theorem 6.5. Assume that $ R_{10} > 1 $ so that the boundary equilibrium $ E_1 $ exists. Then, $ E_1 $ is locally asymptotically stable if $ \widehat{R}'_{20} < 1 $; it becomes unstable if $ \widehat{R}'_{20} > 1 $. In the latter case, there is an interior equilibrium $ E_{\ast} = (S^{\ast}_1, S^{\ast}_2, I^{\ast}_1, I^{\ast}_2) $ (i.e., with $ S^{\ast}_i > 0 $, $ I^{\ast}_i > 0 $, $ i = 1, \; 2 $).

    $ Case \; 2 $: $ R_{10} < 1 $ but $ R'_{20} > 1 $. Going back to (6.11), we know that in this case, all roots of $ \Delta_1(z) = 0 $ and $ \Delta_3(z, \tau) = 0 $ have negative real parts. Thus, the stability of $ E_0 $ is totally determined by $ \widehat{\Delta}_4(z, \tau) $. Note that $ R'_{20} = 1 $ is a critical value for $ \widehat{\Delta}_4(z, \tau) = 0 $: when $ R'_{20} < 1 $, all roots of $ \widehat{\Delta}_4(z, \tau) = 0 $ have negative real parts; at $ R'_{20} = 1 $, $ z = 0 $ is a root of $ \widehat{\Delta}_4(z, \tau) = 0 $ and all other roots have negative real parts; when $ R'_{20} > 1 $, $ \widehat{\Delta}_4(z, \tau) = 0 $ has a positive real root. Thus, when $ R'_{20} $ increases to pass the critical value 1, the DFE $ E_0 $ loses its stability to another non-negative equilibrium. Since there is no boundary equilibrium, this newly bifurcated equilibrium must be an interior one. This analysis leads to the following theorem.

    Theorem 6.6. Assume that $ R_{10} < 1 $ and $ R'_{20} > 1 $. Then there is an interior equilibrium for (6.9).

    In this subsection, we use a similar way as in Subsection 6.2 to discuss the case: $ D^R_{12} = D^R_{21} = D^I_{21} = 0 $, but $ D^I_{12} > 0 $.

    In this case, the matrix $ [b_{ij}(\tau)] = e^{\left(\mathbf{B}\tau\right)} $ is given by

    $ b11(τ)=ed1τ=ϵ1,b22(τ)=ed2τ=ϵ2,b12(τ)=b21(τ)=0. $

    Thus, (3.6) reduces to

    $ {dS1(t)dt=K1d1S1(t)+DS12S2(t)DS21S1(t)λ1S1(t)I1(t),dS2(t)dt=K2d2S2(t)+DS21S1(t)DS12S2(t)λ2S2(t)I2(t),dI1(t)dt=d1I1(t)+λ1S1(t)I1(t)γ1I1(t)+DI12I2(t)+ϵ1γ1I1(tτ),dI2(t)dt=d2I2(t)+λ2S2(t)I2(t)γ2I2(t)DI12I2(t)+ϵ2γ2I2(tτ). $ (6.13)

    The DFE $ E_0 $ is still given by (4.1). A possible boundary equilibrium of the form $ E_1 = (S^{(1)}_1, S^{(1)}_2, I^{(1)}_1, 0) $ is still given by (6.4). Hence, as is seen in Subsection 6.1, $ E_1 $ exists if and only if $ R_{10} > 1 $ where $ R_{10} $ is defined in Subsection 6.1. However, since $ D^I_{12} > 0 $, a boundary equilibrium of the form $ E_2 = (S^{(2)}_1, S^{(2)}_2, 0, I^{(2)}_2) $ becomes impossible.

    Similar to the two composed parameters $ R'_{20} $ and $ \widehat{R'}_{20} $ for (6.9) in Subsection 6.2, the following two new composed parameters play a key role for (6.13):

    $ R20=λ2S(0)2d2+γ2(1ϵ2)+DI12,ˆR20=λ2S(1)2d2+γ2(1ϵ2)+DI12. $

    Parallel to the three theorems for (6.9) in Section 6.2, we can also obtain the following results for (6.13).

    Theorem 6.7. If $ \max\{R_{10}, R''_{20}\} < 1 $, then the disease-free equilibrium $ E_0 $ is globally asymptotically stable for (6.13); it is unstable if $ \max\{R_{10}, R''_{20}\} > 1 $.

    When $ \max\{R_{10}, R''_{20}\} > 1 $, we have the following two theorems, parallel to Theorems 6.5 and 6.6:

    Theorem 6.8. Assume that $ R_{10} > 1 $ so that the boundary equilibrium $ E_1 $ exists. Then, $ E_1 $ is locally asymptotically stable if $ \widehat{R}''_{20} < 1 $; it becomes unstable if $ \widehat{R}''_{20} > 1 $. In the latter case, there is an interior equilibrium $ E_{\ast} = (S^{\ast}_1, S^{\ast}_2, I^{\ast}_1, I^{\ast}_2) $ (i.e., with $ S^{\ast}_i > 0 $, $ I^{\ast}_i > 0 $, $ i = 1, \; 2 $).

    Theorem 6.9. Assume that $ R_{10} < 1 $ and $ R''_{20} > 1 $. Then there is an interior equilibrium for (6.13).

    The proofs for the above three theorems are very much similar to those for Theorems 6.4, 6.5 and 6.9, and thus, we omit them to save space.

    We have derived a new epidemic model in a 2-patch environment to describe the transmission dynamics of a disease for which the infectious individuals, once recovered for a period of fixed length, will relapse back to the infectious class. The derivation makes use of the McKendrick-von Foerster equation with the structure variable being the recovery age (the time elapsed since recovery), incorporated with the dispersals between the patches. By tracking the dispersals of recovered individuals, we have obtained a new model in the form of a system of delay differential equations which, in addition to the linear dispersion terms, contains non-local reverting terms in dynamical equations of the infectious class. The patches can be communities, cities, regions and even countries; and the population dispersals among patches can be interpreted as the movements by which people travel or migrate between patches.

    For this new model (2.12)–(2.13), we have justified the well-posedness by proving the positivity and boundedness of solutions. When all the travel rate matrices are assumed to be irreducible, we have identified concrete conditions for existence and stability/instability of the equilibria for the model. We have shown that if the inequalities (4.6) holds, then the disease dies out and when (4.7) is satisfied, the disease persists globally, (i.e., in these two patches). leading to the existence of an endemic equilibrium. When allowing infection and recovered travel rate matrices to be reducible, we have considered three special cases in Section 6. One important difference is that without the irreducibility of the travel rate matrices, the model may allow boundary equilibrium. For all of these three cases, we have also identified the threshold numbers $ R_{i0}, \, \, i = 1, 2 $, $ R'_{20} $ and $ R''_{20} $ for these three special cases in Sections 6.1, 6.2 and 6.3, respectively.

    Based on the mathematical results, we may discuss the impact of the dispersals on the disease dynamics. To demonstrate, let us take the results in Section 6.1 for (6.2) as an example. Firstly, from Theorems 6.2 and 6.3, we see that $ R_{i0} = 1 $ is the threshold value for the disease to persist in Patch-$ i $. It is thus interesting to compare these two values ($ R_{10} $ and $ R_{20} $) with $ R^{(0)}_{10} $ and $ R^{(0)}_{20} $, the basic reproduction numbers for patch $ 1 $ and patch $ 2 $ respectively when the two patches are disconnected. Indeed, it is easily seen that

    $ R10=λ1d1+γ1(1b11(τ))K1d1d2+DS12+K2K1DS12d2+DS12+d2d1DS21=R(0)10d2+DS12+K2K1DS12d2+DS12+d2d1DS21, $ (7.1)

    and

    $ R20=λ2d1+γ1(1b22(τ))K2d2d1+DS21+K1K2DS21d1+DS21+d1d2DS12=R(0)20d1+DS21+K1K2DS21d1+DS21+d1d2DS12. $ (7.2)

    It is obvious from the above formulas that $ R_{10} $ and $ R_{20} $ reflect the influence of travel of susceptible individuals between the two patches, and hence may be called the travel mediated basic reproduction numbers for patch 1 and patch 2 respectively.

    The following observations are direct consequences of (7.1)–(7.2) and their verifications are straightforward and thus, are omitted.

    $\mathbf{(O1)}$ Assume $ R^{(0)}_{10} < 1 $ and $ R^{(0)}_{20} < 1 $. If $ D^S_{12} > 0 $ and $ D^S_{21} > 0 $ satisfy either

    $ DS12>d2(1R(0)10)+d2d1DS21R(0)10(1+K2K1)1with1>R(0)10>K1K1+K2; $ (7.3)

    or

    $ DS21<d1d2[(R(0)101)(d2+DS12)+R(0)10DS12K2K1]with1>R(0)10>d2+DS12d2+DS12+K2K1DS12, $ (7.4)

    then $ R_{10} > 1 $ and $ R_{20} < 1 $. By symmetry, the conditions parallel to (7.3) or (7.4) can lead to $ R_{10} < 1 $ and $ R_{20} > 1 $. Here and in the sequel in this section, we omit all such parallel statements and the corresponding conditions, as they can be easily obtained by switching the two patches.

    $\mathbf{(O2)}$ Assume $ R^{(0)}_{10} > 1 $ and $ R^{(0)}_{20} > 1 $. If $ D^S_{12} > 0 $ and $ D^S_{21} > 0 $ satisfy either

    $ DS12<d2(1R(0)10)+d2d1DS21R(0)10(1+K2K1)1with1<R(0)10<1+DS21d1; $

    or

    $ DS21>d1d2[(R(0)101)(d2+DS12)+R(0)10DS12K2K1], $

    then $ R_{10} < 1 $ and $ R_{20} > 1 $.

    $\mathbf{(O3)}$ Assume $ R^{(0)}_{10} < 1 $ and $ R^{(0)}_{20} > 1 $. If $ D^S_{12} > 0 $ and $ D^S_{21} > 0 $ satisfy either

    $ d2(1R(0)10)+d2d1DS21R(0)10(1+K2K1)1<DS12<d2d1[(R(0)201)(d1+DS21)+R(0)20DS21K2K1]with1>R(0)10>K1K1+K2; $

    or

    $ d1(1R(0)20)+d1d2DS12R(0)20(1+K1K2)1<DS21<d1d2[(R(0)101)(d2+DS12)+R(0)10DS12K2K1]with1>R(0)10>d2+DS12d2+DS12+K2K1DS12and1<R(0)20<1+DS12d2, $

    then $ R_{10} > 1 $ and $ R_{20} > 1 $.

    $\mathbf{(O4)}$ Assume $ R^{(0)}_{10} < 1 $ and $ R^{(0)}_{20} > 1 $. If $ D^S_{12} > 0 $ and $ D^S_{21} > 0 $ satisfy either

    $ d2d1[(R(0)201)(d1+DS21)+R(0)20DS21K2K1]<DS12<d2(1R(0)10)+d2d1DS21R(0)10(1+K2K1)1with1>R(0)10>K1K1+K2; $

    or

    $ d1d2[(R(0)101)(d2+DS12)+R(0)10DS12K2K1]<DS21<d1(1R(0)20)+d1d2DS12R(0)20(1+K1K2)1with1>R(0)10>d2+DS12d2+DS12+K2K1DS12and1<R(0)20<1+DS12d2, $

    then $ R_{10} > 1 $ but $ R_{20} < 1 $.

    The biological implications of $ \mathbf{(O1)} $–$ \mathbf{(O4)} $ can be explained as follows. $ \mathbf{(O1)} $ implies that travel of the susceptible individuals can help an otherwise dying out disease persist locally. In plain language, a larger inflow of susceptible individuals into a patch will enhance the chance of disease persistence in that patch. $ \mathbf{(O2)} $ implies that travel of the susceptible individuals can also help drive an otherwise globally persistent disease out of a patch. $ \mathbf{(O3)} $ and $ \mathbf{(O4)} $ show that appropriate travel rates may either cause an otherwise partially persistent disease to go to full extinction, or help it persist globally in both patches.

    Similarly, we may explore the impact of the travel of infectious and recovered individuals in the model by using the results, e.g., for the special cases in Sections 6.2 and 6.3. Indeed, from the formulations of $ R'_{20} $ and $ R''_{20} $, we can find that $ R'_{20} $ and $ R''_{20} $ are decreasing functions of $ D^R_{12} $ (the travel rate from Patch 2 to Patch 1 for the recovered individuals) and $ D^I_{12} $ (the travel rate from Patch 2 to Patch 1 for the infected individuals), respectively, so are $ \max \{R_{10}, R'_{20}\} $ and $ \max \{R_{10}, R''_{20}\} $. For example, when we have $ R_{10} < 1 $ and $ R'_{20} > 1 $ which gives $ \max \{R_{10}, R'_{20}\} > 1 $, the increase of $ D^R_{12} $ (the unbalanced travel rate from Patch $ 2 $ to Patch $ 1 $ for the recovered class) will decrease $ R'_{20} $ to a value less than $ 1 $, which results in $ \max \{R_{10}, R'_{20}\} < 1 $. Therefore, $ D^R_{12} $ indeed plays a role of decreasing the threshold number $ \max \{R_{10}, R'_{20}\} $, which is similar to the role of the travel rate of the infected individuals, but differs from the role of the travel rate of the susceptible individuals. More discussions can be expanded, as $ \max \{R_{10}, R'_{20}\} $ also depends on $ D^S_{ij} $ through $ S_2^{(0)} $, however we decide to skip such expansion in this already lengthy paper.

    Finally, we point out that, at the present we are unable to prove the stability of the endemic equilibria when it exists. This seems to be a very challenging mathematical problem due to the presence of the relapse delay and the non-locality in the model. We leave it as a future project.

    The authors are grateful to the two anonymous referees for their valuable comments which have lead to an substantial improvement in the presentation of the paper.

    The first author is supported by China Scholarship Council and NUPTSF (No. NY220093). The corresponding author is supported by NSERC of Canada (No. RGPIN-2016-04665).

    The authors have declared that no competing interests exist.

    [1] G. Ali, V. Furuholt, R. Natalini and I. Torcicollo, A mathematical model of sulphite chemical aggression of limestones with high permeability. Part I. Modeling and qualitative analysis, Transport in Porous Media, 69 (2007), 109-122. doi: 10.1007/s11242-006-9067-2
    [2] G. Ali, V. Furuholt, R. Natalini and I. Torcicollo, A mathematical model of sulphite chemical aggression of limestones with high permeability. Part II: Numerical approximation, Transport in Porous Media, 69 (2007), 175-188. doi: 10.1007/s11242-006-9068-1
    [3] D. Aregba-Driollet, F. Diele and R. Natalini, A Mathematical Model for the SO2 Aggression to Calcium Carbonate Stones: Numerical Approximation and Asymptotic Analysis, SIAM J. APPL. MATH. , 64 (2004), 1636-1667. doi: 10.1137/S003613990342829X
    [4] F. Clareli, A. Fasano and R. Natalini, Mathematics and monument conservation: Free boundary models of marble sulfation, SIAM Journal on Applied Mathematics, 69 (2008), 149-168. doi: 10.1137/070695125
    [5] A. Fasano and R. Natalini, Lost Beauties of the Acropolis: What Mathematics Can Say, SIAM news, 2006.
    [6] T. Fatima, Multiscale Reaction Diffusion Systems Describing Concrete Corrosion: Modelling and Analysis, Ph.D thesis, Technical University of Eindhoven, 2013.
    [7] T. Fatima, N. Arab, E. P. Zemskov and A. Muntean, Homogenization of a reaction - diffusion system modeling sulfate corrosion of concrete in locally periodic perforated domains, Journal of Engineering Mathematics, 69 (2011), 261-276. doi: 10.1007/s10665-010-9396-6
    [8] T. Fatima and A. Muntean, Sulfate attack in sewer pipes: Derivation of a concrete corrosion model via two-scale convergence, Nonlinear Analysis: Real World Applications, 15 (2014), 326-344. doi: 10.1016/j.nonrwa.2012.01.019
    [9] T. Fatima, A. Muntean and T. Aiki, Distributed space scales in a semilinear reaction-diffusion system including a parabolic variational inequality: A well-posedness study, Adv. Math. Sci. Appl., 22 (2012), 295-318.
    [10] T. Fatima, A. Muntean and M. Ptashnyk, Unfolding-based corrector estimates for a reaction - diffusion system predicting concrete corrosion, Applicable Analysis, 91 (2012), 1129-1154. doi: 10.1080/00036811.2011.625016
    [11] F. R. Guarguaglini and R. Natalini, Fast reaction limit and large time behavior of solutions to a nonlinear model of sulphation phenomena, Commun. Partial Differ. Equations, 32 (2007), 163-189. doi: 10.1080/03605300500361438
    [12] F. R. Guarguaglini and R. Natalini, Global existence of solutions to a nonlinear model of sulphation phenomena in calcium carbonate stones, Nonlinear Analysis: Real World Applications, 6 (2005), 477-494. doi: 10.1016/j.nonrwa.2004.09.007
    [13] E. J. Hinch, Perturbation Methods, Cambridge University Press, 1991. doi: 10.1017/CBO9781139172189
    [14] A. A. Lacey and L. A. Herraiz, Macroscopic models for melting derived from averaging microscopic Stefan problems I: Simple geometries with kinetic undercooling or surface tension, Euro. Jnl. of Applied Mathematics, 11 (2002), 153-169. doi: 10.1017/S0956792599004027
    [15] A. A. Lacey and L. A. Herraiz, Macroscopic models for melting derived from averaging microscopic Stefan problems II: Effect of varying geometry and composition, Euro. Jnl. of Applied Mathematics, 13 (2002), 261-282. doi: 10.1017/S0956792501004818
    [16] R. J. Leveque, Finite Volume Methods for Hyperbolic Problems, Caimbridge University Press, 2002. doi: 10.1017/CBO9780511791253
    [17] C. V. Nikolopoulos, A mushy region in concrete corrosion, Applied Mathematical Modelling, 34 (2010), 4012-4030. doi: 10.1016/j.apm.2010.04.005
    [18] C. V. Nikolopoulos, Macroscopic models for a mushy region in concrete corrosion, Journal of Engineering Mathematics, 2014, DOI 10.1007/s10665-014-9743-0.
    [19] J. L. Schnoor, Enviromental Modeling, Fate and transport of pollutants in water, air, and soil, John Willey and Sons, Inc., 1996.
  • This article has been cited by:

    1. Tingting Zheng, Huaiping Zhu, Zhidong Teng, Linfei Nie, Yantao Luo, Patch model for border reopening and control to prevent new outbreaks of COVID-19, 2023, 20, 1551-0018, 7171, 10.3934/mbe.2023310
    2. Yueding Yuan, Zhiming Guo, Global dynamics of a class of delayed differential systems with spatial non-locality, 2023, 349, 00220396, 176, 10.1016/j.jde.2022.12.013
    3. Jian Liu, Qian Ding, Hongpeng Guo, Bo Zheng, DYNAMICS OF AN EPIDEMIC MODEL WITH RELAPSE AND DELAY, 2024, 14, 2156-907X, 2317, 10.11948/20230376
    4. Yuhang Li, Yongzheng Sun, Maoxing Liu, Analysis of a patch epidemic model incorporating population migration and entry–exit screening, 2024, 14, 2158-3226, 10.1063/5.0196679
    5. Dongxue Yan, Yu Cao, Rich dynamics of a delayed SIRS epidemic model with two-age structure and logistic growth, 2023, 2023, 2731-4235, 10.1186/s13662-023-03794-0
    6. Jimmy Calvo-Monge, Jorge Arroyo-Esquivel, Alyssa Gehman, Fabio Sanchez, Source-Sink Dynamics in a Two-Patch SI Epidemic Model with Life Stages and No Recovery from Infection, 2024, 86, 0092-8240, 10.1007/s11538-024-01328-7
  • Reader Comments
  • © 2014 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

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(3879) PDF downloads(64) Cited by(5)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog