Epidemic model of COVID-19 outbreak by inducing behavioural response in population

Abstract

COVID-19 has spread around the world since December 2019, creating one of the greatest pandemics ever witnessed. According to the current reports, this is a situation when people need to be more careful and take the precaution measures more seriously, unless the condition may become even worse. Maintaining social distances and proper hygiene, staying at isolation or adopting the self-quarantine method are some of the common practices that people should use to avoid the infection. And the growing information regarding COVID-19 and its symptoms help the people to take proper precautions. In this present study, we consider an SEIRS epidemiological model on COVID-19 transmission which accounts for the effect of an individual’s behavioural response due to the information regarding proper precautions. Our results indicate that if people respond to the growing information regarding awareness at a higher rate and start to take the protective measures, then the infected population decreases significantly. The disease fatality can be controlled only if a large proportion of individuals become immune, either by natural immunity or by a proper vaccine. In order to apply the latter option, we need to wait until a safe and proper vaccine is developed and it is a time-taking process. Hence, in the latter part of the work, an optimal control problem is considered by implementing control strategies to reduce the disease burden. Numerical figures show that the control denoting behavioural response works with higher intensity immediately after implementation and then gradually decreases with time. Further, the control policy denoting hospitalisation of infected individuals works with its maximum intensity for quite a long time period following a sudden decrease. As, the implementation of the control strategies reduce the infected population and increase the recovered population, so, it may help to reduce the disease transmission at this current epidemic situation.

1 Introduction

The first outbreak caused by a novel betacoronavirus was reported in Wuhan, capital of Hubei Chinese province in December 2019 [8, 9, 28]. Initially, most of the cases were found around the wholesale Huanan seafood market, Wuhan where live animals are also traded [24]. But within one and a half month COVID-19 spread to all Chinese province and to the rest of the world. World Health Organisation (WHO) officially declared COVID-19 as pandemic instead of an epidemic on 17 March 2020. It is an RNA virus which belongs to the Coronaviridae family and of order Nidovirales, known as SARS-CoV-2 [11, 33] and it is reported that the main symptoms of the disease include viral pneumonia, fever, dry cough, tiredness, aches and pains, nasal congestion, breathing problems or even a variety of unspecific symptoms [5, 8, 18, 43]. The severity of the infection is high enough with an estimated case fatality ratio of order 1% [13, 14, 40, 41] and hence, the virus has made a public health priority issue given the anticipated size of the pandemic due to the absence of pre-existing immunity. According to the reports of the dashboard provided by CSSE of John Hopkins University, at April 26, the number of reported infected cases, the number of documented death and the number of documented recovery reached almost 2,912,421; 203,432 and 825,886, respectively, throughout the world [6]. Most of the countries take the epidemic outbreak seriously from the first day by implementing proper public health measures including non-pharmaceutical interventions also. Among 185 countries, United States (939,249 cases), Spain (223,759 cases), Italy (195,351 cases), France (161,644 cases), Germany (156,513 cases), United Kingdom (149,569 cases) are facing worse epidemic situations as the activated infected cases exceed 100,000 there. In particular, the number of infected cases in the United States has grown very fast, the number of reported infected cases increases from 15 to 939,249 till April, 26, though the death case is highest in Italy with number 26,384. China was the first country where quarantine strategy was implemented in Wuhan on 23 January 2020. In the current situation, China has 83,909 confirmed cases with the number of documented death (in Hubei) and recovery are 4642 and 77,346, respectively. Other countries also apply the same strategy, i.e. implementing national lockdown to reduce the infection transmission as France has done on March 17 or United Kingdom has done on March 23 or even India has done on March 25. Comparing the data and strategies with other countries, it looks like there is a large number of cases of Covid-19 in India which is not registered as India still does not has the sufficient number of test kits. The undocumented infected individuals obviously facilitate the rapid spread of COVID-19 [28]. It is the third time when zoonotic human coronavirus has spread in this century. Before this, in 2002, severe acute respiratory syndrome coronavirus (SARS-CoV) spread among 37 countries and also in 2012, Middle East respiratory syndrome coronavirus (MERS-CoV) spread among 27 countries.

The first case of COVID-19 was confirmed at Kerala in India on 30 January 2020. According to NIC, India, there are total 20,177 confirmed cases, 5914 recoveries and 826 deaths are reported in the country till 26 April 2020 [22, 29]. The Indian government has announced to maintain social distance or to adopt self-quarantine strategy as precaution measures to avoid large-scale disease transmission among the population. In fact, on 22 March 2020, central government implemented a 14-hours long “Janata curfew” when the number of affected cases crossed 500. Moreover, the Government of India also announced for a 21-days national lockdown from March 25 in order to reduce the spread of COVID-19. But later, realising the importance of the current problem, the government has increased the lockdown period up to 3 May 2020. So far no vaccine has been discovered for novel coronavirus and so, maintaining social distances or applying self-isolation are taken as common ways for prevention of disease transmission [15]. The lockdown includes the ban on people from stepping out of their homes, closure of all shops except medical stores, hospitals, banks, grocery shops, etc., suspension of all educational institutions and offices (only work-from-home is allowed), suspension of all public and private transport and also the prohibition of all political, cultural, sports, entertainment, religious activities. It has no doubt that this current outbreak has seriously affected life both economically and healthwise. According to the World Bank and RBI, after 1991, this will be the first time when the economic growth rate in India will be decreased by 1.5–2.8% due to novel coronavirus outbreak. So, it has become a matter of concern for all of us how long our social life will go through this calamity.

Till now there are some studies revealing some interesting statistical results about COVID-19 outbreak [12, 19, 30, 32, 35, 37, 38]. Based on the data from 31 December 2019 to 28 January 2020, Wu et al. proposed an SEIR model to analyse the disease transmission dynamics on national and global basis [42]. Tang et al. proposed a compartmental model for COVID-19 where clinical development of the disease, current status of infected patients and control measures are combined. The results reveal that the control reproduction number may reach up to 6.47, and the control policies including social distancing, quarantine and isolation can minimise overall COVID cases [38]. According to a report submitted by Cambridge University, India’s strategy of announcing 21-days lockdown may not be sufficient enough to prevent the large-scale outbreak of coronavirus epidemic as it can bounce back in months and cause infection at a higher rate [36]. They suggested for an extension of two or three lockdowns with 5-days breaks in between or a single 49-days lockdown.

In this manuscript, we have proposed an SEIRS epidemic model to analyse coronavirus transmission where it is considered that the susceptible population can protect themselves from getting infected by taking proper precautions by inducing behavioural changes. India has a population of almost 135 crores and so, it is not possible to lockdown or apply home quarantine on all susceptible population. A proportion of the susceptible population may take the precaution measures successfully but the rest of the people become infective (asymptomatically or symptomatically). Again the recovered people may become susceptible later if they do not maintain the precaution carefully. Rest of the paper is categorised as follows: Sect. 2 contains the proposed epidemiological model which accounts for the information induced behavioural response of susceptible individuals. Section 3 proves that the model is well-posed while in Sect. 4, equilibrium points are obtained with basic reproduction number $(R_{0})$. Sensitivity analysis for different parameters is performed in Sect. 5. Local and global stability conditions of the equilibria are found in Sect. 6. The consequent section shows that the system undergoes a forward bifurcation at $R_{0}=1$. Section 8 contains the pictorial scenarios of the system dynamics without applying any control interventions. Later, a corresponding optimal control problem is formulated. Section 10 contains the numerical simulations by implementing the control strategies and the last section includes a brief conclusion.

2 Mathematical model: basic equations

We have elaborated a compartmental epidemic model here to analyse how the outbreak of COVID-19 affect the population worldwide. The total population N(t) at time t, in this work, is divided into six subclasses such as susceptible (S), exposed (E), asymptomatically infected (A), symptomatically infected (I), hospitalised and under treatment (H) and recovered class (R). The susceptible population become exposed when they come in contact with asymptomatically or symptomatically infected people or even with hospitalised individuals through the term $(\beta _{1}A+\beta _{2}I+\beta _{3}H)S$ where $\beta _{1}, \beta _{2},\ \beta _{3}$ are the rate of disease transmission per contact by an asymptomatic infected, symptomatic infected and hospitalised people, respectively. The constant recruitment rate is denoted as $\Lambda$ which is introduced in susceptible class. The term d denotes the natural death rate in all population, whereas $\delta _{1},\ \delta _{2},\ \delta _{3}$ are disease-related death rates in asymptomatically infected, symptomatically infected and hospitalised individuals, respectively. The people in the exposed class can move into either asymptomatic or symptomatic stage with probabilities $\phi$ and $(1-\phi )$, respectively, depending on whether any physical symptom of the disease has been observed in the infected people [7]. The terms $\eta ,\ \omega$ and $\gamma$ represent the progression rates from asymptomatic to symptomatic, symptomatic to hospitalised and hospitalised to recovered stages, respectively. Also, recovery from the disease does not guarantee permanent recovery and hence some of the recovered people move back to susceptible class further with rate constant $\xi$ [34]. The degradation rate of information with time, caused by natural fading of memory about the consequences, is denoted by $a_{0}$.

Now, when a disease outbreak at a higher rate within a short time period, various media platforms like TV, newspaper try to spread awareness among individuals. The Government also chooses social and educational campaigns to demonstrate the precaution measures. Due to the awareness programs, people become cautious at a higher rate and the disease transmission rate becomes lower and people move to recovered class directly. Now, this density of information and awareness is directly proportional to symptomatically infected individuals and it depends on how many people become infected rapidly in a short time interval. Let, Z(t) denotes the density of awareness due to information in susceptible population so that $Z(t)=0$ in absence of symptomatic infected people. This awareness leads to the behavioural changes among the susceptible population to protect themselves from infection. The importance of information in spreading of Ebola in Senegal has been described in some studies [20]. Though the government tries to spread necessary information. Everyone does not become careful enough all the time: insufficient resources, poor financial condition and heedless nature are some of the reasons in this case. So, it is considered that a proportion of susceptible adopts the changes in their behaviour by responding to the awareness programs and move into recovered class. Now the changes depend not only on susceptible but on information too and so, is a function of S and Z. Moreover, the rate at which information spread depends only on symptomatically infected population.

Moreover, $u_{1}k$ is the rate of behavioural changes of the susceptible population by taking proper protective measures such as maintaining hygiene, social distances, self-isolation, etc., to avoid the disease prevalence. Here, $u_{1}$ and k, respectively, denote the response rate and the information interaction rate through which individuals adopt new behaviour and change their old habits. It is obvious that this type of response is not fully effective because of financial problems, heedless nature, etc., and hence, we have considered $0\le u_{1}\le 1$. Further, the information grows according to a saturation rate function $\displaystyle \frac{pI}{1+qI}$ where p and q represent the ‘growth rate of information’ and the ‘level for unresponsiveness towards information’, respectively [3]. It depends only on the symptomatic infective population. At the time of an epidemic outbreak, government, different health agencies and media platforms become more active to spread awareness among people regarding the protective measures to avoid the disease prevalence. It is assumed that at early stages, the growth of information increases with the increase in symptomatic infective population density but it ultimately comes to saturation with time.

The proposed SEIRS model on COVID-19 mainly analyses the dynamics when the susceptible are provided with necessary information regarding the disease and its precautionary measures. The overall infected class is divided into asymptomatic and symptomatic compartments. There are some reports which reveal that a person may become COVID-19 positive without showing any symptoms. Also, in some cases, the symptoms occur at a later stage and so, incorporation of the asymptomatically infected compartment is justifiable. Moreover, we have considered that the virus can be transmitted to the susceptible population when they get in touch with asymptomatically infected, symptomatically infected and also even hospitalised individuals. The system emphasises how the information regarding the disease, its precautions and also people’s awareness affect the disease propagation in this current pandemic situation.

So, the proposed model with positive parametric values takes the following form:

$$\begin{aligned} \begin{aligned}&\frac{\hbox {d}S}{\hbox {d}t}= \Lambda -(\beta _{1}A+\beta _{2}I+\beta _{3}H)S-dS \\&\quad +\xi R-u_{1}kSZ, \quad S(0)>0, \\&\frac{\hbox {d}E}{\hbox {d}t}= (\beta _{1}A+\beta _{2}I+\beta _{3}H)S-(\kappa +d)E, \quad E(0)\ge 0, \\&\frac{\hbox {d}A}{\hbox {d}t}= \kappa \phi E-\eta A-(d+\delta _{1})A, \ \ A(0)\ge 0, \\&\frac{\hbox {d}I}{\hbox {d}t}= \kappa (1-\phi )E+\eta A-(d+\delta _{2})I-\omega I, \quad I(0)\ge 0, \\&\frac{\hbox {d}H}{\hbox {d}t}=\omega I-\gamma H-(d+\delta _{3})H, \quad H(0)\ge 0, \\&\frac{\hbox {d}R}{\hbox {d}t}= \gamma H-dR-\xi R+u_{1}kSZ, \quad R(0)\ge 0, \\&\frac{\hbox {d}Z}{\hbox {d}t}=\frac{pI}{1+qI}-a_{0}Z, \quad Z(0)\ge 0, \\ \end{aligned} \end{aligned}$$

(1)

A schematic diagram is provided in Fig. 1 to get a better insight of the proposed system.

Fig. 1

Schematic diagram of system (1)

3 Positivity and boundedness

For system (1): the following two theorems prove that the system variables are positive and bounded for all time.

Theorem 3.1

All solutions of system (1) starting from ${\mathbb {R}}_{+}^{7}$ remain positive for all the time.

Proof

Continuity and locally Lipschitzian functions of right-hand side of model (1) on C result in occurrence of a unique solution $(S(t),E(t),A(t),I(t),H(t),R(t),Z(t))$ on $[0,\tau ),$ where $0<\tau \le +\infty$ [21]. First we need to show that, $S(t)>0, \ \forall \, t\in [0,\tau )$. If it does not hold, then $\exists \ t_{1} \in (0,\tau )$ such that $S(t_{1})=0,\ {\dot{S}}(t_{1})\le 0$ and $S(t)>0, \ \forall \, t\in [0,t_{1}).$ Hence we must have $E(t)\ge 0,\ \forall \, t\in [0,t_{1}).$ If it is not true, then $\exists \ t_{2} \in (0,t_{1})$ such that $E(t_{2})=0,\ {\dot{E}}(t_{2})< 0$ and $E(t)\ge 0, \ \forall \, t\in [0,t_{2}).$ Next we claim $A(t)\ge 0,\ \forall \, t\in [0,t_{2}).$ Suppose it is not true. Then $\exists \ t_{3} \in (0,t_{2})$ such that $A(t_{3})=0,\ {\dot{A}}(t_{3})< 0$ and $A(t)\ge 0, \ \forall \, t\in [0,t_{3}).$ From third equation of (1), we have

$$\begin{aligned} \frac{\hbox {d}A}{\hbox {d}t}\bigg |_{t=t_{3}}= & {} \kappa \phi E(t_{3})-(\eta +d+\delta _{1}) A(t_{3}) \\= & {} \kappa \phi E(t_{3})\ge 0, \end{aligned}$$

which is a contradiction to ${\dot{A}}(t_{3})< 0.$ So, $A(t)\ge 0, \ \forall \, t\in [0,t_{2}).$

Our next claim is $I(t)\ge 0, \ \forall \, t\in [0,t_{2}).$ If it is not true, then $\exists \ t_{4} \in (0,t_{2})$ such that $I(t_{4})=0,\ {\dot{I}}(t_{4})< 0$ and $I(t)\ge 0, \ \forall \, t\in [0,t_{4}).$ Now from the fourth equation of (1):

$$\begin{aligned} \frac{\hbox {d}I}{\hbox {d}t}\bigg |_{t=t_{4}}= & {} \kappa (1-\phi )I(t_{4})+\eta I(t_{4}) \\&-(\omega +d+\delta _{2})I(t_{4}) \\= & {} \kappa (1-\phi )I(t_{4})+\eta I(t_{4})\ge 0, \end{aligned}$$

which is a contradiction to ${\dot{I}}(t_{4})< 0.$ Hence, $I(t)\ge 0, \ \forall \, t\in [0,t_{2}).$

Next claim is $H(t)\ge 0, \ \forall \, t\in [0,t_{2}).$ Let the statement is not true. Then $\exists \ t_{5} \in (0,t_{2})$ such that $H(t_{5})=0,\ {\dot{H}}(t_{5})< 0$ and $H(t)\ge 0, \ \forall \, t\in [0,t_{5}).$ From the fifth equation of (1), we have

$$\begin{aligned} \begin{aligned} \frac{\hbox {d}H}{\hbox {d}t}\bigg |_{t=t_{5}}&= \omega I(t_{5})-(\gamma +d+\delta _{3})H(t_{5}) \\&=\omega I(t_{5}) \ge 0, \end{aligned} \end{aligned}$$

which is a contradiction to ${\dot{H}}(t_{5})< 0.$ So, $H(t)\ge 0, \ \forall \, t\in [0,t_{2}).$ From the second equation of (1):

$$\begin{aligned} \begin{aligned} \frac{\hbox {d}E}{\hbox {d}t}\bigg |_{t=t_{2}}&= (\beta _{1}A(t_{2})+\beta _{2}I(t_{2})+\beta _{3}H(t_{2}))S(t_{2}) \\&\ge 0 \ [\because E(t_{2})= 0]. \end{aligned} \end{aligned}$$

It is a contradiction to ${\dot{E}}(t_{2})< 0.$ So, $E(t)\ge 0, \ \forall \, t\in [0,t_{1}).$ Hence $A(t)\ge 0,\ I(t)\ge 0,\ H(t)\ge 0,\ \ \forall \, t\in [0,t_{1}).$

Proceeding as before, it can be proved that $Z(t)\ge 0, \ R(t)\ge 0, \ \forall \, t\in [0,t_{1}).$

From the first equation of (1), we have

$$\begin{aligned} \begin{aligned} \frac{\hbox {d}S}{\hbox {d}t}\bigg |_{t=t_{1}}&= \Lambda + \xi R(t_{1}) > 0 \ [\because S(t_{1})= 0], \end{aligned} \end{aligned}$$

which is a contradiction to ${\dot{S}}(t_{1})\le 0.$ So, $S(t)> 0, \ \forall \, t\in [0,\tau ),$ where $0<\tau \le +\infty$. Also, by the previous steps we have $E(t)\ge 0, \ A(t) \ge 0,\ I(t)\ge 0, \ H(t) \ge 0,\ R(t)\ge 0$ and $Z(t)\ge 0,\ \forall \, t\in [0,\tau ),$ where $0<\tau \le +\infty .$ Hence proved. $\square$

Theorem 3.2

Solutions of system (1) which start from ${\mathbb {R}}_{+}^{7}$ remain uniformly bounded for all the time.

Proof

$$\begin{aligned} \begin{aligned} \text {Let},\ N(t)&=S(t){+}E(t)+A(t){+}I(t)+H(t)+R(t)\\ \therefore \frac{\hbox {d}N}{\hbox {d}t}&=\Lambda -dN \\ \Rightarrow 0&<N(t)\le N(0)e^{-dt} +\frac{\Lambda }{d}\left( 1-e^{-dt}\right) \end{aligned} \end{aligned}$$

where $N(0)=S(0)+E(0)+A(0)+I(0)+H(0)+R(0)$.

As $t\rightarrow \infty ,\ \displaystyle 0<N(t)\le \frac{\Lambda }{d}.$ From the last equation of system (1), we have

$$\begin{aligned} \begin{aligned} \frac{\hbox {d}Z}{\hbox {d}t}&= \frac{pI}{1+qI}-a_{0}Z \le pI-a_{0}Z \\ \Rightarrow \,&\frac{\hbox {d}Z}{\hbox {d}t}+a_{0}Z \le \frac{p\Lambda }{d}\ \ (\text {for large time} \ t). \\ \Rightarrow 0&<Z(t)\le Z(0)e^{-a_{0}t} +\frac{p\Lambda }{a_{0}d}\left( 1-e^{-a_{0}t}\right) \end{aligned} \end{aligned}$$

As $t\rightarrow \infty ,\ \displaystyle 0<Z(t)\le \frac{p\Lambda }{a_{0}d}.$

Hence, all solutions of system (1) enter in the region:

$$\begin{aligned} \Omega= & {} \left\{ (S,E,A,I,H,R,Z)\in {\mathbb {R}}_{+}^{7}: 0\le N(t)\le \frac{\Lambda }{d}; \right. \\&\quad \left. 0<Z(t)\le \frac{p\Lambda }{a_{0}d} \right\} . \end{aligned}$$

$\square$

4 Equilibrium analysis

Here, we obtain the equilibrium points of system (1) by solving the isoclines and determine the basic reproduction number $R_{0}$ to determine existence of endemic equilibrium point. System (1) has following equilibrium points:

• (1) Disease-free equilibrium point (DFE): $\displaystyle E_{0}(S_{0},0,0,0,0,0,0)=\left( \frac{\Lambda }{d}, 0, 0, 0, 0, 0, 0\right)$.

• (2) Endemic equilibrium point: $E^{*}(S^{*}, E^{*}, A^{*}, I^{*}, H^{*}, R^{*}, Z^{*})$.

4.1 Basic reproduction number $(R_{0})$

Basic reproduction number $R_{0}$ basically is a threshold value which denotes the average number of newly infected individuals after coming in contact with a single infected individual in a susceptible population. A method was developed by van den Driessche and Watmough [39] which is used here to determine the basic reproduction number. In system (1), people move into exposed class when the disease is introduced and also the infected classes are asymptomatic (A), symptomatic stages (I) and hospitalised stage (H). Let us take $x\equiv (E,A,I,H)$. Second, third, fourth and fifth equations of model (1) is written as:

$$\begin{aligned} \frac{\hbox {d}x}{\hbox {d}t}= & {} {\mathfrak {F}}(x)-\nu (x), \\ {\mathfrak {F}}(x)= & {} \left( \begin{array}{ccccc} (\beta _{1}A+\beta _{2}I+\beta _{3}H)S \\ 0 \\ 0\\ 0 \end{array} \right) ,\\ \nu (x)= & {} \left( \begin{array}{ccccc} (\kappa +d)E \\ -\kappa \phi E+(\eta +d+\delta _{1})A \\ -\kappa (1-\phi )E-\eta A+(\omega +d+\delta _{2})I \\ -\omega I+(\gamma +d+\delta _{3})H \end{array} \right) , \end{aligned}$$

where ${\mathfrak {F}}(x)$ contains only that compartment where new infection term is introduced and $\nu (x)$ contains rest of the terms. So, corresponding linearised matrices of ${\mathfrak {F}}(x)$ and $\nu (x)$ at disease-free equilibrium $E_{0}=(S_{0},0,0,0,0,0,0)$ are, respectively

$$\begin{aligned} F= & {} \left( D{\mathfrak {F}}(x)\right) (E_{0})=\left( \begin{array}{cccc} 0 &{} \beta _{1}S_{0} &{} \beta _{2}S_{0} &{} \beta _{3}S_{0} \\ 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 \end{array} \right) , \\ V= & {} \left( D\nu (x)\right) (E_{0}) \\= & {} \left( \begin{array}{cccc} \kappa +d &{} 0 &{} 0 &{} 0 \\ -\kappa \phi &{} \eta +d+\delta _{1} &{} 0 &{} 0 \\ -\kappa (1-\phi ) &{} -\eta &{} \omega +d+\delta _{2} &{} 0 \\ 0 &{} 0 &{} -\omega &{} \gamma +d+\delta _{3} \end{array} \right) \end{aligned}$$

Let, $\alpha _{0}=\kappa +d,\ \alpha _{1}=\eta +d+\delta _{1},\ \alpha _{2}=\omega +d+\delta _{2},\ \alpha _{3}=\gamma +d+\delta _{3}$.

Here, $|V|=\det (V) =(\kappa +d)(\eta +d+\delta _{1})(\omega +d+\delta _{2})(\gamma +d+\delta _{3})= \alpha _{0}\alpha _{1}\alpha _{2}\alpha _{3}$.

As the spectral radius of the next generation matrix $FV^{-1}$ is $R_{0}$, then the basic reproduction number of system (1) is given by [39]:

$$\begin{aligned} \begin{aligned} R_{0}&=\frac{\kappa S_{0}}{|V|} \\&\quad [\beta _{1}\phi \alpha _{2}\alpha _{3} +(\beta _{2}\alpha _{3}+\beta _{3}\omega )\{\phi \eta +\alpha _{1}(1-\phi )\}] \\&= \frac{\kappa \Lambda }{d\alpha _{0}\alpha _{1}\alpha _{2}\alpha _{3}}\\&\quad [\beta _{1}\phi \alpha _{2}\alpha _{3}+(\beta _{2}\alpha _{3}+ \beta _{3}\omega )\{\phi \eta +\alpha _{1}(1-\phi )\}] \end{aligned}\nonumber \\ \end{aligned}$$

(2)

Existence of endemic equilibrium point $E^{*}(S^{*}, E^{*}, A^{*}, I^{*}, H^{*}, R^{*}, Z^{*})$

From system (1):

$$\begin{aligned} \begin{aligned}&\Lambda -(\beta _{1}A^{*}+\beta _{2}I^{*}+\beta _{3}H^{*})S^{*}-dS^{*} \\&\quad +\xi R^{*}-u_{1}kS^{*}Z^{*} = 0, \\&(\beta _{1}A^{*}+\beta _{2}I^{*}+\beta _{3}H^{*})S-(\kappa +d)E^{*} = 0, \\&\kappa \phi E^{*}-(\eta +d+\delta _{1})A^{*} = 0, \\&\kappa (1-\phi )E^{*}+\eta A^{*}-(d+\delta _{2}+\omega )I^{*} = 0, \\&\omega I^{*}-(\gamma +d+\delta _{3})H^{*} = 0, \\&\gamma H^{*}-(d+\xi )R^{*}+u_{1}kS^{*}Z^{*} = 0, \\&\frac{pI^{*}}{1+qI^{*}}-a_{0}Z^{*}&=0. \end{aligned} \end{aligned}$$

(3)

Solving these equations, we get

$$\begin{aligned} S^{*}= & {} \frac{\Lambda }{dR_{0}},\ Z^{*}=\frac{pI^{*}}{a_{0}(1+qI^{*})}, \\ H^{*}= & {} \frac{\omega I^{*}}{\alpha _{3}},\ E^{*}=\frac{\alpha _{1}\alpha _{2}I^{*}}{\kappa \{\phi \eta +\alpha _{1}(1-\phi )\}},\\ A^{*}= & {} \frac{\phi \alpha _{2}I^{*}}{\{\phi \eta +\alpha _{1}(1-\phi )\}},\\ R^{*}= & {} \frac{[\gamma \omega da_{0}R_{0}(1+qI^{*})+u_{1}kp\Lambda \alpha _{3}] I^{*}}{a_{0}dR_{0}\alpha _{3}(d+\xi )(1+qI^{*})} \end{aligned}$$

and $I^{*}$ is a positive root of the equation:

$$\begin{aligned} P_{1}I^{2}+P_{2}I+P_{3}=0, \end{aligned}$$

where

$$\begin{aligned} \begin{aligned}&P_{1}=q\left[ \frac{\gamma \xi \omega }{\alpha _{3}(d+\xi )}-\frac{\alpha _{0}\alpha _{1}\alpha _{2}}{\kappa \{\phi \eta +\alpha _{1}(1-\phi )\}}\right] \\&=-\frac{q[\alpha _{1}\{d\alpha _{0}\alpha _{1}\alpha _{2}{+}\xi d\alpha _{2}\alpha _{3}{+}\kappa \xi \alpha _{3}(d{+}\delta _{2}){+}\xi \kappa \omega (d{+}\delta _{3})\}{+}\phi \gamma \xi \omega \kappa (d{+}\delta _{1})]}{\kappa \alpha _{3}(d{+}\xi )\{\phi \eta {+}\alpha _{1}(1{-}\phi )\}} <0 \\&P_{2}=q\Lambda \left( 1-\frac{1}{R_{0}}\right) +\left[ \frac{\gamma \xi \omega }{\alpha _{3}(d+\xi )}- \frac{\alpha _{0}\alpha _{1}\alpha _{2}}{\kappa \{\phi \eta +\alpha _{1}(1-\phi )\}} \right] -\frac{u_{1}kp\Lambda }{a_{0}R_{0}(d+\xi )}, \\&P_{3}=\Lambda \left( 1-\frac{1}{R_{0}}\right) \end{aligned} \end{aligned}$$

Here, $P_{3}>0$ when $R_{0}>1$. As $P_{1}$ is always negative, so, for $R_{0}>1,$ we get a unique endemic equilibrium point. But when $R_{0}<1$, both of $P_{2}$ and $P_{3}$ become negative resulting in non-occurrence of an endemic equilibrium. Hence we have the following theorem.

Theorem 4.1

System (1) has a disease-free equilibrium point $E_{0}\left( S_{0}, 0,0,0,0,0,0\right)$ for any parametric values. And for $R_{0}>1,$ the system possesses a unique endemic equilibrium $E^{*}(S^{*}, E^{*}, A^{*}, I^{*}, H^{*}, R^{*}, Z^{*}).$

5 Sensitivity analysis

For system (1), $R_{0}$ depends on some of the vital parameters such as recruitment rate $(\Lambda )$, transmission rates $(\beta _{1},\ \beta _{2},\ \beta _{3}),$ disease related death rates $(\delta _{1}, \delta _{2},\ \delta _{3})$, natural death rate (d), progression rate of exposed people into infected classes $(\kappa )$, probability at which exposed people move into infected classes $(\phi )$, progression rate of asymptomatic infected individual into symptomatic class $(\eta )$, progression rate of symptomatic infected population into treatment class $(\omega )$, progression rate of hospitalised population into recovered class $(\gamma )$. But among all these parameters, we cannot control only $\beta _{1},\ \beta _{2},\ \beta _{3},\ \kappa ,\ \omega .$

Now, $\displaystyle R_{0}= \frac{\kappa \Lambda }{d\alpha _{0}\alpha _{1}\alpha _{2}\alpha _{3}}[\beta _{1}\phi \alpha _{2}\alpha _{3}+(\beta _{2}\alpha _{3}+\beta _{3}\omega )\{\phi \eta +\alpha _{1}(1-\phi )\}],$ where $\alpha _{0}=\kappa +d,\ \alpha _{1}=\eta +d+\delta _{1},\ \alpha _{2}=\omega +d+\delta _{2}$ and $\alpha _{3}=\gamma +d+\delta _{3}$. From the expression of $R_{0}$:

$$\begin{aligned} \begin{aligned} \frac{\partial R_0}{\partial \beta _{1}}&=\frac{\Lambda \kappa \phi }{d\alpha _{0}\alpha _{1}}>0 \\ \frac{\partial R_0}{\partial \beta _{2}}&=\frac{\Lambda \kappa \{\phi \eta +\alpha _{1}(1-\phi )\}}{d\alpha _{0}\alpha _{1}\alpha _{2}}>0 \\ \frac{\partial R_0}{\partial \beta _{3}}&=\frac{\Lambda \kappa \omega \{\phi \eta +\alpha _{1}(1-\phi )\}}{d\alpha _{0}\alpha _{1}\alpha _{2}\alpha _{3}}>0 \\ \frac{\partial R_0}{\partial \kappa }&=\frac{\Lambda [\beta _{1}\phi \alpha _{2}\alpha _{3} {+}(\beta _{2}\alpha _{3}{+}\beta _{3}\omega )\{\phi \eta {+} \alpha _{1}(1{-}\phi )\}]}{\alpha _{0}^{2}\alpha _{1}\alpha _{2}\alpha _{3}} \\&>0 \\ \frac{\partial R_0}{\partial \omega }&=\frac{\Lambda \kappa \{\phi \eta {+}\alpha _{1}(1{-}\phi )\}}{d\alpha _{0}\alpha _{1}\alpha _{2}^{2} \alpha _{3}} \\&\quad [\beta _{3}(d+\delta _{2})-\beta _{2}(\gamma +d+\delta _{3})] \end{aligned} \end{aligned}$$

So, $\displaystyle \frac{\partial R_0}{\partial \omega }$ is monotonic decreasing when $\beta _{3}(d+\delta _{2})<\beta _{2}(\gamma +d+\delta _{3})$ holds. Next we compute normalised forward sensitivity index with respect to each of the parameters $\beta _{1},\ \beta _{2},\ \beta _{3},\ \kappa$ and $\omega$ to analyse the sensitivity of $R_0$ (to each of the parameters) by the method of Arriola and Hyman [1]:

$$\begin{aligned} \Gamma _{\beta _{1}}= & {} \left| \frac{\frac{\partial R_0}{R_0}}{\frac{\partial \beta _{1}}{\beta _{1}}}\right| =\left| \frac{\beta _{1}}{R_0}\frac{\partial R_0}{\partial \beta _{1}}\right| \\= & {} \left| \frac{\beta _{1}\phi \alpha _{2}\alpha _{3}}{[\beta _{1}\phi \alpha _{2}\alpha _{3}+(\beta _{2}\alpha _{3}+\beta _{3}\omega ) \{\phi \eta +\alpha _{1}(1-\phi )\}]}\right|<1 \\ \Gamma _{\beta _{2}}= & {} \left| \frac{\frac{\partial R_0}{R_0}}{\frac{\partial \beta _{2}}{\beta _{2}}}\right| =\left| \frac{\beta _{2}}{R_0}\frac{\partial R_0}{\partial \beta _{2}}\right| \\= & {} \left| \frac{\beta _{2}\alpha _{3} \{\phi \eta +\alpha _{1}(1-\phi )\}}{[\beta _{1}\phi \alpha _{2} \alpha _{3}+(\beta _{2}\alpha _{3}+\beta _{3}\omega ) \{\phi \eta +\alpha _{1}(1-\phi )\}]}\right|<1 \\ \Gamma _{\beta _{3}}= & {} \left| \frac{\frac{\partial R_0}{R_0}}{\frac{\partial \beta _{3}}{\beta _{3}}}\right| =\left| \frac{\beta _{3}}{R_0}\frac{\partial R_0}{\partial \beta _{3}}\right| \\= & {} \left| \frac{\beta _{3}\omega \{\phi \eta +\alpha _{1}(1-\phi )\}}{[\beta _{1}\phi \alpha _{2}\alpha _{3} +(\beta _{2}\alpha _{3}+\beta _{3}\omega )\{\phi \eta +\alpha _{1}(1-\phi )\}]}\right|<1 \\ \Gamma _{\kappa }= & {} \left| \frac{\frac{\partial R_0}{R_0}}{\frac{\partial \kappa }{\kappa }}\right| =\left| \frac{\hbox {d}}{\kappa +d}\right| <1 \\ \Gamma _{\omega }= & {} \left| \frac{\frac{\partial R_0}{R_0}}{\frac{\partial \omega }{\omega }}\right| =\left| \frac{\omega }{R_0}\frac{\partial R_0}{\partial \omega }\right| \\= & {} \left| \frac{\omega \{\phi \eta +\alpha _{1}(1-\phi )\} \{\beta _{3}(d+\delta _{2})- \beta _{2}(\gamma +d+\delta _{3})\}}{\alpha _{2} [\beta _{1}\phi \alpha _{2}\alpha _{3} +(\beta _{2} \alpha _{3}+\beta _{3}\omega )\{\phi \eta +\alpha _{1}(1-\phi )\}]}\right| \end{aligned}$$

From the calculation, it is observed that the disease transmission rates from asymptomatically infected $(\beta _{1})$, symptomatically infected $(\beta _{2})$ and hospitalised individuals $(\beta _{3})$ maintain direct proportional relation with basic reproduction number which is biologically acceptable too. Increasing disease transmission rates can stimulate the probability of occurrence of an epidemic outbreak. Moreover, if more people from the exposed class move into infected classes whether in asymptomatic stage or symptomatic stage, then the chance of getting an infected system increase. The calculation also shows that $R_{0}$ increases for increasing value of $\kappa$. On the other hand, $\omega$ denotes the rate at which the symptomatically infected individuals are admitted to hospitals for treatment. So, it is evident that increasing $\omega$ helps to reduce the disease prevalence and hence $\omega$ maintains an inversely proportional relation with $R_{0}$, i.e. increment in this parameter leads to a decrease in $R_{0}$. Here, $\displaystyle \frac{\partial R_{0}}{\partial \omega }<0$ only when $\beta _{3}(d+\delta _{2})<\beta _{2}(\gamma +d+\delta _{3})$ holds. So, maintaining this inequality can reduce the disease outbreak for increasing hospitalisation rate. The calculations and numerical simulations reveal that $R_{0}$ is more sensitive to changes in $\beta _{i}$ for $i=1,2,3$ than $\kappa$ and $\omega$. So, if we try to reduce the transmission rates by maintaining social distances and taking proper precautions, then this epidemic situation may be handled.

6 Stability analysis

We discuss the local and global stability conditions for the disease-free equilibrium point as well as endemic equilibrium point in this section. Let, $\alpha _{0}=d+\kappa ,\ \alpha _{1}=\eta +d+\delta _{1},\ \alpha _{2}=\omega +d+\delta _{2}$ and $\alpha _{3}=\gamma +d+\delta _{3}$. The Jacobian matrix of system (1) is given as:

$$\begin{aligned}&{\overline{J}}\nonumber \\&=\left( \begin{array}{ccccccc} a_{11} &{} 0 &{} -\beta _{1}S &{} -\beta _{2}S &{} -\beta _{3}S &{} \xi &{} -u_{1}kS \\ \beta _{1}A+\beta _{2}I+\beta _{3}H &{} -\alpha _{0} &{} \beta _{1}S &{} \beta _{2}S &{} \beta _{3}S &{} 0 &{} 0 \\ 0 &{} \kappa \phi &{} -\alpha _{1} &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} \kappa (1-\phi ) &{} \eta &{} -\alpha _{2} &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} \omega &{} -\alpha _{3} &{} 0 &{} 0 \\ u_{1}kZ &{} 0 &{} 0 &{} 0 &{} \gamma &{} -(d+\xi ) &{} u_{1}kS \\ 0 &{} 0 &{} 0 &{} \frac{p}{(1+qI)^{2}} &{} 0 &{} 0 &{} -a_{0} \end{array} \right) \nonumber \\ \end{aligned}$$

(4)

where $a_{11}=-(\beta _{1}A+\beta _{2}I+\beta _{3}H)-d-u_{1}kZ.$

6.1 Local stability of $E_{0}$

Jacobian matrix at the disease-free equilibrium point $\displaystyle E_{0}{=}\left( \frac{\Lambda }{d},0,0,0,0,0,0\right)$ is given as follows $\displaystyle \left( S_{0}{=}\frac{\Lambda }{d}\right)$:

$$\begin{aligned}\ {\overline{J}}|_{E_{0}}= \left( \begin{array}{ccccccc} -d &{} 0 &{} -\beta _{1}S_{0} &{} -\beta _{2}S_{0} &{} -\beta _{3}S_{0} &{} \xi &{} -u_{1}kS_{0} \\ 0 &{} -\alpha _{0} &{} \beta _{1}S_{0} &{} \beta _{2}S_{0} &{} \beta _{3}S_{0} &{} 0 &{} 0 \\ 0 &{} \kappa \phi &{} -\alpha _{1} &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} \kappa (1-\phi ) &{} \eta &{} -\alpha _{2} &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} \omega &{} -\alpha _{3} &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} \gamma &{} -(d+\xi ) &{} u_{1}kS_{0} \\ 0 &{} 0 &{} 0 &{} p &{} 0 &{} 0 &{} -a_{0} \end{array} \right) \end{aligned}$$

Some of the eigenvalues of the Jacobian matrix ${\overline{J}}|_{E_{0}}$ are $\lambda _{1}=-d,\ \lambda _{2}=-(d+\xi ),\ \lambda _{3}=-a_{0}$ and other four eigenvalues are obtained from the roots of the following equation: $\lambda ^{4}+Q_{1}\lambda ^{3}+Q_{2}\lambda ^{2}+Q_{3}\lambda +Q_{4}=0,$

where, $Q_{1}=\alpha _{0}+\alpha _{1}+\alpha _{2}+\alpha _{3}$, $Q_{2}=\alpha _{0}\alpha _{1}+\alpha _{2}\alpha _{3}+(\alpha _{0}+\alpha _{1})(\alpha _{2}+\alpha _{3})-\kappa S_{0}\{\beta _{1}\phi +\beta _{2}(1-\phi )\}$, $Q_{3}=\alpha _{0}\alpha _{1}(\alpha _{2}+\alpha _{3})+\alpha _{2}\alpha _{3}(\alpha _{0}+\alpha _{1})-S_{0}\kappa [\beta _{1}\phi (\alpha _{2}+\alpha _{3})+\beta _{2}\{\phi \eta +\alpha _{1}(1-\phi )\}+(\beta _{2}\alpha _{3}+\beta _{3}\omega )(1-\phi )]$ and $Q_{4}=\alpha _{0}\alpha _{1}\alpha _{2}\alpha _{3}(1-R_{0})$. So, $Q_{4}>0$ for $R_{0}<1$. So for $R_{0}<1$, the characteristic equation has roots with negative real parts only when $Q_{2},\ Q_{3}>0$ and we have the following theorem

Theorem 6.1

The DFE $(E_{0})$ is locally asymptomatically stable (LAS) for $R_{0}<1$ when $Q_{2}>0$ and $Q_{3}>0$ hold.

6.2 Global stability of $E_{0}$

Theorem 6.2

DFE $E_{0}$ of system (1) is globally asymptotically stable (GAS) when $p<d$ and $u_{1}k\Lambda <a_{0}d$ hold.

Proof

Let us consider the Lyapunov function as

$$\begin{aligned} V_{1}= & {} \left( S-S_{0}-S_{0}ln\left( \frac{S}{S_{0}}\right) \right) \\&+E+A+I+H+R+Z \end{aligned}$$

Here, $V_{1}$ is a positive definite function for all $(S,E,A,I,H,R,Z)$ other than the DFE. Time derivative of $V_1$ computed along the solutions of system (1) is given as:

$$\begin{aligned} \begin{aligned} \frac{\hbox {d}V_{1}}{\hbox {d}t}&=\left( 1-\frac{S_{0}}{S}\right) \\&\quad [\Lambda {-}(\beta _{1}A{+}\beta _{2}I{+}\beta _{3}H)S{-} dS{+}\xi R{-}u_{1}kSZ] \\&\quad +(\beta _{1}A+\beta _{2}I+\beta _{3}H)S \\&\quad -d(E+A+I+H+R)-\delta _{1} A-\delta _{2}I\\&\quad -\delta _{3}H-\xi R+u_{1}kSZ+\frac{pI}{1+qI}-a_{0}Z \\&=\Lambda - d(E+A+I+H+R)-\delta _{1} A-\delta _{2}I\\&\quad -\delta _{3}H -\frac{\Lambda S_{0}}{S} +dS_{0}-\frac{\xi S_{0}R}{S}\\&\quad +u_{1}kS_{0}Z + \frac{pI}{1+qI}-a_{0}Z\\&\le \left( \frac{\Lambda }{S}-d\right) (S-S_{0})+(p-d)I \\&\quad +\frac{1}{d}(u_{1}k\Lambda -a_{0}d)Z \\&=-\frac{\hbox {d}}{S}(S-S_{0})^{2}+(p-d)I \\&\quad +\frac{1}{d}(u_{1}k\Lambda -a_{0}d)Z \end{aligned} \end{aligned}$$

Hence, $\displaystyle \frac{\hbox {d}V_{1}}{\hbox {d}t}<0$ when $p<d$ and $u_{1}k\Lambda <a_{0}d$ hold. Also, $\frac{\hbox {d}V_{1}}{\hbox {d}t}=0$ when $S=S_{0}$ and $E=A=I=H=R=Z=0.$ By LaSalle’s invariance principle [27], $E_{0}$ is globally asymptotically stable when $R_{0}<1$ when the mentioned parametric restrictions are fulfilled. $\square$

6.3 Local stability of $E^{*}$

Theorem 6.3

The endemic equilibrium point $E^{*}$ is LAS for $R_{0}>1$ when the conditions if $R_{i}>0$ for $i=1,7$ and $\Delta _i >0$ for $i=2,3,4,5,6$ hold.

Proof

Proof is given in the “Appendix”. $\square$

6.4 Global stability of $E^{*}$

Theorem 6.4

The endemic equilibrium point $E^{*}$ of system (1) is globally asymptomatically stable (GAS) in

$$\begin{aligned} \begin{aligned} \Phi&=\left\{ (S,E,A,I,H,R,A)\in {\mathbb {R}}^{7}_{+}: \right. \\&\quad (m_{2}\beta _{1}S)^{2}<\frac{4}{9}m_{2}\alpha _{0} \\&\quad (\delta _{1}+m_{3}\alpha _{1}),\ (m_{2}\beta _{3}S-\delta _{3})^{2} \\&<\frac{1}{6}m_{2}\delta _{3}\alpha _{3}, \\&\quad (m_{2}\beta _{2}S+m_{4}\kappa \phi )^{2} \\&<\frac{m_{2}\alpha _{0}}{3}(\delta _{2}+m_{4}\alpha _{2}), \\&\quad \left. \text {with}\ (i) - (iii) \ \text {as stated at the end in the proof}\right\} , \end{aligned} \end{aligned}$$

where $\alpha _{i}$ for $i=0,1,2,3$ and $m_{j}$ for $j=1,2,3,4$ are mentioned in the proof.

Proof

Consider a Lyapunov function $V_{2}$ as:

$$\begin{aligned} \begin{aligned} V_{2}(t)&=\frac{1}{2}[(S-S^{*})+(E-E^{*})+(A-A^{*}) \\&\qquad +(I-I^{*})+(H-H^{*})+(R-R^{*})]^{2}+\\&\quad \frac{m_{1}}{2}(S-S^{*})^{2}+\frac{m_{2}}{2}(E-E^{*})^{2} \\&\quad +\frac{m_{3}}{2}(A-A^{*})^{2} +\frac{m_{4}}{2}(I-I^{*})^{2} \\&\quad +\frac{1}{2}(Z-Z^{*})^{2} \end{aligned} \end{aligned}$$

Time derivative of $V_{2}$ along the solutions of system (1) is given as:

$$\begin{aligned} \begin{aligned} \frac{\hbox {d}V_{2}}{\hbox {d}t}&=[(S-S^{*})+(E-E^{*})+(A-A^{*}) \\&\qquad +(I-I^{*})+(H-H^{*})+(R-R^{*})]\\&\quad \frac{\hbox {d}}{\hbox {d}t}(S+E+A+I+H+R)+m_{1}(S-S^{*})\frac{\hbox {d}S}{\hbox {d}t} \\&\quad +m_{2}(E-E^{*})\frac{\hbox {d}E}{\hbox {d}t}\\&\quad +m_{3}(A-A^{*})\frac{\hbox {d}A}{\hbox {d}t}+m_{4}(I-I^{*})\frac{\hbox {d}I}{\hbox {d}t} \\&\quad +(Z-Z^{*})\frac{\hbox {d}Z}{\hbox {d}t} \end{aligned} \end{aligned}$$

Now steady state of system (1) at $E^{*}$ gives

$$\begin{aligned}&\Lambda -(\beta _{1}A^{*}+\beta _{2}I^{*}+\beta _{3}H^{*})S^{*} -dS^{*} \\&\quad +\xi R^{*}-u_{1}kS^{*}Z^{*} = 0, \\&(\beta _{1}A^{*}+\beta _{2}I^{*}+\beta _{3}H^{*})S-(\kappa +d)E^{*} = 0, \\&\kappa \phi E^{*}-(\eta +d+\delta _{1})A^{*} = 0, \\&\kappa (1-\phi )E^{*}+\eta A^{*}-(d+\delta _{2}+\omega )I^{*} = 0, \\&\omega I^{*}-(\gamma +d+\delta _{3})H^{*} = 0, \\&\gamma H^{*}-(d+\xi )R^{*}+u_{1}kS^{*}Z^{*} = 0, \\&\frac{pI^{*}}{1+qI^{*}}-a_{0}Z^{*} =0, \\&\Lambda -d(S^{*}+E^{*}+A^{*}+I^{*}+H^{*}+R^{*}) \\&\quad -\delta _{1}A^{*}-\delta _{2}I^{*}-\delta _{3}H^{*} =0 \end{aligned}$$

Consider $\alpha _{0}=\kappa +d,\ \alpha _{1}=\eta +d+\delta _{1},\ \alpha _{2}=d+\delta _{2}+\omega ,\ \alpha _{3}=\gamma +d+\delta _{3}$. So, we have

$$\begin{aligned} \frac{\hbox {d}V_{2}}{\hbox {d}t}&=[(S-S^{*})+(E-E^{*})+(A-A^{*}) \\&\qquad +(I-I^{*})+(H-H^{*})+(R-R^{*})]\\&\quad [\Lambda -d(S+E+A+I+H+R)-\delta _{1}A\\&\qquad -\delta _{2}I-\delta _{3}H] \\&\quad +m_{1}(S-S^{*})[\Lambda -dS+\xi R\\&\quad -(\beta _{1}A+\beta _{2}I+\beta _{3}H)S-u_{1}kSZ] \\&\quad +m_{2}(E-E^{*})[(\beta _{1}A+\beta _{2}I+\beta _{3}H)S-\alpha _{0}E]\\&\quad +m_{3}(A-A^{*})[\kappa \phi E-\alpha _{1}A] \\&\quad +m_{4}(I-I^{*})[\kappa (1-\phi )E+\eta A-\alpha _{2}I]\\&\quad +(Z-Z^{*})\left( \frac{pI}{1+qI}-a_{0}Z\right) \\&=[(S-S^{*})+(E-E^{*})+(A-A^{*})+(I-I^{*})\\&\qquad +(H-H^{*})+(R-R^{*})] \\&\quad [-d(S-S^{*})-d(E-E^{*})-(d+\delta _{1})(A-A^{*}) \\&\quad -(d+\delta _{2})(I-I^{*})-(d+\delta _{3})\\&\quad (H-H^{*})-d(R-R^{*})]\\&\quad +m_{1}(S-S^{*}) \\&\quad [\Lambda -(\beta _{1}A+\beta _{2}I+\beta _{3}H)S-dS+\xi R-u_{1}kSZ]\\&\quad +m_{2}(E-E^{*})\\&\quad [(\beta _{1}A+\beta _{2}I+\beta _{3}H)S-\alpha _{0}E]\\&\quad +m_{3}(A-A^{*})[\kappa \phi E-\alpha _{1}A]+m_{4}(I-I^{*})\\&\quad [\kappa (1-\phi )E+\eta A-\alpha _{2}I]+(Z-Z^{*})\\&\quad \left( \frac{pI}{1+qI}-a_{0}Z\right) \\ \frac{\hbox {d}V_{2}}{\hbox {d}t}&=[(S-S^{*})+(E-E^{*})+(A-A^{*})\\&\qquad +(I-I^{*})+(H-H^{*})+(R-R^{*})]\\&\quad [-d(S-S^{*})-d(E-E^{*})-(d+\delta _{1})(A-A^{*})\\&\qquad -(d+\delta _{2})(I-I^{*})\\&\qquad -(d+\delta _{3})(H-H^{*})-d(R-R^{*})]\\&\quad +m_{1}(S-S^{*})\\&\quad [-\beta _{1}(AS-A^{*}S^{*})-\beta _{2}(IS-I^{*}S^{*})\\&\qquad -\beta _{3}(HS-H^{*}S^{*})-d(S-S^{*})\\&\qquad +\xi (R-R^{*})-u_{1}k(SZ-S^{*}Z^{*})] \\&\quad +m_{2}(E-E^{*})\\&\quad [\beta _{1}(AS-A^{*}S^{*})+\beta _{2}(IS-I^{*}S^{*})\\&\qquad +\beta _{3}(HS-H^{*}S^{*})-\alpha _{0}(E-E^{*})]\\&\quad +m_{3}(A-A^{*})\\&\quad [\kappa \phi (E-E^{*})-\alpha _{1}(A-A^{*})]\\&\quad +m_{4}(I-I^{*})\\&\quad [\kappa (1-\phi )(E-E^{*}+\eta (A-A^{*})-\alpha _{2}(I-I^{*})]\\&\quad +(Z-Z^{*})\\&\quad \left[ \frac{p(I-I^{*})}{(1+qI)(1+qI^{*})}-a_{0}(Z-Z^{*})\right] \\&=-d(S-S^{*})^{2}\\&\quad -d[(E-E^{*})+(A-A^{*})+(I-I^{*})+(H-H^{*})]^{2}\\&\quad -m_{2}\alpha _{0}(E-E^{*})^{2}\\&\quad -(\delta _{1}+m_{3}\alpha _{1})(A-A^{*})^{2}\\&\quad -(\delta _{2}+m_{4}\alpha _{2})(I-I^{*})^{2}\\&\quad -\delta _{3}(H-H^{*})^{2}-d(R-R^{*})^{2}\\&\quad +[-2d+m_{2}(\beta _{1}A^{*}+\beta _{2}I^{*}+\beta _{3}H^{*})]\\&\quad (S-S^{*})(E-E^{*})\\&\quad -[2d+\delta _{1}+m_{1}\beta _{1}S^{*}](S-S^{*})(A-A^{*})\\&\quad -[2d+\delta _{2}+m_{1}\beta _{2}S^{*}](S-S^{*})(I-I^{*})\\&\quad -[2d+\delta _{3}+m_{1}\beta _{3}S^{*}](S-S^{*})(H-H^{*})\\&\quad +(-2d+m_{1}\xi )(S-S^{*})(R-R^{*})\\&\quad +[-\delta _{1}+m_{2}\beta _{1}S+m_{3}\kappa \phi ](E-E^{*})(A-A^{*})\\&\quad +[-\delta _{2}+m_{2}\beta _{2}S+m_{4}\kappa \phi ](E-E^{*})(I-I^{*})\\&\quad +[-\delta _{3}+m_{2}\beta _{3}S](E-E^{*})(H-H^{*})\\&\quad -2d(E-E^{*})(R-R^{*})\\&\quad +[-(\delta _{1}+\delta _{2})+m_{4}\eta ](A-A^{*})(I-I^{*})\\&\quad -(\delta _{1}+\delta _{3})(A-A^{*})(H-H^{*})\\&\quad -\delta _{1}(A-A^{*})(R-R^{*})\\&\quad -(\delta _{2}+\delta _{3})(I-I^{*})(H-H^{*})\\&\quad -\delta _{2}(I-I^{*})(R-R^{*})\\&\quad -\delta _{3}(H-H^{*})(R-R^{*})\\&\quad +\frac{p(I-I^{*})(Z-Z^{*})}{(1+qI)(1+qI^{*})}-a_{0}(Z-Z^{*})^{2} \end{aligned}$$

Let, $\displaystyle m_{1}=\frac{2d}{\xi },\ m_{2}=\frac{2d}{\beta _{1}A^{*}+\beta _{2}I^{*}+\beta _{3}H^{*}},\ m_{3}=\frac{\delta _{1}}{\kappa \phi }$ and $\displaystyle m_{4}=\frac{\delta _{1}+\delta _{2}}{\eta }$. Then we have

$$\begin{aligned} \begin{aligned} \frac{\hbox {d}V_{2}}{\hbox {d}t}&=-d(S-S^{*})^{2}\\&\quad -d[(E-E^{*})+(A-A^{*})+(I-I^{*})+(H-H^{*})]^{2}\\&\quad -m_{2}\alpha _{0}(E-E^{*})^{2}\\&\quad -(\delta _{1}+m_{3}\alpha _{1})(A-A^{*})^{2}\\&\quad -(\delta _{2}+m_{4}\alpha _{2})(I-I^{*})^{2}\\&\quad -\delta _{3}(H-H^{*})^{2}-d(R-R^{*})^{2}\\&\quad -[2d+\delta _{1}+m_{1}\beta _{1}S^{*}](S-S^{*})(A-A^{*})\\&\quad -[2d+\delta _{2}+m_{1}\beta _{2}S^{*}](S-S^{*})(I-I^{*})\\&\quad -[2d+\delta _{3}+m_{1}\beta _{3}S^{*}](S-S^{*})(H-H^{*})\\&\quad +m_{2}\beta _{1}S(E-E^{*})(A-A^{*})\\&\quad +[-\delta _{2}+m_{2}\beta _{2}S+m_{4}\kappa \phi ](E-E^{*})(I-I^{*})\\&\quad +[-\delta _{3}+m_{2}\beta _{3}S](E-E^{*})(H-H^{*})\\&\quad -2d(E-E^{*})(R-R^{*})\\&\quad -(\delta _{1}+\delta _{3})(A-A^{*})(H-H^{*})\\&\quad -\delta _{1}(A-A^{*})(R-R^{*})\\&\quad -(\delta _{2}+\delta _{3})(I-I^{*})(H-H^{*})\\&\quad -\delta _{2}(I-I^{*})(R-R^{*})-\delta _{3}(H-H^{*})(R-R^{*})\\&\quad +\frac{p(I-I^{*})(Z-Z^{*})}{(1+qI)(1+qI^{*})}-a_{0}(Z-Z^{*})^{2} \end{aligned} \end{aligned}$$

So, $\displaystyle \frac{\hbox {d}V_{2}}{\hbox {d}t}<0$ in $\Phi$ when the following conditions hold:

• (i) $\displaystyle (\delta _{1}{+}m_{3}\alpha _{1}){>}\max \left\{ \frac{3}{2d}(2d{+}\delta _{1}{+}m_{1}\beta _{1}S^{*})^{2}, \frac{9\delta _{1}^{2}}{4d}\right\}$,

• (ii) $\displaystyle (\delta _{2}+m_{4}\alpha _{2})>\max \left\{ \frac{4}{d}(2d+\delta _{2}+m_{1}\beta _{2}S^{*})^{2},\frac{12\delta _{2}^{2}}{d},\ \frac{8(\delta _{2}+\delta _{3})^{2}}{\delta _{3}},\ \frac{2p^{4}d^{4}}{a_{0}(d+q\Lambda )^{4}}\right\}$,

• (iii) $\displaystyle \frac{48}{d}(2d+\delta _{3}+m_{1}\beta _{3}S^{*})^{2}<12\delta _{3}<d<\frac{m_{2}\alpha _{0}}{18}$.

Moreover, $\displaystyle \frac{\hbox {d}V_{2}}{\hbox {d}t}\bigg |_{E^{*}}=0$. So, by Lyapunov LaSalle’s theorem [27], ${E^{*}}$ is GAS in the interior of $\Phi$ when the mentioned restrictions are fulfilled. $\square$

7 Forward bifurcation at $R_{0} =1$

A unique endemic equilibrium point $E^{*}$ of system (1) exists when basic reproduction number ($R_{0}$) exceeds 1. Also, for $R_{0}$ less than unity, there is no endemic equilibrium point. Hence, there exists of a transcritical bifurcation around the disease-free equilibrium point $E_{0}$ when $R_{0}=1$.

Theorem 7.1

The system undergoes a transcritical bifurcation with respect to the bifurcation parameter $\beta _{1}$ around $\displaystyle E_{0}(S_{0},0,0,0,0,0,0)\equiv \left( \frac{\Lambda }{d},0,0,0,0,0,0\right)$ for $R_{0}=1$.

Proof

Consider the RHS of system (1) as: $g=(g_{1},g_{2},g_{3},g_{4},g_{5},g_{6},g_{7})^{T}$ where

$$\begin{aligned} g_{1}= & {} \Lambda -(\beta _{1}A+\beta _{2}I+\beta _{3}H)S-dS+\xi R-u_{1}kSZ; \\ g_{2}= & {} (\beta _{1}A+\beta _{2}I+\beta _{3}H)S-(\kappa +d)E;\\ g_{3}= & {} \kappa \phi E-(\eta +d+\delta _{1})A;\\ g_{4}= & {} \kappa (1-\phi )E+\eta A-(d+\delta _{2}+\omega )I;\\ g_{5}= & {} \omega I-(\gamma +d+\delta _{3})H;\\ g_{6}= & {} \gamma H-(d+\xi )R+u_{1}kSZ;\\ g_{7}= & {} \frac{pI}{1+qI}-a_{0}Z. \end{aligned}$$

For system (1), $\displaystyle R_{0}=\frac{\kappa \Lambda }{d\alpha _{0}\alpha _{1}\alpha _{2}\alpha _{3}}[\beta _{1}\phi \alpha _{2}\alpha _{3}+(\beta _{2}\alpha _{3}+\beta _{3}\omega )\{\phi \eta +\alpha _{1}(1-\phi )\}]$, where $\alpha _{0}=d+\kappa ,\ \alpha _{1}=\eta +d+\delta _{1},\ \alpha _{2}=\omega +d+\delta _{2},\ \alpha _{3}=\gamma +d+\delta _{3}$.

$$\begin{aligned} {\overline{J}}|_{E_{0}}= \left( \begin{array}{ccccccc} -d &{} 0 &{} -\beta _{1}S_{0} &{} -\beta _{2}S_{0} &{} -\beta _{3}S_{0} &{} \xi &{} -u_{1}kS_{0} \\ 0 &{} -\alpha _{0} &{} \beta _{1}S_{0} &{} \beta _{2}S_{0} &{} \beta _{3}S_{0} &{} 0 &{} 0 \\ 0 &{} \kappa \phi &{} -\alpha _{1} &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} \kappa (1-\phi ) &{} \eta &{} -\alpha _{2} &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} \omega &{} -\alpha _{3} &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} \gamma &{} -(d+\xi ) &{} u_{1}kS_{0} \\ 0 &{} 0 &{} 0 &{} p &{} 0 &{} 0 &{} -a_{0} \end{array} \right) \end{aligned}$$

The eigenvalues are $\lambda _{1}=-d,\ \lambda _{2}=-(d+\xi ),\ \lambda _{3}=-a_{0}$ and other four eigenvalues are roots of the equation: $\lambda ^{4}+Q_{1}\lambda ^{3}+Q_{2}\lambda ^{2}+Q_{3}\lambda +Q_{4}=0$,

where, $Q_{1}=\alpha _{0}+\alpha _{1}+\alpha _{2}+\alpha _{3}$, $Q_{2}=\alpha _{0}\alpha _{1}+\alpha _{2}\alpha _{3}+(\alpha _{0}+\alpha _{1})(\alpha _{2}+\alpha _{3})-\kappa S_{0}\{\beta _{1}\phi +\beta _{2}(1-\phi )\}$, $Q_{3}=\alpha _{0}\alpha _{1}(\alpha _{2}+\alpha _{3})+\alpha _{2}\alpha _{3}(\alpha _{0}+\alpha _{1})-S_{0}\kappa [\beta _{1}\phi (\alpha _{2}+\alpha _{3})+\beta _{2}\{\phi \eta +\alpha _{1}(1-\phi )\}+(\beta _{2}\alpha _{3}+\beta _{3}\omega )(1-\phi )]$ and $Q_{4}=\alpha _{0}\alpha _{1}\alpha _{2}\alpha _{3}(1-R_{0})$.

Let $\displaystyle \beta _{1[TC]}=\frac{1}{\phi \alpha _{2}\alpha _{3}}\left[ \frac{\hbox {d}\alpha _{0}\alpha _{1}\alpha _{2}\alpha _{3}}{\Lambda \kappa }-(\beta _{2}\alpha _{3}+\beta _{3}\omega )\{\phi \eta +\alpha _{1}(1-\phi )\}\right]$ be the value of $\beta _{1}$ such that $R_{0}=1$ which implies ${\overline{J}}|_{E_{0}}$ has a simple zero eigenvalue at $\beta _{1}=\beta _{1[TC]}.$ Also, $V= (v_{1},v_{2},v_{3},v_{4},v_{5},v_{6},v_{7})^{T}$ and $W= (w_{1},w_{2},w_{3},w_{4},w_{5},w_{6},w_{7})^{T}$, respectively, be the eigenvectors of $\displaystyle {\overline{J}}|_{E_{0}}$ and $\displaystyle {\overline{J}}|^{T}_{E_{0}}$ corresponding to the zero eigenvalue.

Calculations give $\displaystyle v_{1}=\kappa (\xi \gamma \omega a_{0}-u_{1}k\Lambda p\alpha _{3})\{\phi \eta +\alpha _{1}(1-\phi )\}-a_{0}\alpha _{0}\alpha _{1}\alpha _{2}\alpha _{3}(d+\xi ),\ v_{2}=da_{0}\alpha _{0}\alpha _{1}\alpha _{2}\alpha _{3}(d+\xi ),\ v_{3}=d\phi a_{0}\kappa \alpha _{2}\alpha _{3}(d+\xi ),\ v_{4}=da_{0}\alpha _{3}\kappa (d+\xi )\{\phi \eta +\alpha _{1}(1-\phi )\},\ v_{5}=\omega da_{0}\kappa (d+\xi )\{\phi \eta +\alpha _{1}(1-\phi )\},\ v_{6}=\kappa \{\phi \eta +\alpha _{1}(1-\phi )\}(\gamma \omega a_{0}d+u_{1}k\Lambda p\alpha _{3})$ and $v_{7}=pd \alpha _{3}\kappa (d+\xi )\{\phi \eta +\alpha _{1}(1-\phi )\}$. Also, $w_{1}=w_{6}=w_{7}=0,\ w_{2}=d\alpha _{1}\alpha _{2}\alpha _{3},\ w_{3}=\Lambda \{\beta _{1}\alpha _{2}\alpha _{3}+\eta (\beta _{2}\alpha _{3}+\beta _{3}\omega )\},\ w_{4}=\Lambda \alpha _{1}(\beta _{2}\alpha _{3}+\beta _{3}\omega )$ and $w_{5}=\Lambda \alpha _{1}\alpha _{2}\beta _{3}$. Therefore,

$$\begin{aligned} \begin{aligned} \Omega _{1}&= W^{T}.g_{\beta _{1}}(E_{0}, \beta _{1[TC]}) \\&=w_{2}(AS)\big |_{E_{0}}=0, \\ \Omega _{2}&= W^{T}\left[ Dg_{\beta _{1}}(E_{0}, \beta _{1[TC]})V\right] \\&=S_{0}v_{3}w_{2} \ne 0, \\ \Omega _{3}&= W^{T}\left[ D^{2}g(E_{0}, \beta _{1[TC]})(V,V)\right] \\&=2(\beta _{1}v_{3}+\beta _{2}v_{4}+\beta _{3}v_{5})v_{1}w_{2} \ne 0. \end{aligned} \end{aligned}$$

Thus, in view of Sotomayor’s Theorem, system (1) possesses a transcritical bifurcation around $E_{0}$ at $R_{0}=1$ taking $\beta _{1}$ as bifurcating parameter. $\square$

8 Numerical simulation without any control policy

Pictorial scenarios help us to understand system dynamics more clearly. The human population in India in April 2020 is about 136 million, the annual birth rate is 18.7 births/1000 people and the annual death rate is 7.3 deaths/1000 people. So, we are taking $S(0)=1.36\times 10^{9}$ and $\Lambda =7\times 10^{4}$ by applying unit conversion from year to day. And the death rate per day we get is near about 0.00002. For the sake of calculation, we are taking $d=0.000055$. From the data provided in the dashboard by the centre for system science and engineering (CSSE) at John Hopkins University on 13th April 2020, India has 9352 corona activated cases [6]. And till the date, total death cases are 324 and recovered cases are 980. Hence, unit conversion to day gives $\delta _{3}$ as $8\times 10^{-4}$ and $\xi$ as $3\times 10^{-3}$. As total active cases till 13th is 8048 among 9352 cases; so, we get $\omega$ as 0.02 approximately [29]. According to the current epidemic situation of coronavirus, the new human cases infected per unit time is denoted by $\beta SI\equiv (\beta _{1}A+\beta _{2}I+\beta _{3}H)S$. Human cases infected by COVID-19 in March was (I) is about 1393, the population in India (S) in April is approximately by $1.35\times 10^{9}$, the new human cases $(\beta _{2} SI)$ till 13th April is about 9352 [29], hence we have $\beta _{2} \approx 1.5\times 10^{-10}$ by doing the unit conversion from month to day. As per the data of March provided by Ministry of Health and Family Welfare, Government of India, the infected cases by COVID-19 is about 1393, so, for sake of simplicity I(0) is taken as 40. Now, all the assumed and estimated parameters are listed in Table 1. By CDC reports, mostly 25% of infected may not show any symptoms, i.e. remain as asymptomatic [7]. So, we have assumed $A(0)=10$. Also, $E(0)=10^{3},\ H(0)=15,\ R(0)=5$ and $Z(0)=0.5$ have been assumed.

Table 1 Parameter values used for numerical simulation of system (1)

Figure 2 shows that for the parametric values in Table 1 and $d=0.6$, the trajectory starting from mentioned initial point ultimately converges to DFE $E_{0}(1.273\times 10^{9},0,0,0,0,0,0)$ and as we get the basic reproduction number $R_{0}$ as 0.000013 here which lies below unity, so, the disease cannot invade in the system in this case.

Fig. 2

Stability of the populations around disease-free equilibrium $E_{0}$

Now if we start to decrease the value of d, then for $d=0.000055$ along with parametric values in Table 1, the trajectory starting from mentioned initial point approaches to unique endemic equilibrium point $E^{*}(1.3788\times 10^{8}, 2.2419\times 10^{6}, 895850.6, 6.5435\times 10^{7}, 2.6129\times 10^{6}, 4.2916\times 10^{8}, 0.1667)$ with time (see Fig. 3). For these parametric values we get $R_{0}=9.2305>1$ indicating the presence of infection in the system.

Fig. 3

Stability of the populations around endemic equilibrium $E^{*}$

Now $E_{0}$ changes its stability when d comes below of a threshold value $d_{[TC]}$ and becomes stable for $d>d_{[TC]}$. So, the system undergoes a transcritical bifurcation at $d=d_{[TC]}=0.0004965$ around DFE $(E_{0})$ (see Fig. 4). Besides of the mortality rate, disease transmission rate from symptomatically infected individuals to susceptible $(\beta _{2})$ and hospitalisation rate of symptomatically infected $(\omega )$ also play key roles to control the system dynamics. Figure 5a, b show the system possesses transcritical bifurcations around $E_{0}$ at $\beta _{2}=\beta _{2[TC]}=1.50478\times 10^{-11}$ and $\omega =\omega _{[TC]}=0.2018$, respectively.

Fig. 4

Trancritical bifurcation around $E_{0}$ taking d as bifurcation parameter

Fig. 5

Trancritical bifurcation around $E_{0}$ taking a $\beta _{2}$ and b $\omega$ as bifurcation parameters

Figure 6 demonstrates the sensibility of some of the vital parameters on disease transmission. Figure 6a shows that $\beta _{2}$ is most sensitive to control the transmission of the disease than $\beta _{1},\ \beta _{3}$ and $\kappa$. A small increase in $\beta _{2}$ can increase the value of $R_{0}$ significantly. On the other hand, $\omega$ is inversely proportional with $R_{0}$, i.e. if more people enter the observation centre or hospital, then the disease fatality starts to decrease with time.

Fig. 6

Relationship between basic reproduction number $R_{0}$ with $\beta _{1},\ \beta _{2},\ \beta _{3},\ \kappa$ and $\omega$

Here, k denotes the information interaction rate which brings behavioural changes in susceptible individuals. People become more cautious when the disease starts to outbreak at higher rate. At early stage, change in information density does not make any significant impact but later, increasing value of k decreases the number of symptomatically infected individuals and it is observed in Fig. 7.

Fig. 7

Trajectory profiles of symptomatically infected population (I) for different values of k

In Fig. 8, the impact of the growth of information and the ‘level for unresponsiveness towards information’ on the infected population have been observed. Here, p and q represent the ‘growth rate of information’ and the ‘level for unresponsiveness towards information’, respectively. Figure 8a depicts that increase in information can lower the infected population with time as people can successfully save themselves from getting infected by induced behavioural changes. On the other hand, decreasing the level of people unresponsiveness towards information (smaller value of q) ultimately decreases the symptomatically infected individuals with time (see Fig. 8b). It is known that increasing information reduce the infected population but if the response rate $(u_{1})$ also starts to increase along with information, then the rate of decrease is higher. Figure 9 depicts that symptomatically infected population decrease significantly for increasing value of $u_{1}$ for increasing p.

9 Optimal control problem

We formulate the corresponding optimal control problem here to observe how suitable control interventions reduce the disease burden on the population. (a) The awareness programs among susceptible and symptomatically infected people regarding the information about COVID-19 and its symptoms so that the susceptible can take precautions and infected can admit to hospitals without neglecting the symptoms and (b) better treatment policies on hospitalised people have been considered as the control policies. We have analysed analytically and also numerically how these control policies make their impact on disease transmission and try to optimise the cost burden for their implementations.

(i) Increase the awareness among susceptible individuals and symptomatically infected individuals through information: Susceptible individuals start to become aware of a disease and its prevention when they are provided with enough information and this results in behavioural changes in population. The awareness programs to spread the information regarding COVID-19 outbreak has been considered as a possible tool to activate the sensibility of those individuals who live in a susceptible environment. These days Government and media sources have spread the news about this disease fatality regularly. And due to the regular broadcast, the people have started to take protective measures at a higher rate by maintaining social distances and proper hygiene, staying at isolation and even adopting the self-quarantine strategy. In system (1), $u_{1}$ represents the intensity of response through information with the restriction $0\le u_{1}\le 1.$ Here, 0 and 1, respectively, denotes no response and the full response of the informed population. As a consequence, $u_1(t)$ changes according to the individual’s behavioural response and we have taken this response intensity as one of the control variables. Incurred cost is involved as a nonlinear function of $u_{1}$ to stimulate the response of individuals and their behavioural changes. If the information starts to spread at a higher rate, then it may help to find the optimal response of susceptible individuals. Moreover, these awareness programs are conducted to aware not only the susceptible population but infected people too. It is a person’s responsibility to consult a medical person or to admit to a hospital if the slightest symptom is shown in his body. Anyone should not ignore by assuming it as a simple cough and cold case. By admitting to the hospital at an early stage can also decrease the disease burden and so, a saturated hospitalisation rate function $\displaystyle \frac{\epsilon _{1}u_{2}I}{1+\zeta _{1}I}$ is incorporated in the system where $\epsilon _{1}$ denotes the rate at which symptomatic infected individuals move to hospitals without ignoring the symptoms with intensity $u_{2}$ and saturation constant is $\zeta _{1}^{-1}$. The costs incurred in hospitalisation, medicines, etc., during the time a person admitted to a hospital or isolation centre is taken into consideration. The awareness intensity $u_{2}$ is taken as another control variable with restriction $0\le u_{2}\le 1$ where 1 denotes when a person admits into a hospital without ignoring the symptoms whenever feel sick. And 0 denotes the case when a person becomes ignorant about his sickness and does not consult a doctor.

Fig. 8

Trajectory profiles of symptomatically infected population (I) for different values of a p and b q

Fig. 9

Variation of symptomatically infected population (I) due to change in growth of information, p for different values of $u_{1}$

(ii) Better treatment policy for hospitalised individual: Providing proper and better antidote of a virus to the hospitalised people at an early stage of infection can lower the disease fatality. It affects disease progression too. It is considered that the treatment which is available and provided to the individuals admitted in the hospitals is of limited quantity. The resource availabilities depend on medical diagnosis, financial stability, treatment, etc., and all these things are limited. Considering this fact, a saturated treatment rate function $\displaystyle \frac{\epsilon _{2}u_{3}H}{1+\zeta _{2}H}$ is incorporated in the system where $\epsilon _{2}$ denotes treatment rate with intensity $u_{3}$ and saturation constant is $\zeta _{2}^{-1}$. The cost incurred in vaccination, medicines, diagnosis, health care, hospitalisation, etc., at the time of treatment period is taken into consideration. The treatment intensity $u_{3}$ is taken as another control variable with restriction $0\le u_{3}\le 1$ where 0 and 1 denote the no response and full response to the given treatment, respectively.

The main work is to determine optimal response intensity and optimal treatment with minimum cost by the help of provided information. So, the region for the control interventions $u_{1}(t),\ u_{2}(t)$ and $u_{3}(t)$ is given as:

$$\begin{aligned} \Pi= & {} \left\{ \left( u_{1}(t),u_{2}(t),u_{3}(t)\right) |\left( u_{1}(t), u_{2}(t),u_{3}(t)\right) \right. \\&\quad \left. \in [0,1]\times [0,1]\times [0,1],t\in [0,T_{f}]\right\} , \end{aligned}$$

where $T_{f}$ is the final time up to which the control policies are executed, and also $u_{i}(t)$ for $i=1,2,3$ are measurable and bounded functions.

9.1 Determination of total cost

We determine the incurred cost which needs to be minimised in order to apply control interventions.

(i) Cost involved to spread awareness among susceptible and symptomatic infected people: The total cost incurred during awareness spreading programs among people is given as:

$$\begin{aligned} \int _{0}^{T_{f}} \left[ w_{1}I(t)+w_{3}u^{2}_{1}(t)+w_{4}u^{2}_{2}(t)\right] \mathrm{d}t \end{aligned}$$

The cost for spreading awareness among susceptible regarding the disease, necessary precautions and prevention by maintaining social distance and hygiene via social campaigns, newspapers, television, social networks, etc., is denoted by $w_{3}u^{2}_{1}(t)$. The term considers the cost of associated efforts to realise the individual about the importance of maintaining social distances and it is observed that this cost is high enough. Moreover, the cost incurred when symptomatically infected individual consults a medical person and admits to a hospital without ignoring the symptoms is represented by the term $w_{1}I(t)+w_{4}u^{2}_{2}(t)$. This cost includes the expenditure of hospitalisation, medicines, etc., and $w_{4}u^{2}_{2}(t)$ considers the productivity loss due to illness. There are some studies already exist which reveal the impact of the cost associated with awareness programs, screening and self-protective measures and nonlinearity up to order two have been taken [2, 4, 25]. We, in this work, now emphasise on the fact how the awareness programs and practical cautions reduce the disease burden at the time of this epidemic outbreak.

(ii) Cost involved during treatment at hospitals: Total cost associated with treatment for symptomatic infected individual is given as:

$$\begin{aligned} \int _{0}^{T_{f}} \left[ w_{2}H(t)+w_{5}u^{2}_{3}(t)\right] \mathrm{d}t, \end{aligned}$$

where $w_{2}H(t)$ and $w_{5}u^{2}_{3}(t)$, respectively, denote the cost associated with hospitalised population for losing man power [17, 23, 25] and the cost at the time of treatment due to diagnosis charges, the expenditure of hospitalisation, etc. The later term considers the opportunity loss including productivity loss due to admittance to the hospital. So, the nonlinearity of $u_{3}(t)$ is taken up to order two for treatment policy [17, 23, 25].

The following control problem is considered based on previous discussions along with the cost functional J to be minimised:

$$\begin{aligned}&J[u_{1}(t),u_{2}(t),u_{3}(t)] \nonumber \\&\quad =\int _{0}^{T_{f}} \left[ w_{1}I(t)+w_{2}H(t)+w_{3}u^{2}_{1}(t) \right. \nonumber \\&\qquad \left. +w_{4}u^{2}_{2}(t)+w_{5}u^{2}_{3}(t)\right] \mathrm{d}t \end{aligned}$$

(5)

subject to the model system:

$$\begin{aligned} \frac{\hbox {d}S}{\hbox {d}t}= & {} \Lambda -(\beta _{1}A+\beta _{2}I+\beta _{3}H)S-dS\nonumber \\&\quad +\xi R-u_{1}(t)kSZ, \nonumber \\ \frac{\hbox {d}E}{\hbox {d}t}= & {} (\beta _{1}A+\beta _{2}I+\beta _{3}H)S-(\kappa +d)E, \nonumber \\ \frac{\hbox {d}A}{\hbox {d}t}= & {} \kappa \phi E-\eta A-(d+\delta _{1})A, \nonumber \\ \frac{\hbox {d}I}{\hbox {d}t}= & {} \kappa (1-\phi )E+\eta A-(\omega +d+\delta _{2})I \nonumber \\&\quad -\frac{\epsilon _{1} u_{2}(t)I}{1+\zeta _{1}I}, \nonumber \\ \frac{\hbox {d}H}{\hbox {d}t}= & {} \omega I-\gamma H-(d+\delta _{3})H+\frac{\epsilon _{1} u_{2}(t)I}{1+\zeta _{1}I} \nonumber \\&\quad - \frac{\epsilon _{2}u_{3}(t)H}{1+\zeta _{2}H}, \nonumber \\ \frac{\hbox {d}R}{\hbox {d}t}= & {} \gamma H-dR-\xi R+u_{1}(t)kSZ \nonumber \\&\quad +\frac{\epsilon _{2}u_{3}(t)H}{1+\zeta _{2}H}, \nonumber \\ \frac{\hbox {d}Z}{\hbox {d}t}= & {} \frac{pI}{1+qI}-a_{0}Z, \end{aligned}$$

(6)

with initial conditions $S(0)>0,\ E(0)\ge 0,\ A(0)\ge 0,\ I(0)\ge 0,\ H(0)\ge 0,\ R(0)\ge 0$ and $Z(0)\ge 0.$ Here, the functional J denotes the total incurred cost as stated and the integrand:

$$\begin{aligned}&L(S,E,A,I,H,R,Z,u_{1}(t),u_{2}(t),u_{3}(t)) \\&\quad =w_{1}I(t)+w_{2}H(t)+ w_{3}u^{2}_{1}(t) \\&\qquad +w_{4}u^{2}_{2}(t)+w_{5}u^{2}_{3}(t) \end{aligned}$$

denotes the cost at time t. Positive parameters $w_{1},w_{2},w_{3},w_{4}\ \text {and}\ w_{5}$ are weight constants balancing the units of the integrand [17, 25]. The optimal control interventions $u^{*}_{1},\ u^{*}_{2}$ and $u^{*}_{3}$, exist in $\Pi$, mainly minimise the cost functional J.

Theorem 9.1

The optimal control interventions $u^{*}_{1},\ u^{*}_{2}$ and $u^{*}_{3}$ in $\Pi$ of the control system (5)–(6) exist such that $J(u^{*}_{1},u^{*}_{2},u^{*}_{3})=\min [J(u_{1},u_{2},u_{3})]$.

Proof

Proof is done in “Appendix”. $\square$

Pontryagin’s Maximum Principle helps to obtain optimal controls $u^{*}_{1},\ u^{*}_{2}$ and $u^{*}_{3}$ of system (5)–(6).

Theorem 9.2

If $u^{*}_{1},\ u^{*}_{2}$ and $u^{*}_{3}$ are the optimal control variables and $S^{*}, E^{*}, A^{*}, I^{*}, H^{*}, R^{*}, Z^{*}$ are corresponding optimal state variables of the control system (5)–(6), then there exist adjoint variables $\lambda =\left( \lambda _{1},\lambda _{2},\ldots , \lambda _{7}\right) \in {\mathbb {R}}^{7}$ satisfying the canonical equations:

$$\begin{aligned} \begin{aligned} \frac{\mathrm{d}\lambda _{1}}{\mathrm{d}t}&=\lambda _{1}[\beta _{1}A+\beta _{2}I+\beta _{3}H+d+u_{1}kZ]\\&\quad -\lambda _{2}(\beta _{1}A+\beta _{2}I+\beta _{3}H)-\lambda _{6}(u_{1}kZ)\\ \frac{\mathrm{d}\lambda _{2}}{\mathrm{d}t}&=\lambda _{2}(d+\kappa )-\lambda _{3}(\kappa \phi ) -\lambda _{4}\{\kappa (1-\phi )\} \\ \frac{\mathrm{d}\lambda _{3}}{\mathrm{d}t}&=\lambda _{1}(\beta _{1}S)-\lambda _{2}(\beta _{1}S)\\&\quad +\lambda _{3}(\eta +d+\delta _{1})-\lambda _{4}(\eta ) \\ \frac{\mathrm{d}\lambda _{4}}{\mathrm{d}t}&=-w_{1}+\lambda _{1}(\beta _{2}S)-\lambda _{2}(\beta _{2}S)\\&\quad +\lambda _{4} \left\{ \omega +d+\delta _{2}+\frac{\epsilon _{1}u_{2}}{(1+\zeta _{1}I)^{2}}\right\} \\&\quad -\lambda _{5} \left\{ \omega +\frac{\epsilon _{1}u_{2}}{(1+\zeta _{1}I)^{2}}\right\} \\&-\lambda _{7}\left\{ \frac{p}{(1+qI)^{2}}\right\} \\ \frac{\mathrm{d}\lambda _{5}}{\mathrm{d}t}&=-w_{2}+\lambda _{1}(\beta _{3}S)-\lambda _{2}(\beta _{3}S)\\&\quad +\lambda _{5} \left\{ \gamma +d+\delta _{3}+\frac{\epsilon _{2}u_{3}}{(1+\zeta _{2}H)^{2}}\right\} \\&\quad -\lambda _{6} \left\{ \gamma +\frac{\epsilon _{2}u_{3}}{(1+\zeta _{2}H)^{2}}\right\} \\ \frac{\mathrm{d}\lambda _{6}}{\mathrm{d}t}&=-\lambda _{1}(\xi )+\lambda _{6}(d+\xi ) \\ \frac{\mathrm{d}\lambda _{7}}{\mathrm{d}t}&=\lambda _{1}(u_{1}kS)-\lambda _{6}(u_{1}kS)+\lambda _{7}(a_{0}) \end{aligned} \end{aligned}$$

(7)

with transversality conditions $\lambda _{i}(T_{f})=0$ for $i=1,2,\ldots ,7$ and corresponding optimal controls $u^{*}_{1},\ u^{*}_{2}$ and $u^{*}_{3}$ are given as:

$$\begin{aligned} \begin{aligned} u^{*}_{1}&= \min \left\{ \max \left\{ 0, \left( \frac{kS^{*}Z^{*}}{2w_{3}}\left( \lambda _{1}-\lambda _{6}\right) \right) \right\} ,1\right\} , \\ u^{*}_{2}&= \min \left\{ \max \left\{ 0, \left( \frac{\epsilon _{1}I^{*}}{2w_{4}\left( 1+\zeta _{1}I^{*}\right) } \left( \lambda _{4}-\lambda _{5}\right) \right) \right\} ,1\right\} , \\ and \qquad u^{*}_{3}&= \min \left\{ \max \left\{ 0, \left( \frac{\epsilon _{2}H^{*}}{2w_{5}\left( 1+\zeta _{2}H^{*}\right) } \left( \lambda _{5}-\lambda _{6}\right) \right) \right\} ,1\right\} . \end{aligned}\nonumber \\ \end{aligned}$$

(8)

Proof

Proof is given in “Appendix”. $\square$

10 Numerical results with control policies

In system (6), different control strategies have been applied to reduce the disease burden and to minimise the total cost by finding the optimal control paths. The growth of information varies with time as it depends on disease fatality and behavioural response. So, $u_{1}$ is taken as a control variable. Also, there exist two other control variables $u_{2}$ and $u_{3}$ representing hospitalisation of symptomatically infected without neglecting symptoms and better treatment of hospitalised people, respectively. The positive weights are taken as $w_{1}=1,\ w_{2}=1,\ w_{3}=2500,\ w_{4}=10$ and $w_{5}=100$ [17, 25]. In order to draw the numerical figures, we slightly adjust $\beta _{i}$ for $i=1,2,3$ and $\phi$ and all the parameters are listed in Table 2. The effects of implementation of one or all control strategies to find the minimal cost have been analysed one by one here. Corresponding control system in Eqs. (5)–(6) is solved here with the initial population size: $S(0)=135\times 10^{4},\ E(0)=500,\ A(0)=10,\ I(0)=40,\ H(0)=15,\ R(0)=5$ and $Z(0)=1$. Four different cases are considered: (i) when only $u_{1}$ is applied, (ii) when $u_{1}$ and $u_{2}$ are applied, (iii) when $u_{1}$ and $u_{3}$ are applied and (iv) when all control policies are applied. The numerical simulation for all cases is obtained by MATLAB. The optimal control variables are found by Forward-backward sweep method where the optimal state system and the adjoint state system are solved by forward and backward in time, respectively. In the next step, the steepest descent method is used to update the optimal controls by Hamiltonian for the optimality of the system [26] and the process continues until the convergence. It is assumed that the control policies are applied for approximately $T_{f}=90$ days. In India, the first case of COVID-19 was registered on 30 January 2020 and still now the number of affected cases are increasing day by day. The first lockdown is announced in March and observing the severity, the Government has extended its duration. The people are strictly advised to maintain social distancing and proper hygiene in order to keep themselves safe. As there is no vaccine discovered still now, so, natural immunity and physical distancing are the only way-outs to avoid from being infected. We have taken (i) awareness among susceptible individuals $(u_{1})$, (ii) awareness among infected individuals $(u_{2})$ and (iii) better treatment $(u_{3})$ as our control policies and looking at the current situation, it is observed that people need to maintain these precautions for quite a long time.

Table 2 Parametric values used in model system (6)

Fig. 10

Profiles of populations with applied optimal control $u^{*}_{1}$ only and $u_{2}=u_{3}=0$

Fig. 11

Profiles of populations with applied optimal control $u^{*}_{1}$ only and $u_{2}=u_{3}=0$

First, we consider the case when people only adopt behavioural changes in a susceptible environment. Figure 10 depicts the population profiles when $u_{1}=u_{1}^{*}$ and $u_{2}=u_{3}=0$. At $T_{f}=90$, the population become (7,647,536.63, 6.2058, 6.0621, 203.3930, 13.79, 561.30, 0.169704). The susceptible population increases with time. Both asymptomatically and symptomatically infected population increase steeply within 7–10 days and reach their maximum values but after that, the slope of the trajectories start to decrease. It is noted that the rate of declination is higher for asymptomatically infected people. Also, the number of asymptomatically infected people is lower than the number of symptomatically infected people. Hospitalised people also decreases with time almost after the first 2 weeks as recovered people increases. The corresponding optimal control path is given in Fig. 11. The intensity of the control variable works with higher intensity at earlier days but later it decreases with time. The declination of the graph may be caused because of people ignorance, etc.

Fig. 12

Profiles of populations with applied optimal controls $u^{*}_{1}$ and $u^{*}_{2}$ only and $u_{3}=0$

Fig. 13

Optimal controls $u^{*}_{1}$ and $u^{*}_{2}$ when $u_{3}=0$

Next, we consider the case when susceptible people take precautions to avoid being infected $(u_{1})$ and also symptomatically infected individuals enter into hospitals even when they are shown slightest symptoms $(u_{2})$. Figure 12 shows the population trajectories when $u_{1}=u_{1}^{*}$ and $u_{2}=u_{2}^{*}$ but $u_{3}=0$. At $T_{f}=90$, the population becomes (7,647,620.98, 3.2498, 3.3487, 106.0726, 8.60, 587.64, 0.169182). It is observed that the number of symptomatically infected people decreases significantly in this case. Also, the recovered population increases here. Both the infected population first increase and then decrease after almost 1 week. Figure 13 shows the optimal control paths of $u_{1}$ and $u_{2}$. From this figure, it is observed that $u_{1}$ works with lesser intensity than $u_{2}$. Moreover, the intensity of $u_{1}$ gradually decreases with time. But $u_{2}$ works with its maximum intensity for quite a long time and then decreases suddenly.

Now we consider the case when people induce behavioural changes to protect themselves in susceptible stage $(u_{1})$ and also hospitalised people are provided with better treatment $(u_{3})$. Figure 14 depicts the population profiles for $u_{1}=u_{1}^{*}$ and $u_{3}=u_{3}^{*}$ but $u_{2}=0$. At $T_{f}=90$, the population becomes (7,647,547.32, 6.1961, 6.0521, 203.2271, 11.7627, 553.23, 0.169703). Implementation of $u_{3}$ increases susceptible population as more people are recovered in this case and so moved to susceptible class further. Corresponding optimal control paths are depicted in Fig. 15. The intensity of the control variable representing behavioural response decreases continuously with time and the intensity itself is not much higher. On the other hand, the intensity of $u_{3}$ first increases for almost 3 weeks and then decreases with a slower rate for next 2 months and then suddenly decreases. It represents that the treatment works with higher intensity at the earlier state and with time the intensity decreases when people become aware of the disease and its precautions.

Fig. 14

Profiles of populations with applied optimal controls $u^{*}_{1}$ and $u^{*}_{3}$ only and $u_{2}=0$

Fig. 15

Optimal controls $u^{*}_{1}$ and $u^{*}_{3}$ when $u_{2}=0$

It is obvious that implementing all the control policies is beneficial for the proposed system. So, we consider the combination of all three control policies in the system, i.e. a system where people induce behavioural changes with time to protect themselves from infection, symptomatically infected people move to a hospital without neglecting symptoms and better treatment is applied to hospitalised people. Figure 16 depicts the population trajectories in the presence of all control policies. At $T_{f}=90$, population becomes (7,647,630.40, 3.2391, 3.3381, 105.8244, 7.2543, 580.2986, 0.169180). The susceptible population increases at the highest level in this case. Also, both the infected population decreases at its lowest population here. Moreover, the recovered population is higher than the case when only $u_{1}$ is applied. Figure 17 shows the paths of optimal control strategies. The control policy on behavioural response works with higher intensity for few days but after that, it starts to decrease. It is justifiable as people become curious when epidemic or pandemic first outbreaks and they try to adopt changes in behaviour based on the information received but later it turns into their usual habit. Moreover, the infected people, if take the slightest symptoms seriously and consult medical person immediately, then it can reduce the disease also. Again, if the hospitalised people are provided with a proper antidote and better treatment, it helps to increase the number of recovered people. Here, $u_{2}$ acts with maximum intensity for a long time and then it decreases. $u_{3}$ also acts with its higher intensity for almost one and a half month following a declination at a very slow rate.

Fig. 16

Profiles of populations with both optimal control policies $u^{*}_{1},\ u^{*}_{2}$ and $u^{*}_{3}$

Fig. 17

Profiles of optimal controls $u^*_1,\ u^*_2$ and $u^*_3$

Fig. 18

a Cost distribution in presence and absence of control policies. b Profiles of symptomatic infected population under different control policies

Fig. 19

a Profiles of symptomatic infective population for various $\zeta _{1}$ with $u^{*}_{1},\ u^{*}_{2}$ and $u^{*}_{3}$. b Profiles of cost for various $\zeta _{1}$ with $u^{*}_{1},\ u^{*}_{2}$ and $u^{*}_{3}$

Fig. 20

a Plots of control $u^{*}_{1}$ for various $\zeta _{1}$. b Plots of control $u^{*}_{2}$ for various $\zeta _{1}$. c Plots of control $u^{*}_{3}$ for various $\zeta _{1}$

Fig. 21

a Profiles of infective population for various $\epsilon _{1}$ with $u^{*}_{1},\ u^{*}_{2}$ and $u^{*}_{3}$. b Profiles of cost for various $\epsilon _{1}$ with $u^{*}_{1},\ u^{*}_{2}$ and $u^{*}_{3}$

Fig. 22

a Plots of control $u^*_1$ for various $\epsilon _{1}$. b Plots of control $u^*_2$ for various $\epsilon _{1}$. c Plots of control $u^*_3$ for various $\epsilon _{1}$

In Fig. 18, cost design analysis has been performed in the absence and presence of $u_{2}$ and $u_{3}$ but we have considered $u_{1}$ in both the cases. There is one case where all three control policies are applied and the next case is when only the behavioural response is considered. Optimal cost profiles are shown in Fig. 18a for the cases and trajectory profiles for symptomatically infected individuals are depicted in Fig. 18b. In the absence of $u_{2}$ and $u_{3}$, cost occurs due to productivity loss by infected only. So, the opportunity loss is higher due to an epidemic outbreak and overall infected population increases in this case. On the other hand, the optimal cost is lower when all control interventions are applied. And as the infected population is lower in this case, it reduces the cost incurred because of opportunity loss.

10.1 Effect of hospitalisation rate and saturation constant on optimal control policies

If more symptomatically infected individuals admit into a hospital, then we observe how the hospitalisation rate and saturation constant have effects on the disease dynamics in the presence of the control policies with optimal intensities. Saturation rate $(\zeta _{1}^{-1})$ and the hospitalisation rate $(\epsilon _{1})$ have been varied in the following figures. Figure 19 shows the graphs of symptomatic infective population (I) and corresponding cost for different values of $\zeta _{1}$. For smaller saturation rate (larger values of $\zeta _{1}$), infective individuals and associated cost increase due to increase in productivity loss. Corresponding optimal control paths are drawn in Fig. 20 from which it is observed that for increasing value of $\zeta _{1}$, the control policy denoting behavioural response in susceptible individuals acts with a higher intensity which means smaller saturation rate increases the time period during which the control policy is implemented with higher intensity successfully. On the other hand, decreasing value of $\zeta _{1}$ increases the intensity level of the other two control interventions. Thus, if the saturation constant is high enough, then it is economically viable and hence, requires comparatively lesser efforts while implementing the control policies.

Further with the increase in treatment rate $(\epsilon _{1})$ from 0.01 to 0.1, the both symptomatic infective population and associated cost decrease with time (Fig. 21). From Fig. 22, it is observed that for increasing value of $\epsilon _{1}$, the control policy denoting behavioural response in susceptible people works for a smaller time period. Also, increasing value of $\epsilon _{1}$ increases the intensity level of the other two control variables. In Fig. 22b, c it is observed that the higher hospitalisation rate increases the lengths of maximum intensities of these two optimal control implementation periods. It means lesser efforts on these two applied controls are sufficient to reduce the overall infective population.

11 Conclusion

Coronavirus or Covid-19 first appeared in China in December 2019, but today it has spread all over the world in the form of a pandemic. Reports of the current situation reveal that almost 30 lakhs people are infected with the virus worldwide. Though the Governments and medical persons of each and every country are trying to provide protective measures to people, the infection rate is still high enough as the proper antidote for this virus is still unknown. According to the data till 26th April 2020, provided by the dashboard of CSSE at John Hopkins University, US has the highest confirmed cases with the number 939, 249 [6]. Though according to the official reports, almost 26, 384 people have died in Italy which is the highest in number among 185 countries or religions. If we consider the current situation in India, then from the reports of NIC, India, there are 20, 177 confirmed cases, 826 death cases and 5, 914 recovered cases are reported till 26th April [29]. Keeping this pandemic situation in mind, in this work, we have formulated a compartmental SEIRS model of Covid-19 where a separate equation is incorporated to reflect the information which induces behavioural changes. The growth rate of information is based on symptomatically infected individuals, awareness programs, social activities, etc. Positivity and boundedness of system variables guarantees that the proposed system is well-defined. Feasibility conditions of equilibrium points show that DFE exists for all parametric values where the unique endemic equilibrium point exists only when basic reproduction number $R_{0}$ exceeds unity. Both local and global stability conditions have been derived in Sect. 6. As we have only one endemic equilibrium point for $R_{0}>1$ and there does not exist any endemic point for $R_{0}<1$, hence, it is concluded that the system undergoes a forward (transcritical) bifurcation around the disease-free equilibrium. If the people start to take the information regarding the virus seriously and take its precautions, then it is obvious that overall infected population is lesser and it is shown in Fig. 7. Moreover, increasing values of the growth of information and decreasing ‘level for unresponsiveness towards information’ also help to reduce the infected population. Also, symptomatically infected population with respect to the growth rate of information decreases for increasing value of $u_{1}$, i.e. if people respond to the growing information regarding the awareness at a higher rate, then infected population decreases significantly. As the behavioural changes include keeping social distances, maintaining proper hygiene, staying at isolation and adopting the self-quarantine method, so, these precautions can really be useful to prevent the disease transmission at a higher rate.

In the later part, a corresponding optimal control problem is considered. Implementation of control interventions helps to reduce the disease burden. The behavioural changes in susceptible population changes with time and so, it is considered as one of the control policy. Further, the symptomatically infected people can also become cautious by the current disease fatality and may consult doctors or admit to hospitals if slightest symptoms are shown. Again, during the treatment period, better and proper medicines or diagnosis can be provided to a hospitalised person. So, all these things can be considered as control strategies. Numerical figures show that the control presenting behavioural response $(u_{1})$ works with higher intensity immediately after implementation but gradually it decreases with time. On the other hand, the control policy denoting hospitalisation of infected individuals $(u_{2})$ works with its maximum intensity for quite a long time period and then it decreases. Again, the control presenting better treatment of hospitalised people works with higher intensity for almost 2 months following a declination at a later stage, though this intensity is lower than the intensity of $u_{2}$. Implementation of a single strategy is useful but all the three control policies together can reduce the infected population and increase the recovered population at a higher rate. Hence, applying all the control policies together may help to reduce disease transmission at this current epidemic situation.

Appendix

1.1 Local stability of the endemic equilibrium point

Proof of Theorem 6.3

The Jacobian matrix at endemic equilibrium point $E^{*}$ is given as

$$\begin{aligned} {\overline{J}}=\left( \begin{array}{ccccccc} a_{11} &{} 0 &{} a_{13} &{} a_{14} &{} a_{15} &{} a_{16} &{} a_{17} \\ a_{21} &{} a_{22} &{} a_{23} &{} a_{24} &{} a_{25} &{} 0 &{} 0 \\ 0 &{} a_{32} &{} a_{33} &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} a_{42} &{} a_{43} &{} a_{44} &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} a_{54} &{} a_{55} &{} 0 &{} 0 \\ a_{61} &{} 0 &{} 0 &{} 0 &{} a_{65} &{} a_{66} &{} a_{67} \\ 0 &{} 0 &{} 0 &{} a_{74} &{} 0 &{} 0 &{} a_{77} \end{array} \right) \end{aligned}$$

where $a_{11}=-(\beta _{1}A^{*}+\beta _{2}I^{*}+\beta _{3}H^{*})-d-u_{1}kZ^{*},\ a_{13}=-\beta _{1}S^{*},\ a_{14}=-\beta _{2}S^{*},\ a_{15}=-\beta _{3}S^{*},\ a_{16}=\xi ,\ a_{17}=u_{1}kS^{*},\ a_{21}=\beta _{1}A^{*}+\beta _{2}I^{*}+\beta _{3}H^{*},\ a_{22}=-\alpha _{0},\ a_{23}=\beta _{1}S^{*},\ a_{24}=\beta _{2}S^{*},\ a_{25}=\beta _{3}S^{*},\ a_{32}=\kappa \phi ,\ a_{33}=-\alpha _{1},\ a_{42}=\kappa (1-\phi ),\ a_{43}=\eta ,\ a_{44}=-\alpha _{2},\ a_{54}=\omega ,\ a_{55}=-\alpha _{3},\ a_{61}=u_{1}kZ^{*},\ a_{65}=\gamma ,\ a_{66}=-(d+\xi ),\ a_{67}=u_{1}kS^{*},\ a_{74}=\frac{p}{(1+qI^{*})^{2}}$ and $a_{77}=-a_{0}$.

Characteristic equation of ${\overline{J}}|_{E^{*}}$ is $\lambda ^{7}+R_{1}\lambda ^{6}+R_{2}\lambda ^{5}+R_{3}\lambda ^{4}+R_{4}\lambda ^{3}+R_{5}\lambda ^{2}+R_{6}\lambda +R_{7}=0,$

$$\begin{aligned} R_{1}&=-(a_{11}+a_{22}+a_{33}+a_{44}+a_{55}+a_{66}+a_{77}),\\ R_{2}&=a_{11}(a_{22}+a_{33}+a_{44}+a_{55}+a_{66}+a_{77}) \\&\quad +a_{22}(a_{33}+a_{44}+a_{55}+a_{66}+a_{77})\\&\quad +a_{33}(a_{44}+a_{55}+a_{66}+a_{77})\\&\quad +a_{44}(a_{55}+a_{66}+a_{77})+a_{55}(a_{66}+a_{77})+a_{66}a_{77}\\&\quad -a_{16}a_{61}-a_{23}a_{32}-a_{24}a_{42} \\ R_{3}&=-[a_{11}a_{22}(a_{33}+a_{44}+a_{55}+a_{66}+a_{77})\\&\qquad -a_{11}a_{23}a_{32}-a_{11}a_{24}a_{42}+a_{11}a_{55}(a_{66}+a_{77})\\&\qquad +a_{11}a_{33}(a_{44}+a_{55}+a_{66}+a_{77})\\&\qquad +a_{11}a_{44}(a_{55} +a_{66}+a_{77})\\&\qquad +a_{11}a_{66}a_{77}+a_{13}a_{21}a_{32}+a_{14}a_{42}a_{21}\\&\qquad -a_{16}a_{61}(a_{22}+a_{33}+a_{44}+a_{55}+a_{77})\\&\qquad +a_{22}a_{33} (a_{44}+a_{55}+a_{66}+a_{77})\\&\qquad +a_{22}a_{44}(a_{55}+a_{66}+a_{77})\\&\qquad +a_{22}a_{55}(a_{66}+a_{77})+a_{22}a_{66}a_{77}\\&\qquad -a_{23}a_{32} (a_{44}+a_{55}+a_{66}+a_{77})+a_{43}a_{32}a_{24}\\&\qquad -a_{24}a_{42}(a_{33}+a_{55}+a_{66}+a_{77})+a_{54}a_{42}a_{25}\\&\qquad +a_{33}a_{44}(a_{55}+a_{66}+a_{77})+a_{33}a_{55}(a_{66}+a_{77})\\&\qquad +a_{33}a_{66}a_{77}+a_{44}a_{55}(a_{66}+a_{77})\\&\qquad +a_{66}a_{77}(a_{44}+a_{55})] \end{aligned}$$

$$\begin{aligned} R_{4}&= -[-a_{11}a_{22}a_{33}(a_{44}+a_{55}+a_{66}+a_{77})\\&\qquad -a_{11}a_{22} a_{44}(a_{55}+a_{66}+a_{77})\\&\qquad -a_{11}a_{22}a_{55}(a_{66}+a_{77}) -a_{11}a_{22}a_{66}a_{77}\\&\qquad +a_{11}a_{23}a_{32}(a_{44}+a_{55}+a_{66}+a_{77})\\&\qquad -a_{11}a_{24}a_{32}a_{43}+a_{11}a_{24}a_{42}\\&\qquad (a_{33}+a_{55}+a_{66} +a_{77})-a_{11}a_{25}a_{42}a_{54}\\&\qquad -a_{11}a_{33}a_{44}(a_{55}+a_{66}+a_{77})\\&\qquad -a_{11}a_{33}a_{55}(a_{66}+a_{77})-a_{11}a_{33}a_{66}a_{77}\\&\qquad -a_{11}a_{44}a_{55}(a_{66}+a_{77})\\&\qquad -a_{11}a_{66}a_{77}(a_{44}+a_{55})\\&\qquad -a_{13}a_{32}a_{21}(a_{44}+a_{55}+a_{66}+a_{77})\\&\qquad +a_{14}a_{43}a_{32}a_{21}\\&\qquad -a_{14}a_{42}a_{21}(a_{33}+a_{55}+a_{66}+a_{77})\\&\qquad +a_{15}a_{54}a_{42}a_{21}\\&\qquad +a_{16}a_{22}a_{61}(a_{33}+a_{44}+a_{55}+a_{77})\\&\qquad -a_{16}a_{61}(a_{23}a_{32}+a_{24}a_{42})\\&\qquad +a_{16}a_{33}a_{61}(a_{44}+a_{55} +a_{77})\\&\qquad +a_{16}a_{44}a_{61}(a_{55}+a_{77})\\&\qquad +a_{16}a_{55}a_{61}a_{77}\\&\qquad +a_{74}a_{42}a_{21}a_{17}\\&\qquad -a_{22}a_{33}a_{44}(a_{55}+a_{66}+a_{77})\\&\qquad -a_{22}a_{33}a_{55}(a_{66}+a_{77})-a_{22}a_{33}a_{66}a_{77}\\&\qquad -a_{22}a_{44}a_{55}(a_{66}+a_{77})\\&\qquad -a_{22}a_{66}a_{77}(a_{44}+a_{55})\\&\qquad +a_{23}a_{32}a_{44}(a_{55}+a_{66}+a_{77})\\&\qquad +a_{23}a_{32}a_{55}(a_{66}+a_{77})+a_{22}a_{32}a_{66}a_{77}\\&\qquad -a_{24} a_{32}a_{43}(a_{55}+a_{66}+a_{77})\\&\qquad +a_{24}a_{33}a_{42}(a_{55}+a_{66}+a_{77})\\&\qquad +a_{24}a_{42}a_{55}(a_{66}+a_{77})+a_{24}a_{42}a_{66}a_{77}\\&\qquad +a_{25} a_{32}a_{43}a_{54}-a_{24}a_{42}a_{54}(a_{33}+a_{66}+a_{77})\\&\qquad -a_{33}a_{44}a_{55}(a_{66}+a_{77})-a_{33}a_{66}\\&\qquad a_{77}(a_{44}+a_{55})-a_{44}a_{55}a_{66}a_{77}] \end{aligned}$$

$$\begin{aligned} R_{5}&=-[a_{11}a_{22}a_{33}a_{44}(a_{55}+a_{66}+a_{77})\\&\qquad +a_{11}a_{22}a_{33} a_{55}(a_{66}+a_{77})\\&\qquad +a_{11}a_{32}a_{33}a_{66}a_{77}\\&\qquad +a_{11}a_{22}a_{44}a_{55}(a_{66}+a_{77})\\&\qquad +a_{11}a_{22}a_{66}a_{77} (a_{44}+a_{55})\\&\qquad -a_{11}a_{23}a_{32}a_{44}(a_{55}+a_{66}+a_{77})\\&\qquad -a_{11}a_{23}a_{32}a_{55}(a_{66}+a_{77})-a_{11}a_{23}a_{32}a_{66}a_{77}\\&\qquad +a_{11}a_{24}a_{32}a_{43}(a_{55}+a_{66}+a_{77})\\&\qquad -a_{11}a_{24}a_{33}a_{42}(a_{55} +a_{66}+a_{77})\\&\qquad -a_{11}a_{24}a_{42}a_{55}(a_{66}+a_{77})\\&\qquad -a_{11}a_{24}a_{42}a_{66} a_{77}-a_{11}a_{25}a_{32}\\&\qquad a_{43}a_{54}+a_{11}a_{25}a_{33}a_{42}a_{54}\\&\qquad +a_{11}a_{25}a_{42}a_{54}(a_{66}+a_{77})\\&\qquad +a_{11}a_{33}a_{44}a_{55}(a_{66}+a_{77})\\&\qquad +a_{11}a_{33}a_{66}a_{77}(a_{44}+a_{55})\\&\qquad +a_{11}a_{44}a_{55}a_{66}a_{77}\\&\qquad +a_{13}a_{21}a_{32}a_{44}(a_{55}+a_{66}+a_{77}) \\&\qquad +a_{13}a_{21}a_{32}a_{55}(a_{66} +a_{77})\\&\qquad +a_{13}a_{21}a_{32}a_{66}a_{77}\\&\qquad -a_{14}a_{21}a_{32}a_{43}(a_{55} +a_{66}+a_{77}) \\&\qquad +a_{14}a_{21}a_{33}a_{42}(a_{55}+a_{66}+a_{77}) \\&\qquad +a_{44}a_{21}a_{42}a_{55}(a_{66}+a_{77})\\&\qquad +a_{14}a_{21}a_{42}a_{66}a_{77} \\&\qquad + a_{15}a_{21}a_{32}a_{43}a_{54}\\&\qquad -a_{15}a_{21}a_{42}a_{54}(a_{33}+a_{66}+a_{77})\\&\qquad +a_{65}a_{54}a_{42}a_{21}a_{16}+a_{16}a_{67}a_{74}a_{42}a_{21} \\&\qquad - a_{16}a_{22}a_{33}a_{61}(a_{44}+a_{55}+a_{77})\\&\qquad -a_{16}a_{22}a_{44}a_{61} (a_{55}+a_{77})\\&\qquad -a_{16}a_{22}a_{55}a_{61}a_{77}\\&\qquad +a_{16}a_{23}a_{32}a_{61}(a_{44} +a_{55}+a_{77})\\&\qquad -a_{16}a_{24}a_{32}a_{43}a_{61}\\&\qquad +a_{16}a_{24}a_{42}a_{61} (a_{33}+a_{55}+a_{77})\\&\qquad -a_{16}a_{25}a_{42}a_{54}a_{61}\\&\qquad -a_{16}a_{33}a_{44}a_{61}(a_{55}+a_{77})\\&\qquad -a_{16}a_{33}a_{55}a_{61}a_{77}\\&\qquad -a_{16}a_{44}a_{55}a_{61}a_{77}+a_{17}a_{21}a_{32}a_{43}a_{74} \\&\qquad - a_{74}a_{42}a_{21}a_{17}(a_{33}a_{55}+a_{66})\\&\qquad +a_{22}a_{33}a_{44}a_{55} (a_{66}+a_{77})\\&\qquad +a_{22}a_{33}(a_{44}+a_{55})a_{66}a_{77}\\&\qquad +a_{22}a_{44}a_{55}a_{66}a_{77}\\&\qquad -a_{23}a_{32}a_{44}a_{55}(a_{66}+a_{77})\\&\qquad -a_{23}a_{32}(a_{44}+a_{55})a_{66}a_{77}\\&\qquad +a_{24}a_{32}a_{43}a_{55}(a_{66}+a_{77})\\&\qquad +a_{24}a_{32}a_{43}a_{66}a_{77}\\&\qquad -a_{24}a_{33}a_{42}a_{55}(a_{66}+a_{77})\\&\qquad -a_{24}a_{42}a_{66}a_{77} (a_{33}+a_{55})\\&\qquad -a_{25}a_{32}a_{43}a_{54}(a_{66}+a_{77})\\&\qquad +a_{25}a_{33}a_{42}a_{54}(a_{66}+a_{77})\\&\qquad +a_{25}a_{42}a_{54}a_{66}a_{77} +a_{33}a_{44}a_{55}a_{66}a_{77}] \end{aligned}$$

$$\begin{aligned} R_{6}&=-[-a_{11}a_{22}a_{33}a_{44}a_{55}(a_{66}+a_{77}) \\&\qquad -a_{11}a_{22}a_{33}(a_{44}+a_{55})a_{66}a_{77}\\&\qquad -a_{11}a_{22}a_{44}a_{55}a_{66}a_{77} \\&\quad + a_{11}a_{23}a_{32}a_{44}a_{55}(a_{66}+a_{77})\\&\qquad +a_{11}a_{23}a_{32}(a_{44}+a_{55})a_{66}\\&\qquad -a_{11}a_{24}a_{32}a_{43}a_{55}(a_{66}+a_{77}) \\&\qquad - a_{11}a_{24}a_{32}a_{43}a_{66}a_{77}\\&\qquad +a_{11}a_{24}a_{33}a_{42}a_{55}(a_{66}+a_{77})\\&\qquad +a_{11}a_{24}a_{42}(a_{55}+a_{33})a_{66}a_{77}\\&\qquad +a_{11}a_{25}a_{32}a_{43}a_{54} (a_{66}+a_{77})\\&\qquad -a_{11}a_{25}a_{33}a_{42}a_{54}(a_{66}+a_{77})\\&\qquad -a_{11}a_{25}a_{42}a_{54}a_{66}a_{77}\\&\qquad -a_{11}a_{33}a_{44}a_{55}a_{66}a_{77} \\&\qquad - a_{13}a_{21}a_{32}a_{44}a_{55}(a_{66}+a_{77})\\&\qquad -a_{13}a_{21}a_{32}(a_{44}+a_{55})a_{66}a_{77}\\&\qquad +a_{14}a_{21}a_{32}a_{43}a_{55}(a_{66}+a_{77})\\&\qquad +a_{14}a_{21}a_{32}a_{43}a_{66}a_{77}\\&\qquad -a_{14}a_{21}a_{33}a_{42}a_{55}(a_{66}+a_{77})\\&\qquad -a_{14}a_{21}(a_{33}+a_{55})a_{42}a_{66}a_{77}\\&\qquad -a_{15}a_{21}a_{32}a_{43}a_{54} (a_{66}+a_{77})\\&\qquad +a_{15}a_{21}a_{33}a_{42}a_{54}(a_{66}+a_{77})\\&\qquad +a_{15}a_{21}a_{42}a_{54}a_{66}a_{77}\\&\qquad +a_{16}a_{21}a_{32}a_{43}(a_{54}a_{65}+a_{67}a_{74})\\&\qquad -a_{16}a_{21}a_{33}a_{42}(a_{54}a_{65}+a_{67}a_{74})\\&\qquad -a_{16}a_{21}a_{42}a_{54}a_{65}a_{77}\\&\qquad -a_{16}a_{21}a_{42}a_{55}a_{67}a_{74}\\&\qquad +a_{16}a_{22}a_{33}a_{44}a_{61} (a_{55}+a_{77}) \\&\qquad +a_{16}a_{22}a_{77}a_{55}a_{61}(a_{33}+a_{44})\\&\qquad -a_{16}a_{23} a_{32}a_{44}a_{61}(a_{55}+a_{77}) \\&\qquad -a_{16}a_{23}a_{32}a_{55}a_{61}a_{77} \\&\qquad + a_{16}a_{24}a_{32}a_{43}a_{61}(a_{55}+a_{77}) \\&\qquad -a_{16}a_{24}a_{33} a_{42}a_{61}(a_{55}+a_{77}) \\&\qquad +a_{16}a_{24}a_{42}a_{55}a_{61}a_{77} \\&\qquad - a_{16}a_{25}a_{32}a_{43}a_{54}a_{61}\\&\qquad +a_{16}a_{25}a_{42}a_{54}a_{61}(a_{33}+a_{77})\\&\qquad +a_{16}a_{33}a_{44}a_{55}a_{61}a_{77}\\&\qquad -a_{17}a_{21}a_{32}a_{43}a_{74} (a_{55}+a_{66})\\&\qquad +a_{17}a_{21}a_{33}a_{42}a_{74}(a_{55}+a_{66})\\&\qquad +a_{17}a_{21}a_{42}a_{55} a_{66}a_{74}\\&\qquad -a_{22}a_{33}a_{44}a_{55}a_{66}a_{77} \\&\qquad + a_{23}a_{32}a_{44}a_{55}a_{66}a_{77}\\&\qquad -a_{24}a_{43}a_{55}a_{66}a_{77}\\&\qquad +a_{24} a_{33}a_{42}a_{55}a_{66}a_{77}\\&\qquad +a_{25}a_{32}a_{43}a_{54}a_{66}a_{77}\\&\qquad -a_{25}a_{33}a_{42}a_{54}a_{66}a_{77}] \\ \end{aligned}$$

$$\begin{aligned} R_{7}&=a_{11}a_{23}a_{32}a_{44}a_{55}a_{66}a_{77}\\&\quad -a_{11}a_{24}a_{32}a_{43}a_{55}a_{66}a_{77}\\&\quad +a_{11}a_{24}a_{33}a_{42}a_{55}a_{66}a_{77}\\&\quad +a_{11}a_{25}a_{32}a_{43}a_{54}a_{66}a_{77}\\&\quad -a_{11}a_{25}a_{33}a_{42}a_{54}a_{66}a_{77}\\&\quad -a_{13}a_{21}a_{32}a_{44}a_{55}a_{66}a_{77}\\&\quad +a_{14}a_{21}a_{32}a_{43}a_{55}a_{66}a_{77}\\&\quad -a_{14}a_{21}a_{33}a_{42}a_{55}a_{66}a_{77}\\&\quad -a_{15}a_{21}a_{32}a_{43}a_{54}a_{66}a_{77}\\&\quad +a_{15}a_{21}a_{33}a_{42}a_{54}a_{66}a_{77}\\&\quad +a_{16}a_{21}a_{32}a_{43}a_{54}a_{65}a_{77}\\&\quad +a_{16}a_{21}a_{32}a_{43}a_{55}a_{67}a_{74}\\&\quad -a_{16}a_{21}a_{33}a_{42}a_{54}a_{65}a_{77}\\&\quad -a_{16}a_{21}a_{33}a_{42}a_{55}a_{67}a_{74}\\&\quad +a_{16}a_{22}a_{33}a_{44}a_{55}a_{61}a_{77}\\&\quad -a_{16}a_{23}a_{32}a_{44}a_{55}a_{61}a_{77}\\&\quad +a_{16}a_{24}(a_{32}a_{43}-a_{33}a_{42})a_{55}a_{61}a_{77}\\&\quad +a_{16}a_{25}(a_{33}a_{42}-a_{32}a_{43})a_{54}a_{61}a_{77}\\&\quad +a_{17}a_{21}(a_{33}a_{42}-a_{32}a_{43})a_{55}a_{66}a_{74}\\&\quad -a_{11}a_{22}a_{33}a_{44}a_{55}a_{66}a_{77} \end{aligned}$$

Let us consider

$$\begin{aligned} \begin{aligned} \Delta _{1}&=R_{1},\ \Delta _{2}=\left| \begin{array}{cc} R_{1} &{} 1 \\ R_{3} &{} R_{2} \end{array} \right| ,\ \Delta _{3}=\left| \begin{array}{ccc} R_{1} &{} 1 &{} 0\\ R_{3} &{} R_{2} &{} R_{1} \\ R_{5} &{} R_{4} &{} R_{3} \end{array} \right| , \\ \Delta _{4}&=\left| \begin{array}{cccc} R_{1} &{} 1 &{} 0 &{} 0\\ R_{3} &{} R_{2} &{} R_{1} &{} 1 \\ R_{5} &{} R_{4} &{} R_{3} &{} R_{2} \\ R_{7} &{} R_{6} &{} R_{5} &{} R_{4} \end{array} \right| , \\ \Delta _{5}&=\left| \begin{array}{ccccc} R_{1} &{} 1 &{} 0 &{} 0 &{} 0\\ R_{3} &{} R_{2} &{} R_{1} &{} 1 &{} 0 \\ R_{5} &{} R_{4} &{} R_{3} &{} R_{2} &{} R_{1} \\ R_{7} &{} R_{6} &{} R_{5} &{} R_{4} &{} R_{3} \\ 0 &{} 0 &{} R_{7} &{} R_{6} &{} R_{5} \end{array} \right| , \\ \Delta _{6}&=\left| \begin{array}{cccccc} R_{1} &{} 1 &{} 0 &{} 0 &{} 0 &{} 0 \\ R_{3} &{} R_{2} &{} R_{1} &{} 1 &{} 0 &{} 0 \\ R_{5} &{} R_{4} &{} R_{3} &{} R_{2} &{} R_{1} &{} 1 \\ R_{7} &{} R_{6} &{} R_{5} &{} R_{4} &{} R_{3} &{} R_{2} \\ 0 &{} 0 &{} R_{7} &{} R_{6} &{} R_{5} &{} R_{4} \\ 0 &{} 0 &{} 0 &{} 0 &{} R_{7} &{} R_{6} \end{array} \right| , \\ \Delta _{7}&=Det({\overline{J}}|_{E^{*}}). \end{aligned} \end{aligned}$$

By Routh–Hurwitz criterion, $E^{*}$ is locally asymptomatically stable (LAS) if and only if $\Delta _{i}>0$ for $i=1,2,3,4,5,6,7$, i.e. equivalently

• (i) $R_{i}>0$ for $i=1,7;$

• (ii) $\Delta _i >0$ for $i=2,3,4,5,6$.

$\square$

1.2 Existence of optimal control functions

Now we derive the conditions for existence of optimal control interventions which also minimise the cost function J in a finite time period.

Proof of Theorem 9.1

The optimal control variables, when exist, satisfy the following conditions:

• (i) Solutions of system (6) with control variables $u_{1},\ u_{2}$ and $u_{3}$ in $\Pi \ne \phi$.

• (ii) The mentioned set $\Pi$ is closed, convex and the state system is represented with linear function of control variables where coefficients depend on time and also on state variables.

• (iii) Integrand of (5): L is convex on $\Pi$ and $L(S,E,A,I,H,R,u_{1},u_{2},u_{3})\ge f(u_{1},u_{2},u_{3})$ where $f(u_{1},u_{2},u_{3})$ is continuous and $||(u_{1},u_{2},u_{3})||^{-1}f(u_{1},u_{2},u_{3})\rightarrow \infty$ when $||(u_{1},u_{2},u_{3})||\rightarrow \infty$; ||.|| represents the $L^{3}(0,T_{f})$ norm.

From (6), the total population $N=S+E+A+I+H+R$.

$$\begin{aligned} \begin{aligned} \text {So,}\ \frac{\hbox {d}N}{\hbox {d}t}&= \Lambda -dN-\delta _{1}A-\delta _{2}I-\delta _{3}H \le \Lambda -dN \\ \Rightarrow 0&<N(t)\le N(0)e^{-dt} +\frac{\Lambda }{d}\left( 1-e^{-dt}\right) , \end{aligned} \end{aligned}$$

where $N(0)=S(0)+E(0)+A(0)+I(0)+H(0)+R(0)$.

As $\displaystyle t\rightarrow \infty ,\ 0<N(t)\le \frac{\Lambda }{d}$.

$$\begin{aligned} \begin{aligned} \text {Also},\ \frac{\hbox {d}Z}{\hbox {d}t}&= \frac{pI}{1+qI}-a_{0}Z\le pI-a_{0}Z \\ \Rightarrow \,&\frac{\hbox {d}Z}{\hbox {d}t}+a_{0}Z\le \frac{p\Lambda }{d}\ (\text {for large time} \ t)\\ \Rightarrow&0<Z(t)\le Z(0)e^{-a_{0}t} +\frac{p\Lambda }{da_{0}}\left( 1-e^{-a_{0}t}\right) \end{aligned} \end{aligned}$$

As $\displaystyle t\rightarrow \infty ,\ 0<Z(t)\le \frac{p\Lambda }{a_{0}d}$.

For each of the control variable in $\Pi$, solution of (6) is bounded and right-hand side functions are locally Lipschitzian too. $Picard-Lindel\ddot{o}f$ theorem shows that condition (i) is satisfied [10].

The control set $\Pi$ is closed and convex by definition. Again all the equations of system (6) are written as linear equations in $u_{1},\ u_{2}$ and $u_{3}$ where state variables depend on coefficients and hence condition (ii) is satisfied also. Moreover, the quadratic nature of all control variables guarantee the convex property of integrand $L(S,E,A,I,H,R,Z,u_{1},u_{2},u_{3})$.

$$\begin{aligned} \begin{aligned}&\text {Also,}\ L(S,E,A,I,H,R,Z,u_{1},u_{2},u_{3}) \\&\quad = w_{1}I+w_{2}H+w_{3}u^{2}_{1}+w_{4}u^{2}_{2}+w_{5}u^{2}_{3} \\&\quad \ge w_{3}u^{2}_{1}+w_{4}u^{2}_{2}+w_{5}u^{2}_{3} \end{aligned} \end{aligned}$$

Let, ${\overline{u}} = \min (w_{3},w_{4},w_{5})>0$ and $f(u_{1},u_{2},u_{3})= {\overline{u}}(u^{2}_{1}+ u^{2}_{2}+u^{2}_{3})$.

Then $L(S,E,A,I,H,R,Z,u_{1},u_{2},u_{3}){\ge } f(u_{1},u_{2},u_{3}).$

Here, f is continuous and $||(u_{1},u_{2},u_{3})||^{-1}f(u_{1},u_{2},u_{3})\rightarrow \infty$ whenever $||(u_{1},u_{2},u_{3})||\rightarrow \infty .$ Hence, condition (iii) is also satisfied. So, it is concluded that there exist control variables $u^{*}_{1},\ u^{*}_{2}$ and $u^{*}_{3}$ with the condition $J[u^{*}_{1},u^{*}_{2},u^{*}_{3}]=\min [J[u_{1},u_{2},u_{3}]]$ [16, 17]. $\square$

1.3 Characterisation of optimal control functions

By Pontryagin’s Maximum Principle, we have derived here the necessary conditions for optimal control functions for system (5)–(6) [16, 31]. Let us define the Hamiltonian function as:

$$\begin{aligned} \begin{aligned}&{\overline{H}}\left( S,E,A,I,H,R,Z,u_{1},u_{2},u_{3},\lambda \right) \\&\quad =L(S,E,A,I,H,R,Z,u_{1},u_{2},u_{3}) \\&\qquad +\lambda _{1}\frac{\hbox {d}S}{\hbox {d}t}+\lambda _{2}\frac{\hbox {d}E}{\hbox {d}t}+\lambda _{3}\frac{\hbox {d}A}{\hbox {d}t}\\&\qquad +\lambda _{4}\frac{\hbox {d}I}{\hbox {d}t}+\lambda _{5}\frac{\hbox {d}H}{\hbox {d}t}+\lambda _{6}\frac{\hbox {d}R}{\hbox {d}t} +\lambda _{7}\frac{\hbox {d}Z}{\hbox {d}t} \\ \text {So},\,&{\overline{H}} =w_{1}I+w_{2}H+w_{3}u^{2}_{1}+w_{4}u^{2}_{2} +w_{5}u^{2}_{3} \\&\qquad +\lambda _{1}[\Lambda -(\beta _{1}A+\beta _{2}I+\beta _{3}H)S\\&\quad \qquad -dS+\xi R -u_{1}(t)kSZ] \\&\qquad +\lambda _{2}[(\beta _{1}A+\beta _{2}I+\beta _{3}H)S-(\kappa +d)E]\\&\qquad +\lambda _{3}[\kappa \phi E-\eta A-(d+\delta _{1})A]\\&\qquad +\lambda _{4}\left[ \kappa (1-\phi )E+\eta A-(\omega +d+\delta _{2})I \right. \\&\quad \qquad \left. -\frac{\epsilon _{1} u_{2}(t)I}{1+\zeta _{1}I}\right] \\&\qquad +\lambda _{5}\left[ \omega I-\gamma H-(d+\delta _{3})H +\frac{\epsilon _{1} u_{2}(t)I}{1+\zeta _{1}I}\right. \\&\quad \qquad \left. - \frac{\epsilon _{2}u_{3}(t)H}{1+\zeta _{2}H}\right] \\&\qquad + \lambda _{6}\left[ \gamma H-dR-\xi R+u_{1}(t)kSZ \right. \\&\quad \qquad \left. + \frac{\epsilon _{2}u_{3}(t)H}{1+\zeta _{2}H}\right] \\&\qquad +\lambda _{7}\left[ \frac{pI}{1+qI}-a_{0}Z\right] \end{aligned} \end{aligned}$$

(9)

Here, $\lambda =\left( \lambda _{1},\lambda _{2},\lambda _{3}, \lambda _{4},\lambda _{5},\lambda _{6},\lambda _{7}\right)$ are the adjoint variables. We get minimised Hamiltonian by Pontryagin’s Maximum Principle to minimise the cost functional. Pontryagin’s Maximum Principle mainly adjoin the cost functional with the state equations by introducing adjoint variables.

Proof of Theorem 9.2

Let $u^{*}_{1},\ u^{*}_{2}$ and $u^{*}_{3}$ be optimal control variables and $S^{*}, E^{*}, A^{*}, I^{*}, H^{*}, R^{*}, Z^{*}$ are corresponding optimal state variables of the control system (6) which minimise the cost functional (5). So, by Pontryagin’s Maximum Principle, there exist adjoint variables $\lambda _{1},\ \lambda _{2},\ \lambda _{3},\ \lambda _{4},\ \lambda _{5},\ \lambda _{6},\ \lambda _{7}$ which satisfy the following canonical equations:

$$\begin{aligned} \begin{aligned} \frac{\hbox {d}\lambda _{1}}{\hbox {d}t}&=-\frac{\partial {\overline{H}}}{\partial S}, \qquad \frac{\hbox {d}\lambda _{2}}{\hbox {d}t}=-\frac{\partial {\overline{H}}}{\partial E}, \\ \frac{\hbox {d}\lambda _{3}}{\hbox {d}t}&=-\frac{\partial {\overline{H}}}{\partial A}, \qquad \frac{\hbox {d}\lambda _{4}}{\hbox {d}t}=-\frac{\partial {\overline{H}}}{\partial I}, \\ \frac{\hbox {d}\lambda _{5}}{\hbox {d}t}&=-\frac{\partial {\overline{H}}}{\partial H}, \\ \frac{\hbox {d}\lambda _{6}}{\hbox {d}t}&=-\frac{\partial {\overline{H}}}{\partial R}, \qquad \frac{\hbox {d}\lambda _{7}}{\hbox {d}t}=-\frac{\partial {\overline{H}}}{\partial Z}. \end{aligned} \end{aligned}$$

So, we have

$$\begin{aligned} \begin{aligned} \frac{\hbox {d}\lambda _{1}}{\hbox {d}t}&=\lambda _{1}[\beta _{1}A+\beta _{2}I+\beta _{3}H+d+u_{1}kZ] \\&\quad -\lambda _{2}(\beta _{1}A+\beta _{2}I+\beta _{3}H)-\lambda _{6}(u_{1}kZ)\\ \frac{\hbox {d}\lambda _{2}}{\hbox {d}t}&=\lambda _{2}(d+\kappa )-\lambda _{3}(\kappa \phi )\\&\quad -\lambda _{4}\{\kappa (1-\phi )\} \\ \frac{\hbox {d}\lambda _{3}}{\hbox {d}t}&=\lambda _{1}(\beta _{1}S)-\lambda _{2}(\beta _{1}S)\\&\quad +\lambda _{3}(\eta +d+\delta _{1})-\lambda _{4}(\eta ) \\ \frac{\hbox {d}\lambda _{4}}{\hbox {d}t}&=-w_{1}+\lambda _{1}(\beta _{2}S)-\lambda _{2}(\beta _{2}S)\\&\quad +\lambda _{4}\left\{ \omega +d+\delta _{2} +\frac{\epsilon _{1}u_{2}}{(1+\zeta _{1}I)^{2}}\right\} \\&\quad -\lambda _{5}\left\{ \omega +\frac{\epsilon _{1}u_{2}}{(1+\zeta _{1}I)^{2}}\right\} \\&-\lambda _{7}\left\{ \frac{p}{(1+qI)^{2}}\right\} \\ \frac{\hbox {d}\lambda _{5}}{\hbox {d}t}&=-w_{2}+\lambda _{1}(\beta _{3}S)-\lambda _{2}(\beta _{3}S)\\&\quad +\lambda _{5}\left\{ \gamma +d+\delta _{3} +\frac{\epsilon _{2}u_{3}}{(1+\zeta _{2}H)^{2}}\right\} \\&\quad -\lambda _{6}\left\{ \gamma +\frac{\epsilon _{2}u_{3}}{(1+\zeta _{2}H)^{2}}\right\} \\ \frac{\hbox {d}\lambda _{6}}{\hbox {d}t}&=-\lambda _{1}(\xi )+\lambda _{6}(d+\xi ) \\ \frac{\hbox {d}\lambda _{7}}{\hbox {d}t}&=\lambda _{1}(u_{1}kS)-\lambda _{6}(u_{1}kS)+\lambda _{7}(a_{0}) \end{aligned} \end{aligned}$$

(10)

with the transversality conditions $\lambda _{i}(T_{f})=0$, for $i=1,2,3,4,5,6,7.$

$$\begin{aligned}&\text {From optimality conditions}: \\&\frac{\partial H}{\partial u_{1}}\bigg |_{u_{1}=u^{*}_{1}}=0,\ \frac{\partial H}{\partial u_{2}}\bigg |_{u_{2}=u^{*}_{2}}=0 \ \text {and} \ \frac{\partial H}{\partial u_{3}}\bigg |_{u_{3}=u^{*}_{3}}=0. \end{aligned}$$

So, $\displaystyle u^{*}_{1}= \frac{kS^{*}Z^{*}}{2w_{3}}\left( \lambda _{1}-\lambda _{6}\right)$, $u^{*}_{2}=\frac{\epsilon _{1}I^{*}}{2w_{4}\left( 1+\zeta _{1}I^{*}\right) } \left( \lambda _{4}-\lambda _{5}\right)$ and $\displaystyle u^{*}_{3}=\frac{\epsilon _{2}H^{*}}{2w_{5} \left( 1+\zeta _{2}H^{*}\right) } \left( \lambda _{5}-\lambda _{6}\right) .$ Now from these findings along with the characteristics of control set $\Pi ,$ we have

$$\begin{aligned} \begin{aligned} u^{*}_{1}&= {\left\{ \begin{array}{ll} 0, &{} \text {if}\ \frac{kS^{*}Z^{*}}{2w_{3}}\left( \lambda _{1}-\lambda _{6}\right)< 0 \\ \\ \frac{kS^{*}Z^{*}}{2w_{3}}\left( \lambda _{1}-\lambda _{6}\right) , &{} \text {if}\ 0\le \frac{kS^{*}Z^{*}}{2w_{3}}\left( \lambda _{1}-\lambda _{6}\right) \le 1 \\ \\ 1, &{} \text {if}\ \frac{kS^{*}Z^{*}}{2w_{3}}\left( \lambda _{1}-\lambda _{6}\right)> 1 \end{array}\right. } \\ u^{*}_{2}&= {\left\{ \begin{array}{ll} 0, &{} \text {if}\ \frac{\epsilon _{1}I^{*}}{2w_{4}\left( 1+\zeta _{1}I^{*}\right) }\left( \lambda _{4}-\lambda _{5}\right)< 0 \\ \\ \frac{\epsilon _{1}I^{*}}{2w_{4}\left( 1+\zeta _{1}I^{*}\right) }\left( \lambda _{4}-\lambda _{5}\right) , &{} \text {if}\ 0\le \frac{\epsilon _{1}I^{*}}{2w_{4}\left( 1+\zeta _{1}I^{*}\right) }\left( \lambda _{4}-\lambda _{5}\right) \le 1 \\ \\ 1, &{} \text {if}\ \frac{\epsilon _{1}I^{*}}{2w_{4}\left( 1+\zeta _{1}I^{*}\right) }\left( \lambda _{4}-\lambda _{5}\right)> 1 \end{array}\right. }\\ u^{*}_{3}&= {\left\{ \begin{array}{ll} 0, &{} \text {if}\ \frac{\epsilon _{2}H^{*}}{2w_{5}\left( 1+\zeta _{2}H^{*}\right) }\left( \lambda _{5}-\lambda _{6}\right) < 0 \\ \\ \frac{\epsilon _{2}H^{*}}{2w_{5}\left( 1+\zeta _{2}H^{*}\right) }\left( \lambda _{5}-\lambda _{6}\right) , &{} \text {if}\ 0\le \frac{\epsilon _{2}H^{*}}{2w_{5}\left( 1+\zeta _{2}H^{*}\right) }\left( \lambda _{5}-\lambda _{6}\right) \le 1 \\ \\ 1, &{} \text {if}\ \frac{\epsilon _{2}H^{*}}{2w_{5}\left( 1+\zeta _{2}H^{*}\right) }\left( \lambda _{5}-\lambda _{6}\right) > 1 \end{array}\right. } \end{aligned} \end{aligned}$$

which is equivalent as (8). $\square$

1.4 Optimal system

We state the optimal system with optimal control variables $u^{*}_{1},\ u^{*}_{2}$ and $u^{*}_{3}$ below. The optimal system with minimised Hamiltonian ${\overline{H}}^{*}$ at $(S^{*},E^{*}, A^{*}, I^{*}, H^{*}, R^{*}, Z^{*}$, $\lambda _{1}, \lambda _{2}, \lambda _{3}, \lambda _{4}, \lambda _{5}$, $\lambda _{6}, \lambda _{7})$ is as follows:

$$\begin{aligned} \frac{\hbox {d}S^{*}}{\hbox {d}t}&= \Lambda -(\beta _{1}A^{*}+\beta _{2}I^{*}+\beta _{3}H^{*})S^{*}\nonumber \\&\quad -dS^{*}+\xi R^{*}-u_{1}^{*}kS^{*}Z^{*}, \nonumber \\ \frac{\hbox {d}E^{*}}{\hbox {d}t}&= (\beta _{1}A^{*}+\beta _{2}I^{*}+\beta _{3}H^{*})S^{*}-(\kappa +d)E^{*}, \nonumber \\ \frac{\hbox {d}A^{*}}{\hbox {d}t}&= \kappa \phi E^{*}-\eta A^{*}-(d+\delta _{1})A^{*}, \nonumber \\ \frac{\hbox {d}I}{\hbox {d}t}&= \kappa (1-\phi )E^{*}+\eta A^{*}-(\omega +d+\delta _{2})I^{*}\nonumber \\&\quad -\frac{\epsilon _{1} u_{2}^{*}I^{*}}{1+\zeta _{1}I^{*}}, \nonumber \\ \frac{\hbox {d}H^{*}}{\hbox {d}t}&=\omega I^{*}-\gamma H^{*}-(d+\delta _{3})H^{*}\nonumber \\&\quad +\frac{\epsilon _{1} u_{2}^{*}I^{*}}{1+\zeta _{1}I^{*}} -\frac{\epsilon _{2}u_{3}^{*}H^{*}}{1+\zeta _{2}H^{*}}, \nonumber \\ \frac{\hbox {d}R^{*}}{\hbox {d}t}&= \gamma H^{*}-dR^{*}-\xi R^{*}+u_{1}^{*}kS^{*}Z^{*}\nonumber \\&\quad +\frac{\epsilon _{2}u_{3}^{*}H^{*}}{1+\zeta _{2}H^{*}}, \nonumber \\ \frac{\hbox {d}Z^{*}}{\hbox {d}t}&=\frac{pI^{*}}{1+qI^{*}}-a_{0}Z^{*}, \end{aligned}$$

(11)

with initial conditions: $S^{*}(0)>0,\ E^{*}(0)\ge 0,\ A^{*}(0)\ge 0,\ I^{*}(0)\ge 0, H^{*}\ge 0,\ R^{*}(0)\ge 0$ and $Z^{*}\ge 0$. The corresponding adjoint system is given as:

$$\begin{aligned} \frac{\hbox {d}\lambda _{1}}{\hbox {d}t}&=\lambda _{1}[\beta _{1}A^{*}+\beta _{2}I^{*} +\beta _{3}H^{*}+d+u_{1}^{*}kZ^{*}]\nonumber \\&\quad -\lambda _{2}(\beta _{1}A^{*}+\beta _{2}I^{*}+\beta _{3}H^{*})\nonumber \\&\quad -\lambda _{6}(u_{1}^{*}kZ^{*})\nonumber \\ \frac{\hbox {d}\lambda _{2}}{\hbox {d}t}&=\lambda _{2}(d+\kappa )-\lambda _{3}(\kappa \phi )\nonumber \\&\quad -\lambda _{4}\{\kappa (1-\phi )\} \nonumber \\ \frac{\hbox {d}\lambda _{3}}{\hbox {d}t}&=\lambda _{1}(\beta _{1}S^{*})-\lambda _{2}(\beta _{1}S^{*})\nonumber \\&\quad +\lambda _{3}(\eta +d+\delta _{1})-\lambda _{4}(\eta ) \nonumber \\ \frac{\hbox {d}\lambda _{4}}{\hbox {d}t}&=-w_{1}+\lambda _{1}(\beta _{2}S^{*})-\lambda _{2}(\beta _{2}S^{*})\nonumber \\&\quad +\lambda _{4}\left\{ \omega +d+\delta _{2} +\frac{\epsilon _{1}u_{2}^{*}}{(1+\zeta _{1}I^{*})^{2}}\right\} \nonumber \\&\quad -\lambda _{5}\left\{ \omega +\frac{\epsilon _{1}u_{2}^{*}}{(1+\zeta _{1}I^{*})^{2}}\right\} \nonumber \\&-\lambda _{7}\left\{ \frac{p}{(1+qI^{*})^{2}}\right\} \nonumber \\ \frac{\hbox {d}\lambda _{5}}{\hbox {d}t}&=-w_{2}+\lambda _{1}(\beta _{3}S^{*})-\lambda _{2}(\beta _{3}S^{*})\nonumber \\&\quad +\lambda _{5}\left\{ \gamma +d+\delta _{3} +\frac{\epsilon _{2}u_{3}^{*}}{(1+\zeta _{2}H^{*})^{2}}\right\} \nonumber \\&\quad -\lambda _{6}\left\{ \gamma +\frac{\epsilon _{2}u_{3}^{*}}{(1+\zeta _{2}H^{*})^{2}}\right\} \nonumber \\ \frac{\hbox {d}\lambda _{6}}{\hbox {d}t}&=-\lambda _{1}(\xi )+\lambda _{6}(d+\xi ) \nonumber \\ \frac{\hbox {d}\lambda _{7}}{\hbox {d}t}&=\lambda _{1}(u_{1}^{*}kS^{*})\nonumber \\&\quad -\lambda _{6}(u_{1}^{*}kS^{*})+\lambda _{7}(a_{0}), \end{aligned}$$

(12)

with transversality conditions $\lambda _{i}(T_{f})=0$, for $i=1,2,\ldots ,7$ and $u_{1}^{*},\ u_{2}^{*}$ and $u_{3}^{*}$ are same as (8).