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

1.   Introduction

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

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

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) $:

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

and accordingly, (1.3) reduces to

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

and accordingly, (1.3) splits to

and

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.

2.   Model formulation

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

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:

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

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

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:

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

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

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

where

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

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

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

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

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:

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:

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

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

where

Hence the model becomes

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

Then

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

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

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

Thus,

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,

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.

3.   Well-posedness

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:

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

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

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

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:

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

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

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

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

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

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

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

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

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):

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

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

4.   Disease-free equilibrium and its stability

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

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

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

The characteristic equation of (4.2) is given by

where

and

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

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

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

with

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 $:

Hence, if

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

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.

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

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

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

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):

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.

5.   Disease persistence and endemic equilibrium

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

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

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

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

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

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

Let us consider the following linear system

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

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

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 $,

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

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

and hence

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.

6.   Allowing reducible travel rate matrices

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

Thus, (3.6) is decoupled to

for patch 1, and

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

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

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.

6.1.   Sub-case 1: Only susceptible individuals travel

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:

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

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

The characteristic equation of (6.3) is given by

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

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

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

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

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

Similarly, for $ E_2 $ we have

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

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,

Solving these equations for positive components leads to

Define

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:

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

Remark 6.1. Direct computations show that

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) $,

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)

Also, by Theorem 3.1, we know that

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

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

In a similar way, we can establish

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.

6.2.   Sub-case 2: Travel of recovered individuals is unidirectional

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

Thus, (3.6) reduces to

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

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

The characteristic equation of (6.10)

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:

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

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).

6.3.   Sub-case 3: Travel of infected individuals is unidirectional

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

Thus, (3.6) reduces to

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):

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.

7.   Discussions

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

and

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

or

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

or

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

or

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

or

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.

Acknowledgments

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).

Conflict of interest

The authors have declared that no competing interests exist.