A New Mathematical Modeling of the COVID-19 Pandemic Including the Vaccination Campaign

Mehmet Yavuz^*, Fatma Özlem Coşar, Fatma Günay, Feyza Nur Özdemir

Department of Mathematics and Computer Sciences, Necmettin Erbakan University, Konya, Turkey.

DOI: 10.4236/ojmsi.2021.93020                

Abstract

In a short time, many illustrative studies have been conducted on the mathematical modeling and analysis of COVID-19. There are not enough studies taking into account the vaccine campaign among these studies. In this context, a mathematical model is developed to reveal the effects of vaccine treatment, which has been performed recently, on COVID-19 in this study. In the proposed model, as well as the vaccinated individuals, a five-dimensional compartment system including the susceptible, infected, exposed and recovered population is constructed. Moreover, besides the positivity, existence and uniqueness of the solution, biologically feasible region are provided. The basic reproduction number known as expected secondary infection which means that expected infection among the susceptible populations caused by this infection is evaluated. In the numerical simulations, the parameter values taken from the literature and estimated are used to perform the solutions of the proposed model. Fourth-order Runge-Kutta numerical scheme is applied to obtain the results.

Keywords

COVID-19 Model, Vaccination Campaign, Stability Analysis, Runge-Kutta Numerical Scheme

1. Introduction

In today’s world, more and more attention is being paid over time to research of epidemics such as HIV, HBV, Ebola, H1N1 and malaria, and controlling the spread of epidemics in the population is a major challenge. On the other hand, the world continues to fight existing infectious diseases, while on the other hand, changing world conditions lead to the emergence of different types of viruses. The newest of these viruses, and the most effective of recent years, is the new type of coronavirus, called COVID-19, which appeared in early 2020 and is still not fully controlled. So far the biological origin of the disease has not been fully determined, while the first cases were traced in Wuhan, China, on December 31, 2019. Lung disease and if left untreated, can cause the rapid spread of viruses that cause diseases such as severe acute respiratory syndrome, has a high rate of death and be seen throughout the world due to the World Health Organization (WHO) has been declared as a pandemic by the international.

Coronaviruses have three subgroups known as alpha, beta, and gamma. There is also a fourth new group called delta coronaviruses: SARS-CoV. Human coronaviruses were first identified in the mid-1960s. By 2020, the virus, which only appeared in Saudi Arabia, Qatar and Jordan, has killed three people. The case reported in Wuhan city of 11 million people in China’s Hubei province on December 31, 2019, was found to be infected with a new coronavirus that has never been seen before. According to who reports [1] , this virus is thought to be transmitted from animals to humans, such as SARS-CoV and MERS-CoV. Today, the disease is transmitted from person to person, on April 27 2021, the number of infected cases confirmed by who reached 147.539.302, and by this date there were 3.116.444 deaths related to the disease worldwide [2] . A study of patients who died found that the majority of them were elderly or patients diagnosed with chronic heart, lung, kidney, Parkinson’s and diabetes. Viruses that come into contact with the mouth and nose can be easily transmitted, such as flu, by removing them after touching the infected substance. In recent years, the concept of mathematical modeling has become increasingly important. Today’s information society now associates mathematics directly with everyday life, and it seems that weight is given to the embodiment of this abstract science by establishing its relationship with real life. One of the most important methods used to embody mathematics, which is an abstract science in its structure, and apply it to real life, is modeling. Mathematical modeling in the most general sense; an attempt to mathematically express an event, phenomenon, and relations between events other than mathematics or mathematics is the process of revealing mathematical patterns within these events and phenomena [3] . In fact, mathematics is a systematic way of thinking that generates solutions to real-world events and problems through modeling. When mathematics is associated with the real world, it is seen that the roots of all mathematical concepts exist and correspond in the real world.

On the other hand, the processes of creating a mathematical model of a problem, transforming it, developing it, solving it, interpreting the solution, validating the model and finally making the model ready for use are not processes that can be overcome immediately. Especially nowadays, many of the problems are complex, non-linear, have memory effect (memory effect) or have a stochastic structure, so modeling ways and solution methods specific to these problems should be developed.

While mathematical models do not provide a cure for a particular infectious disease, they can be used to present and analyze possible scenarios of the dynamics at hand. Currently, medical institutions, politicians, armies, law enforcement, industrialists, chemists, physicists, engineers, etc. sectors around the world are making efforts to stop the spread of COVID-19. Mathematicians have also developed new mathematical models that can be used for simulation to flatten the infection and death curve, predict the future behavior of its spread, and make decisions to control it, and are trying to create new optimal models every day. Many studies have been conducted in the literature on the mathematical model of the COVID-19 pandemic in a short time.

Zu et al. (2020) [4] in their data and model-oriented study proposed and examined the contagion models of COVID-19 on the Chinese mainland and the effectiveness of different control strategies. February 10 January-February 17, 2020 based on data reported by the National Health Commission of the People’s Republic of China in this study; they estimated the parameters of the model. They then successfully estimated the epidemic trend and risk of transmission of COVID-19 and predicted the effectiveness of various intervention strategies with the help of the sensitivity analysis method.

Tang et al. (2020) [5] examined the effectiveness of quarantine and isolation, which determined the trend of the COVID-19 outbreak in the last phase of the epidemic in China. In their study, they found that uncertainty analyzes revealed that the outbreak is still uncertain and concluded that it is important to continue to develop the quarantine and isolation strategy and improve detection rate in mainland China, which helped the country combat the COVID-19 outbreak.

Ahmed et al. (2020) [6] , in their study, where they analyzed a mathematical model of COVID-19 using numerical approaches and logistic models, they analyzed and introduced some COVID-19 models that contain important questions about global health services and suggest important notes. They proposed three well-known numerical techniques for solving the proposed equations; these are Euler’s method, second order (RK2) and fourth order (RK4) Runge-Kutta method. Recommended numerical techniques based on the results of the approximate solution are obtained, this is a global issue that is important to provide key answers have analyzed the results for the two countries, including Turkey and Iraq. More interestingly, these two countries, the basic reproductive number (R0), according to a survey conducted on April 9, 2020 for the estimated 7.4 Turkey is estimated to be 3.4 to Iraq. But these numbers have increased noticeably since the beginning of the pandemic. Therefore, using the model proposed by researching the logistic model to estimate the size of the outbreak of the epidemic in Turkey and Iraq in modeling and estimating results showed that a reasonable revealed.

In another study provided by Okuonghae and Omame (2020) [7] , using numerical simulations, they examined the impact of control measures, particularly common social distancing, face mask use, and case detection (through contact tracking and subsequent testing) on COVID-19 dynamics. They also made significant estimates for the cumulative values of reported cases and active cases for different levels of control measures implemented. Numerical simulations of the model suggested that if at least 55% of the population complied with the social distance regulation and about 55% of the population used face masks effectively in society, the disease would eventually disappear in the population. They also found that if the case detection rate for symptomatic individuals could be increased to about 0.8 per day with about 55% of the population complying with social distance regulations, they could achieve a large decrease in the incidence (and prevalence) of COVID-19.

Foy et al. (2020) [8] , they used an extended SEIR model, age-structured with social communication matrices, to evaluate age-specific vaccine distribution strategies in India. Their findings in this model support global recommendations for prioritizing Covid-19 vaccine delivery for older age groups. As the rate of vaccine administration increased, relative differences between allocation strategies were observed to decrease. Optimal vaccine delivery strategies; they concluded that it would depend on vaccine characteristics, the strength of simultaneous non-pharmaceutical interventions, and site-specific targets.

Farooq and Bazaz (2020) [9] , they used Deep Learning to propose an Artificial Neural Network (ANN)-based and data-stream driven real-time learning algorithm to predict the parameters of the non-competitive, intelligent adaptive and online analytical model of COVID-19 disease. Moreover, there have been a number of modeling studies related to COVID-19 and other important infectious diseases, such as [10] - [21] . Especially, modeling the HIV infection [22] , computer viruses [23] [24] [25] , pricing of different types of options [26] [27] [28] , chaotic systems [29] [30] , neuronal dynamics [31] , nerve impulses [32] , fluid flow dynamics [33] , wave problems [34] [35] , nonlinear evolution dynamics [36] can be considered as prominent illustrative modeling applications.

On the other hand, effective vaccination studies on COVID-19 continue at full speed. There are some academic studies conducted in this context. Kaur and Gupta (2020) [37] , one of them, showed in their study that different types of vaccine strategies for COVID-19 were developed and presented a review on the studies of finding an effective vaccine for this new coronavirus, which affects the world in terms of economy, human health and life.

The rest of the article is organized as follows: After the introduction given in Section 1, the formulation of the first-order mathematical model is proposed in the next Section 2. The assumptions considered in the formulation and the description of the model are given in Section 2.1. In Section 3, the existence-uniqueness of the solution of the constructed model has been investigated. In Section 4, the basic reproductive number of the model and equilibrium points are calculated and the stability conditions of the system are proposed. Mathematical investigations are provided for the solution of the model in Section 5, and numerical simulations are performed to verify the results. The behavior of the solutions obtained is also discussed in this section. Finally, Section 6 concludes all the important findings of this research study and the future direction.

2. Model Formulation

The aim of this article is to create a new COVID-19 model, taking into account the vaccination process. The model is based on the dynamics of contact with the virus, the infection of individuals who come into contact with the virus or asymptomatic recovery, and the exposure or absence of the virus at the end of the vaccination process in individuals who are not exposed to the virus.

$\begin{array}{l}\frac{\text{d}S}{\text{d}t}=\Lambda -\left(\alpha E+m+\mu \right)S,\\ \frac{\text{d}E}{\text{d}t}=\alpha SE+pVE-\left(fI+c+\mu \right)E,\\ \frac{\text{d}I}{\text{d}t}=fEI-\left(z+\mu +\sigma \right)I,\\ \frac{\text{d}V}{\text{d}t}=mS-\left(pE+\mu \right)V,\\ \frac{\text{d}R}{\text{d}t}=zI+cE-\mu R,\end{array}$ (1)

initial conditions,

$S\left(0\right)=S_0,E\left(0\right)=E_0,I\left(0\right)=I_0,V\left(0\right)=V_0,R\left(0\right)=R_0,$ (2)

where $\left(S\left(t\right)\mathrm{,}E\left(t\right)\mathrm{,}I\left(t\right)\mathrm{,}V\left(t\right)\mathrm{,}R\left(t\right)\right)\in {\mathcal{R}}_{+}^5$ . Here, it is assumed that the functions $S\left(t\right)\mathrm{,}E\left(t\right)\mathrm{,}I\left(t\right)\mathrm{,}V\left(t\right)\mathrm{,}R\left(t\right)$ and their derivatives are continuous at $t\ge 0$ . To start, non-negativity, existence and uniqueness of the solution of system (1) are analyzed.

Biological Assumptions of the Model

The biological assumptions of the model are as follows:

Ÿ We assume that individuals to be vaccinated are selected from individuals who have not been exposed to or immunized against the virus.

Ÿ In individuals in the vaccine compartment, the vaccine may not be completely protective (failed vaccine). In this case, we assume that the vaccinated individuals become infected when exposed to the virus.

COVID-19 is a type of epidemic disease that is effective today. On the other hand, different types of vaccines have been developed to end the epidemic disease in the world and are being used these days. Although the abundance of mathematical models related to COVID-19 that do not contain vaccines draw attention in the literature, there are almost no studies taking into account the vaccine parameter. In this context, the structure of the model (SEIVR) consists of the following compartments: $S\left(t\right)$ : Individuals susceptible to the disease but not susceptible to the disease (Susceptible), $E\left(t\right)$ : Individuals exposed to the disease but not showing symptoms (Exposed), $I\left(t\right)$ : Infected and individuals likely to infect susceptible individuals (Infected), $V\left(t\right)$ : Vaccinated individuals from susceptible individuals (Vaccinated), $R\left(t\right)$ : Recovered individuals (Figure 1).

Figure 1. Dynamic flow diagram of the COVID-19 model.

In the model, transition from all people ( $\Lambda$ ) to individuals sensitive to the disease is provided. The transition from susceptible individuals to individuals who are exposed to the disease as much as the rate ( $\alpha$ ), exposed, and those who have not been exposed to the disease as much as the rate (m), are passed on to individuals to be vaccinated. Considering the possibility that the vaccine developed against the COVID-19 virus may be successful or unsuccessful, in case the vaccine fails, the vaccine may be exposed to the disease as much as (p) ratio of individuals vaccinated. Individuals exposed to the virus are likely to transmit the virus to other individuals. These individuals may or may not show symptoms of the disease. Symptomatic individuals pass to the active patient part as much as (f). Individuals who recover without showing symptoms pass to the healing part as much as (c). Among active patients, deaths occur at the rate of ( $\sigma$ ) depending on the disease, and also recover at rate (z). Natural deaths occur at the rate of ( $\mu$ ) in all these compartments.

3. Mathematical Investigation of the Proposed Model

Stability analysis of a system in epidemiology and immunology determines the behavior of the system in disease dynamics. In this section, we give the positivity and boundedness of the solution for the proposed system (1). We then determine the equilibria and their stability results. Finally, the existence and uniqueness conditions for the solution are provided.

3.1. Positivity of Solutions and Determining the Biologically Invariant Region

In order for the COVID-19 outbreak model, which includes the vaccination strategy given in System (1), (2), to be epidemiologically realistic, it is necessary to show that all state variables always remain positive. In this context, this subsection addresses the positivity of the solution and limited the biologically invariant region for the proposed system (1) to the initial conditions (2).

Theorem 1. Let ${ℝ}_{+}^5=\left\{\zeta \left(t\right)\in {ℝ}^5\mathrm{<mo>\; :\; </mo>}\zeta \left(t\right)\ge 0\right\}$ and $\zeta \left(t\right)={\left[S\left(t\right)\mathrm{,}E\left(t\right)\mathrm{,}I\left(t\right)\mathrm{,}V\left(t\right)\mathrm{,}R\left(t\right)\right]}^{\text{T}}$ . The solution set $\zeta \left(t\right)$ of the proposed system (1) along initial conditions (2) is non-negative for all $t>0$ in ${ℝ}_{+}^5$ .

Proof. As proposed in the study [38] , by taking into account the nonlinear system of equations (1), we consider the first equation

$\frac{\text{d}S}{\text{d}t}=\Lambda -\left(\alpha E+m+\mu \right)S\mathrm{,}$ (3)

which means that

$\frac{\text{d}S}{\text{d}t}\ge -\left(\alpha E+m+\mu \right)S\mathrm{.}$ (4)

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

$S\left(t\right)\ge S\left(0\right){\text{e}}^{-\left(\alpha E+m+\mu \right)t}\mathrm{,}$ (5)

which implies

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

By following the similar steps with the condition in $S\left(t\right)$ , it can easily be shown that $E\left(t\right)\ge 0$ , $I\left(t\right)\ge 0$ , $V\left(t\right)\ge 0$ and $R\left(t\right)\ge 0$ for all $t>0$ . Therefore, the solution set cannot escape from the hyperplanes $S=E=I=V=R=0$ . Moreover, on each hyperplane the vector field falls into ${ℝ}_{+}^5$ which means that the feasible region ${ℝ}_{+}^5$ is positively invariant.

In the next theorem, we provide a biological feasible region for the proposed model (1).

Theorem 2. The solutions of system (1) with the initial condition (2) are stated in the region $\mathcal{B}\subset {ℝ}_{+}^5$ , given by

$\mathcal{B}=\left\{\left(S\left(t\right)\mathrm{,}E\left(t\right)\mathrm{,}I\left(t\right)\mathrm{,}V\left(t\right)\mathrm{,}R\left(t\right)\right)\in {\mathcal{R}}_{+}^5\mathrm{<mo>\; |\; </mo>}N\left(t\right)\le \frac{\Lambda }{\mu }\right\}\mathrm{.}$ (7)

Proof. Considering the summation of all equations in the system yields

$\frac{\text{d}N\left(t\right)}{\text{d}t}=\frac{\text{d}S\left(t\right)}{\text{d}t}+\frac{\text{d}E\left(t\right)}{\text{d}t}+\frac{\text{d}I\left(t\right)}{\text{d}t}+\frac{\text{d}V\left(t\right)}{\text{d}t}+\frac{\text{d}R\left(t\right)}{\text{d}t}.$ (8)

Then we have the following for the whole population

$\frac{\text{d}N}{\text{d}t}=\Lambda -\mu N-\sigma I\le \Lambda -\mu N.$ (9)

The solution of Equation (9) is given as

$N\left(t\right)\le \frac{\Lambda }{\mu }-\left(\frac{\Lambda }{\mu }-N_0\right){\text{e}}^{-\mu t}\mathrm{,}$ (10)

where the initial population is defined as $N_0=N\left(0\right)$ . Benefiting the Birkhoff-Rota theorem [39] , we can say that, if $N_0<\frac{\Lambda }{\mu }$ , then as $t\to \infty$ , asymptotically $N\left(t\right)\to \frac{\Lambda }{\mu }$ in Equation (7) and the total population size becomes $N\left(t\right)\to \frac{\Lambda }{\mu }$ which means that $0\le N\le \frac{\Lambda }{\mu }$ . Thus, all the feasible solutions of the model converge in the region $\mathcal{B}$ [40] .

3.2. Existence and Uniqueness

In this part, we give the qualitative properties related to the solutions of the COVID-19 model which is given in system (1) - (2). Firstly, we start by taking the classical integral of both sides the mentioned system and we get the following Volterra-type integral equations:

$\begin{array}{l}S\left(t\right)-S\left(0\right)={\displaystyle {\int }_0^t}\left[\Lambda -\left(\alpha E\left(\tau \right)+m+\mu \right)S\left(\tau \right)\right]\text{d}\tau ,\\ E\left(t\right)-E\left(0\right)={\displaystyle {\int }_0^t}\left[E\left(\tau \right)\left(\alpha S+pV-fI-c-\mu \right)\right]\text{d}\tau ,\\ I\left(t\right)-I\left(0\right)={\displaystyle {\int }_0^t}\left[fE\left(\tau \right)I\left(\tau \right)-I\left(\tau \right)\left(z+\mu +\sigma \right)\right]\text{d}\tau ,\\ V\left(t\right)-V\left(0\right)={\displaystyle {\int }_0^t}\left[mS\left(\tau \right)-pV\left(\tau \right)E\left(\tau \right)-\mu V\left(\tau \right)\right]\text{d}\tau ,\\ R\left(t\right)-R\left(0\right)={\displaystyle {\int }_0^t}\left[zI\left(\tau \right)+cE\left(\tau \right)-\mu R\left(\tau \right)\right]\text{d}\tau .\end{array}$ (11)

Below, let’s define the kernels as follows:

$\begin{array}{l}{\phi }_1\left(t,S\right)=\Lambda -\left(\alpha E\left(t\right)+m+\mu \right)S\left(t\right),\\ {\phi }_2\left(t,E\right)=E\left(t\right)\left(\alpha S+pV-fI-c-\mu \right),\\ {\phi }_3\left(t,I\right)=fE\left(t\right)I\left(t\right)-\left(z+\mu +\sigma \right)I\left(t\right),\\ {\phi }_4\left(t,V\right)=mS\left(t\right)-pV\left(t\right)E\left(t\right)-\mu V\left(t\right),\\ {\phi }_5\left(t,R\right)=zI\left(t\right)+cE\left(t\right)-\mu R\left(t\right).\end{array}$ (12)

Then the following theorem arises:

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

$0\le q_1,q_2,q_3,q_4,q_5<1.$ (13)

where, $\Vert S\Vert \le m_1$ , $\Vert E\Vert \le m_2$ , $\Vert I\Vert \le m_3$ , $\Vert V\Vert \le m_4$ , $\Vert R\Vert \le m_5$ , $q_1=\alpha m_2+m+\mu$ , $q_2=\alpha m_1+p{m}_4+f{m}_3+c+\mu$ , $q_3=f{m}_2+z+\mu +\sigma$ , $q_4=p{m}_2+\mu$ and $q_5=\mu$ .

Proof. Let $S_1$ and $S_2$ be two functions for the kernel ${\phi }_1$ ; $E_1$ and $E_2$ be two functions for the kernel ${\phi }_2$ ; $I_1$ and $I_2$ be two functions for the kernel ${\phi }_3$ ; $V_1$ and $V_2$ be two functions for the kernel ${\phi }_4$ ; and $R_1$ and $R_2$ be two functions for the kernel ${\phi }_5$ . Then we have

$\begin{array}{l}\Vert {\phi }_1\left(t,S_1\right)-{\phi }_1\left(t,S_2\right)\Vert \\ =\Vert \Lambda -\left(\alpha E+m+\mu \right)S_1-\left[\Lambda -\left(\alpha E+m+\mu \right)S_2\right]\Vert \\ \le \left(\alpha m_2+m+\mu \right)\Vert S_1-S_2\Vert \\ =q_1\Vert S_1-S_2\Vert ,\end{array}$

$\begin{array}{l}\Vert {\phi }_2\left(t,E_1\right)-{\phi }_2\left(t,E_2\right)\Vert \\ =\Vert E_1\left(\alpha S+pV-fI-c-\mu \right)-E_2\left(\alpha S+pV-fI-c-\mu \right)\Vert \\ \le \left(\alpha m_1+p{m}_4+f{m}_3+c+\mu \right)\Vert E_1-E_2\Vert \\ =q_2\Vert E_1-E_2\Vert ,\end{array}$

$\begin{array}{l}\Vert {\phi }_3\left(t,I_1\right)-{\phi }_3\left(t,I_2\right)\Vert \\ =\Vert fE{I}_1-\left(z+\mu +\sigma \right)I_1-\left[fE{I}_2-\left(z+\mu +\sigma \right)I_2\right]\Vert \\ \le \left(f{m}_2+z+\mu +\sigma \right)\Vert I_1-I_2\Vert \\ =q_3\Vert I_1-I_2\Vert ,\end{array}$ (14)

$\begin{array}{l}\Vert {\phi }_4\left(t,V_1\right)-{\phi }_4\left(t,V_2\right)\Vert \\ =\Vert mS-pE{V}_1-\mu V_1-\left[mS-pE{V}_2-\mu V_2\right]\Vert \\ \le \left(p{m}_2+\mu \right)\Vert V_1-V_2\Vert \\ =q_4\Vert V_1-V_2\Vert \end{array}$

and

$\begin{array}{l}\Vert {\phi }_5\left(t,R_1\right)-{\phi }_5\left(t,R_2\right)\Vert \\ =\Vert zI+cE-\mu R_1-\left[zI+cE-\mu R_2\right]\Vert \\ \le \mu \Vert R_1-R_2\Vert \\ =q_5\Vert R_1-R_2\Vert .\end{array}$

Therefore, the Lipschitz conditions are satisfied for kernels ${\phi }_1,{\phi }_2,{\phi }_3,{\phi }_4,{\phi }_5$ and if $0\le q_1,q_2,q_3,q_4,q_5<1$ then $q_1,q_2,q_3,q_4,q_5$ are also contractions for ${\phi }_1,{\phi }_2,{\phi }_3,{\phi }_4,{\phi }_5$ respectively. This proves the theorem.

By considering the kernels ${\phi }_1,{\phi }_2,{\phi }_3,{\phi }_4$ and ${\phi }_5$ , we can rewrite the system which is given in Equation (11) as follows:

$\begin{array}{l}S\left(t\right)=S\left(0\right)+{\displaystyle {\int }_0^t}{\phi }_1\left(\tau ,S\right)\text{d}\tau ,\\ E\left(t\right)=E\left(0\right)+{\displaystyle {\int }_0^t}{\phi }_2\left(\tau ,E\right)\text{d}\tau ,\\ I\left(t\right)=I\left(0\right)+{\displaystyle {\int }_0^t}{\phi }_3\left(\tau ,I\right)\text{d}\tau ,\\ V\left(t\right)=V\left(0\right)+{\displaystyle {\int }_0^t}{\phi }_4\left(\tau ,V\right)\text{d}\tau ,\\ R\left(t\right)=R\left(0\right)+{\displaystyle {\int }_0^t}{\phi }_5\left(\tau ,R\right)\text{d}\tau .\end{array}$ (15)

We can proceed with the following recursive formula

$\begin{array}{l}{S}_n\left(t\right)=S\left(0\right)+{\displaystyle {\int }_0^t}{\phi }_1\left(\tau ,S_{n-1}\right)\text{d}\tau ,\\ E_n\left(t\right)=E\left(0\right)+{\displaystyle {\int }_0^t}{\phi }_2\left(\tau ,E_{n-1}\right)\text{d}\tau ,\\ I_n\left(t\right)=I\left(0\right)+{\displaystyle {\int }_0^t}{\phi }_3\left(\tau ,I_{n-1}\right)\text{d}\tau ,\\ V_n\left(t\right)=V\left(0\right)+{\displaystyle {\int }_0^t}{\phi }_4\left(\tau ,V_{n-1}\right)\text{d}\tau ,\\ R_n\left(t\right)=R\left(0\right)+{\displaystyle {\int }_0^t}{\phi }_5\left(\tau ,R_{n-1}\right)\text{d}\tau ,\end{array}$ (16)

where $S_0\left(t\right)=S\left(0\right)$ , $E_0\left(t\right)=E\left(0\right)$ , $I_0\left(t\right)=I\left(0\right)$ , $V_0\left(t\right)=V\left(0\right)$ and $R_0\left(t\right)=R\left(0\right)$ . Then we can write

$\begin{array}{l}{\Psi }_n\left(t\right)=S_n\left(t\right)-S_{n-1}\left(t\right)={\displaystyle {\int }_0^t}\left[{\phi }_1\left(\tau ,S_{n-1}\right)-{\phi }_1\left(\tau ,S_{n-2}\right)\right]\text{d}\tau ,\\ {\Delta }_n\left(t\right)=E_n\left(t\right)-E_{n-1}\left(t\right)={\displaystyle {\int }_0^t}\left[{\phi }_2\left(\tau ,E_{n-1}\right)-{\phi }_2\left(\tau ,E_{n-2}\right)\right]\text{d}\tau ,\\ {\Phi }_n\left(t\right)=I_n\left(t\right)-I_{n-1}\left(t\right)={\displaystyle {\int }_0^t}\left[{\phi }_3\left(\tau ,I_{n-1}\right)-{\phi }_3\left(\tau ,I_{n-2}\right)\right]\text{d}\tau ,\\ {\Upsilon }_n\left(t\right)=V_n\left(t\right)-V_{n-1}\left(t\right)={\displaystyle {\int }_0^t}\left[{\phi }_4\left(\tau ,V_{n-1}\right)-{\phi }_4\left(\tau ,V_{n-2}\right)\right]\text{d}\tau ,\\ {\partial }_n\left(t\right)=R_n\left(t\right)-R_{n-1}\left(t\right)={\displaystyle {\int }_0^t}\left[{\phi }_5\left(\tau ,R_{n-1}\right)-{\phi }_5\left(\tau ,R_{n-2}\right)\right]\text{d}\tau ,\end{array}$ (17)

where $S_n\left(t\right)={\displaystyle {\sum }_{j=1}^n}{\Psi }_n\left(t\right)$ , $E_n\left(t\right)={\displaystyle {\sum }_{j=1}^n}{\Delta }_n\left(t\right)$ , $I_n\left(t\right)={\displaystyle {\sum }_{j=1}^n}{\Phi }_n\left(t\right)$ , $V_n\left(t\right)={\displaystyle {\sum }_{j=1}^n}{\Upsilon }_n\left(t\right)$ and $R_n\left(t\right)={\displaystyle {\sum }_{j=1}^n}{\partial }_n\left(t\right)$ . By taking the norm of both sides of Equation (17), we have

$\begin{array}{l}\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(\tau ,S_{n-1}\right)-{\phi }_1\left(\tau ,S_{n-2}\right)\right]\text{d}\tau \Vert ,\\ \Vert {\Delta }_n\left(t\right)\Vert =\Vert E_n\left(t\right)-E_{n-1}\left(t\right)\Vert \le \Vert {\displaystyle {\int }_0^t}\left[{\phi }_2\left(\tau ,E_{n-1}\right)-{\phi }_2\left(\tau ,E_{n-2}\right)\right]\text{d}\tau \Vert ,\\ \Vert {\Phi }_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 }_3\left(\tau ,I_{n-1}\right)-{\phi }_3\left(\tau ,I_{n-2}\right)\right]\text{d}\tau \Vert ,\\ \Vert {\Upsilon }_n\left(t\right)\Vert =\Vert V_n\left(t\right)-V_{n-1}\left(t\right)\Vert \le \Vert {\displaystyle {\int }_0^t}\left[{\phi }_4\left(\tau ,V_{n-1}\right)-{\phi }_4\left(\tau ,V_{n-2}\right)\right]\text{d}\tau \Vert ,\\ \Vert {\partial }_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 }_5\left(\tau ,R_{n-1}\right)-{\phi }_5\left(\tau ,R_{n-2}\right)\right]\text{d}\tau \Vert .\end{array}$ (18)

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

$\begin{array}{l}\Vert S_n\left(t\right)-S_{n-1}\left(t\right)\Vert \le q_1{\displaystyle {\int }_0^t}\Vert S_{n-1}-S_{n-2}\Vert \text{d}\tau \mathrm{,}\\ \Vert E_n\left(t\right)-E_{n-1}\left(t\right)\Vert \le q_2{\displaystyle {\int }_0^t}\Vert E_{n-1}-E_{n-2}\Vert \text{d}\tau \mathrm{,}\\ \Vert I_n\left(t\right)-I_{n-1}\left(t\right)\Vert \le q_3{\displaystyle {\int }_0^t}\Vert I_{n-1}-I_{n-2}\Vert \text{d}\tau \mathrm{,}\\ \Vert V_n\left(t\right)-V_{n-1}\left(t\right)\Vert \le q_4{\displaystyle {\int }_0^t}\Vert V_{n-1}-V_{n-2}\Vert \text{d}\tau \mathrm{,}\\ \Vert R_n\left(t\right)-R_{n-1}\left(t\right)\Vert \le q_5{\displaystyle {\int }_0^t}\Vert R_{n-1}-R_{n-2}\Vert \text{d}\tau \mathrm{.}\end{array}$ (19)

Then we achieve from the last inequality

$\begin{array}{l}\Vert {\Psi }_n\left(t\right)\Vert \le q_1{\displaystyle {\int }_0^t}\Vert {\Psi }_{n-1}\left(\tau \right)\Vert \text{d}\tau \mathrm{,}\\ \Vert {\Delta }_n\left(t\right)\Vert \le q_2{\displaystyle {\int }_0^t}\Vert {\Delta }_{n-1}\left(\tau \right)\Vert \text{d}\tau \mathrm{,}\\ \Vert {\Phi }_n\left(t\right)\Vert \le q_3{\displaystyle {\int }_0^t}\Vert {\Phi }_{n-1}\left(\tau \right)\Vert \text{d}\tau \mathrm{,}\\ \Vert {\Upsilon }_n\left(t\right)\Vert \le q_4{\displaystyle {\int }_0^t}\Vert {\Upsilon }_{n-1}\left(\tau \right)\Vert \text{d}\tau \mathrm{,}\\ \Vert {\partial }_n\left(t\right)\Vert \le q_5{\displaystyle {\int }_0^t}\Vert {\partial }_{n-1}\left(\tau \right)\Vert \text{d}\tau \mathrm{.}\end{array}$ (20)

These results give us the following theorem.

Theorem 4. The proposed COVID-19 model has a solution under the condition that can be founded $t_{\max }$ holding:

$q_i{t}_{\max }<1,i=1,2,3,4,5.$ (21)

Proof. Considering the functions $S\left(t\right)\mathrm{,}E\left(t\right)\mathrm{,}I\left(t\right)\mathrm{,}V\left(t\right)$ and $R\left(t\right)$ are bounded and their kernels ${\phi }_1,{\phi }_2,{\phi }_3,{\phi }_4$ and ${\phi }_5$ hold the Lipschitz condition, we can give the following by taking Equation (20) into account

$\begin{array}{l}\Vert {\Psi }_n\left(t\right)\Vert \le \Vert S_0\left(t\right)\Vert {\left\{q_1{t}_{\max }\right\}}^n,\\ \Vert {\Delta }_n\left(t\right)\Vert \le \Vert E_0\left(t\right)\Vert {\left\{q_2{t}_{\max }\right\}}^n,\\ \Vert {\Phi }_n\left(t\right)\Vert \le \Vert I_0\left(t\right)\Vert {\left\{q_3{t}_{\max }\right\}}^n,\\ \Vert {\Upsilon }_n\left(t\right)\Vert \le \Vert V_0\left(t\right)\Vert {\left\{q_4{t}_{\max }\right\}}^n,\\ \Vert {\partial }_n\left(t\right)\Vert \le \Vert R_0\left(t\right)\Vert {\left\{q_5{t}_{\max }\right\}}^n.\end{array}$ (22)

Now we show that the functions in Equation (22) are the solutions of the given COVID-19 model. We suppose

$\begin{array}{l}S\left(t\right)-S\left(0\right)=S_n\left(t\right)-{\omega }_n\left(t\right),\\ E\left(t\right)-E\left(0\right)=E_n\left(t\right)-{\vartheta }_n\left(t\right),\\ I\left(t\right)-I\left(0\right)=I_n\left(t\right)-{\varkappa }_n\left(t\right),\\ V\left(t\right)-V\left(0\right)=V_n\left(t\right)-{\rho }_n\left(t\right),\\ R\left(t\right)-R\left(0\right)=R_n\left(t\right)-{\kappa }_n\left(t\right).\end{array}$ (23)

Then we will demonstrate that the terms which are given in Equation (23) hold that $\Vert {\omega }_{\infty }\left(t\right)\Vert \to 0$ , $\Vert {\vartheta }_{\infty }\left(t\right)\Vert \to 0$ , $\Vert {\varkappa }_{\infty }\left(t\right)\Vert \to 0$ , $\Vert {\rho }_{\infty }\left(t\right)\Vert \to 0$ and $\Vert {\kappa }_{\infty }\left(t\right)\Vert \to 0$ . Since we have

$\begin{array}{c}\Vert {\omega }_n\left(t\right)\Vert \le \Vert {\displaystyle {\int }_0^t}\left[{\phi }_1\left(\tau ,S\right)-{\phi }_1\left(\tau ,S_{n-1}\right)\right]\text{d}\tau \Vert \\ \le {\displaystyle {\int }_0^t}\Vert {\phi }_1\left(\tau ,S\right)-{\phi }_1\left(\tau ,S_{n-1}\right)\Vert \text{d}\tau \\ \le t{q}_1\Vert S-S_{n-1}\Vert ,\end{array}$

$\begin{array}{c}\Vert {\vartheta }_n\left(t\right)\Vert \le \Vert {\displaystyle {\int }_0^t}\left[{\phi }_2\left(\tau \mathrm{,}E\right)-{\phi }_2\left(\tau \mathrm{,}{E}_{n-1}\right)\right]\text{d}\tau \Vert \\ \le {\displaystyle {\int }_0^t}\Vert {\phi }_2\left(\tau \mathrm{,}E\right)-{\phi }_2\left(\tau \mathrm{,}{E}_{n-1}\right)\Vert \text{d}\tau \\ \le t{q}_2\Vert E-E_{n-1}\Vert \mathrm{,}\end{array}$

$\begin{array}{c}\Vert {\varkappa }_n\left(t\right)\Vert \le \Vert {\displaystyle {\int }_0^t}\left[{\phi }_3\left(\tau \mathrm{,}I\right)-{\phi }_3\left(\tau \mathrm{,}{I}_{n-1}\right)\right]\text{d}\tau \Vert \\ \le {\displaystyle {\int }_0^t}\Vert {\phi }_3\left(\tau \mathrm{,}I\right)-{\phi }_3\left(\tau \mathrm{,}{I}_{n-1}\right)\Vert \text{d}\tau \\ \le t{q}_3\Vert I-I_{n-1}\Vert \mathrm{,}\end{array}$ (24)

$\begin{array}{c}\Vert {\rho }_n\left(t\right)\Vert \le \Vert {\displaystyle {\int }_0^t}\left[{\phi }_4\left(\tau \mathrm{,}V\right)-{\phi }_4\left(\tau \mathrm{,}{V}_{n-1}\right)\right]\text{d}\tau \Vert \\ \le {\displaystyle {\int }_0^t}\Vert {\phi }_4\left(\tau \mathrm{,}V\right)-{\phi }_4\left(\tau \mathrm{,}{V}_{n-1}\right)\Vert \text{d}\tau \\ \le t{q}_4\Vert V-V_{n-1}\Vert \end{array}$

and

$\begin{array}{c}\Vert {\kappa }_n\left(t\right)\Vert \le \Vert {\displaystyle {\int }_0^t}\left[{\phi }_5\left(\tau \mathrm{,}R\right)-{\phi }_5\left(\tau \mathrm{,}{R}_{n-1}\right)\right]\text{d}\tau \Vert \\ \le {\displaystyle {\int }_0^t}\Vert {\phi }_5\left(\tau \mathrm{,}R\right)-{\phi }_5\left(\tau \mathrm{,}{R}_{n-1}\right)\Vert \text{d}\tau \\ \le t{q}_5\Vert R-R_{n-1}\Vert \mathrm{,}\end{array}$

repeating this process recursively, we get

$\Vert {\omega }_n\left(t\right)\Vert \le {\left\{t\right\}}^{n+1}{q}_1^n\Theta \mathrm{,}$

$\Vert {\vartheta }_n\left(t\right)\Vert \le {\left\{t\right\}}^{n+1}{q}_2^n\Theta \mathrm{,}$

$\Vert {\varkappa }_n\left(t\right)\Vert \le {\left\{t\right\}}^{n+1}{q}_3^n\Theta \mathrm{,}$

$\Vert {\rho }_n\left(t\right)\Vert \le {\left\{t\right\}}^{n+1}{q}_4^n\Theta$

and

$\Vert {\kappa }_n\left(t\right)\Vert \le {\left\{t\right\}}^{n+1}{q}_5^n\Theta \mathrm{.}$

Taking the last inequalities into account at $t_{\max }$ point, we get

$\Vert {\omega }_n\left(t\right)\Vert \le {\left\{t_{\max }\right\}}^{n+1}{q}_1^n\Theta ,$

$\Vert {\vartheta }_n\left(t\right)\Vert \le {\left\{t_{\max }\right\}}^{n+1}{q}_2^n\Theta ,$

$\Vert {\varkappa }_n\left(t\right)\Vert \le {\left\{t_{\max }\right\}}^{n+1}{q}_3^n\Theta \mathrm{,}$

$\Vert {\rho }_n\left(t\right)\Vert \le {\left\{t_{\max }\right\}}^{n+1}{q}_4^n\Theta ,$

and

$\Vert {\kappa }_n\left(t\right)\Vert \le {\left\{t_{\max }\right\}}^{n+1}{q}_5^n\Theta .$

For the last step, after applying the limit to both sides of the last inequalities as $n\to \infty$ , and by considering the results of Theorem 3, we have $\Vert {\omega }_{\infty }\left(t\right)\Vert \to 0$ , $\Vert {\vartheta }_{\infty }\left(t\right)\Vert \to 0$ , $\Vert {\varkappa }_{\infty }\left(t\right)\Vert \to 0$ , $\Vert {\rho }_{\infty }\left(t\right)\Vert \to 0$ and $\Vert {\kappa }_{\infty }\left(t\right)\Vert \to 0$ .

Theorem 5. The suggested COVID-19 model has only one solution.

Proof. Let us say there exists another solution of the system namely $S_1\left(t\right)$ , $E_1\left(t\right)$ , $I_1\left(t\right)$ , $V_1\left(t\right)$ and $R_1\left(t\right)$ . Then we can write

$\begin{array}{l}S\left(t\right)-S_1\left(t\right)={\displaystyle {\int }_0^t}\left[{\phi }_1\left(\tau ,S\right)-{\phi }_1\left(\tau ,S_1\right)\right]\text{d}\tau ,\\ E\left(t\right)-E_1\left(t\right)={\displaystyle {\int }_0^t}\left[{\phi }_2\left(\tau ,E\right)-{\phi }_2\left(\tau ,E_1\right)\right]\text{d}\tau ,\\ I\left(t\right)-I_1\left(t\right)={\displaystyle {\int }_0^t}\left[{\phi }_3\left(\tau ,I\right)-{\phi }_3\left(\tau ,I_1\right)\right]\text{d}\tau ,\\ V\left(t\right)-V_1\left(t\right)={\displaystyle {\int }_0^t}\left[{\phi }_4\left(\tau ,V\right)-{\phi }_4\left(\tau ,V_1\right)\right]\text{d}\tau ,\\ R\left(t\right)-R_1\left(t\right)={\displaystyle {\int }_0^t}\left[{\phi }_5\left(\tau ,R\right)-{\phi }_5\left(\tau ,R_1\right)\right]\text{d}\tau .\end{array}$ (25)

If we apply the norm to both sides of Equation (25), we obtain

$\begin{array}{l}\Vert S\left(t\right)-S_1\left(t\right)\Vert \le {\displaystyle {\int }_0^t}\Vert {\phi }_1\left(\tau \mathrm{,}S\right)-{\phi }_1\left(\tau \mathrm{,}{S}_1\right)\Vert \text{d}\tau \mathrm{,}\\ \Vert E\left(t\right)-E_1\left(t\right)\Vert \le {\displaystyle {\int }_0^t}\Vert {\phi }_2\left(\tau \mathrm{,}E\right)-{\phi }_2\left(\tau \mathrm{,}{E}_1\right)\Vert \text{d}\tau \mathrm{,}\\ \Vert I\left(t\right)-I_1\left(t\right)\Vert \le {\displaystyle {\int }_0^t}\Vert {\phi }_3\left(\tau \mathrm{,}I\right)-{\phi }_3\left(\tau \mathrm{,}{I}_1\right)\Vert \text{d}\tau \mathrm{,}\\ \Vert V\left(t\right)-V_1\left(t\right)\Vert \le {\displaystyle {\int }_0^t}\Vert {\phi }_4\left(\tau \mathrm{,}V\right)-{\phi }_4\left(\tau \mathrm{,}{V}_1\right)\Vert \text{d}\tau \mathrm{,}\\ \Vert R\left(t\right)-R_1\left(t\right)\Vert \le {\displaystyle {\int }_0^t}\Vert {\phi }_5\left(\tau \mathrm{,}R\right)-{\phi }_5\left(\tau \mathrm{,}{R}_1\right)\Vert \text{d}\tau \mathrm{.}\end{array}$ (26)

Since the Lipschitz condition is satisfied by the kernels ${\phi }_1,{\phi }_2,{\phi }_3,{\phi }_4$ and ${\phi }_5$ , we can write

$\begin{array}{l}\Vert S\left(t\right)-S_1\left(t\right)\Vert \le q_{1}t\Vert S\left(t\right)-S_1\left(t\right)\Vert \mathrm{,}\\ \Vert E\left(t\right)-E_1\left(t\right)\Vert \le q_{2}t\Vert E\left(t\right)-E_1\left(t\right)\Vert \mathrm{,}\\ \Vert I\left(t\right)-I_1\left(t\right)\Vert \le q_{3}t\Vert I\left(t\right)-I_1\left(t\right)\Vert \mathrm{,}\\ \Vert V\left(t\right)-V_1\left(t\right)\Vert \le q_{4}t\Vert V\left(t\right)-V_1\left(t\right)\Vert \mathrm{,}\\ \Vert R\left(t\right)-R_1\left(t\right)\Vert \le q_{5}t\Vert R\left(t\right)-R_1\left(t\right)\Vert \mathrm{,}\end{array}$ (27)

which gives

$\begin{array}{l}\Vert S\left(t\right)-S_1\left(t\right)\Vert \left(1-q_{1}t\right)\le \mathrm{0,}\\ \Vert E\left(t\right)-E_1\left(t\right)\Vert \left(1-q_{2}t\right)\le \mathrm{0,}\\ \Vert I\left(t\right)-I_1\left(t\right)\Vert \left(1-q_{3}t\right)\le \mathrm{0,}\\ \Vert V\left(t\right)-V_1\left(t\right)\Vert \left(1-q_{4}t\right)\le \mathrm{0,}\\ \Vert R\left(t\right)-R_1\left(t\right)\Vert \left(1-q_{5}t\right)\le 0.\end{array}$ (28)

Therefore, we obtain $\Vert S\left(t\right)-S_1\left(t\right)\Vert =0$ , $\Vert E\left(t\right)-E_1\left(t\right)\Vert =0$ , $\Vert I\left(t\right)-I_1\left(t\right)\Vert =0$ , $\Vert V\left(t\right)-V_1\left(t\right)\Vert =0$ and $\Vert R\left(t\right)-R_1\left(t\right)\Vert =0$ which represents $S\left(t\right)=S_1\left(t\right)$ , $E\left(t\right)=E_1\left(t\right)$ , $I\left(t\right)=I_1\left(t\right)$ , $V\left(t\right)=V_1\left(t\right)$ and $R\left(t\right)=R_1\left(t\right)$ This explains that the model has a unique solution which is the proof of the theorem.

4. Equilibria, Their Stabilities and Basic Reproduction Number

In order to have the equilibrium points of the system (1), we take

$\{\begin{array}{l}\Lambda -\left(\alpha E+m+\mu \right)S=0,\\ \alpha SE+pVE-\left(fI+c+\mu \right)E=0,\\ fEI-\left(z+\mu +\sigma \right)I=0,\\ mS-\left(pE+\mu \right)V=0,\\ zI+cE-\mu R=0.\end{array}$ (29)

The solution of system Eq. (29) together yields five steady states. We give these equilibria and discuss their local behavior according to their biological relevance. The first equilibrium point is disease-free equilibrium (DFE) given as

${\Omega }_1=\left(\frac{\Lambda }{\mu +m}\mathrm{,0,0,}\frac{\Lambda m}{\mu \left(\mu +m\right)}\mathrm{,0}\right)$, which means none of the cell populationexists. The second equilibrium point is the endemic equilibrium given as

$\begin{array}{l}{\Omega }_2=(\frac{f\Lambda }{\alpha \mu +\alpha \sigma +f\mu +fm+\alpha z}\mathrm{,}\frac{\mu +\sigma +z}{f}\mathrm{,}\\ \frac{\Lambda \left(f\left(\alpha \mu +mp\right)+\alpha p\left(\mu +\sigma +z\right)\right)}{\left(f\mu +p\left(\mu +\sigma +z\right)\right)\left(f\left(\mu +m\right)+\alpha \left(\mu +\sigma +z\right)\right)}-\frac{c+\mu }{f}\mathrm{,}\\ \frac{f^2\Lambda m}{\left(f\mu +\mu p+p\sigma +pz\right)\left(\alpha \mu +\alpha \sigma +f\mu +fm+\alpha z\right)}\mathrm{,}\end{array}$

$\begin{array}{l}\frac{f^2\left(m\left(\mu \left(c\left(\mu +\sigma \right)-\mu z\right)+\Lambda pz\right)+\mu \left(c\mu \mathrm{<mo\; stretchy="false">\; (\; </mo>}\mu +\sigma \mathrm{<mo\; stretchy="false">\; )\; </mo>}+z\left(\alpha \Lambda -{\mu }^2\right)\right)\right)}{f\mu \left(f\mu +p\left(\mu +\sigma +z\right)\right)\left(f\left(\mu +m\right)+\alpha \left(\mu +\sigma +z\right)\right)}\\ +\frac{f\left(\mu +\sigma +z\right)p\left(c\left(\mu +\sigma \right)\left(\mu +m\right)+z\left(\alpha \Lambda -\mu \left(\mu +m\right)\right)\right)}{f\mu \left(f\mu +p\left(\mu +\sigma +z\right)\right)\left(f\left(\mu +m\right)+\alpha \left(\mu +\sigma +z\right)\right)}\\ +\frac{\alpha \mu \left(c\left(\mu +\sigma \right)-\mu z\right)+\alpha p{\left(\mu +\sigma +z\right)}^2\left(c\left(\mu +\sigma \right)-\mu z\right)}{f\mu \left(f\mu +p\left(\mu +\sigma +z\right)\right)\left(f\left(\mu +m\right)+\alpha \left(\mu +\sigma +z\right)\right)})\mathrm{.}\end{array}$

Reproduction Number

One of the commonly known scientific tests in epidemics is the estimated number of cases, known by names such as the number of basic reproduction number or basic reproduction ratio, and the concept with ${\mathcal{R}}_0$ that shows how many people on average can be infected with the virus. The model we consider here is assuming a heterogeneous population group whose non-homogeneous individuals can be grouped. In order to evaluate the basic reproduction number of the system (1) by using the next generation matrix method given in , firstly we give the solution set as

$\zeta \left(t\right)={\left[S\left(t\right)\mathrm{,}E\left(t\right)\mathrm{,}I\left(t\right)\mathrm{,}V\left(t\right)\mathrm{,}R\left(t\right)\right]}^{\text{T}}\mathrm{.}$ (30)

Then, it is essential to distinguish new infections from all other cases in population. Therefore, we define the solution set given by Equation (30) as the difference two matrices:

$\zeta \left(t\right)=\mathcal{F}\left(\zeta \right)-\mathcal{V}\left(\zeta \right)\mathrm{,}$

where

$\mathcal{F}\left(\zeta \right)=\left(\begin{array}{c}0\\ \alpha SE+pVE\\ fEI\\ 0\\ 0\end{array}\right)\mathrm{,}\mathcal{V}\left(\zeta \right)=\left(\begin{array}{c}S\left(\alpha E+\mu +m\right)\\ E\left(fI+c+\mu \right)\\ I\left(z+\mu +\sigma \right)\\ -mS+V\left(pE+\mu \right)\\ -zI-cE-\mu R\end{array}\right)\mathrm{.}$

In what follows, F and V matrices of uninfected division at the equilibrium point ${\Omega }_1$ are defined by the matrix generation method as $F=\left[\frac{\partial {\mathcal{F}}_i\left({\Omega }_1\right)}{\partial t_j}\right]$ and $V=\left[\frac{\partial {\mathcal{V}}_i\left({\Omega }_1\right)}{\partial t_j}\right]$ , $1\le i\mathrm{,}j\le 2$ . As more explicitly the can be written as the following forms:

$F=\left(\begin{array}{cc}\frac{\alpha \Lambda }{m+\mu }+\frac{pm\Lambda }{\mu \left(m+\mu \right)}& 0\\ 0& 0\end{array}\right)\mathrm{,}V=\left(\begin{array}{cc}c+\mu & 0\\ 0& z+\mu +\sigma \end{array}\right)\mathrm{,}$

where the matrix F is non-negative, and the matrix V is non-singular. The basic reproduction number of the disease is calculated with the spectral radius of the matrix $F{V}^{-1}$ at the equilibrium point ${\Omega }_1$ , basic reproductive coefficient denoted by:

${\mathcal{R}}_0=\frac{\Lambda \left(\alpha \mu +mp\right)}{\mu \left(c+\mu \right)\left(m+\mu \right)}\mathrm{.}$ (31)

If ${\mathcal{R}}_0<1$ , on average, an infected individual, less than one newly infected individual occurs during the contagion period, and the rate of spread of the infection decreases. This shows that the infection will gradually disappear and the epidemic will stop. If ${\mathcal{R}}_0>1$ , on average each infected individual will have more than one newly infected individuals and the rate of spread of infection increases. This means that the epidemic will not end, maintaining the presence of the disease. ${\mathcal{R}}_0=1$ means that each patient with an infection leads to an infection on average. In this case, the disease will maintain its existence, but it will neither turn into an epidemic nor end. ${\mathcal{R}}_0$ basic reproduction number of reproduction indicates that the systems determine their movements.

Theorem 6. The disease-free equilibrium point ${\Omega }_1$ of the epidemic model is locally asymptotically stable if ${\mathcal{R}}_0<1$ , otherwise unstable.

Proof 6. In order to point out the stability conditions at the DFE point given as ${\Omega }_1$ let us consider the following Jakobian matrix:

$\mathcal{J}\left({\Omega }_1\right)=\left(\begin{array}{ccccc}-\left(m+\mu \right)& \frac{-\alpha \Lambda }{m+\mu }& 0& 0& 0\\ 0& \left({\mathcal{R}}_0-1\right)\left(c+\mu \right)& 0& 0& 0\\ 0& 0& -\left(z+\mu +\sigma \right)& 0& 0\\ m& \frac{-pm\Lambda }{\left(m+\mu \right)\mu }& 0& -\mu & 0\\ 0& c& z& 0& -\mu \end{array}\right)\mathrm{.}$

Then, one can obtain the characteristic equation of the matrix $\mathcal{J}\left({\Omega }_1\right)$ as:

${\Omega }_1\left(\lambda \right)={\left(\lambda +\mu \right)}^2\left(\lambda +m+\mu \right)\left(\lambda +z+\mu +\sigma \right)\left(\lambda -\left({\mathcal{R}}_0-1\right)\left(c+\mu \right)\right)=0$

where

${\lambda }_{1,2}=-\mu ,$

${\lambda }_3=-\left(m+\mu \right),$

${\lambda }_4=-\left(z+\mu +\sigma \right)\mathrm{,}$

${\lambda }_5=\left({\mathcal{R}}_0-1\right)\left(c+\mu \right)\mathrm{,}$

are the solutions of the characteristic equation. It is clear that ${\lambda }_{\mathrm{1,2,3,4}}$ are negative. Also, we If ${\mathcal{R}}_0<1$ , then we can obtain that ${\lambda }_5$ is negative, which means that DFE is locally asymptotically stable. This completes the proof.

5. Numerical Results and Discussion

In this section, we discuss the importance of numerical results for the model we have constructed by considering the vaccination strategies. For this simulation we have taken few parameter values from the literature [41] and we have estimated several of them as given in Table 1. In Figure 2, we have depicted the

Figure 2. Population densities of the COVID-19 model for the estimated values.

Table 1. Parameters used for the Covid 19 model.

behaviors and densities of the population dynamics in the model. In Figure 3, we have represented the effect of the parameter $\alpha$ on the Exposed and Infected individuals. According to Figure 3, it can be understood that as the parameter $\alpha$ increases from 0.002 to 0.01, the I and E populations are getting increases too, accordingly.

In Figure 4, we have showed the effectiveness of the parameter f on the Exposed and Infected individuals. By taking into account Figure 4, one can conclude that as the parameter increases from 0.01 to 0.018, when the populations

Figure 3. Population densities of infected and exposed individuals for various values of $\alpha =0.002,0.004,0.006,0.008,0.01$ bottom to top.

Figure 4. Population densities of infected (I) and exposed (E) individuals for various values of $f=0.01,0.012,0.014,0.016,0.018$ bottom to top for I and top to bottom for E.

in I increases, the population of E decreases. In another Figure 5 and Figure 6, we have explained the dynamics of vaccinated population by considering the values of p and m, respectively. According to these figures, we can say that as the parameter p increases from 0.1 to 0.2, the number of vaccinated individuals decreases, while the parameter m increases from 0.1 to 0.5, the number of vaccinated individuals increases, too.

In Figure 7, we have depicted the population density of the recovered individuals with respect to the parameter c. It is clear from the figure that as the values of the parameter increases from 0.01 to 0.05, the density of the recovered population increases as well.

Figure 5. Vaccinated population density (V) for various values of $p=0.1,0.12,0.14,\cdots ,0.2$ top to bottom.

Figure 6. Vaccinated population density (V) for various values of $m=0.1,0.2,0.3,0.4,0.5$ bottom to top.

Figure 7. Recovered population density (R) for various values of $c=0.01,0.02,0.03,0.04,0.05$ bottom to top.

6. Concluding Remarks

In this study, by taking into account an SEIR-type model, an SEIVR model that contains an effective vaccination strategy has been constructed. We have pointed out the effectiveness and different aspects of the vaccination by the mentioned model. We have determined the biological feasible region, positiveness of the solutions and the boundedness. By using the Lipschitz conditions we have presented the existence and uniqueness of the solutions.

Moreover, we have provided both the disease-free equilibrium and the endemic equilibrium points. Basic reproduction number has been evaluated by the benefiting the matrix generation method. In what follows, we have investigated the stability analyses of the mentioned equilibrium points and have proved that the DFE is asymptotically stable when ${\mathcal{R}}_0<1$ , and unstable when ${\mathcal{R}}_0>1$ . In numerical simulations we have used ode45 symbolic code of MATLAB R2020b to obtain the solutions to the system of differential equations.

According to the results that have been represented in the figures, we have provided the effectiveness of the parameters on the population dynamics. In Figure 2, we have depicted the behaviors and densities of the population dynamics in the model. We have represented the effect of the parameter $\alpha$ and f on the exposed and infected individuals in Figure 3 and Figure 4, respectively. In Figure 5 and Figure 6, we have explained the dynamics of vaccinated population by considering the values of p and m, respectively. According to these figures, we can say that as the parameter p increases from 0.1 to 0.2, the number of vaccinated individuals decreases, while the parameter m increases from 0.1 to 0.5, the number of vaccinated individuals increases, too. In Figure 7, we have depicted the population density of the recovered individuals with respect to the parameter c.

In the future studies, optimal control strategies can be performed to point out the effects of the treatment and vaccination thoughts. Also, the extension of the integer-order model to the fractional order can be considered.

Acknowledgements

The authors would like to thank the reviewers and editors of this paper for their careful attention to detail and constructive feedback that improved the presentation of the paper greatly.

The authors were supported by TUBITAK (The Scientific and Technological Research Council of Turkey).

Conflicts of Interest

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

References

[1]	World Health Organization Weakly Report. https://www.who.int/publications/m/item/weekly-epidemiological-update---5-january-2021
[2]	World Health Organization. https://covid19.who.int
[3]	Verschaffel, L., Greer, B. and De Corte, E. (2002) Everyday Knowledge and Mathematical Modeling of School Word Problems. In: Symbolizing, Modeling and Tool Use in Mathematics Education, Springer, Berlin, 257-276. https://doi.org/10.1007/978-94-017-3194-2_16
[4]	Zu, J., Li, M.L., Li, Z.F., Shen, M.W., Xiao, Y.N. and Ji, F.P. (2020) Transmission Patterns of COVID-19 in the Mainland of China and the Efficacy of Different Control Strategies: A Data-and Model-Driven Study. Infectious Diseases of Poverty, 9, 1-14. https://doi.org/10.1186/s40249-020-00709-z
[5]	Tang, B., Xia, F., Tang, S., Bragazzi, N.L., Li, Q., Sun, X., Liang, J., Xiao, Y. and Wu, J. (2020) The Effectiveness of Quarantine and İsolation Determine the Trend of the COVID-19 Epidemic in the Final Phase of the Current Outbreak in China. International Journal of Infectious Diseases, 96, 636-647. https://doi.org/10.1016/j.ijid.2020.05.113
[6]	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
[7]	Okuonghae, D. and Omame, A. (2020) 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
[8]	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.1101/2020.11.22.20236091
[9]	Farooq, J. and Bazaz, M.A. (2020) A Novel Adaptive Deep Learning Model of COVID-19 with Focus on Mortality Reduction Strategies. Chaos, Solitons and Fractals, 138, Article ID: 110148. https://doi.org/10.1016/j.chaos.2020.110148
[10]	Naik, P.A., Yavuz, M., Qureshi, S., Zu, J. and Townley, S. (2020) Modeling and Analysis of COVID-19 Epidemics with Treatment in Fractional Derivatives Using Real Data from Pakistan. The European Physical Journal Plus, 135, Article No. 795. https://doi.org/10.1140/epjp/s13360-020-00819-5
[11]	Uçar, E., Uçar, S., Evirgen, F. and ÖZdemir, N. (2021) Investigation of E-Cigarette Smoking Model with Mittag-Leffler Kernel. Foundations of Computing and Decision Sciences, 46, 97-109. https://doi.org/10.2478/fcds-2021-0007
[12]	Naik, P.A., Owolabi, K.M., Yavuz, M. and Zu, J. (2020) Chaotic Dynamics of a Fractional Order HIV-1 Model Involving AIDS-Related Cancer Cells. Chaos, Solitons and Fractals, 140, Article ID: 110272. https://doi.org/10.1016/j.chaos.2020.110272
[13]	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
[14]	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
[15]	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
[16]	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
[17]	Gao, W., Veeresha, P., Baskonus, H.M., Prakasha, D.G. and Kumar, P. (2020) A New Study of Unreported Cases of 2019-nCOV Epidemic Outbreaks. Chaos, Solitons and Fractals, 138, Article ID: 109929. https://doi.org/10.1016/j.chaos.2020.109929
[18]	Yavuz, M. and Özdemir, N. (2020) Analysis of an Epidemic Spreading Model with Exponential Decay Law. Mathematical Sciences and Applications E-Notes, 8, 142-154. https://doi.org/10.36753/mathenot.691638
[19]	Gao, W., Veeresha, P., Prakasha, D.G. and Baskonus, H.M. (2020) Novel Dynamic Structures of 2019-nCoV with Nonlocal Operator via Powerful Computational Technique. Biology, 9, 107. https://doi.org/10.3390/biology9050107
[20]	Jena, R.M., Chakraverty, S., Yavuz, M. and Abdeljawad, T. (2021) A New Modeling and Existence-Uniqueness Analysis for Babesiosis Disease of Fractional Order. Modern Physics Letters B.
[21]	Yavuz, M. and Bonyah, E. (2019) New Approaches to the Fractional Dynamics of Schistosomiasis Disease Model. Physica A: Statistical Mechanics and Its Applications, 525, 373-393. https://doi.org/10.1016/j.physa.2019.03.069
[22]	Bulut, H., Kumar, D., Singh, J., Swroop, R. and Baskonus, H.M. (2018) Analytic Study for a Fractional Model of HIV Infection of CD4+ T Lymphocyte Cells. Mathematics in Natural Science, 2, 33-43. https://doi.org/10.22436/mns.02.01.04
[23]	Uçar, E., Uçar, S., Evirgen, F. and Özdemir, N. (2021) A Fractional SAIDR Model in the Frame of Atangana-Baleanu Derivative. Fractal and Fractional, 5, 32. https://doi.org/10.3390/fractalfract5020032
[24]	Uçar, S. (2020) Analysis of a Basic SEIRA Model with Atangana-Baleanu Derivative. AIMS Mathematics, 5, 1411-1424. https://doi.org/10.3934/math.2020097
[25]	Singh, J., Kumar, D., Hammouch, Z. and Atangana, A. (2018) A Fractional Epidemiological Model for Computer Viruses Pertaining to a New Fractional Derivative. Applied Mathematics and Computation, 316, 504-515. https://doi.org/10.1016/j.amc.2017.08.048
[26]	Yavuz, M. (2021) Characterizations of Two Different Fractional Operators without Singular Kernel. Mathematical Modelling of Natural Phenomena, 14, 13.
[27]	Akgül, E.K., Akgül, A. and Yavuz, M. (2021) New Illustrative Applications of Integral Transforms to Financial Models with Different Fractional Derivatives. Chaos, Solitons and Fractals, 146, Article ID: 110877. https://doi.org/10.1016/j.chaos.2021.110877
[28]	Yavuz, M. (2021) European option pricing models described by fractional operators with classical and generalized Mittag-Leffler kernels. Numerical Methods for Partial Differential Equations. https://doi.org/10.1002/num.22645
[29]	Sene, N. (2021) Introduction to the Fractional-Order Chaotic System under Fractional Operator in Caputo Sense. Alexandria Engineering Journal, 60, 3997-4014. https://doi.org/10.1016/j.aej.2021.02.056
[30]	Sene, N. (2021) Analysis of a Fractional-Order Chaotic System in the Context of the Caputo Fractional Derivative via Bifurcation and Lyapunov Exponents. Journal of King Saud University—Science, 33, Article ID: 101275. https://doi.org/10.1016/j.jksus.2020.101275
[31]	Yavuz, M. and Yaşkıran, B. (2018) Conformable Derivative Operator in Modelling Neuronal Dynamics. Applications and Applied Mathematics, an International Journal, 13, 803-817.
[32]	Yavuz, M. and Yokus, A. (2020) Analytical and Numerical Approaches to Nerve Impulse Model of Fractional-Order. Numerical Methods for Partial Differential Equations, 36, 1348-1368. https://doi.org/10.1002/num.22476
[33]	Yavuz, M. and Sene, N. (2020) Approximate Solutions of the Model Describing Fluid Flow Using Generalized ρ-Laplace Transform Method and Heat Balance Integral Method. Axioms, 9, 123. https://doi.org/10.3390/axioms9040123
[34]	Yavuz, M. (2019) Dynamical Behaviors of Separated Homotopy Method Defined by Conformable Operator. Konuralp Journal of Mathematics, 7, 1-6.
[35]	Yokus, A. and Yavuz, M. (2021) Novel Comparison of Numerical and Analytical Methods for Fractional Burger-Fisher Equation. Discrete and Continuous Dynamical Systems-S. https://doi.org/10.3934/dcdss.2020258
[36]	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
[37]	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
[38]	Mandal, M., Jana, S., Nandi, S.K., Khatua, A., Adak, S. and Kar, T.K. (2020) A Model Based Study on the Dynamics of COVID-19: Prediction and Control. Chaos, Solitons and Fractals, 136, Article ID: 109889. https://doi.org/10.1016/j.chaos.2020.109889
[39]	Obasi, C. and Mbah, G.C.E. (2019) On the Stability Analysis of a Mathematical Model of Lassa Fever Disease Dynamics. Journal of the Nigerian Society for Mathematical Biology, 2, 135-144.
[40]	Birkhoff, G. and Rota, G.C. (1989) Ordinary Differential Equations. Wiley, Hoboken.
[41]	Hethcote, H.W. (2000) The Mathematics of Infectious Diseases. SIAM Review, 42, 599-653. https://doi.org/10.1137/S0036144500371907