Analysis the Dynamics of SIHR Model: Covid-19 Case in Djibouti

Yahyeh Souleiman^*, Abdoulrazack Mohamed, Liban Ismail

Mathematics and Computer Science Laboratory, University of Djibouti, Campus Balbala, Djibouti.

DOI: 10.4236/am.2021.1210057                

Abstract

The Covid-19 epidemic is an emerging infectious disease of the viral zoonosis type caused by the coronavirus strain SARS-CoV-2, it is classified as a human-to-human communicable disease and is currently a pandemic worldwide. In this paper, we propose conceptual mathematical models of the epidemic dynamics of four compartments. We have collected data from the Djibouti health ministry. We study the positivity, boundedness, existence and uniqueness of the weak solution. Next, we define the Basic reproduction number by the method of the DFE and EEP. Then, we study the local and global stability and the bifurcation analysis of equilibrium to examine its epidemiological relevance. Finally, we analyze the fit of the data in comparison with the result of our mathematical results, to validate the model and estimate the important model parameters and prediction about the disease. We consider the real cases of Djibouti from 15th March to 15th May 2021.

Keywords

SIHR Model, Covid-19, Basic Reproduction Number, Stability Analysis, Mod-el Validation

1. Introduction

Covid-19 is an infectious disease caused by a coronavirus, Sars-CoV-2, which primarily affects the respiratory tract. Transmission occurs through close contact with a person carrying the virus (mainly through direct contact or through respiratory secretions emitted into the air when coughing, sneezing or talking). This disease erupted in China at the end of 2019 and the virus was identified in early January 2020. The World Health Organization (WHO) declared the disease as a pandemic and was named SARS-CoV-2 virus (March 11, 2020). In August 2021, the number of infected cases was confirmed by who reached 202,608,306, and by this date, there were 4,293,591 deaths related to the disease worldwide [1] . On the other hand, the Republic of Djibouti declared its first case of Covid-19 on March 23 (2020) and recorded a very low mortality rate.

Disease models play an important role in understanding and managing the transmission dynamics of various pathogens. We can use them to describe the spatial and temporal patterns of disease prevalence, as well as to explore or better understand the factors that influence infection incidence. Modeling is a key step in understanding what treatments and interventions can be the most effective, how cost-effective these approaches may be, and what specific factors need to be considered when trying to eradicate disease. Some recent studies provided different guidelines by introducing the basic reproduction number, education and socio-economic index and lock-down strategies [2] [3] [4] [5] [6] .

New variants of the virus are circulating, in particular, the Alpha, Beta, Gamma and Delta variants, initially discovered in the United Kingdom, South Africa, Brazil and India respectively. In France, the Delta variant is very predominant. This variant is more contagious, but the vaccines generally retain their efficacy against this variant provided that the vaccination is complete. Djibouti has not been spared by the new variant of Covid-19. The most contagious variant spread is hit hard compared to the first wave. The highest average number of daily contamination cases was reported on April 4. There have been 11,663 cases of contamination and 156 deaths linked to the coronavirus recorded in Djibouti since the start of the epidemic with 53,016 people vaccinated. Some academic studies are conducted in this context, Foy, B.H., etal. [7] , Yavuz, M., etal. [8] , Kaur and Gupta [9] , have shown different types of vaccine strategies for Covid-19.

The simplest way to model epidemic spread in populations is to classify people into different population groups or compartments. Compartmental models are governed by a system of differential equations that track the population as a function of time, stratifying it into different groups based on the risk or infection status, several extensions have been proposed to explain the evolution of Covid-19 disease [10] [11] . These models have been widely used by scientists around the world to model more complex compartments [12] [13] [14] .

Compartmental models are deterministic, that is, given the same inputs, they produce the same results every time. They are able to predict the various properties of pathogen spread, can estimate the duration of epidemics, and can be used to understand how different situations or interventions can impact the outcome of pathogen spread. Several extensions have been proposed to explain the evolution of diseases Covid-19 [15] [16] .

Therefore, the model with multiple compartments is a useful tool to predict the nature of Covid-19. Thereby, to predict the nature of the Covid-19 Djibouti, in this works we develop the SIHR model that stratifies a population based on their infection status as a function of time such as susceptible class (S), infected class (I), persons Hospitalized (H) and removed class (R). The key objectives of this study are as follows:

1) We analyze the stability of the equilibrium of the model using the basic reproduction number to understand its severity.

2) Theoretical results are established using local and global analysis of the model.

3) Numerical demonstrations validated the analytic outcomes.

The paper is organized as follows: Mathematical Model is elaborately discussed in section 2. Positivity and boundedness of solution including auxiliary results are described in section 3. In section 4, we present the result of existence and uniqueness of the solution. Then, in section 5, we study the basic reproduction number DFE and EEP. Local-global stability analyses are prescribed in section 6. Finally, Section 7 provides the data analysis in comparison with the model solution with further prediction to control the epidemic, as a case study in Djibouti.

In the following section, we will thoroughly discuss our mathematical model and formulation of the model elaborately.

2. Model

To study the epidemic of Covid-19 spread at Djibouti, we extend the classical deterministic susceptible-infectious-removed (SIR) epidemic model by adding a hospitalized compartment. In the proposed model, let the total population is divided into four compartments; those are the susceptible population (S), infected population (I), hospitalized population (H) and recovered population (R).

The differential system that describes this SIHR model is

$\{\begin{array}{l}\frac{\text{d}S}{\text{d}t}=\tau -\mu S-\beta SI\hfill \\ \frac{\text{d}I}{\text{d}t}=\beta SI-\left(\nu +\mu +\gamma +\alpha \right)I\hfill \\ \frac{\text{d}H}{\text{d}t}=\alpha I-\left(\nu +\mu +\lambda \right)H\hfill \\ \frac{\text{d}R}{\text{d}t}=\gamma I+\lambda H-\mu R.\hfill \end{array}$ (1)

with the initial conditions $S\left(0\right)>0,I\left(0\right)>0,H\left(0\right)\ge 0,R\left(0\right)\ge 0$ where the interpretation of the parameters is presented in Table 1 below.

The flow diagram of the proposed model in this paper is given in Figure 1.

3. Positivity and Boundedness of Solution

In this section, we present the study of positivity and limitation of the solution. Taking into account the nonlinear system of Equation (1), we consider the first equation

$\frac{\text{d}S}{\text{d}t}=\tau -\mu S-\beta SI,$ (2)

Table 1. Model parameters and their descriptions.

Figure 1. SIHR model transmission diagram.

which means that

$\frac{\text{d}S}{\text{d}t}\ge -\left(\mu +\beta I\right)S\mathrm{,}$ (3)

By using the exponential growth criterion and integrating (3) gives

$S\left(t\right)\ge S\left(0\right){\text{e}}^{-\left(\mu +\beta I\right)t}\mathrm{,}$ (4)

which implies

$S\left(t\right)\ge 0.$

By following the similar steps with the condition in $S\left(t\right)$, it can easily be shown that $I\left(t\right)\ge \mathrm{0,}H\left(t\right)\ge 0$ and $R\left(t\right)\ge 0$ $\forall t>0$.

Next, the population (N) is then divided into four classes; the susceptible (S), the infected (I), the hospitalized (H), and the recovered (R) at any time $t\ge 0$,

$N\left(t\right)=S\left(t\right)+I\left(t\right)+H\left(t\right)+R(\; t\; )$

with the initial conditions $S\left(0\right)>0,I\left(0\right)>0,H\left(0\right)\ge 0,R\left(0\right)\ge 0$.

Now we use model (1) we have

$\frac{\text{d}N}{\text{d}t}=\frac{\text{d}S}{\text{d}t}+\frac{\text{d}I}{\text{d}t}+\frac{\text{d}H}{\text{d}t}+\frac{\text{d}R}{\text{d}t}=\tau -\mu N-\nu I-\nu H$.

Which give

$\frac{\text{d}N}{\text{d}t}\le \tau -\mu N$ (5)

Integrating the inequality (5), using initial condition, we obtain

$N\left(t\right)\le N\left(0\right){\text{e}}^{-\mu t}+\frac{\tau }{\mu }\left(1-{\text{e}}^{-\mu t}\right)$.

Letting t tends to infinity, asymptotically we get $N\left(t\right)<\frac{\tau }{\mu }$.

Thus we can summarize the above results in the following theorem.

Theorem 1. The closed region $\Omega =\left\{\left(S,I,H,R\right)\in {ℝ}_{+}^4:0<N<\frac{\tau }{\mu }\right\}$ is positively invariant set for the system (1).

4. Existence and Uniqueness

In this section, we present the result of existence and uniqueness of the solution of the COVID-19 model which is given in system (1). First, by integrating the system (1) over the interval $\left[\mathrm{0,}t\right]$, we obtain

$S\left(t\right)-S\left(0\right)={\displaystyle {\int }_0^t}{\Phi }_1\left(z,S\right)\text{d}z,I\left(t\right)-I\left(0\right)={\displaystyle {\int }_0^t}{\Phi }_2\left(z,I\right)\text{d}z,$

$H\left(t\right)-H\left(0\right)={\displaystyle {\int }_0^t}{\Phi }_3\left(z,H\right)\text{d}z,R\left(t\right)-R\left(0\right)={\displaystyle {\int }_0^t}{\Phi }_4\left(z,R\right)\text{d}z.$

With

${\Phi }_1\left(t\right)=\tau -\left(\mu +\beta I\left(t\right)\right)S\left(t\right),{\Phi }_2\left(t\right)=\left(\beta S\left(t\right)-\left(\nu +\mu +\gamma +\alpha \right)\right)I\left(t\right),$

${\Phi }_3\left(t\right)=\alpha I\left(t\right)-\left(\nu +\mu +\lambda \right)H\left(t\right),{\Phi }_4\left(t\right)=\gamma I\left(t\right)+\lambda H\left(t\right)-\mu R\left(t\right),$

We have the result of the following theorem:

Theorem 2. The kernels ${\Phi }_1\mathrm{,}{\Phi }_2\mathrm{,}{\Phi }_3$ and ${\Phi }_4$ satisfy the Lipschitz assumptions and contractions if the following inequality is verified:

$0\le k_1,k_2,k_3,k_4<1.$ (6)

where, $\Vert S\Vert \le c_1\mathrm{,}\Vert I\Vert \le c_2\mathrm{,}\Vert H\Vert \le c_3\mathrm{,}\Vert R\Vert \le c_4$, $k_1=\mu +\beta c_2,k_2=\beta c_1+\nu +\mu +\gamma +\alpha ,k_3=\nu +\mu +\lambda$ and $k_4=\mu$.

Proof 3. We consider the two functions $S_1$ and $S_2$ for the kernel ${\Phi }_1$ ; $I_1$ and $I_2$ be two functions for the kernel ${\Phi }_2$ ; $H_1$ and $H_2$ be two functions for the kernel ${\Phi }_3$ ; and $R_1$ and $R_2$ be two functions for the kernel ${\Phi }_4$, Then we have

$\Vert {\Phi }_1\left(t\mathrm{,}{S}_1\right)-{\Phi }_1\left(t\mathrm{,}{S}_2\right)\Vert =\Vert \left(\mu +\beta I\right)\left(S_1-S_2\right)\Vert \le k_1\Vert S_1-S_2\Vert$,

$\Vert {\Phi }_2\left(t\mathrm{,}{I}_1\right)-{\Phi }_2\left(t\mathrm{,}{I}_2\right)\Vert =\Vert \left(\beta S-\left(\nu +\mu +\gamma +\alpha \right)\right)\left(I_1-I_2\right)\Vert \le k_2\Vert I_1-I_2\Vert$,

$\Vert {\Phi }_3\left(t\mathrm{,}{H}_1\right)-{\Phi }_2\left(t\mathrm{,}{H}_2\right)\Vert =\Vert \left(\nu +\mu +\lambda \right)\left(H_1-H_2\right)\Vert \le k_3\Vert H_1-H_2\Vert$,

and

$\Vert {\Phi }_4\left(t\mathrm{,}{R}_1\right)-{\Phi }_4\left(t\mathrm{,}{R}_2\right)\Vert =\Vert \mu \left(R_1-R_2\right)\Vert \le k_4\Vert R_1-R_2\Vert$.

Therefore, the Lipschitz conditions are satisfied for kernels ${\Phi }_1\mathrm{,}{\Phi }_2\mathrm{,}{\Phi }_3\mathrm{,}{\Phi }_4$ and if and if condition (6) is verified then $k_1\mathrm{,}{k}_2\mathrm{,}{k}_3\mathrm{,}{k}_4$ are also contractions for ${\Phi }_1\mathrm{,}{\Phi }_2\mathrm{,}{\Phi }_3\mathrm{,}{\Phi }_4$ respectively. $\square$

By considering the kernels ${\Phi }_1\mathrm{,}{\Phi }_2\mathrm{,}{\Phi }_3$ and ${\Phi }_4$ previously defined, we can proceed with the following recursive formula

$S_n\left(t\right)=S\left(0\right)+{\displaystyle {\int }_0^t}{\Phi }_1\left(z,S_{n-1}\right)\text{d}z,I_n\left(t\right)=I\left(0\right)+{\displaystyle {\int }_0^t}{\Phi }_2\left(z,I_{n-1}\right)\text{d}z,$

$H_n\left(t\right)=H\left(0\right)+{\displaystyle {\int }_0^t}{\Phi }_3\left(z,H_{n-1}\right)\text{d}z,R_n\left(t\right)=R\left(0\right)+{\displaystyle {\int }_0^t}{\Phi }_4\left(z,R_{n-1}\right)\text{d}z,$

where $S_0\left(t\right)=S\left(0\right),I_0\left(t\right)=I\left(0\right),H_0\left(t\right)=H\left(0\right)$ and $R_0\left(t\right)=R\left(0\right)$. Then we can write

$\Vert {\Psi }_n\left(t\right)\Vert =\Vert S_n\left(t\right)-S_{n-1}\left(t\right)\Vert \le \Vert {\displaystyle {\int }_0^t}\left({\Phi }_1\left(z,S_{n-1}\right)-{\Phi }_1\left(z,S_{n-2}\right)\right)\text{d}z\Vert ,$

$\Vert {\Theta }_n\left(t\right)\Vert =\Vert I_n\left(t\right)-I_{n-1}\left(t\right)\Vert \le \Vert {\displaystyle {\int }_0^t}\left({\Phi }_2\left(z\mathrm{,}{I}_{n-1}\right)-{\Phi }_2\left(z\mathrm{,}{I}_{n-2}\right)\right)\text{d}z\Vert \mathrm{,}$

$\Vert {\Gamma }_n\left(t\right)\Vert =\Vert H_n\left(t\right)-H_{n-1}\left(t\right)\Vert \le \Vert {\displaystyle {\int }_0^t}\left({\Phi }_3\left(z,H_{n-1}\right)-{\Phi }_3\left(z,H_{n-2}\right)\right)\text{d}z\Vert ,$

$\Vert {\Sigma }_n\left(t\right)\Vert =\Vert R_n\left(t\right)-R_{n-1}\left(t\right)\Vert \le \Vert {\displaystyle {\int }_0^t}\left({\Phi }_4\left(z,R_{n-1}\right)-{\Phi }_4\left(z,R_{n-2}\right)\right)\text{d}z\Vert .$

Since the kernels satisfy the Lipschitz condition (see Theorem 2), we get

$\{\begin{array}{l}\Vert {\Psi }_n\left(t\right)\Vert \le k_1{\displaystyle {\int }_0^t}\Vert S_{n-1}-S_{n-2}\Vert \text{d}z\le k_1{\displaystyle {\int }_0^t}\Vert {\Psi }_{n-1}\left(z\right)\Vert \text{d}z\hfill \\ \Vert {\Theta }_n\left(t\right)\Vert \le k_2{\displaystyle {\int }_0^t}\Vert I_{n-1}-I_{n-2}\Vert \text{d}z\le k_2{\displaystyle {\int }_0^t}\Vert {\Theta }_{n-1}\left(z\right)\Vert \text{d}z\mathrm{,}\hfill \\ \Vert {\Gamma }_n\left(t\right)\Vert \le k_3{\displaystyle {\int }_0^t}\Vert H_{n-1}-H_{n-2}\Vert \text{d}z\le k_3{\displaystyle {\int }_0^t}\Vert {\Gamma }_{n-1}\left(z\right)\Vert \text{d}z\mathrm{,}\hfill \\ \Vert {\Sigma }_n\left(t\right)\Vert \le k_4{\displaystyle {\int }_0^t}\Vert R_{n-1}-R_{n-2}\Vert \text{d}z\le k_4{\displaystyle {\int }_0^t}\Vert {\Sigma }_{n-1}\left(z\right)\Vert \text{d}z\mathrm{.}\hfill \end{array}$ (7)

Therefore, we obtain results of the following theorem.

Theorem 4. The model of the system (1) that we have proposed has a solution under the condition:

$k_i{t}_{\max }<1,i=1,2,3,4.$ (8)

Proof 5. Assuming that the hypothesis of Theorem 2 for the functions 1, 2, 3, and 4 hold, we can give the following by taking Equation (7) into account

$\{\begin{array}{l}\Vert {\Psi }_n\left(t\right)\Vert \le \Vert S_0\left(t\right)\Vert {\left(k_1{t}_{\max }\right)}^n,\Vert {\Theta }_n\left(t\right)\Vert \le \Vert I_0\left(t\right)\Vert {\left(k_2{t}_{\max }\right)}^n,\hfill \\ \Vert {\Gamma }_n\left(t\right)\Vert \le \Vert H_0\left(t\right)\Vert {\left(k_3{t}_{\max }\right)}^n,\Vert {\Sigma }_n\left(t\right)\Vert \le \Vert R_0\left(t\right)\Vert {\left(k_4{t}_{\max }\right)}^n.\hfill \end{array}$ (9)

Now we assume

$\{\begin{array}{l}S\left(t\right)-S\left(0\right)=S_n\left(t\right)-f_n\left(t\right),I\left(t\right)-I\left(0\right)=I_n\left(t\right)-g_n\left(t\right),\hfill \\ H\left(t\right)-H\left(0\right)=H_n\left(t\right)-h_n\left(t\right),R\left(t\right)-R\left(0\right)=R_n\left(t\right)-{\rho }_n\left(t\right).\hfill \end{array}$ (10)

Then we get

$\Vert f_n\left(t\right)\Vert \le \Vert {\displaystyle {\int }_0^t}\left({\Phi }_1\left(z\mathrm{,}S\right)-{\Phi }_1\left(z\mathrm{,}{S}_{n-1}\right)\right)\text{d}z\Vert \le t{k}_1\Vert S-S_{n-1}\Vert \mathrm{,}$

$\Vert g_n\left(t\right)\Vert \le \Vert {\displaystyle {\int }_0^t}\left({\Phi }_2\left(z\mathrm{,}I\right)-{\Phi }_2\left(z\mathrm{,}{I}_{n-1}\right)\right)\text{d}z\Vert \le t{k}_2\Vert I-I_{n-1}\Vert \mathrm{,}$

$\Vert h_n\left(t\right)\Vert \le \Vert {\displaystyle {\int }_0^t}\left({\Phi }_3\left(z\mathrm{,}H\right)-{\Phi }_3\left(z\mathrm{,}{H}_{n-1}\right)\right)\text{d}z\Vert \le t{k}_3\Vert H-H_{n-1}\Vert \mathrm{,}$

$\Vert {\rho }_n\left(t\right)\Vert \le \Vert {\displaystyle {\int }_0^t}\left({\Phi }_4\left(z\mathrm{,}R\right)-{\Phi }_4\left(z\mathrm{,}{R}_{n-1}\right)\right)\text{d}z\Vert \le t{k}_4\Vert R-R_{n-1}\Vert \mathrm{,}$

repeating this process recursively and considering $t=t_{\max }$, taking the upper bound of the time interval into account, we get

$\{\begin{array}{l}\Vert f_n\left(t\right)\Vert \le {\left(t_{\max }\right)}^{n+1}{k}_1^n\sigma ,\Vert g_n\left(t\right)\Vert \le {\left(t_{\max }\right)}^{n+1}{k}_2^n\sigma ,\hfill \\ \Vert h_n\left(t\right)\Vert \le {\left(t_{\max }\right)}^{n+1}{k}_3^n\sigma ,\Vert {\rho }_n\left(t\right)\Vert \le {\left(t_{\max }\right)}^{n+1}{k}_4^n\sigma .\hfill \end{array}$ (11)

Let n tends to infinity, we then obtain

$\Vert f_{\infty }\left(t\right)\Vert \to \mathrm{0,}\Vert g_{\infty }\left(t\right)\Vert \to \mathrm{0,}\Vert h_{\infty }\left(t\right)\Vert \to 0\text{and}\Vert {\rho }_{\infty }\left(t\right)\Vert \to 0.$

$\square$

Theorem 6. System model (1) has only one solution.

Proof 7. We assume that system (1) admits two solutions, $S_1\mathrm{,}{S}_2\mathrm{,}{I}_1\mathrm{,}{I}_2\mathrm{,}{H}_1\mathrm{,}{H}_2$ and $R_1\mathrm{,}{R}_2$. Using the Lipschitz condition satisfied by the functions ${\Phi }_1\mathrm{,}{\Phi }_2\mathrm{,}{\Phi }_3$ and ${\Phi }_4$, we obtain

$\Vert S_1\left(t\right)-S_2\left(t\right)\Vert =\Vert {\displaystyle {\int }_0^t}\left({\Phi }_1\left(z,S_1\right)-{\Phi }_1\left(z,S_2\right)\right)\text{d}z\Vert \le k_{1}t\Vert S_1\left(t\right)-S_2\left(t\right)\Vert ,$

$\Vert I_1\left(t\right)-I_2\left(t\right)\Vert =\Vert {\displaystyle {\int }_0^t}\left({\Phi }_2\left(z,I_1\right)-{\Phi }_2\left(z,I_2\right)\right)\text{d}z\Vert \le k_{2}t\Vert I_1\left(t\right)-I_2\left(t\right)\Vert ,$

$\Vert H_1\left(t\right)-H_2\left(t\right)\Vert =\Vert {\displaystyle {\int }_0^t}\left({\Phi }_3\left(z,H_1\right)-{\Phi }_3\left(z,H_2\right)\right)\text{d}z\Vert \le k_{3}t\Vert H_1\left(t\right)-H_2\left(t\right)\Vert ,$

$\Vert R_1\left(t\right)-R_2\left(t\right)\Vert =\Vert {\displaystyle {\int }_0^t}\left({\Phi }_4\left(z,R_1\right)-{\Phi }_4\left(z,R_2\right)\right)\text{d}z\Vert \le k_{4}t\Vert R_1\left(t\right)-R_2\left(t\right)\Vert ,$

we deduce that

$\Vert S_1\left(t\right)-S_2\left(t\right)\Vert \left(1-k_{1}t\right)\le \mathrm{0,}\Vert I_1\left(t\right)-I_2\left(t\right)\Vert \left(1-k_{2}t\right)\le \mathrm{0,}$

$\Vert H_1\left(t\right)-H_2\left(t\right)\Vert \left(1-k_{3}t\right)\le \mathrm{0,}\Vert R_1\left(t\right)-R_2\left(t\right)\Vert \left(1-k_{4}t\right)\le 0.$

Therefore, we obtain $\Vert S_1\left(t\right)-S_2\left(t\right)\Vert =0$, $\Vert I_1\left(t\right)-I_2\left(t\right)\Vert =0$, $\Vert H_1\left(t\right)-H_2\left(t\right)\Vert =0$ and $\Vert R_1\left(t\right)-R_2\left(t\right)\Vert =0$ which means that $S_1\left(t\right)=S_2\left(t\right)$, $I_1\left(t\right)=I_2\left(t\right)$, $H_1\left(t\right)=H_2\left(t\right)$ and $R_1\left(t\right)=R_2\left(t\right)$. This explains that the model has a unique solution which is the proof of the theorem. $\square$

5. Basic Reproduction Number DFE and EEP

Compartmental models are deterministic, that is, given the same inputs, they produce the same results every time. They are able to predict the various properties of the spread of the virus, can estimate the duration of epidemics and can be used to understand how different situations or interventions can affect the results of the spread. To do this, the $R_0$ parameter, describing the average number of new infections due to a sick individual, plays a crucial role. As you can imagine, if this number is less than 1 then the epidemic will tend to die out. In this case, the disease-free equilibrium (DFE) will be locally asymptotically stable and the disease cannot persist in the population. However, if $R_0>1$, it may persist or even spread to the whole population. This implies that the disease-

free equilibrium (DFE) is unstable. Using next generation matrix [17] the basic reproduction of (1) is found here. Since the DFE is $E_0={\left(\frac{\tau }{\mu }\mathrm{,0,0,0}\right)}^t$ and hence the basic reproduction number can be found using the analytical approach.

Let $G=\left(\begin{array}{c}\tau \\ \beta SI\end{array}\right)$ represents the rate of new infection matrix and $U=\left(\begin{array}{c}\mu S+\beta SI\\ \left(\nu +\mu +\gamma +\alpha \right)I\end{array}\right)$ denotes the transfer rate matrix of the individuals.

Let us we define $F=\frac{\partial G}{\partial x_j}\left(E_0\right)$ and $V=\frac{\partial U}{\partial x_j}\left(E_0\right)$, the reproduction number for the Covid-19 model given by (1) can be calculated from the relation $R_0=\rho \left(F{V}^{-1}\right)$, the spectral radius of $F{V}^{-1}$ is given below

$R_0=\frac{\beta \tau }{\mu \left(\nu +\mu +\gamma +\alpha \right)}$ (12)

To find the endemic equilibrium state of the model we set

$\frac{\text{d}S}{\text{d}t}=0,\frac{\text{d}I}{\text{d}t}=0,\frac{\text{d}H}{\text{d}t}=0,\frac{\text{d}R}{\text{d}t}=0$.

Solving the above system, we get the epidemic equilibrium (EEF) state

$E_1=\left(S_1\mathrm{,}{I}_1\mathrm{,}{H}_1\mathrm{,}{R}_1\right)$,

where

$S_1=\frac{\nu +\mu +\gamma +\alpha }{\beta },I_1=\frac{\tau \beta -\mu \left(\nu +\mu +\gamma +\alpha \right)}{\beta \left(\nu +\mu +\gamma +\alpha \right)},$

$H_1=\frac{\alpha \tau \beta -\alpha \mu \left(\nu +\mu +\gamma +\alpha \right)}{\beta \left(\nu +\mu +\lambda \right)\left(\nu +\mu +\gamma +\alpha \right)}\mathrm{,}$

$R_1=\frac{\left(\gamma \left(\nu +\mu +\lambda \right)+\lambda \alpha \right)\left(\tau \beta -\mu \left(\nu +\mu +\gamma +\alpha \right)\right)}{\beta \mu \left(\nu +\mu +\lambda \right)\left(\nu +\mu +\gamma +\alpha \right)}$.

6. Stability and Bifurcation the Equilibrium States

In this section, we shall establish the stability and bifurcation condition of the equilibrium points. In Theorem 8, we shall establish nature of the $E_0$ and in Theorem 10 nature of $E_1$.

6.1. Stability and Bifurcation of Disease-Free Equilibrium State (E₀)

Theorem 8. The DFE will be locally asymptotically stable if $R_0<1$ and unstable if $R_0>1$.

Proof 9. The Jacobian matrix corresponding to the system (1) at DFE point $E_0=\left(\frac{\tau }{\mu }\mathrm{,0,0,0}\right)$,

$J\left(E_0\right)\mathrm{<mo>\; =\; </mo>}\left(\begin{array}{cccc}-\mu & -\frac{\beta \tau }{\mu }& 0& 0\\ 0& \frac{\beta \tau }{\mu }-\left(\nu +\mu +\gamma +\alpha \right)& 0& 0\\ 0& \alpha & -\left(\nu +\mu +\lambda \right)& 0\\ 0& \gamma & \lambda & -\mu \end{array}\right)$.

The characteristic roots of the Jacobian matrix at $J\left(E_0\right)$ are $-\mu \mathrm{,}-\lambda -\nu -\mu \mathrm{,}-\mu$ and $\left(\nu +\mu +\gamma +\alpha \right)\left(R_0-1\right)$. Since the first three roots are negative the other is negative if $R_0<1$ and positive if $R_0>1$. Therefore, the disease-free equilibrium state $\left(E_0\right)$ is locally asymptotically stable if $R_0<1$ and unstable if $R_0>1$. $\square$

Since when $R_0=1$ i.e. when

$\beta =\tilde{\beta }=\frac{\mu \left(\nu +\mu +\gamma +\alpha \right)}{\tau }$

then one of the eigenvalues of the Jacobin matrix corresponding to the system (1) at DFE is zero. Using the Theorem by Castillo-Chavez and Song [18] and Martcheva Maia [19] , we investigate the nature of the disease free equilibrium points. Let $V_d$ and $V_g$ be the eigenvector corresponding to the zero eigenvalue of $J\left(E_0\right)$ and ${\left[J\left(E_0\right)\right]}^t$ respectively then

$V_d=\left(\begin{array}{c}\frac{-\left(\nu +\mu +\gamma +\alpha \right)}{\mu }\\ 1\\ \frac{\alpha }{\nu +\mu +\lambda }\\ \frac{\alpha \lambda +\gamma \left(\nu +\mu +\lambda \right)}{\mu \left(\nu +\mu +\lambda \right)}\end{array}\right)\text{and}{V}_g=\left(\begin{array}{c}0\\ 1\\ 0\\ 0\end{array}\right)$

and

let $T=\left(\begin{array}{c}\tau -\mu S-\beta SI\\ \beta SI-\left(\nu +\mu +\gamma +\alpha \right)I\\ \alpha I-\left(\nu +\mu +\lambda \right)H\\ \gamma I+\lambda H-\mu R\end{array}\right)$ then

${V_g^t{T}_{\beta }|}_{E_0,\beta =\tilde{\beta }}=0,{V_g^{t}D{T}_{\beta }|}_{E_0,\beta =\tilde{\beta }}{V}_d=\frac{\tau }{\mu }\ne 0$

$\begin{array}{l}{V_g^t{D}^{2}T|}_{E_0,\beta =\tilde{\beta }}\left(V_d,V_d\right)\\ =\frac{-2{\left(\nu +\mu +\gamma +\alpha \right)}^2}{\tau }\ne 0.\end{array}$

Hence the system experiences transcritical bifurcation when the rate of infection crosses the critical value $\beta =\tilde{\beta }$. Thus to spreading the disease the rate of transmission plays an important role. There are critical values of the rate of infection which the disease is easy to control but above of which the society will experience endemic disease spreading.

6.2. Stability of Endemic Equilibrium State (E₁)

Theorem 10. The endemic equilibrium state $E_1=\left(S_1\mathrm{,}{I}_1\mathrm{,}{H}_1\mathrm{,}{R}_1\right)$ is stable if $R_0>1$.

Proof 11. The Jacobian corresponding to the endemic equilibrium point $E_1$ is $J\left(E_1\right)$, is given below

$J\left(E_1\right)=\left(\begin{array}{cccc}-\mu -\beta I_1& -\beta S_1& 0& 0\\ \beta I_1& \beta S_1-\left(\nu +\mu +\gamma +\alpha \right)& 0& 0\\ 0& \alpha & -\left(\nu +\mu +\lambda \right)& 0\\ 0& \gamma & \lambda & -\mu \end{array}\right)$.

The characteristic polynomial of the Jacobian matrix at $J\left(E_1\right)$ is

$\det \left(J\left(E_1\right)-X{I}_4\right)=\left(-\mu -X\right)\left(-\nu -\mu -\lambda -X\right)\left(X^2-aX+b\right)$.

The roots of the characteristic equations are $X_1=-\mu <0$, $X_2=-\nu -\mu -\lambda <0$ and other two satisfies the following quadratic equation

$X^2-aX+b=0,$ (13)

the two roots of (13) are negative because

$a=\frac{-\tau \beta }{\nu +\mu +\gamma +\alpha }<0$

and

$b=\left(\mu \left(\nu +\mu +\gamma +\alpha \right)\left(R_0-1\right)\right)>0$

if $R_0>1$. $\square$

7. Parameter Estimation, Model Validation and Prediction

In this study, we considered the country Djibouti less infected with the Covid-19 virus and collected the data available online in chronological order from World Health Organization.

7.1. Parameter Estimation

Using Matlab minimization technique and fminsearch software package, we have estimated the important model parameters using the Djibouti infection cases from 15th March to 15th May, 2021 (total days 60) which are given in Table 2. To estimate the important model parameters we consider the cumulative number of infected persons from the real data source and the model predicted cumulative number of infected persons. The fitness of the model compared to the real data can be verified by computing the residual. The residuals are defined as

Table 2. Parameters estimation for Djibouti.

$\text{residuals}=\left\{Y_j-I\left(t_j\right)|j=1,2,3,\cdots ,n\right\}$

where $Y_j$ is the $j^{th}$ day cumulative infection data and $I\left(t_j\right)$ model predicted cumulative infected data of the same day. If the residuals are randomly distributed then we can say that the fitness is reasonably good [20] .

7.2. Model Validation and Prediction

To study the particular case: the spreading of Covid-19 in Djibouti we have considered the real cases in cumulative number of infection (from March 15 to May 15, 2021) and using the above described formula we have estimated the model parameters, using the global sensitivity index method we found the sensitivity index of the parameters. To execute the Matlab package, we have considered the initial population size as

$S\left(0\right)=986363,I\left(0\right)=2,H\left(0\right)=0\text{and}R\left(0\right)=0$.

The estimated model parameters and their sensitivity induces are given in Table 2. To validate the model we consider the real case of Covid-19 infection of Djibouti from 15 March to 15 May, 2021 i.e. the real cases for 60 days. In Figure 2(a) and Figure 2(b) respectively we have presented the best fitted cumulative infected cases and the curves word missing new infection of real cases and the model predictions.

Figure 3 shows that the peak of death corresponding to the peak of number of confirmed cases in Figure 2(b). As a result, the more the Covid-19 epidemic spreads, the more the number increases.

Figure 4(a) shows that the infected gradually decreases compared to the day before while we see a strong increase in deaths between the 10th day and 20th day. for whoever is in Figure 4(b), the percentage of the number of deaths and infected evolves in a regressive manner. As a result, the more the number of infected increases, the more the number of deaths from the epidemic increases.

We observe in Figure 5(a) and Figure 5(b) that the susceptible population is decreasing with time, which means more and more people are getting exposed. Since infected population is increasing, therefore they are getting a large number of infected individuals over time which can lead to an outbreak in a very short time.

If we analyze the rest of the population dynamics, we see that the infected population grew faster over the last 40 days, which shows that the spread of the epidemic over time and this corresponds well to the actual data of Figure 2(b). We also see in Figure 5(a) a lower rate infected initially meanwhile the recovered population is not growing as much as infected population in Figure 5(b). Using the parameters that are estimated in Table 2, we can determine the value of the basic reproduction number in Equation (12), with $R_0=1.8235$ and $R_0=2.5493$ of Figure 5(a) and Figure 5(b) respectively. Finally, the more the hospitalized population increases the more the number of deaths over the last 40 days, which corresponds to the period of containment.

(a) (b)

Figure 2. (a) Fitting model to cumulative cases in Djibouti; (b) The time series of new cases of Covid-19 in Djibouti in 15 March to 15 May, 2021.

Figure 3. Number of death cases in Djibouti.

(a) (b)

Figure 4. Ratio of infections and deaths from March 15 to May 15, 2021.

(a) (b)

Figure 5. Validation of the SIHR model from March 15 to May 15, 2021.

8. Conclusion

In this study, we have formulated a SIHR epidemic model for pandemic Covid-19. All the properties necessary for epidemiological relevance have also been proved. Theoretically, it is proven that the dynamics depend on the basic reproduction number to examine the stability of the system. All the properties necessary for epidemiological relevance have also been proved. We have estimated the parametric values for Djibouti using the data of the ministry of health. The maximum number of reported cases was observed on 4 April 2021, which means that the number of infections was on an upward quickly trend for the next 10 days. After that day, the number of daily reported cases was observed to decrease asymptotically. In future studies, vaccination strategies may be implemented to control the COVID-19 epidemic. Also, the extension of the compartment model can be considered.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	World Health Organization. https://covid19.who.int
[2]	Ahmed, A., Salam, B., Mohammad, M., Akgul, A. and Khoshnaw, S.H.A. (2020) Analysis Coronavirus Disease (COVID-19) Model Using Numerical Approaches and Logistic Model. AIMS Bioengineering, 7, 130-146. https://doi.org/10.3934/bioeng.2020013
[3]	Kyrychko, Y.N., Blyuss, K.B. and Brovchenko, I. (2020) Mathematical Modelling of the Dynamics and Containment of COVID-19 in Ukraine. Scientific Reports, 10, Article No. 19662. https://doi.org/10.1038/s41598-020-76710-1
[4]	Piu, S., Jayanta, M. and Subhas, K. (2020) A Mathematical Model for COVID-19 Transmission Dynamics with a Case Study of India. Chaos, Solitons and Fractals, 140, Article ID: 110173. https://doi.org/10.1016/j.chaos.2020.110173
[5]	Nguyen, H.T., Hakimeh, M. and Shahram, R. (2020) A Mathematical Model for COVID-19 Transmission by Using the Caputo Fractional Derivative. Chaos, Solitons and Fractals, 140, Article ID: 110107. https://doi.org/10.1016/j.chaos.2020.110107
[6]	Okuonghae, D. and Omame, A. (2020) A Analysis of a Mathematical Model for COVID-19 Population Dynamics in Lagos, Nigeria. Chaos, Solitons and Fractals, 139, Article ID: 110032. https://doi.org/10.1016/j.chaos.2020.110032
[7]	Foy, B.H., Wahl, B., Mehta, K., Shet, A., Menon, G.I. and Britto, C. (2020) Comparing COVID-19 Vaccine Allocation Strategies in India: A Mathematical Modelling Study. International Journal of Infectious Diseases, 103, 431-438. https://doi.org/10.1016/j.ijid.2020.12.075
[8]	Yavuz, M., Coar, F., Gnay, F. and Zdemir, F. (2021) A New Mathematical Modeling of the COVID-19 Pandemic Including the Vaccination Campaign. Open Journal of Modelling and Simulation, 9, 299-321. https://doi.org/10.4236/ojmsi.2021.93020
[9]	Kaur, S.P. and Gupta, V. (2020) COVID-19 Vaccine: A Comprehensive Status Report. Virus Research, 288, Article ID: 198114. https://doi.org/10.1016/j.virusres.2020.198114
[10]	Yokus, A. and Yavuz, M. (2021) Novel Comparison of Numerical and Analytical Methods for Fractional Burger-Fisher Equation. Discrete and Continuous Dynamical Systems, 14, 2591-2606. https://doi.org/10.3934/dcdss.2020258
[11]	Yoku, A., Durur, H. and Abro, K.A. (2021) Symbolic Computation of Caudrey-Dodd-Gibbon Equation Subject to Periodic Trigonometric and Hyperbolic Symmetries. The European Physical Journal Plus, 136, 1-16. https://doi.org/10.1140/epjp/s13360-021-01350-x
[12]	Yavuz, M. and Sene, N. (2020) Stability Analysis and Numerical Computation of the Fractional Predator-Prey Model with the Harvesting Rate. Fractal and Fractional, 4, 35. https://doi.org/10.3390/fractalfract4030035
[13]	Inc, M., Acay, B., Berhe, H.W., Yusuf, A., Khan, A. and Yao, S.W. (2021) Analysis of Novel Fractional COVID-19 Model with Real-Life Data Application. Results in Physics, 23, Article ID: 103968. https://doi.org/10.1016/j.rinp.2021.103968
[14]	Yusuf, A., Acay, B., Mustapha, U.T., Inc, M. and Baleanu, D. (2021) Mathematical Modeling of Pine Wilt Disease with Caputo Fractional Operator. Chaos, Solitons and Fractals, 143, Article ID: 110569. https://doi.org/10.1016/j.chaos.2020.110569
[15]	Naik, P.A., Yavuz, M. and Zu, J. (2020) The Role of Prostitution on HIV Transmission with Memory: A Modeling Approach. Alexandria Engineering Journal, 59, 2513-2531. https://doi.org/10.1016/j.aej.2020.04.016
[16]	Faal, N., Ivn, A., Juan, J.N. and Delfim, F.M.T. (2020) Mathematical Modeling of COVID-19 Transmission Dynamics with a Case Study of Wuhan. Chaos, Solitons and Fractals, 135, Article ID: 109846. https://doi.org/10.1016/j.chaos.2020.109846
[17]	Diekmann, O., Heesterbeek, J. and Roberts, M. (2009) The Construction of Next Generation Matrices for Compartmental Epidemic Models. Journal of the Royal Society Interface, 7, 873-885. https://doi.org/10.1098/rsif.2009.0386
[18]	Castillo-Chavez, C. and Song, B. (2004) Dynamical Model of Tuberculosis and Their Applications. Mathematical Biosciences and Engineering, 1, 361-404. https://doi.org/10.3934/mbe.2004.1.361
[19]	Martcheva, M. (2015) An Introduction to Mathematical Epidemiology. Springer, New York. https://doi.org/10.1007/978-1-4899-7612-3
[20]	Biswas, S.K., Ghosh, U. and Sarkar, S. (2020) Mathematical Model of Zika Virus Dynamics with Vector Control and Sensitivity Analysis. Infectious Disease Modelling, 5, 23-41. https://doi.org/10.1016/j.idm.2019.12.001