Proper water resource management is essential for maintaining a sustainable supply chain and meeting water demand. The urgent need to preserve river ecosystems by sustaining environmental flow (EF) in the realm of environmental management has been highlighted by the drastic changes to river ecosystems and upstream flow dynamics brought about by careless river exploitation in the last few decades. To optimize EF in river basin management, we present an integrated modeling approach. We focused on the Pir Khezran River basin. Our objective was to estimate EF and generalize the findings to adjacent rivers using modeling techniques, thus providing valuable insights for environmental management applications. The assessment and optimization of EF under uncertain conditions was achieved by combining physical habitat simulation (PHABSIM) modeling with advanced techniques like Adaptive Neuro-Fuzzy Inference Systems (ANFIS) and Multilayer Perceptron (MLP) neural networks. This integrated modeling approach contributes to sustainable solutions for river basin management and environmental conservation by effectively optimizing EF, as demonstrated by the results. This research, therefore, makes valuable contributions to environmental management in various areas such as ecological preservation, modeling and optimizing environmental systems, and policy considerations.
Citation: Seiran Haghgoo, Jamil Amanollahi, Barzan Bahrami Kamangar, Shahryar Sorooshian. Decision models enhancing environmental flow sustainability: A strategic approach to water resource management[J]. AIMS Environmental Science, 2024, 11(6): 900-917. doi: 10.3934/environsci.2024045
[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 |
Proper water resource management is essential for maintaining a sustainable supply chain and meeting water demand. The urgent need to preserve river ecosystems by sustaining environmental flow (EF) in the realm of environmental management has been highlighted by the drastic changes to river ecosystems and upstream flow dynamics brought about by careless river exploitation in the last few decades. To optimize EF in river basin management, we present an integrated modeling approach. We focused on the Pir Khezran River basin. Our objective was to estimate EF and generalize the findings to adjacent rivers using modeling techniques, thus providing valuable insights for environmental management applications. The assessment and optimization of EF under uncertain conditions was achieved by combining physical habitat simulation (PHABSIM) modeling with advanced techniques like Adaptive Neuro-Fuzzy Inference Systems (ANFIS) and Multilayer Perceptron (MLP) neural networks. This integrated modeling approach contributes to sustainable solutions for river basin management and environmental conservation by effectively optimizing EF, as demonstrated by the results. This research, therefore, makes valuable contributions to environmental management in various areas such as ecological preservation, modeling and optimizing environmental systems, and policy considerations.
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 γ 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,∞)→[0,∞) is differentiable except at possibly finite many points where it may have jump discontinuities, non-increasing and satisfies P(0)=1,P(∞)=0 and ∫∞0P(t)dt positive and finite.
Then, the population of the recovered class at time t is given by
R(t)=∫t0γI(ξ)e−d(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(ξ)e−d(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)=K−dS(t)−λS(t)I(t),I′(t)=λS(t)I(t)−dI(t)−γI(t)−∫t0γI(ξ)e−d(t−ξ)dtP(t−ξ)dξ,R′(t)=γI(t)−dR(t)+∫t0γI(ξ)e−d(t−ξ)dtP(t−ξ)dξ. | (1.3) |
Here, a simple demographic dynamics S′(t)=K−dS(t) and a mass action infection mechanism λ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≥0 where r>0 is a constant;
(II) step function, i.e., P(t)=1 for t∈[0,τ) and P(t)=0 for t≥τ, where τ>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 τ>0.
With choice (I), the integral term in (1.2) becomes
∫t0γI(ξ)e−d(t−ξ)dtP(t−ξ)dξ=−rR(t), |
and accordingly, (1.3) reduces to
{S′(t)=K−dS(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 τ whose derivative is the Dirac delta function at τ, i.e., H′(t)=δ(t−τ). Hence dtP(t−ξ)=−δ(t−ξ−τ) and therefore, the integral in (1.2) becomes
∫t0γI(ξ)e−d(t−ξ)dtP(t−ξ)dξ=−γI(t−τ)e−dτ, |
and accordingly, (1.3) splits to
{S′(t)=K−dS(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)=K−dS(t)−λS(t)I(t),I′(t)=λS(t)I(t)−dI(t)−γI(t)+e−dτγI(t−τ),R′(t)=γI(t)−dR(t)−e−dτγ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,⋯,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 Si(t), Ii(t), Ri(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 ri(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=τ. 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 ri(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 DRijrj(t,a) corresponds to the dispersal of the recovered individuals at the recovery age a from patch j to patch i; constant di>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 ri(0,a)=0 for a∈(0,τ].
Obviously, ri(t,0) corresponds to the new recovery individuals in patch i who come from the infectious individuals. Assuming a constant recovery rate γ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 τ 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 Ki>0 in patch i, i=1,2, and a constant natural death rate for each class denoted still by di>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=K1−d1S1(t)+DS12S2(t)−DS21S1(t)−λ1S1(t)I1(t),dS2(t)dt=K2−d2S2(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 DSij is the rate at which susceptible individuals migrate from patch j to patch i, i≠j, and DIij is the rate at which infectious individuals migrate from patch j to patch i, i≠j. Here λi>0, i=1,2 are the transmission rate in patch i. Note that the equations for the recovered class Ri, i=1,2 are decoupled from the equations for Si and Ii, i=1,2. Thus we only need to consider the 4 equations for Si and Ii, i=1,2 in (2.5).
Obviously, ri(t,τ) is the rate at which patch i gains relapse individuals, which can be determined below in terms of Ij(t) for j=1,2.
For fixed ξ≥0, let
Vξi(t)=ri(t,t−ξ),forξ≤t≤ξ+τandi=1,2. |
Then for 1≤i≠j≤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 Vξ(t)=(Vξ1(t),Vξ2(t))⊤, where ⊤ represents the transpose of a vector. Then Vξ(t) satisfies
ddtVξ(t)=BVξ(t), | (2.7) |
where
B=(−d1−DR21DR12DR21−d2−DR12). |
Integrating (2.7) with respect to t from ξ to t, we have
Vξ(t)=exp(B(t−ξ))(Vξ1(ξ),Vξ2(ξ))⊤,ξ≤t≤ξ+τ. | (2.8) |
By using the definition of Vξi(t) and (2.3),
Vξ(t)=exp(B(t−ξ))(r1(ξ,0),r2(ξ,0))⊤,ξ≤t≤ξ+τ=exp(B(t−ξ))(γ1I1(ξ),γ2I2(ξ))⊤. | (2.9) |
For t≥τ (hence t−τ≥0), letting r(t,τ)=(r1(t,τ),r2(t,τ))⊤, we obtain
r(t,τ)=Vt−τ(t)=exp(Bτ)(γ1I1(t−τ),γ2I2(t−τ))⊤. | (2.10) |
Denoting [bij(τ)]2×2:=exp(Bτ), it follows that
ri(t,τ)=2∑j=1bij(τ)γjIj(t−τ). | (2.11) |
Substituting ri(t,τ) back into the Ii equations in (2.5) and taking out the first 4 equations for Si, and Ii, i=1,2, results in the following new model:
{dS1(t)dt=K1−d1S1(t)+DS12S2(t)−DS21S1(t)−λ1S1(t)I1(t),dS2(t)dt=K2−d2S2(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)+2∑j=1b1j(τ)γjIj(t−τ),dI2(t)dt=−d2I1(t)+DI21I1(t)−DI12I2(t)+λ2S2(t)I2(t)−γ2I2(t)+2∑j=1b2j(τ)γjIj(t−τ). | (2.12) |
For 0<t≤τ, 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=K1−d1S1(t)+DS12S2(t)−DS21S1(t)−λ1S1(t)I1(t),dS2(t)dt=K2−d2S2(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 Ii equation in (2.12) accounts for non-local reverting force, reflecting how the individuals recovered τ 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 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 d1=d2=d, then
B=(−d−DR21DR12DR21−d−DR12)=(−d00−d)+(−DR21DR12DR21−DR12), |
and we can obtain [bij(τ)]=exp(Bτ) explicitly as
b11(τ)=e−dτ(1−δ1(τ)),b12(τ)=e−dτδ2(τ),b22(τ)=e−dτ(1−δ2(τ)),b21(τ)=e−dτδ1(τ), | (2.14) |
where
δi(τ)=DRjiDRji+DRij(1−e−(DRji+DRij)τ),for1≤i≠j≤2. | (2.15) |
Hence the model becomes
{dS1(t)dt=K1−d1S1(t)+DS12S2(t)−DS21S1(t)−λ1S1(t)I1(t),dS2(t)dt=K2−d2S2(t)+DS21S1(t)−DS12S2(t)−λ2S2(t)I2(t),dI1(t)dt=−d1I1(t)+DI12I2(t)−DI21I1(t)+λ1S1(t)I1(t)−γ1I1(t)+e−dτ(1−δ1(τ))γ1I1(t−τ)+e−dτδ2(τ)γ2I2(t−τ),dI2(t)dt=−d2I2(t)+DI21I1(t)−DI12I2(t)+λ2S2(t)I2(t)−γ2I2(t)+e−dτ(1−δ2(τ))γ2I2(t−τ)+e−dτδ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τ. Due to the mobility during the recovered period between the two patches, τ time units later, a survived recovered individual may relapse in patch 1 with probability (1−δ1(τ)) or in patch 2 with probability δ1(τ). This explains the term e−dτ[1−δ1(τ)]γ1I1(t−τ) in the I1 equation and the term e−dτδ1(τ)γ1I1(t−τ) in the I2 equation. The terms e−dτ(1−δ2(τ))γ2I2(t−τ) in I2 equation and the term e−dτδ2(τ)γ2I2(t−τ) in I1 equation are explained similarly. Alternatively, we may explain these terms in light of fractions as below: among the individuals recovered in the first patch τ time units ago, a fraction e−dτ can survive the relapse period, and a fraction (1−δ1(τ)) of these survived individuals is now still in patch 1 while the remaining fraction δ1(τ) of them has now moved to patch 2.
For bij(τ), one should expect the relation 0≤bij(τ)≤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
d_=min{d1,d2},andˉd=max{d1,d2}. |
Then
e−ˉdτ≤2∑i=1bij(τ)≤e−d_τ,forj=1,2. | (2.17) |
Proof. Choose a constant H>0 sufficiently large such that
H>max{d1τ+DR21τ,d2τ+DR12τ}. |
Write Bτ as Bτ=−HE+HE+D0+D1, where E is the 2×2 identity matrix and
D0=(−d1τ00−d2τ),D1=(−DR21τDR12τDR21τ−DR12τ). |
Let D_=−d_τE. Then both HE+D0+D1 and HE+D_+D1 are nonnegative matrices and
HE+D0+D1≤HE+D_+D1. |
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 VD1=0, and hence Vexp(D1)=VE. 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∈[0,τ] and a system of DDEs (2.12) for t≥τ. 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≥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 (S01(t),S02(t),I01(t),I02(t)), which exists globally, remains non-negative and is bounded. Consider the restriction of this solution on [0,τ] and denote its components by
ϕi(θ)=S0i(θ),andψi(θ)=I0i(θ),fori=1,2,andθ∈[0,τ]. |
Then, ϕi(θ) and ψi(θ) are continuous and non-negative functions on [0,τ]. By the fundamental theory of delay differential equations (see, e.g., [18]), we know that the DDE system (2.12) with the initial conditions
Si(θ)=ϕi(θ)andIi(θ)=ψi(θ),fori=1,2, |
has a unique solution (S(t,ϕ,ψ),I(t,ϕ,ψ)), which is well-defined on its maximal interval of existence [τ,tmax(ϕ,ψ)), where
(S(t,ϕ,ψ),I(t,ϕ,ψ)):=(S1(t,ϕ,ψ),S2(t,ϕ,ψ),I1(t,ϕ,ψ),I2(t,ϕ,ψ)), |
(ϕ,ψ):=(ϕ1(θ),ϕ2(θ),ψ1(θ),ψ2(θ)). |
Firstly, we show the non-negativity of the solution for t∈[τ,tmax(ϕ,ψ)). 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 S(t)=(S1(t),S2(t))⊤, I(t)=(I1(t),I2(t))⊤ and K=(K1,K2)⊤, and
D(t)=(−d1−DS21−λ1I1(t)DS12DS21−d2−DS12−λ2I2(t)),A=(b11(τ)γ1b12(τ)γ2b21(τ)γ1b22(τ)γ2), |
C(t)=(−d1−DI21+λ1S1(t)−γ1DI12DI21−d2−DI12+λ2S2(t)−γ2). |
Noting that the off-diagonal elements of matrix D(t) are non-negative, we conclude that the entries of the matrix e∫tτD(ξ)dξ are all nonnegative. Indeed, let
G(t)=max{d1+DS21+λ1I1(t)+1,d2+DS12+λ2I2(t)+1} |
and rewrite D(t) as
D(t)=(−G(t)00−G(t))+(G(t)−d1−DS21−λ1I1(t)DS12DS21G(t)−d2−DS12−λ2I2(t))≜−G(t)E+ˉD(t). |
Then all entries of ˉD(t) are nonnegative, and hence, so are the entries of e∫tτˉD(ξ)dξ. It is obvious that
e∫tτ(−G(ξ)E)dξ=(e∫tτ(−G(ξ))dξ00e∫tτ(−G(ξ))dξ). |
Noting that the scalar matrix −G(t)E commutes with any 2×2 matrix (hence with ˉD(t)), we have
e∫tτD(ξ)dξ=e∫tτ(−G(ξ)E)dξ⋅e∫tτˉD(ξ)dξ, |
implying that all entries of e∫tτD(ξ)dξ are nonnegative. Now from (3.2), we have
S(t)=e∫tτD(ξ)dξS(τ)+∫tτKe∫t−sτD(ξ)dξds≥0,fort∈[τ,tmax(ϕ,ψ)). | (3.4) |
Similarly, for any t≥τ, all entries of e∫tτC(ξ)dξ are nonnegative. Moreover, it is obvious that all entries of A are all non-negative. Now, (3.3) leads to
I(t)=e∫tτC(ξ)dξI(τ)+∫tτe∫t−sτC(ξ)dξAI(s−τ)ds,t≥τ, | (3.5) |
implying I(t)≥0 for t∈[τ,2τ] from the initial condition Ii(θ)≥0 for θ∈[0,τ] and i=1,2. This and (3.5) ensure I(t)≥0 for t∈[2τ,3τ]. By induction, we then conclude that I(t)≥0 for t∈[τ,tmax(ϕ,ψ)).
Now, we show that Si(t), Ii(t) and Ri(t) are bounded for t∈[τ,tmax(ϕ,ψ)) and i=1,2. Noting that, by using the method of characteristic lines for the model (2.2), we can derive that ri(t,a)≥0, as well as Ri(t)≥0, i=1,2. Let N(t)=S1(t)+S2(t)+I1(t)+I2(t)+R1(t)+R2(t). Then from System (2.5), we have
ddtN(t)=K1+K2−d1S1(t)−d2S2(t)−d1I1(t)−d2I2(t)−d1R1(t)−d2R2(t)≤K1+K2−min{d1,d2}N(t). |
This implies that N(t) is bounded, and so are Si(t), Ii(t) and Ri(t) for i=1,2 and t∈[τ,tmax(ϕ,ψ)). By the theory of continuation of solutions (see, e.g., [18]), we conclude that tmax(ϕ,ψ)=∞, which means the solution (S(t,ϕ,ψ),I(t,ϕ,ψ)) exists globally. This together with the results on S0i(t) and I0i(t) on t∈[0,τ) implies that all of the above results actually hold for all t≥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 ψ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,τ] 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 τ to 0 and accordingly consider initial functions on the interval [−τ,0], leading to the following equivalent system for (2.12):
{dS1(t)dt=K1−d1S1(t)+DS12S2(t)−DS21S1(t)−λ1S1(t)I1(t),dS2(t)dt=K2−d2S2(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)+2∑j=1b1j(τ)γjIj(t−τ),dI2(t)dt=−d2I1(t)+DI21I1(t)−DI12I2(t)+λ2S2(t)I2(t)−γ2I2(t)+2∑j=1b2j(τ)γjIj(t−τ). | (3.6) |
with the initial conditions specified in [−τ,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 [DSij], [DIij] and [DRij] 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 E0=(S(0)1,S(0)2,0,0) with S(0)=(S(0)1,S(0)2)⊤ satisfying the linear system MS(0)=K, where
M=(d1+DS21−DS12−DS21d2+DS12). |
Note that matrix M is irreducible, has positive column sums and negative off-diagonal entries. Thus M is a non-singular M-matrix (see [17], page 141) with M−1>0, and therefore, the linear system has a unique solution given by S(0)=M−1K>0. Indeed, one can explicitly solve MS(0)=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 E0. To this end, we consider the linearization of (3.6) at E0:
{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+DS21−DS12−DS21z+d2+DS12]=z2+(d1+d2+DS12+DS21)z+d1d2+d1DS12+d2DS21, |
and
Δ2(z,τ)=det[z+d1+DI21+γ1−λ1S(0)1−b11(τ)γ1e−zτ−DI12−b12(τ)γ2e−zτ−DI21−b21(τ)γ1e−zτz+d2+DI12+γ2−λ2S(0)2−b22(τ)γ2e−zτ]. |
It is obvious that all roots of Δ1(z) have negative real parts. Thus, the stability of E0 is fully determined by the roots of Δ2(z,τ). Note that the I equations of the linearization (4.2) is decoupled from the S equations and Δ2(z,τ)=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 A is defined in Theorem 3.1, and
F=(−d1−DI21+λ1S01−γ1DI12DI21−d2−DI12+λ2S02−γ2). |
Note that A and 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−τ) with I(t) in (4.4). This leads to the system
I′(t)=(F+A)I(t)≜WI(t), | (4.5) |
with
W=(−d1−DI21+λ1S(0)1−γ1+b11(τ)γ1DI12+b12(τ)γ2DI21+b21(τ)γ1−d2−DI12+λ2S(0)2−γ2+b22(τ)γ2)≜(w11w12w21w22). |
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 W are non-negative, we conclude that the characteristic equation of (4.5) has two real zeros z2<z1:
z1=w11+w22+√(w11−w22)2+4w12w212,z2=w11+w22−√(w11−w22)2+4w12w212. |
Hence, if
w11+w22+√(w11−w22)2+4w12w21<0, | (4.6) |
the trivial solution of the system (4.4) is stable, and so is E0 for (3.6); and when
w11+w22+√(w11−w22)2+4w12w21>0, | (4.7) |
the trivial solution of the system (4.4) and E0 are both unstable and so is E0 for (3.6).
By estimating the trace and determinant of 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(1−b11(τ))<1, | (4.8) |
λ2S(0)2d2+DI12+γ2(1−b22(τ))<1, | (4.9) |
(−d1−DI21+λ1S(0)1−γ1+b11(τ)γ1)(−d2−DI12+λ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 E0=(S(0)1,S(0)2,0,0) of the system (3.6) is locally asymptotically stable.
The next theorem shows that E0=(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 E0=(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 E0 is locally asymptotically stable if (4.8)–(4.10) satisfied. It merely remains to prove that E0 is globally attractive in the case (4.8)–(4.10) held, that is, for any non-negative solutions (S1(t),S2(t),I1(t),I2(t)) of (3.6), we will prove that limt→+∞(S1(t),S1(t),I1(t),I2(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=K1−d1S1(t)+DS12S2(t)−DS21S1(t)−λ1S1(t)I1(t)≤K1−d1S1(t)+DS12S2(t)−DS21S1(t),dS2(t)dt=K2−d2S2(t)+DS21S1(t)−DS12S2(t)−λ2S2(t)I2(t)≤K2−d2S2(t)+DS21S1(t)−DS12S2(t). | (4.11) |
This suggests the following comparison system for the S-equations of (3.6)
{du1(t)dt=K1−d1u1(t)+DS12u2(t)−DS21u1(t),du2(t)dt=K2−d2u2(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 Δ1(z)=0 where Δ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
S∞i≜lim supt→+∞Si(t)≤limt→+∞ui(t)=S(0)i,i=1,2. | (4.13) |
Thus, for any constant ϵ>0, there is a large enough T such that Si(t)≤S(0)i+ϵ, for all t≥T.
Now, for t≥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→∞. 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→∞. 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 E0 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 [DSij], [DIij] and [DRij] are irreducible, we have shown DFE E0 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([−τ,0],R2) denote the set of all continuous functions from [−τ,0] to R2. As is customary, C+:=C([−τ,0],R2+) denotes the subset of C consisting of all non-negative functions. By Theorem 3.1 and Remark 3.1, for any (ϕ,ψ)∈C+×C+, with ψ(0)>0 there is a unique solution to (3.6), denoted by
(S(t,ϕ,ψ),I(t,ϕ,ψ)):=(S1(t,ϕ,ψ),S2(t,ϕ,ψ),I1(t,ϕ,ψ),I2(t,ϕ,ψ)), |
whose components are all positive and bounded for t>0.
Theorem 5.1. Assume that all three travel rate matrices [DSij], [DIij] and [DRij] are irreducible. Suppose (4.7) hold, then there is an ε>0 such that for every (ϕ,ψ)∈C+×C+ with ψ(0)=(ψ1(0),ψ2(0))>0, meaning that ψi(0)≥0 for i=1,2 and ψ1(0)+ψ2(0)>0, then the solution (S(t),I(t))=(S(t,ϕ,ψ),I(t,ϕ,ψ)) of (3.6) satisfies
lim inft→∞Ii(t,ϕ,ψ)≥ε,i=1,2. |
Moreover, the model (3.6) admits at least one (componentwise) positive equilibrium.
Proof. Define X:={(ϕ,ψ)∈C+×C+}, X0:={(ϕ,ψ)∈X,ψi(0)>0,i=1,2} and ∂X0:=X∖X0. It then suffices to show that (3.6) is uniformly persistent with respect to (X0,∂X0).
Let Φ(t):X→X be the solution semiflow of (3.6)-(3.7), that is, Φ(t)(ϕ,ψ)=(St(ϕ,ψ),It(ϕ,ψ)) holds for t≥0, with S(t)=φ(t), I(t)=ψ(t), t∈[−τ,0]. By the fundamental theory for functional differential equations with bounded delays established in [18], the solution semin-flow Φ(t) is actually compact for t≥τ (consequence of the Arzela-Ascoli Theorem.) By Theorem 3.1 and Remark 3.1, X and X0 are positively invariant for Φ(t). Clearly, ∂X0={(ϕ,ψ)∈X,ψi(0)=0,foratleastonei∈{1,2,}} and it is relatively closed in X. Furthermore, system (3.6) is point dissipative in R2+ since nonnegative solutions of (3.6) are ultimately bounded (see Theorem 3.1).
Define Ω∂={(ϕ,ψ)∈X:(St(ϕ,ψ),It(ϕ,ψ))∈∂X0,∀t≥0}. We now show that
Ω∂={(ϕ,ψ)∈∂X0:I(t,ϕ,ψ)=0,∀t≥0}. | (5.1) |
Assume (ϕ,ψ)∈Ω∂. It suffices to show that I(t,ϕ,ψ)=0, ∀t≥0. For the sake of contradiction, assume that there is a t0≥0 such that I1(t0)>0. Then by the definition of Ω∂, we can derive that I2(t0,ϕ,ψ)=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 ϵ0>0 such that I2(t,ϕ,ψ)>0 and t0<t<t0+ϵ0. Clearly, we can restrict ϵ0>0 small enough such that I1(t,ϕ,ψ)>0 for t0<t<t0+ϵ0. This means that (St(ϕ,ψ),It(ϕ,ψ))∉∂X0 for t0<t<t0+ϵ0, which contradicts the assumption that (ϕ,ψ)∈Ω∂. This proves (5.1).
Choose ξ>0 small enough such that
wξ11+wξ22+√(wξ11−wξ22)2+4w12w21>0, | (5.3) |
where w12 and w21 are as in Section 4 and
wξ11=−d1−DI21+λ1(S(0)1−ξ)−γ1+b11(τ)γ1,wξ22=−d2−DI12+λ2(S(0)2−ξ)−γ2+b22(τ)γ2. |
Let us consider the following linear system
{dS1(t)dt=K1−d1S1(t)+DS12S2(t)−DS21S1(t)−ελ1S1(t)=K1−(d1+DS21+ελ1)S1(t)+DS12S2(t),dS2(t)dt=K2−d2S2(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 ε>0 small enough such that (5.4), just as (4.12), has a unique positive equilibrium (S(0)1(ε),S(0)2(ε)) which is globally asymptotically stable. By the implicit function theorem, it follows that (S(0)1(ε),S(0)2(ε)) is continuous in ε. Thus, we can further restrict ε small enough such that (S(0)1(ε),S(0)2(ε))>(S(0)1−ξ,S(0)2−ξ).
Next for the solution (S(t,ϕ,ψ),I(t,ϕ,ψ)) of (3.6) through (ϕ,ψ), we claim that
lim supt→∞max{I1(t,ϕ,ψ),I2(t,ϕ,ψ)}>ε,forall(ϕ,ψ)∈X0. | (5.5) |
Otherwise, there is a T1>0 such that 0<Ii(t,ϕ,ψ)≤ε, i=1,2, for all t≥T1. Then for t≥T1, we have
{dS1(t)dt≥K1−d1S1(t)+DS12S2(t)−DS21S1(t)−λ1S1(t)ε=K1−(d1+DS21+λ1ε)S1(t)+DS12S2(t),dS2(t)dt≥K2−d2S2(t)+DS21S1(t)−DS12S2(t)−λ2S2(t)ε=K2−(d2+DS12+λ2ε)S2(t)+DS21S1(t). | (5.6) |
Since the equilibrium (S(0)1(ε),S(0)2(ε)) of (5.4) is globally asymptotically stable and (S(0)1(ε),S(0)2(ε))>(S(0)1−ξ,S(0)2−ξ), there is a T2 such that (S1(t),S2(t))>(S(0)1−ξ,S(0)2−ξ) for t≥T1+T2. Consequently, for t≥T1+T2,
{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.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∈{1,2} such that Ii(t)→∞ as t→∞, a contradiction to the boundedness of solutions. Therefore (5.5) holds.
Note that (S(0)1,S(0)2) is globally asymptotically stable in R2+∖{0} for system (4.12). By the aforementioned claim, it then follows that E0 is an isolated invariant set in X, and Ws(E0)∩X0=∅. Clearly, every orbit in Ω∂ converges to E0, and E0 is the only invariant set in Ω∂. By Theorem 4.6 in [25], we conclude that system (3.6) is indeed uniformly persistent with respect to (X0,∂X0). Moreover, by Theorem 2.4 in [27], system (3.6) has an equilibrium (S∗1,S∗2,I∗1,I∗2)∈X0, implying that (S∗1,S∗2)∈R2+ and (I∗1,I∗2)∈int(R2+). We further claim that (S∗1,S∗2)∈R2+∖{0}. Suppose that (S∗1,S∗2)=0. By the I -equations in (3.6), we then obtain
0=−d1I∗1+DI12I∗2−DI21I∗1−γ1I∗1+b11(τ)γ1I∗1+b12(τ)γ2I∗2,0=−d2I∗2+DI21I∗1−DI12I∗2−γ2I∗2+b21(τ)γ1I∗1+b22(τ)γ2I∗2. |
and hence
0=[−d1+(b11(τ)+b21(τ)−1)γ1]I∗1+[−d2+(b22(τ)+b12(τ)−1)γ2]I∗2. |
By Lemma 2.1, we know that ∑2i=1bij(τ)≤e−d_τ≤1,forj=1,2, therefore, I∗1=I∗2=0, a contradiction. By the S-equation in (3.6) and the irreducibility of the cooperative matrix [DSij], it follows that S∗=S(t,S∗,I∗)∈int(R2+) with S∗:=(S∗1,S∗2) and I∗:=(I∗1,I∗2), for ∀t>0. Then (S∗,I∗) 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(τ)=e−d1τ:=ϵ1,b22(τ)=e−d2τ:=ϵ2. | (6.1) |
Thus, (3.6) is decoupled to
{dS1(t)dt=K1−d1S1(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=K2−d2S2(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(0)i0≜Kidiλidi+γi(1−ϵi),i=1,2, |
as summarized below.
Theorem 6.1. If R(0)i0<1, then the disease dies out in Patch i (i = 1, 2) in the sense that the disease-free equilibrium (Kidi,0) is globally asymptotically stable; if R(0)i0>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,I∗i)=(di+γi(1−ϵi)λi,Kiλi−di[di+γi(1−ϵi)]λi[di+γi(1−ϵi)]), |
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 DS12 and DS21are positive, but DI12=DI21=DR12=DR21=0. Accordingly, one can compute to obtain the following:
B=(−d100−d2),and[bij(τ)]=e(Bτ)=(ϵ100ϵ2), |
where ϵi, i=1,2 are defined in (6.1). In such a case, (3.6) reduces to
{dS1(t)dt=K1−d1S1(t)+DS12S2(t)−DS21S1(t)−λ1S1(t)I1(t),dS2(t)dt=K2−d2S2(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 E0 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 [DIij] and [DRij] does not hold. Linearizing (6.2) at E0 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 Δ1(z) is given by (4.3) and
Δ3(z,τ)=z+d1+γ1−λ1S(0)1−ϵ1γ1e−zτ, |
Δ4(z,τ)=z+d2+γ2−λ2S(0)2−ϵ2γ2e−zτ. |
It is obvious that all roots of Δ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 Δi(z,τ)=0 have negative real parts if and only if
Ri0=λiS(0)idi+γi(1−ϵi)<1. |
Therefore, the DFE E0 is asymptotically stable if max{R10,R20}<1 and it is unstable if max{R10,R20}>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 E1=(S(1)1,S(1)2,I(1)1,0) or E2=(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 E1, we need to solve the algebraical equations
{K1−d1S(1)1+DS12S(1)2−DS21S(1)1−λ1S(1)1I(1)1=0,K2−d2S(1)2+DS21S(1)1−DS12S(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)(R10−1). | (6.4) |
Thus, E1 exists (I(1)1>0) if and only if
R10=λ1S(0)1d1+γ1(1−ϵ1)>1. |
Similarly, for E2 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)(R20−1). |
Hence, E2 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∗1,S∗2,I∗1,I∗2) with all components positive, which can be determined from the following equations,
{K1−d1S∗1+DS12S∗2−DS21S∗1−λ1S∗1I∗1=0,K2−d2S∗2+DS21S∗1−DS12S∗2−λ2S∗2I∗2=0,−d1I∗1+λ1S∗1I∗1−γ1I∗1+ϵ1γ1I∗1=0,−d2I∗2+λ2S∗2I∗2−γ2I∗2+ϵ2γ2I∗2=0. |
Solving these equations for positive components leads to
S∗1=d1+γ1(1−ϵ1)λ1,S∗2=d2+γ2(1−ϵ2)λ2,I∗1=1d1+γ1(1−ϵ1)[K1−(d1+DS21)(d1+γ1(1−ϵ1))λ1+DS12d2+γ2(1−ϵ2)λ2],I∗2=1d2+γ2(1−ϵ2)[K2−(d2+DS12)(d2+γ2(1−ϵ2))λ2+DS21d1+γ1(1−ϵ1)λ1]. |
Define
\begin{equation} \nonumber \hat{R}_{10} = \frac{\lambda_1S_1^{(2)}}{d_1+\gamma_1(1-\epsilon_1)}, \; \; \; \; \; \; \hat{R}_{20} = \frac{\lambda_2S_2^{(1)}}{d_2+\gamma_2(1-\epsilon_2)}. \end{equation} |
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:
\begin{equation} \nonumber I^{\ast}_1 = \frac{d_1+D^S_{21}}{\lambda_1}\left(\hat{R}_{10}-1\right), \end{equation} |
\begin{equation} \nonumber I^{\ast}_2 = \frac{d_2+D^S_{12}}{\lambda_2}\left(\hat{R}_{20}-1\right). \end{equation} |
Thus, the interior equilibrium E^* exists if and only if
\begin{equation} \nonumber \hat{R}_{10} \gt 1\; \; \; \; and\; \; \; \; \hat{R}_{20} \gt 1. \end{equation} |
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)
\begin{equation} \begin{array}{l} 0\leq S_{1\infty}\leq S^{\infty}_1\leq S^{(0)}_{1}, \\ 0\leq S_{2\infty}\leq S^{\infty}_2\leq S^{(0)}_{2}. \end{array} \end{equation} | (6.5) |
Also, by Theorem 3.1, we know that
\begin{equation} \begin{array}{l} 0\leq I_{1\infty}\leq I^{\infty}_1 \lt \infty, \\ 0\leq I_{2\infty}\leq I^{\infty}_2 \lt \infty. \end{array} \end{equation} | (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
\begin{equation} [d_1+\gamma_1(1-\epsilon_1)]I_{1}^{\infty} \leq \lambda_1I^{\infty}_1S^{\infty}_1 \lt \lambda_1I^{\infty}_1S^{(0)}_1. \end{equation} | (6.7) |
In a similar way, we can establish
\begin{equation} [d_2+\gamma_2(1-\epsilon_2)]I_{2}^{\infty} \leq \lambda_2I^{\infty}_2S^{\infty}_2 \lt \lambda_2I^{\infty}_2S^{(0)}_2. \end{equation} | (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
\begin{equation} \nonumber b_{11}(\tau) = e^{-d_1\tau} = \epsilon_1, \; \; b_{22}(\tau) = e^{-(d_2+D^R_{12})\tau}, \; \; b_{12}(\tau) = e^{-d_1\tau}\left(1-e^{-D^R_{12}\tau}\right), \; \; b_{21}(\tau) = 0. \end{equation} |
Thus, (3.6) reduces to
\begin{equation} \left\{\begin{array}{l} {\frac{d S_1(t)}{d t}} = K_1-d_1S_1(t)+D^S_{12}S_2(t)-D^S_{21}S_1(t)-\lambda_1S_1(t)I_1(t), \\ {\frac{d S_2(t)}{d t}} = K_2-d_2S_2(t)+D^S_{21}S_1(t)-D^S_{12}S_2(t)-\lambda_2S_2(t)I_2(t), \\ {\frac{d I_1(t)}{d t}} = -d_1I_1(t)+\lambda_1S_1(t)I_1(t)-\gamma_1I_1(t)+b_{11}(\tau)\gamma_1I_1(t-\tau)+b_{12}(\tau)\gamma_2I_2(t-\tau), \\ {\frac{d I_2(t)}{d t}} = -d_2I_2(t)+\lambda_2S_2(t)I_2(t)-\gamma_2I_2(t)+b_{22}(\tau)\gamma_2I_2(t-\tau). \end{array}\right. \end{equation} | (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
\begin{equation} \nonumber R'_{20} = \frac{\lambda_2S_2^{(0)}}{d_2+\gamma_2(1-b_{22}(\tau))}, \; \; \; \; \; \; \widehat{R}'_{20} = \frac{\lambda_2S_2^{(1)}}{d_2+\gamma_2(1-b_{22}(\tau))}. \end{equation} |
Linearizing (6.9) at E_0 = (S^{(0)}_1, S^{(0)}_2, 0, 0) leads to
\begin{equation} \left\{\begin{array}{l} {\frac{d S_1(t)}{d t}} = -(d_1+D^S_{21})S_1(t)+D^S_{12}S_2(t)-\lambda_1S^{(0)}_1I_1(t), \\ {\frac{d S_2(t)}{d t}} = -(d_2+D^S_{12})S_2(t)+D^S_{21}S_1(t)-\lambda_2S^{(0)}_2(t)I_2(t), \\ {\frac{d I_1(t)}{d t}} = -(d_1+\gamma_1)I_1(t)+\lambda_1S^{(0)}_1I_1(t)+b_{11}(\tau)\gamma_1I_1(t-\tau)+b_{12}(\tau)\gamma_2I_2(t-\tau), \\ {\frac{d I_2(t)}{d t}} = -(d_2+\gamma_2)I_2(t)+\lambda_2S^{(0)}_2I_2(t)+b_{22}(\tau)\gamma_2I_2(t-\tau). \end{array}\right. \end{equation} | (6.10) |
The characteristic equation of (6.10)
\begin{equation} \Delta_1(z)\Delta_3(z, \tau)\widehat{\Delta}_4(z, \tau) = 0, \end{equation} | (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
\begin{equation} \left\{\begin{array}{l} {\frac{d S_1(t)}{d t}} = -(d_1+D^S_{21}+\lambda_1I^{(1)}_1)S_1(t)+D^S_{12}S_2(t)-\lambda_1S^{(1)}_1I_1(t), \\ {\frac{d S_2(t)}{d t}} = -(d_2+D^S_{12})S_2(t)+D^S_{21}S_1(t)-\lambda_2S^{(1)}_2I_2(t), \\ {\frac{d I_1(t)}{d t}} = -d_1I_1(t)+\lambda_1S^{(1)}_1I_1(t)+\lambda_1 I_1^{(1)}S_1(t)-\gamma_1I_1(t)+b_{11}(\tau)\gamma_1I_1(t-\tau)+b_{12}(\tau)\gamma_2I_2(t-\tau), \\ {\frac{d I_2(t)}{d t}} = -d_2I_2(t)+\lambda_2S^{(1)}_2I_2(t)-\gamma_2I_2(t)+b_{22}(\tau)\gamma_2I_2(t-\tau). \end{array}\right. \end{equation} | (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
\begin{equation} \nonumber b_{11}(\tau) = e^{-d_1\tau} = \epsilon_1, \; \; b_{22}(\tau) = e^{-d_2\tau} = \epsilon_2, \; \; b_{12}(\tau) = b_{21}(\tau) = 0. \end{equation} |
Thus, (3.6) reduces to
\begin{equation} \left\{\begin{array}{l} {\frac{d S_1(t)}{d t}} = K_1-d_1S_1(t)+D^S_{12}S_2(t)-D^S_{21}S_1(t)-\lambda_1S_1(t)I_1(t), \\ {\frac{d S_2(t)}{d t}} = K_2-d_2S_2(t)+D^S_{21}S_1(t)-D^S_{12}S_2(t)-\lambda_2S_2(t)I_2(t), \\ {\frac{d I_1(t)}{d t}} = -d_1I_1(t)+\lambda_1S_1(t)I_1(t)-\gamma_1I_1(t)+D^I_{12}I_2(t)+\epsilon_1\gamma_1I_1(t-\tau), \\ {\frac{d I_2(t)}{d t}} = -d_2I_2(t)+\lambda_2S_2(t)I_2(t)-\gamma_2I_2(t)-D^I_{12}I_2(t)+\epsilon_2\gamma_2I_2(t-\tau). \end{array}\right. \end{equation} | (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):
\begin{equation} \nonumber {R''}_{20} = \frac{\lambda_2S_2^{(0)}}{d_2+\gamma_2(1-\epsilon_2)+D^I_{12}}, \; \; \; \; \; \; \widehat{R}''_{20} = \frac{\lambda_2S_2^{(1)}}{d_2+\gamma_2(1-\epsilon_2)+D^I_{12}}. \end{equation} |
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
\begin{equation} \begin{split} R_{10}& = \frac{\lambda_1}{d_1+\gamma_1(1-b_{11}(\tau))}\cdot \frac{K_1}{d_1}\cdot \frac{d_2+D^S_{12}+\frac{K_2}{K_1}D^S_{12}}{d_2+D^S_{12}+\frac{d_2}{d_1}D^S_{21}}\\ & = R^{(0)}_{10}\cdot \frac{d_2+D^S_{12}+\frac{K_2}{K_1}D^S_{12}}{d_2+D^S_{12}+\frac{d_2}{d_1}D^S_{21}}, \end{split} \end{equation} | (7.1) |
and
\begin{equation} \begin{split} R_{20}& = \frac{\lambda_2}{d_1+\gamma_1(1-b_{22}(\tau))}\cdot \frac{K_2}{d_2}\cdot \frac{d_1+D^S_{21}+\frac{K_1}{K_2}D^S_{21}}{d_1+D^S_{21}+\frac{d_1}{d_2}D^S_{12}}\\ & = R^{(0)}_{20}\cdot \frac{d_1+D^S_{21}+\frac{K_1}{K_2}D^S_{21}}{d_1+D^S_{21}+\frac{d_1}{d_2}D^S_{12}}. \end{split} \end{equation} | (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
\begin{equation} D^S_{12} \gt \frac{d_2(1-R^{(0)}_{10})+\frac{d_2}{d_1}\cdot D^S_{21}}{R^{(0)}_{10}\cdot(1+\frac{K_2}{K_1})-1} \quad \mbox{with} \, \, 1 \gt R^{(0)}_{10} \gt \frac{K_1}{K_1+K_2}; \end{equation} | (7.3) |
or
\begin{equation} D^S_{21} \lt \frac{d_1}{d_2}\left[(R^{(0)}_{10}-1)(d_2+D^S_{12})+R^{(0)}_{10}D^S_{12}\frac{K_2}{K_1}\right] \quad \mbox{with} \; \; 1 \gt R^{(0)}_{10} \gt \frac{d_2+D^S_{12}}{d_2+D^S_{12}+\frac{K_2}{K_1}D^S_{12}}, \end{equation} | (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
\begin{equation} \nonumber D^S_{12} \lt \frac{d_2(1-R^{(0)}_{10})+\frac{d_2}{d_1}\cdot D^S_{21}}{R^{(0)}_{10}\cdot(1+\frac{K_2}{K_1})-1} \quad \mbox{with} \; \; 1 \lt R^{(0)}_{10} \lt 1+\frac{D^S_{21}}{d_1}; \end{equation} |
or
\begin{equation} \nonumber \begin{array}{l} D^S_{21} \gt \frac{d_1}{d_2}\left[(R^{(0)}_{10}-1)(d_2+D^S_{12})+R^{(0)}_{10}D^S_{12}\frac{K_2}{K_1}\right], \end{array} \end{equation} |
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
\begin{equation} \nonumber \begin{split} &\frac{d_2(1-R^{(0)}_{10})+\frac{d_2}{d_1}\cdot D^S_{21}}{R^{(0)}_{10}\cdot(1+\frac{K_2}{K_1})-1} \lt D^S_{12} \lt \frac{d_2}{d_1}\left[(R^{(0)}_{20}-1)(d_1+D^S_{21})+R^{(0)}_{20}D^S_{21}\frac{K_2}{K_1}\right] \\ & \mbox{with} \; \; 1 \gt R^{(0)}_{10} \gt \frac{K_1}{K_1+K_2}; \end{split} \end{equation} |
or
\begin{equation} \nonumber \begin{split} &\frac{d_1(1-R^{(0)}_{20})+\frac{d_1}{d_2}\cdot D^S_{12}}{R^{(0)}_{20}\cdot(1+\frac{K_1}{K_2})-1} \lt D^S_{21} \lt \frac{d_1}{d_2}\left[(R^{(0)}_{10}-1)(d_2+D^S_{12})+R^{(0)}_{10}D^S_{12}\frac{K_2}{K_1}\right]\\ &\mbox{with} \; \; 1 \gt R^{(0)}_{10} \gt \frac{d_2+D^S_{12}}{d_2+D^S_{12}+\frac{K_2}{K_1}D^S_{12}}\; \; \mbox{and}\; \; 1 \lt R^{(0)}_{20} \lt 1+\frac{D^S_{12}}{d_2}, \end{split} \end{equation} |
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
\begin{equation} \nonumber \begin{split} &\frac{d_2}{d_1}\left[(R^{(0)}_{20}-1)(d_1+D^S_{21})+R^{(0)}_{20}D^S_{21}\frac{K_2}{K_1}\right] \lt D^S_{12} \lt \frac{d_2(1-R^{(0)}_{10})+\frac{d_2}{d_1}\cdot D^S_{21}}{R^{(0)}_{10}\cdot(1+\frac{K_2}{K_1})-1} \\ & \mbox{with} \; \; 1 \gt R^{(0)}_{10} \gt \frac{K_1}{K_1+K_2}; \end{split} \end{equation} |
or
\begin{equation} \nonumber \begin{split} &\frac{d_1}{d_2}\left[(R^{(0)}_{10}-1)(d_2+D^S_{12})+R^{(0)}_{10}D^S_{12}\frac{K_2}{K_1}\right] \lt D^S_{21} \lt \frac{d_1(1-R^{(0)}_{20})+\frac{d_1}{d_2}\cdot D^S_{12}}{R^{(0)}_{20}\cdot(1+\frac{K_1}{K_2})-1} \\ &\mbox{with} \; \; 1 \gt R^{(0)}_{10} \gt \frac{d_2+D^S_{12}}{d_2+D^S_{12}+\frac{K_2}{K_1}D^S_{12}}\; \; \mbox{and}\; \; 1 \lt R^{(0)}_{20} \lt 1+\frac{D^S_{12}}{d_2}, \end{split} \end{equation} |
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] |
Pandey CL (2021) Managing urban water security: challenges and prospects in Nepal. Environ Dev Sustain 23: 241–257. https://doi.org/10.1007/s10668-019-00577-0 doi: 10.1007/s10668-019-00577-0
![]() |
[2] | IUCN (World Conservation Union) (2000). Vision for Water and Nature. A World Strategy for conservation and sustainable management of water resources in the 21st Century. IUCN, Gland, Switzerland and Cambridge, UK, p 52. |
[3] |
Perez-Blanco CD, Gil-Garcia L, Saiz-Santiago (2021) An actionable hydroeconomic descision support system for the assessment of water reallocations in irrigated agriculture. Astudy of minimum environmental flows in the Douro River Basin, Spain. J Environ Managet 298: 113432. https://doi.org/10.1016/j.jenvman.2021.113432 doi: 10.1016/j.jenvman.2021.113432
![]() |
[4] |
Senent-Aparicio J, George C, Srinivasan R (2021) Introducing a new post-processing tool for the SWAT+model to evaluate environmental flows. Environ Modell Softw 136: 104944. https://doi.org/10.1016/j.envsoft.2020.104944 doi: 10.1016/j.envsoft.2020.104944
![]() |
[5] |
Zolfagharpour F, Saghafian B, Delavar M (2021) Adapting reservoir operation rules to hydrological drought state and environmental flow requirements. J Hydrol 660: 126581. https://doi.org/10.1016/j.jhydrol.2021.126581 doi: 10.1016/j.jhydrol.2021.126581
![]() |
[6] |
Sedighkia M, Datta B (2023) Analyzing environmental flow supply in the semi-arid area through integrating drought analysis and optimal operation of reservoir. J Arid Land 15: 1439–1454. https://doi.org/10.1007/s40333-023-0035-2 doi: 10.1007/s40333-023-0035-2
![]() |
[7] |
Aazami J, Motevalli A, Savabieasfahani M (2022) Correction to: Evaluation of three environmental flow techniques in Shoor wetland of Golpayegan, Iran. Int J Environ Sci Technol 19: 7899. https://doi.org/10.1007/s13762-022-04289-3 doi: 10.1007/s13762-022-04289-3
![]() |
[8] | Wang H, Cong P, Zhu Z, et al. (2022) Analysis of environmental dispersion in wetland flows with floating vegetation islands. J Hydrol 606: 127359. https://doi.org/10.1016/j.jhydrol.2021.127359 |
[9] |
Arthington ÁH, Naiman RJ, Mcclain ME, et al. (2010). Preserving the biodiversity and ecological services of rivers: new challenges and research opportunities. Freshwater Biol 55: 1–16. https://doi.org/10.1111/j.1365-2427.2009.02340.x doi: 10.1111/j.1365-2427.2009.02340.x
![]() |
[10] |
Chilton D, Hamilton DP, Nagelkerken I, et al. (2021). Environmental flow requirements of estuaries: providing resilience to current and future climate and direct anthropogenic changes. Front Environ Sci 9: 764218. https://doi.org/10.3389/fenvs.2021.764218 doi: 10.3389/fenvs.2021.764218
![]() |
[11] |
Espinoza T, Burke CL, Carpenter-Bundhoo L, et al. (2021) Quantifying movement of multiple threatened species to inform adaptive management of environmental flows. J Environ Managet 295: 113067. https://doi.org/10.1016/j.jenvman.2021.113067 doi: 10.1016/j.jenvman.2021.113067
![]() |
[12] |
Zhang Y, Yu H, Yu H, et al. (2020) Optimization of environmental variables in habitat suitability modeling for mantis shrimp Oratosquilla oratoria in the Haizhou Bay and adjacent waters. Acta Oceanol Sin 39: 36–47. https://doi.org/10.1007/s13131-020-1546-8 doi: 10.1007/s13131-020-1546-8
![]() |
[13] |
Vivian LM, Godfree RC, Colloff MJ, et al. (2014) Wetland plant growth under contrasting water regimes associated with river regulation and drought: implications for environmental water management. Plant Ecol 215: 997–1011. https://doi.org/10.1007/s11258-014-0357-4 doi: 10.1007/s11258-014-0357-4
![]() |
[14] |
Gain AK, Apel H, Renaud FG, et al. (2013) Thresholds of hydrologic flow regime of a river and investigation of climate change impact—the case of the Lower Brahmaputra river Basin. Climatic Change 120: 463–475. https://doi.org/10.1007/s10584-013-0800-x doi: 10.1007/s10584-013-0800-x
![]() |
[15] |
Belmar O, Velasco J, Martinez-Capel F (2011) Hydrological Classification of Natural Flow Regimes to Support Environmental Flow Assessments in Intensively Regulated Mediterranean Rivers, Segura River Basin (Spain). Environ Manage 47: 992–1004. https://doi.org/10.1007/s00267-011-9661-0 doi: 10.1007/s00267-011-9661-0
![]() |
[16] |
Dumitriu D (2020) Sediment flux during flood events along the Trotuș River channel: hydrogeomorphological approach. J Soils Sediments 20: 4083–4102. https://doi.org/10.1007/s11368-020-02763-4 doi: 10.1007/s11368-020-02763-4
![]() |
[17] |
Bouska K (2020) Regime change in a large-floodplain river ecosystem: patterns in body-size and functional biomass indicate a shift in fish communities. Biol Invasions 22: 3371–3389. https://doi.org/10.1007/s10530-020-02330-5 doi: 10.1007/s10530-020-02330-5
![]() |
[18] |
Poff NL, Allan JD, Bain MB, et al. (1997). The natural flow regime. Bio Science 47: 769–784. https://doi.org/10.2307/1313099 doi: 10.2307/1313099
![]() |
[19] | Nasiri Khiavi A, Mostafazadeh R, Ghanbari Talouki F (2024) Using game theory algorithm to identify critical watersheds based on environmental flow components and hydrological indicators. Environ Dev Sustain https://doi-org.access.semantak.com/10.1007/s10668-023-04390-8 |
[20] |
Sedighkia M, Abdoli A (2023) Design of optimal environmental flow regime at downstream of multireservoir systems by a coupled SWAT-reservoir operation optimization method. Environ Dev Sustain 25: 834–854. https://doi.org/10.1007/s10668-021-02081-w doi: 10.1007/s10668-021-02081-w
![]() |
[21] |
Yu L, Wu X, Wu S, et al. (2021) Multi-objective optimal operation of cascade hydropower plants considering ecological flow under different ecological conditions. J Hydrol 601: 126599. https://doi.org/10.1016/j.jhydrol.2021.126599 doi: 10.1016/j.jhydrol.2021.126599
![]() |
[22] |
Mlynski D, Operacz A, Walega A (2020) Sensitivity of methods for calculating environmental flows based on hydrological characteristics of watercourses regarding te hydropower potential of rivers. J Clean Pro 250: 119527. https://doi.org/10.1016/j.jclepro.2019.119527 doi: 10.1016/j.jclepro.2019.119527
![]() |
[23] |
Bejarano MD, Sord-Ward A, Gabriel-Martin I, et al. (2019) Tradeoff between economic and environmental costs and benefits of hydropower production at run-of-river-diversion schemes under different environmental flows scenarios. J Hydrol 572: 790–804. https://doi.org/10.1016/j.jhydrol.2019.03.048 doi: 10.1016/j.jhydrol.2019.03.048
![]() |
[24] | Bovee KD (1982) Instream flow methodology. US Fish and Wildlife Service. FWS/OBS, 82, 26 |
[25] |
Stucchi L. Bocchiola D (2023) Environmental Flow Assessment using multiple criteria: A case study in the Kumbih river, West Sumatra (Indonesia). Sci Total Environ 901: 166516. https://doi.org/10.1016/j.scitotenv.2023.166516 doi: 10.1016/j.scitotenv.2023.166516
![]() |
[26] |
Mwamila TB, Kimwaga RJ, Mtalo FW (2008) Eco-hydrology of the Pangani River downstream of Nyumba ya Mungu reservoir, Tanzania. Phys Chem Earth 33: 695–700. https://doi.org/10.1016/j.pce.2008.06.054 doi: 10.1016/j.pce.2008.06.054
![]() |
[27] |
Nagaya T, Shiraishi Y, Onitsuka K, et al (2008) Evaluation of suitable hydraulic conditions for spawning of ayu with horizontal 2D numerical simulation and PHABSIM. Ecol Modell 215: 133–143. https://doi.org/10.1016/j.ecolmodel.2008.02.043 doi: 10.1016/j.ecolmodel.2008.02.043
![]() |
[28] |
Nikghalb S, Shokoohi A, Singh VP et al (2016) Ecological Regime versus Minimum Environmental Flow:Comparison of Results for a River in a Semi Mediterranean Region. Water Resour Manage 30: 4969–4984. https://doi.org/10.1007/s11269-016-1488-2 doi: 10.1007/s11269-016-1488-2
![]() |
[29] |
Spence R, Hickley P (2000) The use of PHABSIM in the management of water resources and fisheries in England and Wales. Ecol Eng 16: 153–158. https://doi.org/10.1016/S0925-8574(00)00099-9 doi: 10.1016/S0925-8574(00)00099-9
![]() |
[30] |
Miao Y, Li J, Feng P, et al. (2020) Effects of land use changes on the ecological operation of the Panjiakou-Daheiting Reservoir system, China Ecol Eng 152: 10585. https://doi.org/10.1016/j.ecoleng.2020.105851 doi: 10.1016/j.ecoleng.2020.105851
![]() |
[31] |
Dana K, Jana C, Tomáš V (2017) Analysis of environmental flow requirements for macroinvertebrates in a creek affected by urban drainage (Prague metropolitan area, Czech Republic). Urban Ecosyst 20: 785–797. https://doi.org/10.1007/s11252-017-0649-2 doi: 10.1007/s11252-017-0649-2
![]() |
[32] |
Weng X, Jiang C, Yuan M, et al. (2021) An ecologically dispatch strategy using environmental flows for a cascade multi-sluice system: A case study of the Yongjiang River Basin, China Ecol Indic 121: 107053. https://doi.org/10.1016/j.ecolind.2020.107053 doi: 10.1016/j.ecolind.2020.107053
![]() |
[33] |
Knack IM, Huang F, Shen HT (2020) Modeling fish habitat condition in ice affected rivers. Cold Reg Sci Technol 176: 103086. https://doi.org/10.1016/j.coldregions.2020.103086 doi: 10.1016/j.coldregions.2020.103086
![]() |
[34] |
Peng L, Sun L (2016) Minimum instream flow requirement for the water-reduction section of diversion-type hydropower station: a case study of the Zagunao River, China. Environ Earth Sci 75: 1210. https://doi.org/10.1007/s12665-016-6019-1 doi: 10.1007/s12665-016-6019-1
![]() |
[35] |
Gholami V, Khalili A, Sahour H, et al. (2020) Assessment of environmental water requirement for rivers of the Miankaleh wetland drainage basin. Appl Water Sci 10: 233. https://doi.org/10.1007/s13201-020-01319-8 doi: 10.1007/s13201-020-01319-8
![]() |
[36] |
Im D, Choi SU, Choi B (2018) Physical habitat simulation for a fish community using the ANFIS method. Ecol Inform 43: 73-83. https://doi.org/10.1016/j.ecoinf.2017.09.001 doi: 10.1016/j.ecoinf.2017.09.001
![]() |
[37] |
Zamanzad-Ghavidel S, Fazeli S, Mozaffari S, et al. (2023) Estimating of aqueduct water withdrawal via a wavelet-hybrid soft-computing approach under uniform and non-uniform climatic conditions. Environ Dev Sustain 25: 5283–5314. https://doi.org/10.1007/s10668-022-02265-y doi: 10.1007/s10668-022-02265-y
![]() |
[38] | Emadi A, Sobhani R, Ahmadi H, et al. Multivariate modeling of agricultural river water abstraction via novel integrated-wavelet methods in various climatic conditions. Environ Dev Sustain 24: 4845–4871. https://doi.org/10.1007/s10668-021-01637-0 |
[39] | Amanollahi J, Ausati Sh (2020a) PM2.5 concentration forecasting using ANFIS, EEMD-GRNN, MLP, and MLR models: A case study of Tehran, Iran. Air Qual Atmos Health 13: 161–17. https://doi.org/10.1007/s11869-019-00779-5 |
[40] |
Ghasemi A, Amanollahi J (2019) Integration of ANFIS model and forward selection method for air quality forecasting. Air Qual Atmos Health 12: 59–72. https://doi.org/10.1007/s11869-018-0630-0 doi: 10.1007/s11869-018-0630-0
![]() |
[41] |
Kaboodvandpour S, Amanollahi J, Qhavami S, et al. (2015) Assessing the accuracy of multiple regressions, ANFIS, and ANN models in predicting dust storm occurrences in Sanandaj. Iran Nat Hazards 78: 879–893. https://doi.org/10.1007/s11069-015-1748-0 doi: 10.1007/s11069-015-1748-0
![]() |
[42] |
Bice CM, Gehrig SL, Zampatti BP, et al. (2014) Flow-induced alterations to fish assemblages, habitat and fish–habitat associations in a regulated lowland river. Hydrobiologia 722: 205–222. https://doi.org/10.1007/s10750-013-1701-8 doi: 10.1007/s10750-013-1701-8
![]() |
[43] |
Kim HJ, Kim JH, Ji U, et al. (2020) Effect of Probability Distribution-Based Physical Habitat Suitability Index on Environmental-Flow Estimation. KSCE. J Civ Eng 24: 2393–2402. https://doi.org/10.1007/s12205-020-1923-z doi: 10.1007/s12205-020-1923-z
![]() |
[44] | EPA, 2015. Fish Field Sampling. SESDPROC-512-R4. |
[45] | Standard Methods for the Examination of Water and Wastewater, 1999. 20th Edition, Part 10600 B, Electrofishing, P 2206. |
[46] |
Van der Lee GEM, Van der Molen DT, Van den Boogaard HFP, et al. (2006) Uncertainty analysis of a spatial habitat suitability model and implications for ecological management of water bodies. Landscape Ecol 21: 1019–1032. https://doi.org/10.1007/s10980-006-6587-7 doi: 10.1007/s10980-006-6587-7
![]() |
[47] |
Sekine M, Wang J, Yamamoto K, et al. (2020) Fish habitat evaluation based on width-to-depth ratio and eco-environmental diversity index in small rivers. Environ Sci Pollut Res 27: 34781–34795. https://doi.org/10.1007/s11356-020-08691-7 doi: 10.1007/s11356-020-08691-7
![]() |
[48] |
Wen X, Lv Y, Liu Z, et al. (2021) Operation chart optimization of multi-hydropower system incorporating the long- and short-term fish habitat requirements. J Clen Pro 121: 107053. https://doi.org/10.1016/j.jclepro.2020.125292 doi: 10.1016/j.jclepro.2020.125292
![]() |
[49] |
Jung SH, Choi S-Uk (2015) Prediction of composite suitability index for physical habitatsimulations using the ANFIS method. Appl Soft Comput 34: 502–512. https://doi.org/10.1016/j.asoc.2015.05.028 doi: 10.1016/j.asoc.2015.05.028
![]() |
[50] |
Fraternali P, Castelletti A, Soncini-Sessa R, et al. (2012) Putting humans in the loop: social computing for Water Resources Management. Environ Modell Softw 37: 68–77. https://doi.org/10.1016/j.envsoft.2012.03.002 doi: 10.1016/j.envsoft.2012.03.002
![]() |
[51] |
Rinderknecht SL, Borsuk ME, Reichert P (2012) Bridging uncertain and ambiguous knowledge with imprecise probabilities. Environ Modell Softw 36: 122–130. https://doi.org/10.1016/j.envsoft.2011.07.022 doi: 10.1016/j.envsoft.2011.07.022
![]() |
[52] |
Fukuda S (2009) Consideration of fuzziness: Is it necessary in modelling fish habitat preference of Japanese medaka (Oryzias latipes)? Ecol Model 220: 2877–2884. https://doi.org/10.1016/j.ecolmodel.2008.12.025 doi: 10.1016/j.ecolmodel.2008.12.025
![]() |
[53] |
Fukuda S, Hiramatsu B (2008) Prediction ability and sensitivity of artificial intelligence-based habitat preference models for predicting spatial distribution of Japanese medaka (Oryzias latipes). Ecol Model 215: 301–313. https://doi.org/10.1016/j.ecolmodel.2008.03.022 doi: 10.1016/j.ecolmodel.2008.03.022
![]() |
[54] |
Dawson, CW, Wilby R L (2001) Hydrological modelling using artificial neural networks. Progr Phys Geogr 25: 80–108. https://doi.org/10.1177/030913330102500104 doi: 10.1177/030913330102500104
![]() |
[55] |
Sykes AO (1993) An Introduction to Regression Analysis. Am Stat 61: 101. https://doi.org/10.1198/tas.2007.s74 doi: 10.1198/tas.2007.s74
![]() |
[56] | Kayes I, Shahriar SA, Hasan K, et al. (2019) The relationships between meteorological parameters and air pollutants in an urban environment. Global J Environ Sci Manag 5: 265–278. |
[57] |
Ceylan Ζ, Bulkan S (2018) Forecasting PM10 levels using ANN and MLR: A case study for Sakarya City. Global Nest J 20: 281–290. https://doi.org/10.30955/gnj.002522 doi: 10.30955/gnj.002522
![]() |
[58] |
Karatzas K, Katsifarakis N, Orlowski C, et al. (2018) Revisiting urban air quality forecasting: a regression approach. Vietnam J Comput Sci 5: 177–184. https://doi.org/10.1007/s40595-018-0113-0 doi: 10.1007/s40595-018-0113-0
![]() |
[59] |
Mishra D, Goyal P (2015) Development of artificial intelligence based NO2 forecasting models at Taj Mahal, Agra. Atmos Pollut Res 6: 99–106. https://doi.org/10.5094/APR.2015.012 doi: 10.5094/APR.2015.012
![]() |
[60] |
Mishra D, Goyal P (2016) Neuro-Fuzzy approach to forecasting ozone episodes over the urban area of Delhi, India. Environ Technol Inno 5: 83–94. https://doi.org/10.1016/j.eti.2016.01.001 doi: 10.1016/j.eti.2016.01.001
![]() |
[61] |
Rahimi A (2017) Short-term prediction of NO 2 and NO x concentrations using multilayer perceptron neural network: a case study of Tabriz, Iran. Ecol Process 6: 4. https://doi.org/10.1186/s13717-016-0069-x doi: 10.1186/s13717-016-0069-x
![]() |
[62] |
Abba S, Hadi SJ, Abdullahi J (2017) River water modelling prediction using multi linear regression, artificial neural network, and adaptive neuro-fuzzy inference system techniques. Procedia Comput Sci 120: 75–82. https://doi.org/10.1016/j.procs.2017.11.212 doi: 10.1016/j.procs.2017.11.212
![]() |
[63] |
Abba SI, Elkirn G (2017) Effluent prediction of chemical oxygen demand from the astewater treatment plant using artificial neural network application. Procedia Comput Sci 120: 156–163. https://doi.org/10.1016/j.procs.2017.11.223 doi: 10.1016/j.procs.2017.11.223
![]() |
[64] |
Mohammadi S, Naseri F, Abri R (2019) Simulating soil loss rate in Ekbatan Dam watershed using experimental and statistical approaches. Int J Sediment Res 34: 226–239. https://doi.org/10.1016/j.ijsrc.2018.10.013 doi: 10.1016/j.ijsrc.2018.10.013
![]() |
[65] |
Amanollahi J, Ausati Sh (2020b). Validation of linear, nonlinear, and hybrid models for predicting particulate matter concentration in Tehran. Iran. Theor Appl Climatol 140: 709–717. https://doi.org/10.1007/s00704-020-03115-5 doi: 10.1007/s00704-020-03115-5
![]() |
[66] | Jorde K, Schneider M, Peter A, et al. (2001) Fuzzy based Models for the Evaluation of Fish Habitat Quality and Instream Flow Assessment. Proceedings of the 2001 International Symposium on Environmental Hydraulics Copyright © 2001, ISEH. |
[67] |
Fukuda Sh, Hiramatsu K, Mori M (2006) Fuzzy neural network model for habitat prediction and HEP for habitat quality estimation focusing on 61 Japanese medaka (Oryzias latipes) in agricultural canals. Paddy Water Environ 4: 119–124. https://doi.org/10.1007/s10333-006-0039-5 doi: 10.1007/s10333-006-0039-5
![]() |
[68] |
Marsili-Libelli S, Giusti E, Nocita A (2012) A new instream flow assessment method based on fuzzy habitat suitability and large scale river modelling. Environ Modell Softw 41: 27–38. https://doi.org/10.1016/j.envsoft.2012.10.005 doi: 10.1016/j.envsoft.2012.10.005
![]() |
[69] |
Junga SH, Choib S-Uk (2015) Prediction of composite suitability index for physical habitatsimulations using the ANFIS method. Appl Soft Comput 34: 502–512. https://doi.org/10.1016/j.asoc.2015.05.028 doi: 10.1016/j.asoc.2015.05.028
![]() |
[70] |
Im D, Choi S-U, Choi B (2017) Physical habitat simulation for a fish community using the ANFIS method. Ecol Inform 43: 73–83. https://doi.org/10.1016/j.ecoinf.2017.09.001 doi: 10.1016/j.ecoinf.2017.09.001
![]() |
[71] | Elkiran G, Nourani V, Abba1 SI, et al. (2018) Artificial intelligence-based approaches for multi-station modelling of dissolve oxygen in river. Global J Environ Sci Manage 4: 439–450. |
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 |