
Citation: Misrak Girma, Abebayehu Assefa, Marta Molinas. Feasibility study of a solar photovoltaic water pumping system for rural Ethiopia[J]. AIMS Environmental Science, 2015, 2(3): 697-717. doi: 10.3934/environsci.2015.3.697
[1] | Andreas Widder, Christian Kuehn . Heterogeneous population dynamics and scaling laws near epidemic outbreaks. Mathematical Biosciences and Engineering, 2016, 13(5): 1093-1118. doi: 10.3934/mbe.2016032 |
[2] | Yicang Zhou, Zhien Ma . Global stability of a class of discrete age-structured SIS models with immigration. Mathematical Biosciences and Engineering, 2009, 6(2): 409-425. doi: 10.3934/mbe.2009.6.409 |
[3] | Marcin Choiński . An inherently discrete–time $ SIS $ model based on the mass action law for a heterogeneous population. Mathematical Biosciences and Engineering, 2024, 21(12): 7740-7759. doi: 10.3934/mbe.2024340 |
[4] | Peter Rashkov, Ezio Venturino, Maira Aguiar, Nico Stollenwerk, Bob W. Kooi . On the role of vector modeling in a minimalistic epidemic model. Mathematical Biosciences and Engineering, 2019, 16(5): 4314-4338. doi: 10.3934/mbe.2019215 |
[5] | Joan Ponce, Horst R. Thieme . Can infectious diseases eradicate host species? The effect of infection-age structure. Mathematical Biosciences and Engineering, 2023, 20(10): 18717-18760. doi: 10.3934/mbe.2023830 |
[6] | Xiaowei An, Xianfa Song . A spatial SIS model in heterogeneous environments with vary advective rate. Mathematical Biosciences and Engineering, 2021, 18(5): 5449-5477. doi: 10.3934/mbe.2021276 |
[7] | Marcin Choiński, Mariusz Bodzioch, Urszula Foryś . A non-standard discretized SIS model of epidemics. Mathematical Biosciences and Engineering, 2022, 19(1): 115-133. doi: 10.3934/mbe.2022006 |
[8] | Philip K. Pollett . An SIS epidemic model with individual variation. Mathematical Biosciences and Engineering, 2024, 21(4): 5446-5455. doi: 10.3934/mbe.2024240 |
[9] | Zhen Jin, Guiquan Sun, Huaiping Zhu . Epidemic models for complex networks with demographics. Mathematical Biosciences and Engineering, 2014, 11(6): 1295-1317. doi: 10.3934/mbe.2014.11.1295 |
[10] | Jinzhe Suo, Bo Li . Analysis on a diffusive SIS epidemic system with linear source and frequency-dependent incidence function in a heterogeneous environment. Mathematical Biosciences and Engineering, 2020, 17(1): 418-441. doi: 10.3934/mbe.2020023 |
There is a large number of papers with mathematical models of epidemic dynamics. These papers relate to both homogeneous and heterogeneous populations. By heterogeneous population we understand a population with at least two subpopulations that differ in the risk of getting infected. Because of convenience and lack of appropriate data, only two subpopulations are often distinguished. In recent papers [1] and [2], one can find exemplary models of epidemic spread in such populations with mathematical analysis.
Nowadays, the risk of co-infections is increasing worldwide, including co-infections with COVID-19 [3]. This increase is caused by, for example, increased people's mobility [4], drug resistance [5], impaired immunity after suffering from long COVID-19, and reduction of exposure to microbiota caused by home isolation during the COVID-19 pandemic [6]. Also, getting an infection raises the probability of simultaneous development of another illness [7]. Co-infection also boosts healthcare costs [8] and number of deaths [9]. It is therefore essential to put an effort into reducing co-infection cases. This reduction can be achieved by applying mathematical modeling. Thanks to mathematical models of co-infection spreading and their mathematical analysis, one can predict co-infections and implement proper therapeutic approaches. This is desirable since data concerning co-infections are less accessible than for single epidemics.
Generally, models of co-infection dynamics are based on models that describe the spread of a single epidemic. Hence, they have a more complex form and are, consequently, more difficult to analyze. However, one can find literature dealing with mathematical analysis of co-infection spread. A recent paper [10] conducted a systemic review of mathematical models. A model of general co-infection for an acute and a chronic disease was presented in [11]. The authors in [12] proposed and analyzed the model of COVID-19 and tuberculosis $ (TB) $ infection. The same type of epidemics, together with measuring impact on isolation, was investigated in [13]. Papers [14] and [15] provide mathematical models for the spread of SARS-CoV-2 with hepatitis B virus and SARS-CoV-2 with human T-cell lymphotropic virus type-I, respectively. The analysis of models concerning individuals suffering from COVID-19 and kidney disease was presented in a recent work [16]. An interesting approach is presented in [17], where the authors investigated the co-infection of airborne and vector-host diseases, namely COVID-19 and dengue.
Not only COVID-19 is investigated in co-infection modeling. Since $ HIV $ infection increases the probability of developing $ TB $ [18], papers related to the modeling of $ TB $ and $ HIV $ spread contribute significantly to the literature. Authors in [19] distinguish two $ TB $ infected classes: fast and slow latent. They also considered acute and chronic $ HIV $-infected groups. A recent paper [20] focused on an analysis of a $ TB $/$ HIV $ epidemic spread in Ethiopia, with two infected groups for each disease. Naturally, other co-infections are also analyzed in the literature. Authors in [21] investigated the epidemic of $ HIV $ and hepatitis C virus, whereas paper [22] dealt with co-infection of $ TB $ and pneumonia.
In each paper described above, authors assumed that there is only one class of people that are not infected with any disease; this means they are susceptible to both infections. This leads to situations where the probability of co-infection is raised only by encountering a single infection. That supposition implies no population heterogeneity for healthy individuals. Clearly, this case is not valid. Therefore, there is a need to include heterogeneity in a susceptible class. While modeling an epidemic of a single infection for a heterogeneous population, authors assume that values of parameters reflecting a given subpopulation are the same. This assumption actually leads to the case of a homogeneous population. To adequately describe the dynamics of infection, and consequently of co-infection, in a heterogeneous population, one must consider different values of parameters for each subpopulation. Such consideration makes the mathematical analysis of the model more difficult; thus, this approach is uncommon.
Our paper aims to construct and analyze the mathematical models of co-infection in a heterogeneous population consisting of two subpopulations. We assume that parameter values reflecting a given process in each subpopulation differ. This assumption and the heterogeneity in the population make our paper novel among the literature considering the modeling of co-infections. In our model, we also introduce a case of a spread of each disease not only among one subpopulation but also between subpopulations. This introduction enables our system to be classified as a criss-cross model. We propose a system of ordinary differential equations that is a $ SIS $-type (susceptible – infected – susceptible) model. In models of such type, there is no immunity after recovery and an individual becomes susceptible again. Including heterogeneity constrains computations. For this reason, we do not incorporate a recovered class to maintain explicit results of the analysis. Our model does not relate to particular diseases. Hence, one can apply the obtained results to different infections. Such applications include co-infections that combine sole airborne illnesses, such as $ TB $, COVID-19, or influenza, sexually transmitted diseases, such as gonorrhea, $ HIV $, and hepatitis $ B $, or incorporate illnesses from both types.
This paper is a continuation of our work from [23] and [24], in which we investigated the criss-cross model of epidemic spread of a single disease for a heterogeneous population. In [23], we conducted the mathematical analysis of the model that was proposed in [25]. The aim of that analysis was to confirm, from a mathematical point of view, the medical hypothesis that stated that to control the epidemic spread in the heterogeneous population, one must consider criss-cross illness transmission. Investigating such spread in a single subpopulation does not provide complete insight into the epidemic dynamics. Because of the proposed model's unexpected properties, such as a possible unbounded population growth, we modified the system from [23] by assuming a constant inflow into each subpopulation. We analyzed that modified system in [24] and again obtained consistency with the medical hypothesis. Considering that the second disease in epidemic dynamics for a heterogeneous population is medically driven by an increasing number of co-infections worldwide and different scenarios regarding individuals' susceptibility can provide mathematical results that can help predict a co-infection course.
In [23] and [24], we focused on the stability analysis; we indicated stationary states appearing in the given systems and determined the conditions for their local stability. In this paper, we also apply this approach. The model presented herein relies on the system from [24].
The paper is organized as follows: In Section 2, we describe and introduce our model. Then we find stationary states of the system and describe conditions for their existence. The basic reproduction number of the model is computed in Section 4. The following section deals with local stability of the found stationary state. We summarize our results in Section 7.
Let us first describe the assumptions concerning the proposed model. In a population, we indicate two subpopulations, a low-risk ($ LS $) and a high-risk ($ HS $) subpopulation, relating to the risk of getting infected. $ LS $ and $ HS $ have, respectively, lower and higher susceptibility to each disease. For every variable and parameter, we assign a subscript $ i $ equal to 1 and 2 for $ LS $ and $ HS $ respectively. If $ i $ has no assigned value, then $ i \in \{1, 2 \} $. By $ S_1 $ and $ S_2 $, we denote a density of healthy people in $ LS $ and $ HS $, respectively. The variables $ I_i $ refer to the density of individuals from the given subpopulation that are infected by a pathogen of the disease that we call disease $ A $ ($ DA $). Similarly, we define $ J_i $ as the density of individuals suffering from disease $ B $ ($ DB $). The density of a group infected by pathogens from both diseases is denoted by $ K_i $.
Migrating and newborn individuals join each subpopulation through $ S_i $ class with a recruitment rate $ C_i $. A natural death rate for each subpopulation is equal to $ \mu_i $. For $ DA $, we introduce the transmission rates $ \beta_{11} $, $ \beta_{22} $, $ \beta_{12} $ and $ \beta_{21} $, reflecting transmission among $ LS $, among $ HS $, from $ HS $ to $ LS $, and from $ LS $ to $ HS $, respectively. These four different rates mean that $ DA $ differs in spreading and contracting a pathogen. To get a preliminary insight on co-infection dynamics for the heterogeneous population, for $ DB $ we assume that individuals differ only in contracting a pathogen. For this reason, we take only two transmission coefficients for $ DB $: $ \sigma_1 $ for $ LS $ and $ \sigma_2 $ for $ HS $. By $ \gamma_i $ and $ g_i $, we denote the recovery rate for $ DA $ and $ DB $, respectively. The disease-mortality rate for $ DA $ and $ DB $ is depicted by $ \alpha_i $ and $ a_i $.
The proposed model of co-infection is
$ ˙S1=C1−β11S1I1−β12S1I2+γ1I1−μ1S1−σ1S1(J1+J2)+g1J1, $ | (2.1a) |
$ ˙I1=β11S1I1+β12S1I2−(γ1+α1+μ1)I1−σ1I1(J1+J2)+g1K1, $ | (2.1b) |
$ ˙J1=σ1S1(J1+J2)−(g1+a1+μ1)J1−β11J1I1−β12J1I2+γ1K1, $ | (2.1c) |
$ ˙K1=σ1I1(J1+J2)+β11J1I1+β12J1I2−(g1+a1+γ1+α1+μ1)K1, $ | (2.1d) |
$ ˙S2=C2−β22S2I2−β21S2I1+γ2I2−μ2S2−σ2S2(J1+J2)+g2J2, $ | (2.1e) |
$ ˙I2=β22S2I2+β21S2I1−(γ2+α2+μ2)I2−σ2I2(J1+J2)+g2K2, $ | (2.1f) |
$ ˙J2=σ2S2(J1+J2)−(g2+a2+μ2)J2−β22J2I2−β21J2I1+γ2K2, $ | (2.1g) |
$ ˙K2=σ2I2(J1+J2)+β22J2I2+β21J2I1−(g2+a2+γ2+α2+μ2)K2. $ | (2.1h) |
Each parameter is fixed and positive. In particular, every parameter besides $ C_i $ is in the range $ (0, 1) $. If we assume that $ \sigma_i, g_i, a_i = 0 $, the above system would reduce to the system that we introduced and analyzed in [24].
Figure 1 is a schematic drawing of the proposed model.
In order to make the form of system (2.1) more transparent, we define:
$ ki:=γi+αi+μi,qi:=gi+ai+μi,ri:=gi+ai+γi+αi+μi. $ | (2.2) |
Now we indicate the basic properties of the model. The form of the right-hand side of system (2.1) implies that its solutions exist and are unique and positive for any positive initial condition. Let us introduce a variable $ N_i : = S_i + I_i + J_i + K_i $ that naturally means a density of the whole $ HS $ or $ LS $. After adding both sides of Eqs (2.1a)–(2.1d) or Eqs (2.1e)–(2.1h), we get
$ ˙Ni=˙Si+˙Ii+˙Ji+˙Ki=Ci−μiNi−αiIi−aiJi−(αi+ai)Ki. $ | (2.3) |
See that we can estimate Eq (2.3) from above by the inequality
$ ˙Ni≤Ci−μiNi, $ | (2.4) |
which solution is
$ N_i (t) \leq \left( N_i(0) - \frac{C_i}{\mu_i} \right) e^{- \mu_i t} + \frac{C_i}{\mu_i}. $ |
For $ N_i(0) > \tfrac{C_i}{\mu_i} $, we observe a decrease of the population. Clearly, we have
$ S_i(t), I_i(t), J_i(t), K_i(t) \leq N_i(t) \leq N_i(0). $ |
For $ N_i(0) < \tfrac{C_i}{\mu_i} $, we have a limited growth of the population, since
$ S_i(t), I_i(t), J_i(t), K_i(t) \leq N_i(t) \leq - \tilde{C} e^{- \mu_i t} + \frac{C_i}{\mu_i} \leq \frac{C_i}{\mu_i}, $ |
where $ \tilde{C} $ is any positive constant.
Now we estimate Eq (2.3) by
$ ˙Ni≥Ci−(μi+αi+ai)Ni. $ | (2.5) |
Combining inequalities (2.4) and (2.5), we get
$ C_i - (\mu_i + \alpha_i + a_i) N_i \leq \dot{N}_i \leq C_i - \mu_i N_i, $ |
which produces the invariant set
$ \Omega : \left \{ (S_1, I_1, J_1, K_1, S_2, I_2, J_2, K_2): S_i+I_i+J_i+K_i \in \left[ \frac{C_i}{\mu_i + \alpha_i + a_i}, \frac{C_i}{\mu_i} \right] \right \}. $ |
This set attracts all solutions of system (2.1). We therefore conclude that the variables $ S_i(t) $, $ I_i(t) $, $ J_i(t) $, and $ K_i(t) $ are defined for every $ t > 0 $.
In this section, we indicate stationary states of system (2.1) and determine conditions for their existence.
Firstly, we find stationary states that are not endemic, meaning that at least one of their coordinates equals zero.
We consider separated cases.
1) Firstly, let us assume that $ I_i = J_i = K_i = 0 $. We immediately get the form of the disease-free stationary state:
$ Edf=(ˆS1,0,0,0,ˆS2,0,0,0),whereˆS1=C1μ1,ˆS2=C2μ2. $ | (3.1) |
Clearly, this state always exists.
2) Now, consider the case $ K_i = 0 $. Without loss of generality, we take $ K_1 = 0 $. From Eq (2.1d) for any stationary, state we have $ 0 = \sigma_1 I_1 (J_1 + J_2) + \beta_{11} J_1 I_1 + \beta_{12} J_1 I_2, $ which implies
$ (I1=0∨J1+J2=0)∧(J1=0∨I1=0)∧(J1=0∨I2=0). $ | (3.2) |
Let us consider the first alternative from (3.2).
a) We first take $ I_1 = 0 $. From Eq (2.1b), we have $ \beta_{12} S_1 I_{2} = 0 $, which yields $ S_1 = 0 $ or $ I_2 = 0 $. If $ S_1 = 0 $, then Eq (2.1b) gives the contradiction $ 0 = C_1 + g_1 J_1 $. If $ I_2 = 0 $, then from Eq (2.1h) we have $ K_2 = 0 $ and our system reduces to:
$ ˙Si=0=Ci−μiSi−σiSi(J1+J2)+giJi,˙Ji=0=σiSi(J1+J2)−qiJi. $ | (3.3) |
The above system's form suggests that there is a stationary state with present $ DB $ and absent $ DA $. We will investigate the existence of such postulated state later. System (3.3) fulfills the case $ I_1 = 0 \land I_1 = 0 \land I_2 = 0 $ from (3.2).
Now, assume that $ I_1 = 0 \land J_1 = 0 \land J_1 = 0 $ holds. Then from Eq (2.1b), we get $ \beta_{11} S_1 I_2 + g_1 K_1 = 0 $. It provides $ S_1 = 0 $ or $ I_2 = 0 $. Case $ S_1 = 0 $ linked to Eq (2.1a) yields the contrary, $ 0 = C_1 $, hence we must have $ I_2 = 0 $. Then from Eqs (2.1c) and (2.1h), we get $ J_2 = 0 $ and $ K_2 = 0 $, respectively. We obtain state $ E_{df} $.
It is easy to check that the case $ I_1 = 0 \land J_1 = 0 \land I_2 = 0 $ gives $ E_{df} $ as well. The case $ I_1 = 0 \land I_1 = 0 \land J_1 = 0 $ is equivalent to $ I_1 = 0 \land J_1 = 0 \land J_1 = 0 $.
b) Now assume that $ J_1+J_2 = 0 $. This assumption obviously provides $ J_1 = 0 $ and $ J_2 = 0 $. From Eq (2.1g), we get $ \dot J_2 = \gamma_2 K_2 $, giving $ K_2 = 0 $. We obtain the system
$ ˙Si=0=Ci−βiiSiIi−βijSiIj+γiIi−μiSi,˙Ii=0=βiiSiIi+βijSiIj−kiIi, $ | (3.4) |
where $ j = 3-i $. The above system fulfills condition $ J_1+J_2 = 0 \land J_1 = 0 \land J_1 = 0 $, emerged from (3.2).
If we take $ J_1+J_2 = 0 $ and simultaneously one of the cases $ J_1 = 0 \land I_2 = 0 $, $ I_1 = 0 \land J_1 = 0 $ or $ I_1 = 0 \land I_2 = 0 $, then we obtain $ E_{df} $.
3) Now suppose that $ I_1 = 0 $. From Eq (2.1b), we get $ 0 = \beta_{12} S_1I_2 + g_1 K_1 $. From $ S_1 I_2 = 0 $, we get $ I_2 = 0 $, which gives system (3.3). Choosing $ K_1 = 0 $ leads to system (3.4).
4) Consider $ J_1 = 0 $. Then Eq (2.1c) yields $ 0 = \sigma_1 S_1 J_2 + \gamma_1 K_1 $. Both cases $ J_2 = 0 $ and $ K_1 = 0 $ provide system (3.4).
5) Assuming $ I_1 = 0 \land J_1 = 0 $ leads to state $ E_{df} $.
6) Supposing $ I_1 = 0 \land K_1 = 0 $ and $ J_1 = 0 \land K_1 = 0 $ results in systems (3.3) and (3.4), respectively.
Reasoning from points 1)–6) conducted for the variables with subscript 2 gives the same conclusions.
System (3.4) appears in our previous paper [24]. The solution of this system provides the endemic state $ E^* = (S_1^*, I_1^*, S_2^*, I_2^*) $. This state is the projection of the state
$ EA=(S∗1,I∗1,0,0,S∗2,I∗2,0,0),S∗i,I∗i>0, $ | (3.5) |
onto the non-negative subspace $ (S_1, I_1, S_2, I_2) \in \mathbb{R}^4 $. State $ E_A $ reflects the presence of $ DA $ and the absence of $ DB $. Observe that adding both sides of Eqs (2.1a)–(2.1b) and Eqs (2.1e)–(2.1f) for $ E_A $ yields $ 0 = C_i - \mu_i (I_i + S_i) - \alpha_i I_i $, providing
$ Ii=Ci−μiSiμi+αi. $ | (3.6) |
Relying on results for state $ E^* $ from [24], we formulate conditions for the existence of state $ E_A $.
Proposition 1. In System (2.1), there exists the stationary state $ E_A $ defined in (3.5) that reflects the presence of disease $ A $ and the absence of disease $ B $. This state exists if at least one of three cases holds:
1) $ \beta_{11} C_1 \geq\mu_1 k_1 $;
2) $ \beta_{22} C_2 \geq\mu_2 k_2 $;
3) $ \beta_{ii} C_i < \mu_i k_i $ and
$ (μ1k1−β11C1)(μ2k2−β22C2)≤β12β21C1C2. $ | (3.7) |
For $ E_A $, we have
$ 0<S∗i<kiβii,max(0,βiiCi−μikiμi+αi)<I∗i<βiiCiμi+αi. $ |
Now we focus on system (3.3). Equation (2.3), being the sum of equations from this system, simplifies to $ \dot{N}_i = 0 = C_i - \mu_i S_i - (\mu_i + a_i) J_i $, which gives
$ Si=Ci−(μi+ai)Jiμi. $ | (3.8) |
Computations concerning this system indicate the expected stationary state. Observe that the structure of system (3.3) is analogous to the structure of system (3.4). Based on this analogousness and Proposition 1, we formulate the following proposition concerning the other stationary state.
Proposition 2. In System (2.1), there exists the stationary state
$ EB:=(ˉS1,0,ˉJ1,0,ˉS2,0,ˉJ2,0),ˉSi,ˉJi>0, $ | (3.9) |
with present disease $ B $ and absent disease $ A $. This state exists if at least one of three cases holds:
1) $ \sigma_1 C_1 \geq\mu_1 q_1 $;
2) $ \sigma_2 C_2 \geq\mu_2 q_2 $;
3) $ \sigma_i C_i < \mu_i q_i $ and
$ σ1C1μ1q1+σ2C2μ2q2≥1. $ | (3.10) |
For $ E_B $, we have
$ 0<ˉSi<qiβii,max(0,σiCi−μiqiμi+ai)<ˉIi<σiCiμi+ai. $ |
Now we investigate the endemic stationary state, reflecting the presence of both diseases $ A $ and $ B $. Since it is difficult to obtain its explicit form, we limit our investigation to indicate some properties concerning its existence. Let us define $ Z_i : = I_i + J_I + K_i $. Obviously, variable $ Z_i $ is the density of individuals from the particular subpopulation infected by disease $ A $, $ B $, or both. Summing both sides of Eqs (2.1b)–(2.1d) gives
$ \dot{Z}_1 = \beta_{11} I_1 S_1 + \beta_{12} I_2 S_1 - g_1 J_1 - a_1 (J_1 + K_1) + \sigma_1 S_1 (J_1 + J_2) - \gamma_1 I_1 - \mu_1 Z_1 - \alpha_1 (I_1 + K_1). $ |
From the above equation for the postulated endemic state, we have
$ S1=g1J1+a1(J1+K1)+γ1I1+μ1Z1+α1(I1+K1)β11I1+β12I2+σ1(J1+J2). $ | (3.11) |
Equation (2.3) and definition of $ N_1 $ for this state yield
$ S1=1μ1(C1−μ1(I1+J1+K1)−α1I1−a1J1−(α1+a1)K1). $ | (3.12) |
Since $ S_1 > 0 $, one must fulfill
$ C_1 > \mu_1 (I_1+J_1+K_1) + \alpha_1 I_1 + a_1 J_1 + (\alpha_1 + a_1) K_1. $ |
Combining reorganized Eqs (3.12) and (3.11), we get
$ (C1−μ1(I1+J1)−α1I1−a1J1−s1K1)(β11I1+β12I2+σ1(J1+J2))=μ1(k1I1+w1J1+s1K1), $ | (3.13) |
where $ s_1 : = \alpha_1 + a_1 + \mu_1 $ and $ w_1 : = g_1 + a_1 + \mu_1 $.
From Eqs (2.1a)–(2.1d), we get the formula for coordinates $ S_1, I_1, J_1, K_1 $:
$ S1=C1+g1J1+γ1I1β11I1+β12I2+μ1+σ1(J1+J2), $ | (3.14a) |
$ I1=β12S1I2+g1K1σ1(J1+J2)+k1I1−β11S1, $ | (3.14b) |
$ J1=σ1S1J2+γ1K1q1+β11I1+β12I2−σ1S1 $ | (3.14c) |
$ K1=β11I1J1+β12J1I2+σ1I1(J1+J2)r1. $ | (3.14d) |
Positivity of coordinates $ I_1 $ and $ J_1 $ yields inequalities:
$ S_1 < \frac{ \sigma_1 (J_1 + J_2) + k_1 I_1}{\beta_{11}}, \qquad S_1 < \frac{q_1 + \beta_{11} I_1 + \beta_{12} I_2}{\sigma_1}. $ |
Combining the above dependences with condition (3.14a), we get
$ C1<min(σ1(J1+J2)+k1I1β11,q1+β11I1+β12I2σ1)(σ1(J1+J2)+μ1+β11I1+β12I2)−γ1I1−g1J1. $ |
If we substitute Eq (3.14d) into the above inequality, we get
$ C_1 > \mu_1 (I_1+J_1) + \alpha_1 I_1 + a_1 J_1 + \frac{\alpha_1 + a_1 + \mu_1}{r_1} \Bigg( \beta_{11} I_1 J_1 + \beta_{12} J_1 I_2 + \sigma_1 I_1 (J_1 + J_2) \Bigg) . $ |
Observe that the reasoning presented in this subsection can be applied to variables with subscript 2. This application provides another condition for the existence of the postulated endemic state.
In this section, we find and investigate the basic reproduction number $ \mathcal{R}_0 $ of system (2.1). According to the definition from [26], $ \mathcal{R}_0 $ refers to the number of new infections produced by a single infectious individual in a population at a disease-free stationary state. To compute $ \mathcal{R}_0 $, we will rely on the next generation method described in [26]. In this section, we provide a sketch of computations leading to the formula for $ \mathcal{R}_0 $. We consider the subsystem of system (2.1), including the equations only for infected variables:
$ [˙I1,˙J1,˙K1,˙I2,˙J2,˙K2]^T = \mathcal{F} - \mathcal{V}, $ |
where vector $ \mathcal{F} $ concerns the terms related to new infections, and vector $ \mathcal{V} $ reflects the remaining processes. These vectors read
$ \mathcal{F} = [β11S1I1+β12S1I2σ1S1(J1+J2)σ1I1(J1+J2)+β11J1I1+β12J1I2β22S2I2+β21S2I1σ2S2(J1+J2)b2I2(J1+J2)+β22J2I2+β12J2I1], \qquad \mathcal{V} = [k1I1+σ1I1(J1+J2)−g1K1q1J1+β11J1I1+β12J1I2−γ1K1r1K1k2I2+σ2I1(J1+J2)−g2K2q2J2+β22J2I2+β21J2I1−γ2K2r2K2], $ |
respectively. We construct matrices $ F $ and $ V $ that are the Jacobian matrices of $ \mathcal{F} $ and $ \mathcal{V} $ evaluated at $ E_{df} $. These matrices have the form
$ F = [β11ˆS1β12ˆS10000β21ˆS2β22ˆS2000000σ1ˆS1σ1ˆS10000σ2ˆS2σ2ˆS200000000000000], \qquad V = [k1000−g100k2000−g200q10−γ10000q20−γ20000r1000000r2], $ |
Then we compute $ F V^{-1} $ and obtain
$ FV^{-1} = [β11k1ˆS1β12k2ˆS100g1β11k1r1ˆS1g2β12k2r2ˆS1β21k1S2β22k2ˆS200g1β21k1r1ˆS2g2β22k2r2ˆS200σ1q1ˆS1σ1q2ˆS1γ1σ1q1r1ˆS1γ2σ1q2r2ˆS100σ2q1ˆS2σ2q2ˆS2γ1σ2q1r1ˆS2γ2σ2q2r2ˆS2000000000000]. $ |
The basic reproduction number is the spectral radius of the matrix $ FV^{-1} $. The eigenvalues of this read $ \lambda_{1, 2, 3} = 0, \lambda_4 = \tfrac{\sigma_1}{q_1} \widehat{S}_1 + \tfrac{\sigma_2}{q_2} \widehat{S}_2 $, and
$ \lambda_{5,6} = \frac{1}{2 k_1 k_2} \left( k_2 \beta_{11} \widehat{S}_1 + k_1 \beta_{22} \widehat{S}_2 \mp \sqrt{ ( k_2 \beta_{11} \widehat{S}_1 - k_1 \beta_{22} \widehat{S}_2 )^2 + 4 k_1 k_2 \beta_{12} \beta_{21} \widehat{S}_1 \widehat{S}_2 } \right) $ |
We finally get
$ \mathcal{R}_0 = \max (\lambda_4, \lambda_6). $ |
See that $ \lambda_4 $ consists of terms relying to only infection $ B $, whereas $ \lambda_6 $ refers only to infection $ A $.
Now we formulate the theorem relied on the dependence between $ \lambda_4 $ and $ \lambda_6 $.
Theorem 4.1. $ \mathcal{R}_0 = \lambda_4 $ if
$ ˆS1ˆS2(β12β21−β11β12)<ˆS1ˆS2(σ1u2q1+σ2u1q2)+σ1u1q1ˆS21+σ2u2q2ˆS22, $ | (4.1) |
where
$ ui:=k3−i(kiσiqi−βii), $ | (4.2) |
and
$ ˆS1(2σ1q1−k2β11)+ˆS2(2σ2q2−k1β22)>0. $ | (4.3) |
Proof. Observe that $ \lambda_4 > \lambda_6 $ if
$ √(k2β11ˆS1−k1β22ˆS2)2+4k1k2β12β21ˆS1ˆS2<2λ4k1k2−k2β11ˆS1−k1β22ˆS2. $ | (4.4) |
Under fulfillment of
$ 2λ4k1k2>k2β11ˆS1+k1β22ˆS2 $ | (4.5) |
we raise both sides of inequality (4.4) to the square and transform the result to
$ \widehat{S}_1 \widehat{S}_2 (\beta_{12} \beta_{21} - \beta_{11} \beta_{12} ) < \lambda_4 (\lambda_4 - k_2 \beta_{11} \widehat{S}_1 - k_1 \beta_{22} \widehat{S}_2). $ |
Using the definition of $ \lambda_4 $ in the above inequality and transforming the obtained expression, we get
$ ˆS1ˆS2(β12β21−β11β12)<(σ1ˆS1q1+σ2ˆS2q2)(ˆS1u1+ˆS2u2), $ | (4.6) |
where $ u_i $ is defined by Eq (4.2). We rewrite inequality (4.6) as inequality (4.1). Inequality (4.5) with the definition of $ \lambda_4 $ can be transformed into inequality (4.3).
Let us strengthen the assumptions from Theorem 4.1 so that we obtain more explicit ones. Suppose that $ u_i > 0 $, which can be written as
$ βiiqi<kiσi. $ | (4.7) |
Then condition
$ σ1u2q1+σ2u1q2>β12β21−β11β12 $ | (4.8) |
suffices fulfillment of inequality (4.1). Moreover, if
$ βiiqi<2σik3−i, $ | (4.9) |
then inequality (4.3) is always true. Combining inequalities (4.7) and (4.9) yields
$ βiiqi<max(kiσi,2σik3−i). $ | (4.10) |
We conclude that
Corollary 1. If inequality (4.10), then $ \mathcal{R}_0 = \lambda_4 $.
The above corollary confirms the obvious dependence that if the transmission of infection $ B $, represented by $ \sigma_i $, is sufficiently stronger than the transmission of infection $ A $ between different subpopulations, reflected by $ \beta_{ii} $, then infection $ B $ plays a bigger role in the whole population.
Observe that if in inequality (4.3) we replace sign $ > $ by sign $ < $, then $ \lambda_4 < \lambda_6 $. It yields $ \mathcal{R}_0 = \lambda_6 $, which is the same as $ \mathcal{R}_0 $ for the system with one infection from [24].
From a medical point of view, the desirable situation is when $ \mathcal{R}_0 < 1 $. Let us check when this case holds. We formulate the theorem
Theorem 4.2. For system (2.1) $ \mathcal{R}_0 < 1 $ if
$ σ1q1ˆS1+σ2q2ˆS2<1, $ | (4.11) |
$ ˆSi<kiβii $ | (4.12) |
and
$ β12β21ˆS1ˆS2<(β11ˆS1−k1)(β22S2−k2). $ | (4.13) |
Proof. Condition (4.11) is obvious from the definition of $ \lambda_4 $. The part of the proof related to $ \lambda_6 $ is similar to the proof of Theorem 4.1. The inequality $ \lambda_6 < 1 $ can be transformed into
$ √(k2β11ˆS1−k1β22ˆS2)2+4k1k2β12β21ˆS1ˆS2<2k1k2−(k2β11ˆS1+k1β22ˆS2). $ | (4.14) |
Under condition
$ 2>β11k1ˆS1+β22k2ˆS2 $ | (4.15) |
multiplying both sides of inequality (4.14) yields
$ \beta_{12} \beta_{21} \widehat{S}_1 \widehat{S}_2 < k_1 k_2 - k_2 \beta_{11} \widehat{S}_1 - k_1 \beta_{22} \widehat{S}_2 + \beta_{11} \beta_{12} \widehat{S}_1 \widehat{S}_2, $ |
which can be written as inequality (4.13). The right-hand side of inequality (4.13) must be positive. It is true when $ \widehat{S}_i > \tfrac{k_i}{\beta_{ii}} $ or (4.12). The first case is contrary to inequality (4.15). Hence, inequality (4.12) must hold, which is stronger than inequality (4.15).
If is easy to check that inequality (4.13) is opposite to inequality (3.7) from Preposition 1 discussing the existence of state $ E_A $. Similarly, inequality (4.11) is opposite to inequality (3.10) from Preposition 2 providing conditions for the $ E_B $ existence.
Now we investigate the local stability of states $ E_{df} $, $ E_A $, and $ E_B $. The Jacobian matrix of system (2.1) can be written as $ J = (M1M2) $, where
$ M_1 = (G1−β11S1+γ1−σ1S1+g10F1β11S1−k1−σ1(J1+J2)−σ1I1g1σ1(J1+J2)−β11J1σ1S1−q1−F1γ10σ1(J1+J2)+β11J1σ1I1+F1−r10−β21S2−σ2S200β21S2−σ2I200−β21J2−σ2S200β21J2σ2I20), $ |
and
$ M_2 = (0−β12S1−σ1S100β12S1−σ1I100−β12J1σ1S100β12J1σ1I10G2−β22S2+γ2−σ2S2+g20F2β22S2−k2−σ2(J1+J2)−σ2I2g2σ2(J1+J2)−β22J2σ2S2−q2−F2γ20σ2(J1+J2)+β22J2σ2I2+F2−r2) $ |
with
$ Fi:=Fi(I1,I2)=βiiIi+βijIj,Gi:=Gi(I1,I2,J1,J2)=−Fi−μi−σi(J1+J2),j=3−i. $ |
We start from the local stability of $ E_{df} $.
Theorem 5.1. $ E_{df} $ is locally stable if $ \mathcal{R}_0 < 1 $.
Proof. The Jacobian matrix for state $ E_{df} $ reads
$ (−μ1−β11ˆS1+γ1−σ1ˆS1+g100−β12ˆS1−σ1ˆS100β11ˆS1−k10g10β12ˆS10000σ1ˆS1−q1γ100σ1ˆS10000−r100000−β21ˆS2−σ2ˆS20−μ2−β22ˆS2+γ2−σ2ˆS2+g200β21ˆS2000β22ˆS2−k20g200−σ2ˆS2000σ2ˆS2−q2γ20000000−r2). $ |
Immediately, we get four negative eigenvalues: $ - \mu_1 $, $ -\mu_2 $, $ -r_1 $, $ -r_2 $. The remaining eigenvalues are the zeros of the characteristic polynomial of the matrix:
$ (β11ˆS1−k10β12ˆS100σ1ˆS1−q10σ1ˆS1β21ˆS20β22ˆS2−k200−σ2ˆS20σ2ˆS2−q2). $ |
This polynomial reads $ P(\lambda) = P_1(\lambda) P_2(\lambda) $, where
$ P1(λ)=λ2−(β11ˆS1−k1+β22ˆS2−k2)λ+(β11ˆS1−k1)(β22S2−k2)−β12β21ˆS1ˆS2,P2(λ)=λ2−(σ1ˆS1−q1+σ2ˆS2−q2)λ+(σ1ˆS1−q1)(σ2ˆS2−q2)−σ1σ2ˆS1ˆS2. $ |
It is easy to check that their discriminants are positive. The zeros of $ P_1 $ are negative if inequalities (4.12) and (4.13) hold. Analogously, $ P_2 $ has negative zeros if
$ σ1ˆS1−q1+σ2S2−q2<0 $ | (5.1) |
and
$ (σ1ˆS1−q1)(σ2ˆS2−q2)>σ1σ2ˆS1ˆS2. $ | (5.2) |
Inequality (5.2) can be transformed into inequality (4.11). According to Theorem 4.2, merging inequality (4.11)–(4.13) yields $ \mathcal{R}_0 < 1 $. Inequality (5.2) yields two exclusive possibilities: $ \sigma_i \widehat{S}_i > q_i $ or
$ σiˆSi<qi. $ | (5.3) |
The first case is contrary to inequality (5.1), whereas inequality (5.3), under condition (5.2), yields inequality (5.1). Hence, we replace inequality (5.1) by inequality (5.3). We rewrite inequality (5.3) as
$ \widehat{S}_i < \frac{q_i}{\sigma_i}. $ |
Observe that the above inequality is weaker than inequality (4.11); hence, it is omitted in the thesis of Theorem 5.1.
Theorem 5.1 is in line with the analogical result from [24], where we also obtained the local stability of the given disease-free stationary for the basic reproduction number smaller than one.
Now we provide a theorem guaranteeing the local stability of state $ E_A $.
Theorem 5.2. $ E_A $ is locally stable if
$ γiβii<S∗i<kiβii, $ | (5.4) |
$ (k1−β11S∗1)(k2−β22S∗2)>β12β21S∗1S∗2, $ | (5.5) |
$ βiiI∗i+βijI∗j+qi>σiSi,j=3−i, $ | (5.6) |
$ S∗1S∗2<r1r2σ1σ2, $ | (5.7) |
$ ri(βiiI∗i+βijI∗j−σiS∗i+qi)>γi(βiiI∗i+βijI∗j+σiI∗i), $ | (5.8) |
$ (rjγj(βjjI∗j+βjiI∗i−σjS∗j+qj)−βjjI∗j−βjiI∗i−σjI∗j)⋅(βiiI∗i+βijI∗j−σiS∗i+qi+ri)>σiS∗i(rjσjγjS∗j+σjI∗j),j=3−i, $ | (5.9) |
and
$ 2∏m=1(rm(βmmI∗m+βmjI∗j−σmS∗m+qm)−γm(βmmI∗m+βmjI∗j+σmI∗m))>2∏m=1(rmσmS∗m+γmσmI∗m),j=3−m. $ | (5.10) |
Proof. Let us rearrange matrix $ J(E_A) $ so that the characteristic polynomial of the new matrix remains the same. We get $ M^* = (M∗1 M∗20M∗3), $ where
$ M_1^* = \left(−β11I∗1−β12I∗2−μ1−β11S∗1+γ10−β12S∗1β11I∗1+β12I∗2β11S∗1−k10β12S∗10−β21S∗2−β22I∗2−β21I∗1−μ2−β22S∗2+γ20β21S∗2β22I∗2+β21I∗1β22S∗2−k2\right) $ |
and
$ M_3^* = \left(−β11I∗1−β12I∗2+σ1S∗1−q1γ1σ1S∗10β11I∗1+β12I∗2+σ1I∗1−r1σ1I∗10σ2S∗20−β22I∗2−β21I∗1+σ2S∗2−q2γ2σ2I∗20β22I∗2+β21I∗1+σ2I∗2−r2\right). $ |
The form of $ M_2^* \in M_4(\mathbb{R}) $ is not needed for further computations.
We start by investigating matrix $ M_1^* $. To simplify computations, we define auxiliary notations:
$ ai=βiiI∗i+βijI∗j,ui=βiiS∗i−γi,ti=ki−βiiS∗i,zi=βijS∗i,wherej=3−i. $ | (5.11) |
Thanks to this simplification, matrix $ M_1^* $ reads
$ M_1^* = \left(−a1−μ1−u10−z1a1−t10z10−z2−a2−μ2−u20z2a2−t2\right). $ |
The characteristic polynomial of $ M_1^* $ has the form
$ P_1 (\lambda) : = \lambda^4 + c_3 \lambda^3 + c_2 \lambda^2 + c_1 \lambda + c_0, $ |
where
$ c3=a1+a2+t1+t2+μ1+μ2,c2=(a1+μ1)(a2+μ2)+(t1+t2)(a1+a2+μ1+μ2)+a1u1+a2u2+t1t2−z1z2,c1=a1a2(t1+t2+u1+u2)+μ1μ2(t1+t2)+t1t2(a1+a2)+a1μ2(t1+t2+u1)+a2μ1(t1+t2+u2)+a1u1t2+a2u2t1+(μ1+μ2)(t1t2−z1z2),c0=a1a2(t1+u1)(t2+u2)+a1t2μ2(t1+u1)+a2t1μ1(t2+u2)+μ1μ2(t1t2−z1z2). $ |
Observe that if
$ ui>0,ti>0 $ | (5.12) |
and
$ t1t2−z1z2>0, $ | (5.13) |
then each coefficient of $ P_1 $ is positive. Hence, from Descartes' rule of signs, we get that $ P_4 $ has real negative roots or complex roots with negative real parts. Substituting definition (5.11) into inequalities (5.12) and (5.13) provides inequalities (5.4) and (5.5), respectively.
Let us focus now on matrix $ M_3^* $. After using notations:
$ ti=βiiI∗i+βijI∗j−σiS∗i+qi,ai=βiiI∗i+βijI∗j+σiI∗i,si=σiS∗i,yi=σiI∗i,wherej=3−i, $ | (5.14) |
we rewrite $ M_3^* $ as
$ M∗3=(−t1γ1s10a1−r1y10s20−t2γ2y20a2−r2). $ | (5.15) |
The characteristic polynomial of the matrix reads $ P_3 (\lambda) : = \lambda^4 + c_3 \lambda^3 + c_2 \lambda^2 + c_1 \lambda + c_0 $, where
$ c3=t1+t2+r1+r2>0,c2=(r1+t1)(r2+t2)−s1s2+r1t1−γ1a1+r2t2−γ2a2,c1=2∑j=1((tj+rj)(t3−jr3−j−γ3−ja3−j)−sj(r3−js3−j+γ3−jy3−j)),c0=(γ1a1−r1t1)(γ2a2−r2t2)−(r1s1+γ1y1)(r2s2+γ2y2). $ |
Let us investigate the signs of the coefficients of polynomial $ P_3 $. Observe that if inequality (5.6) holds, then $ t_1, t_2 > 0 $, which yields $ c_3 > 0 $. Conditions $ r_1 r_2 > s_1 s_2 $ and $ r_i t_i > \gamma_i a_i $, which can be written as inequalities (5.7) and (5.8), respectively, yield $ c_2 > 0 $. See that $ c_1 > 0 $ if $ (t_i + r_i) (t_j r_j - \gamma_j a_j) > s_i (r_j s_j + \gamma_j y_j) $ for $ j = 3-i $, which we transform to
$ (t_i + r_i) \left(\frac{r_j }{\gamma_j } t_j - a_j \right) > s_i \left( \frac{r_j}{\gamma_j} s_j + y_j \right). $ |
Using expressions from (5.14), we rewrite the above inequality as inequality (5.9). Condition $ c_0 > 0 $ can be written as
$ \prod\limits_{m = 1}^2 (\gamma_m a_m - r_m t_m) > \prod\limits_{m = 1}^2 (r_m s_m + \gamma_m y_m) $ |
With (5.14), the above inequality transforms into inequality (5.10).
Now let us strengthen conditions from Theorem 5.2 providing the local stability of $ E_A $ so that they have a more explicit form. Observe that the condition
$ S∗i<qiσi. $ | (5.16) |
yields fulfillment of inequality (5.6). Moreover, from definitions (2.2), clearly we have $ r_i > q_i $. Hence, inequality (5.16) implies inequality (5.7).
Again relying on (2.2), we get $ r_i > \gamma_i $. Instead of inequality (5.8), it is therefore enough to investigate the inequality
$ r_i ( q_i - \sigma_i S_i^* ) > \gamma_i \sigma_i I_i^*. $ |
If inequality (5.16) holds, then the left-hand side of the above inequality is always positive. Using Eq (3.6), we transform this inequality into
$ S∗i<riqi−Ciγiσiαi+μiriσi−γiσiμiαi+μi. $ | (5.17) |
From the definition of $ r_i $, the denominator of the right-hand side of inequality (5.17) is positive if
$ (g_i + a_i + \alpha_i + \mu_i) (\mu_i + \alpha_i) + \gamma_i \alpha_i > 0 , $ |
which is always true. The positivity of the numerator of the right-hand side of inequality (5.17) is maintained if
$ Ci<riqi(μi+αi)γiσi. $ | (5.18) |
Now observe that since $ r_i > \gamma_i $, inequality (5.9) can be, under fulfillment of inequality (5.16), strengthened to
$ (qj−σj(S∗j+I∗j))⋅(βiiI∗i+βijI∗j−σiS∗i+qi+ri)>σiS∗i(rjσjγjS∗j+σjI∗j). $ | (5.19) |
If inequality (5.16) holds, then one must fulfill
$ S∗i+I∗i<qiσi. $ | (5.20) |
so that inequality (5.19) makes sense. Using Eq (3.6), we transform inequality (5.20) into
$ S∗i<qi(μi+αi)−Ciσiσiαi, $ | (5.21) |
which is reasonable if
$ Ci<qi(μi+αi)σi. $ | (5.22) |
Observe that inequality (5.22) is stricter than inequality (5.18). Rewriting inequality (5.19), we get
$ (βiiI∗i+βijI∗j)(qj−σj(S∗j+I∗j))+(qi+ri)qj>(rjγj−1)σ1σ2S∗1S∗2+(σiqjS∗i+(qi+ri)σj(S∗j+I∗j)). $ |
Using again Eq (3.6), we transform the above inequality into
$ −βijσjμ2j(μj+αj)2(S∗j)2−βiiμiμi+αi⋅σjμjμj+αjS∗1S∗2+(rjγj−1)σ1σ2S∗1S∗2+(βiiCiμi+αi⋅βijCjμj+αj)σjμjμj+αjS∗j−βijμjμj+αj(qj−σjCjμj+αj)S∗j+αj(qi+ri)μj+αjS∗j+(qj−σjCjμj+αj)(βii(Ci−μiS∗i)μi+αi+βijCjμj+αj)+σiqjS∗i+(qi+ri)σjCjμj+αj+(qi+ri)qj>0,j=3−i. $ | (5.23) |
Now we use the dependence $ r_i > \gamma_i $ and transform inequality (5.10) into
$ \prod\limits_{m = 1}^2 \Big( r_m (q_m - \sigma_m S_m^*) - \gamma_m \sigma_m I_m^* \Big) > \prod\limits_{m = 1}^2 (r_m \sigma_m S_m^* + \gamma_m \sigma_m I_m^*), $ |
which can be simplified to
$ σ1q1(S∗1+I∗1)+σ2q2(S∗2+I∗2)<1. $ | (5.24) |
Clearly, inequality (5.24) is stricter than inequalities (5.16) and (5.21).
Using Eq (3.6), we rewrite inequality (5.24) as
$ 2∑j=1(σj(S∗jαj+Cj)qj(μj+αj))<1. $ | (5.25) |
We finally conclude that
Corollary 2. If (5.4), (5.5), (5.17), (5.22), (5.23) and (5.25), then $ E_A $ is locally stable.
Let us treat the left-hand side of inequality (5.23) as a quadratic trinomial $ P(S_j^*) $. If inequality (5.22) holds, then the condition $ C_i > \sigma_i S_i^* $ suffices for the constant of $ P(S_j^*) $ to be positive. Hence, the form of $ P(S_j^*) $ implies that inequality (5.23) is true for $ 0 < S_j^* < \mathcal{S} $, where $ \mathcal{S} $ is the positive zero of $ P(S_j^*) $.
Now we provide the theorem indicating the conditions for the local stability of state $ E_B $:
Theorem 5.3. $ E_B $ is locally stable if
$ ˉSi<kiβii, $ | (5.26) |
$ ˉS1ˉS2<r1r2β12β21, $ | (5.27) |
$ ri(σi(ˉJ1+ˉJ2)−βiiˉSi+ki)>gi(σi(ˉJ1+ˉJ2)+βiiˉJi), $ | (5.28) |
$ (rjgj(σj(ˉJ1+ˉJ2)−βjjˉSj+kj)−σj(ˉJ1+ˉJ2)−βjjˉJj)⋅(σi(ˉJ1+ˉJ2)−βiiˉSi+ki+ri)>β12β21ˉSi(rjgjˉSj+ˉJj),j=3−i, $ | (5.29) |
$ 2∏m=1(gm(σm(ˉJ1+ˉJ2)+βmmˉJm)−rm(σm(ˉJ1+ˉJ2)−βmmˉSm+km))>2∏m=1(rmβmjˉSm+gmβmjˉJm),j=3−m. $ | (5.30) |
Proof. Similarly as in the proof of Theorem 5.2, we transform matrix $ J(E_B) $ into $ M_B = (ˉM1ˉM20ˉM3), $ where
$ \bar{M}_1 = \left(−μ1−σ1(ˉJ1+ˉJ2)−σ1ˉS1+g10−σ1ˉS1σ1(ˉJ1+ˉJ2)σ1ˉS1−q10σ1ˉS10−σ2ˉS2−μ2−σ2(ˉJ1+ˉJ2)−σ2ˉS2+g20σ2ˉS2σ2(ˉJ1+ˉJ2)σ2ˉS2−q2\right), $ |
$ \bar{M}_3 = \left(β11ˉS1−k1−σ1(ˉJ1+ˉJ2)g1β12ˉS10σ1(ˉJ1+ˉJ2)+β11ˉJ1−r1β12ˉJ10β21ˉS20β22ˉS2−k2+σ2(ˉJ1+ˉJ2)g2β21ˉJ20σ2(ˉJ1+ˉJ2)+β22ˉJ2−r2\right). $ |
We start from matrix $ \bar{M}_1 $. After using notations:
$ ai=σi(ˉJ1+ˉJ2),si=σiˉSi,ti=qi−σiˉSi, $ | (5.31) |
we transform $ \bar{M}_1 $ to
$ \bar{M}_1 = \left(−a1−μ1g1−s10−s1a1−t10s10−s2−a2−μ2g2−s20s2a2−t2\right). $ |
The characteristic polynomial of matrix $ \bar{M}_1 $ reads $ P_1 (\lambda) = \lambda^4 + c_3 \lambda^3 + c_2 \lambda^2 + c_1 \lambda + c_0 $, where
$ c3=a1+a2+t1+t2+γ1+γ2,c2=t1t2−s1s2+a1(s1−g1)+a2(s2−g2)+(t1+t2)(a1+μ1+a2+μ2)+(a1+μ1)(a2+μ2),c1=a1a2(s1+s2−g1−g2)+(μ1+μ2)(t1t2−s1s2)+a1(t2+μ2)(s1−g1)+a2(t1+μ1)(s2−g2),+a1t1(t2+μ2)+a2t2(t1+μ1)+(a1a2+μ1μ2)(t1+t2)+a1t2μ2+a2t1μ1,c0=(t1t2−s1s2)(μ1μ2+a1a2)+a1t2(a2+μ2)(s2−g2)+t1t2(a1μ2+a2μ1)+a2t1(a1+μ1)(s1−g1)+a1a2(s1−g1)(s2−g2)+a1a2s1s2. $ |
Observe that if
$ s1>g1,s2>g2, $ | (5.32) |
and
$ t1t2>s1s2 $ | (5.33) |
then $ c_2, c_1, c_0 > 0 $. Hence, from Descartes' rule of signs, we get that $ P_1 $ has real negative roots or complex roots with negative real parts.
Using definitions (5.31), we rewrite inequality (5.32) as inequality (5.6) and inequality (5.33) as
$ σ1q1ˉS1+σ2q2ˉS2<1. $ | (5.34) |
Now we focus on matrix $ \bar{M}_3 $. With the use of expressions,
$ ti=σi(ˉJ1+ˉJ2)−βiiˉSi+ki,ai=σi(ˉJ1+ˉJ2)+βiiˉJi,si=βijˉSi,yi=βijˉJi,wherej=3−i, $ | (5.35) |
we rewrite $ \bar{M}_3 $ as
$ \bar{M}_3 = \left(−t1g1s10a1−r1y10s20−t2g2y20a2−r2\right), $ |
which has a similar form as matrix $ M_3^* $ from (5.15). This similarity allows us to apply reasoning for $ M_3^* $ to $ \bar{M}_3 $. Analogically to conditions (5.6)–(5.10), we obtain inequalities (5.26)–(5.30).
Similarly as for Theorem 5.2, let us strengthen conditions from Theorem 5.3 providing the local stability of $ E_B $. Since $ r_i > g_i $, we replace inequality (5.28) by
$ ri(ki−βiiˉSi)>giβiiˉJi, $ | (5.36) |
which obviously requires fulfillment of inequality (5.26). Using Eq (3.8), we express the above inequality as
$ ˉSi<ri(μi+αi)ki−Ciβiigiβii((ai+γi+αi+μi)μi+riαi), $ | (5.37) |
under fulfillment of
$ Ci<ri(μi+αi)kiβiigi. $ | (5.38) |
Again using the dependence $ r_i > g_i $, we simplify inequality (5.29) to
$ (rjgj(kj−βjjˉSj)−βjjˉJj)(σi(ˉJ1+ˉJ2)−βiiˉSi+ki+ri)>β12β21ˉSi(rjgjˉSj+ˉJj),j=3−i. $ | (5.39) |
Since inequality (5.26) holds, it is enough that inequality (5.36) holds so that inequality (5.39) makes sense.
Inequality (5.30) can be strengthened by
$ \prod\limits_{m = 1}^2 \Big( g_m \beta_{mm} J_m + r_m ( \beta_{mm} \bar{S}_m - k_m ) \Big) > \prod\limits_{m = 1}^2 \Big( r_m \beta_{mj} \bar{S}_m + g_m \beta_{mj} \bar{J}_m \Big), $ |
which can be expressed as
$ (β11β22−β12β21)(g1ˉJ1+r1ˉS1)(g2ˉJ2+r2ˉS2)+β11k2r2(r1ˉS1−g1ˉJ1)+β22k2r1(r2ˉS2−g2ˉJ2)+r1r2k1k2>0. $ | (5.40) |
The above inequality is always true if
$ β11β22>β12β21 $ | (5.41) |
and
$ riˉSi>giˉJi. $ | (5.42) |
Using Eq (3.8), we rewrite inequality (5.42) as
$ ˉSi>Cigiμigi+μiri+αiri. $ | (5.43) |
Let us compare inequalities (5.26) and (5.27). Observe that $ r_i > k_i $. Moreover, if inequality (5.41) holds, then inequality (5.26) is stronger than inequality (5.27).
Finally, we conclude that
Corollary 3. If inequalities (5.26), (5.37), (5.38), (5.39), (5.40), (5.41), and (5.43) hold, then $ E_B $ is locally stable.
In Subsection 3.2, we provided only a slight analysis of the existence of the endemic stationary state $ E_E $, with two diseases present. We are therefore not certain if this state exists. Furthermore, the complexity of the proper Jacobian matrix does not allow us to obtain the explicit conditions for local stability of such a postulated equilibrium. However, the conditions from Theorems 5.2 and 5.3 restrict a set of parameters' values guaranteeing the local stability of existing states $ E_A $ and $ E_B $. Such restriction suggests that there should be ranges of the values for the $ E_E $ local stability under its existence. Indicating these ranges is difficult, even numerically, because of the system's intricacy. For this reason, for each parameter, we only give specific values that yield the desirable local stability. We rely on the values from [10] that concern $ TB $ and COVID-19. In our system, these diseases correspond to diseases $ A $ and $ B $, respectively. Since paper [10] relates to co-infection dynamics in a homogeneous population, we arbitrarily choose the values of incompatible parameters from our system. These values are chosen so that we reach the $ E_E $ local stability. In Table 1, one can find the taken numbers.
Symbol | Value |
$ C_1 $ | 130 |
$ C_2 $ | 10 |
$ \beta_{11} $ | $ 2 \cdot 10^{-6} $ |
$ \beta_{12} $ | $ 8 \cdot 10^{-6} $ |
$ \beta_{21} $ | $ 3 \cdot 10^{-6} $ |
$ \beta_{22} $ | $ 6.5 \cdot 10^{-6} $ |
$ \sigma_1 $, $ \sigma_2 $ | $ 5.5 \cdot 10^{-6} $ |
$ \gamma_1 $, $ \gamma_2 $ | $ 0.02 $ |
$ \mu_1 $, $ \mu_2 $ | $ \tfrac{1}{59.365} $ |
$ g_1 $, $ g_2 $ | $ 0.015 $ |
$ \alpha_1, \alpha_2 $ | $ 0.004 $ |
$ a_1 $, $ a_2 $ | $ 0.0018 $ |
We illustrate the local stability of $ E_E $ on the plots showing the dependence of the system's solution on time for each particular variable. For the illustration, we use the Matlab software, which provides a built-in $ ode45 $ function. This function numerically solves a given differential equation system for a specified initial condition [27]. For the simulation, we take the parameters' values from Table 1 and the arbitrarily chosen initial condition
$ \Big( S_1(0), I_1(0), J_1(0), K_1(0), S_2(0), I_2(0), J_2(0), K_2(0) \Big) = (5000,100,300, 10,500, 60, 70, 5). $ |
Figures 2–3 depict the result of the simulation. For the figures' transparency, we show plots for the infected variables for each subpopulation and the non-infected variables in separate picture graphs. The obtained figures suggest the existence of the stationary state $ E_e $ that is locally stable for the parameter values from Table 1.
Now for the illustrated example of epidemic, we depict the relative sizes of the infections for time $ t $ that we define by ratios:
$ R_1(t): = \frac{I_1(t)+K_1(t)}{J_1(t)+K_1(t)}, \qquad R_2(t): = \frac{I_2(t)+K_2(t)}{J_2(t)+K_2(t)}, \qquad R(t) : = \frac{I_1(t)+I_2(t)+K_1(t)+K_2(t)}{J_1(t)+J_2(t)+K_1(t)+K_2(t)}. $ |
These ratios correspond to $ LS $, $ HS $, and the whole population, respectively, and in our case represent the number of $ TB $-infected individuals relative to the number of COVID-infected ones. Figure 4 shows the dependence of the relative sizes on time.
For the last point of the simulation timescale, i.e., $ \bar{t} = 1600 $, we get $ R_1(\bar{t}) = 0.1176 $ and $ R_2(\bar{t}) = 0.3432 $. Plots from the figure suggest the convergence of each relative size. We therefore state that for the stabilized co-infection epidemic, for one $ TB $-infected person, there are approximately nine and three COVID-infected people in $ LS $ and $ HS $, respectively.
In this paper, we proposed and analyzed the continuous-time model (2.1) describing co-infection in a heterogeneous population, in which we distinguish two subpopulations. These subpopulations, low-risk $ LS $ and high-risk $ HS $, differ in the risk of getting infected by any of two diseases, called disease $ A $ ($ DA $) and disease $ B $ ($ DA $). The values of the parameters for every subpopulation are different, which guarantees complete population heterogeneity. System (2.1) has three stationary states: disease-free ($ E_{df} $), with sole $ DA $ or $ DB $ ($ E_A $ and $ E_B $). We also suspect that the endemic state, with two diseases present, exists, but we did not manage to prove it because of complicated computations. State $ E_{df} $ exists unconditionally, while provided conditions determine the existence of $ E_A $ and $ E_B $. For state $ E_e $, we only gave insight into its existence because of the complexity of the computations. For system (2.1), we computed the basic reproduction number $ \mathcal{R}_0 $. This number is the maximum of two terms, whose forms depend on parameters corresponding to the particular sole infection. Later, we investigated the local stability of the stationary state. State $ E_{df} $ is locally stable if $ \mathcal{R}_0 < 1 $, which is expected. Analysis of the local stability for $ E_A $ and $ E_B $ provided the list of conditions. Importantly, the parameters from both diseases affect the local stability of both states.
The proposed model expands the system from [24], where we investigated the epidemic dynamics of one disease in heterogeneous populations. In that system, there are only two stationary states: the disease-free state and the endemic state, which is a counterpart of state $ E_A $ of system (2.1). The results for their local stabilities are analogical to those for state $ E_{df} $ and $ E_A $ in this paper.
The next step of our work will be an analysis of the proposed model with some assumptions simplifying this form. We hope to obtain the endemic state of the simplified model and get explicit results concerning its local stability.
The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.
The author declares there is no conflict of interest.
[1] | UNICEF Report, Eastern and Southern Africa. Available from: https://www.unicef.org. |
[2] | Dessalegn R (1999) Water resource and irrigation development in Ethiopia, Forum for Social Studies. Available from: http://www.ethiopians.com/Main_FSS_Paper1.htm. |
[3] |
Ramos JS, Helena M, Ramos (2009) Solar powered pumps to supply water for rural or isolated zones: a case study. Energ Sust Develop 13: 151-158. doi: 10.1016/j.esd.2009.06.006
![]() |
[4] | Asefa K, Abha R, Chaubey UC (2013) Solar pump application in rural water supply - a case study from Ethiopia. Int J Energ Eng 3: 176-182. |
[5] | Argaw N, Foster R, Ellis A (2003) Renewable energy for water pumping applications in rural villages, New Mexico State University Lass Cruces, New Mexico. |
[6] | Meah K, Fletcher S, Ula S (2006) Solar photovoltaic water pumping for remote locations. Renew Sust Energ Rev 12: 472-487. |
[7] |
Raturi A (2011) Feasibility study of a solar water pumping system, university of the south pacific, Fiji islands. Appl Sol Energy 47: 11-13. doi: 10.3103/S0003701X11010129
![]() |
[8] |
Gopal C, Mohanraj M, Chandramohan P, et al. (2013) Renewable energy source water pumping systems—A literature review. Renew Sustain Energ Rev 25: 351-370. doi: 10.1016/j.rser.2013.04.012
![]() |
[9] |
Pietro E, Hailong L, Jinyue Y. (2013) Dynamic modeling of a PV pumping system with special consideration on water demand. Appl Energ 112: 635-645. doi: 10.1016/j.apenergy.2012.12.073
![]() |
[10] | Panigrahi CK, Parida PR, Das M, et al. (2014) Design and modeling of photovoltaic water pumping system. IJLTEMAS 3: 2278-2540. |
[11] | Ghoneim AA (2005) Design optimization of photovoltaic powered water pumping systems. Energ Convers Manage 47: 1449-1463. |
[12] | Mansaray KG (2014) Optimum design of solar photovoltaic pumping systems by computer simulation. Int J Emer Techn Adv Eng 4: 2250-2459. |
[13] | Robert F, Cota A (2013) Solar water pumping advances and comparative economics, solar world congress. International Solar Energy Society (ISES), Cancun, Quintana Roo, Mexico. |
[14] | Robert F, Cisneros G, Hanley C (1998) Life cycle cost analysis for photovoltaic water pumping system in Mexico, 2nd world conference on photovoltaic solar energy conversion, Vienna, Austria. |
[15] | Solar Electric Light Fund, A cost and reliability comparison between solar and diesel powered pumps, solar electric light fund. 2008. Available from: self.org/SELF_White_Paper_-_Solar_vs_Diesel.pdf. |
[16] |
Getachew B, Palm B (2010) Feasibility study for a standalone solar-wind-based hybrid energy system for application in Ethiopia. Appl Energ 87: 487-495. doi: 10.1016/j.apenergy.2009.06.006
![]() |
[17] |
Wolde-Ghiorgis W (2002) Renewable energy for rural development in Ethiopia: the case for new energy policies and institutional reform. Energ Policy 30: 1095-1105. doi: 10.1016/S0301-4215(02)00061-7
![]() |
[18] | Master plan report of wind and solar energy in the federal democratic republic of Ethiopia, final version, HYDROCHINA Corporation, July 2012. Available from: http://www.mowie.gov.et/documents/714785/1953720/MP+Report+of+Wind+and+Solar+Energy/c53ac65a-a84d-448c-86d8-2e6f44aa3e39?version=1.0 |
[19] | Renewable energy in Norway. Available from: http://en.wikipedia.org/wiki/Renewable energy in Norway. |
[20] | Chapter 6 photovoltaic pumping system. Available from: http://link.springer.com/content/pdf/. |
[21] | User's guide, PVsyst contextual help, PVsyst SA 1994-2012. Available from: http://files.pvsyst.com/pvsyst5.pdf. |
[22] | Muluken Z, Tassew T, Abdulkadir A (2014) Optimal sizing of solar water pumping system for small, scale irrigation: case study of Dangila. Int J Renew Sust Energ 3: 99-107. |
[23] | Irfan G, Nevzat O, (2010) Cost Calculation Algorithm for Photovoltaic Systems Paths to Sustainable Energy, Paths to Sustainable Energy. Available from: http://www.intechopen.com/books/paths-to-sustainable-energy/cost-calculation-algorithm-for-photovoltaic-systems. |
[24] | Pandey R, Gaur MK, Malvi CS (2012) Estimation of Cost Analysis for 4 kw Grids Connected Solar Photovoltaic Plant. Int J Mod Eng Res 2: 4292-4294. |
[25] | Avinash K, Purva S, Apurva V (2014) Design and Cost Analysis of PV System Using Nano Solar Cell. Int J Sci Res Publ 4: 2250-3153. |
[26] | Chan S. Park (2012) Fundamentals of engineering economics, 3 Eds, Prentice Hall. |
Symbol | Value |
$ C_1 $ | 130 |
$ C_2 $ | 10 |
$ \beta_{11} $ | $ 2 \cdot 10^{-6} $ |
$ \beta_{12} $ | $ 8 \cdot 10^{-6} $ |
$ \beta_{21} $ | $ 3 \cdot 10^{-6} $ |
$ \beta_{22} $ | $ 6.5 \cdot 10^{-6} $ |
$ \sigma_1 $, $ \sigma_2 $ | $ 5.5 \cdot 10^{-6} $ |
$ \gamma_1 $, $ \gamma_2 $ | $ 0.02 $ |
$ \mu_1 $, $ \mu_2 $ | $ \tfrac{1}{59.365} $ |
$ g_1 $, $ g_2 $ | $ 0.015 $ |
$ \alpha_1, \alpha_2 $ | $ 0.004 $ |
$ a_1 $, $ a_2 $ | $ 0.0018 $ |
Symbol | Value |
$ C_1 $ | 130 |
$ C_2 $ | 10 |
$ \beta_{11} $ | $ 2 \cdot 10^{-6} $ |
$ \beta_{12} $ | $ 8 \cdot 10^{-6} $ |
$ \beta_{21} $ | $ 3 \cdot 10^{-6} $ |
$ \beta_{22} $ | $ 6.5 \cdot 10^{-6} $ |
$ \sigma_1 $, $ \sigma_2 $ | $ 5.5 \cdot 10^{-6} $ |
$ \gamma_1 $, $ \gamma_2 $ | $ 0.02 $ |
$ \mu_1 $, $ \mu_2 $ | $ \tfrac{1}{59.365} $ |
$ g_1 $, $ g_2 $ | $ 0.015 $ |
$ \alpha_1, \alpha_2 $ | $ 0.004 $ |
$ a_1 $, $ a_2 $ | $ 0.0018 $ |