Introduction

Many types of viruses in the coronavirus family can infect people and cause illnesses like the common cold and severe acute respiratory syndrome (SARS). Two epidemics of coronaviruses, SARS and Middle East Respiratory Syndrome Coronavirus (MERS), have been reported in the past two decades with more than 10 thousand confirmed cases and 1600 deaths1,2,3. Initially, the infection of Middle East Respiratory Syndrome Coronavirus (MERS) was identified in Saudi Arabia while spreading in many other countries4. In late 2019, a severe respiratory illness identified with a causative agent isolated from a single patient reported in Wuhan, China, known as the novel coronavirus (2019-nVoV)5. The scientific data clarify that animals were the virus’s primary transmission source, but most cases increase due to interaction between vulnerable and infected people. The signs of a new coronavirus infection include fever, cough, exhaustion, breathing issues, etc. The novel coronavirus was a public health emergency of international concern, according to the World Health Organization (WHO), which declared it a global pandemic due to a pressing issue.

Due to the new nature of the disease, the coronavirus has been identified as a global danger and has drawn the interest of numerous researchers. Infectious disease epidemiology has a rich literature (see for detail6,7,8,9). Various epidemiological models investigated the dynamics and control analysis of different infectious diseases10,11,12. Similarly, the new coronavirus has a rich literature, and numerous articles have been reported to study the future spread of the disease dynamics. Particularly, a model was introduced to investigate the new coronavirus transmission by Wu et al. in13. To estimate the outbreak based on human interaction, a computational study has been reported by Imai et al.14. Further, the infectivity of the newly reported virus has been investigated by Zhu et al.15. Likewise, various other studies have been performed to investigate the dynamics of 2019-nCoV virus transmission16,17,18,19,20.

The biological interpretation of the novel coronavirus reveals that the pandemic rises due to close contact of human interactions but the initial transmission source is an animal. Therefore environmental reservoirs play an essential role in the spreading of novel coronavirus transmission. Thus, we assume both transmission sources according to the characteristics of the novel infection and follow the classical susceptible-infected-recovered (SIR) model to develop a new mathematical model while studying the temporal dynamics of the novel coronavirus disease propagation. More specifically, we assume humans and reservoirs as a source of disease transmission and propose an epidemiological model by considering the extending version of the classical susceptible-infected-recovered (SIR) model. In addition, we formulate an optimal control mechanism to minimize novel coronavirus spreading in the community. We first show the model validity to discuss the mathematical and biological properties of the epidemic problem that is under consideration. We calculate the threshold parameter and show its sensitivity to estimate the role of every epidemic parameter involved in the modeling process of the proposed problem. It will clarify the importance and role of every model parameter. We also investigate the stability of the model to find the stability conditions using linear stability analysis. Further, we develop a control mechanism based on the temporal dynamics of the model and sensitivity analysis to identify the best control policy for eradicating the disease. A large-scale numerical simulation will be also performed to validate the theoretical results of the model by providing numerical experiments.

Mathematical formulation of the model

We examine the dynamics of the newly reported disease in this section by keeping in view the properties of 2019-nCoV to propose the model, whose schematic process is given in the flowchart as depicted in Fig. 1. For this purpose, the different compartments of humans and reservoirs are taken. We also take multiple transmission routes, i.e. from humans and reservoirs. We impose some assumptions on the model stated as:

  1. 1.

    All the state variables and parameters are taken to be positive or non-negative.

  2. 2.

    The total human population is classified into three different compartments and one class of reservoir.

  3. 3.

    Humans and reservoirs play a crucial role in the rise of the pandemic, so both are taken as a source of disease transmission.

  4. 4.

    All newborns will lead to the susceptible group of individuals because there is still no evidence regarding vertical transmission of the disease.

Thus, the epidemic problem can take the combination of the coupled nonlinear differential equation which looks like:

Figure 1
figure 1

The schematic diagram for the transmission of 2019-nCoV.

Full size image
$$\begin{aligned} {\left\{ \begin{array}{ll} \begin{aligned} \frac{dS_{h}(t)}{dt}& =\Pi -\left\{ \beta _1I_h(t)+\beta _2W(t)+d\right\} S_h(t),\\ \frac{dI_{h}(t)}{dt}& =\left\{ \beta _1I_h(t)+\beta _2W(t)\right\} S_h(t)-\left\{ d+\sigma +d_1\right\} I_h(t),\\ \frac{dR_{h}(t)}{dt}& =\sigma I_h(t)-dR_h(t),\\ \frac{dW(t)}{dt}& =\alpha I_h(t)-\eta W(t), \end{aligned} \end{array}\right. } \end{aligned}$$
(1)

with non-negative initial size of population

$$\begin{aligned} S_h(0)>0,\quad I_h(0),~R_h(0),~W(0)\ge 0, \end{aligned}$$
(2)

where \(\Pi\) is the proportion of newborns. \(\beta _1\) and \(\beta _2\) are the disease transmission rates occurring from humans and reservoirs, respectively. d demonstrates the natural mortality of humans while \(d_1\) is the death rate of infected humans. The recovery/removal rate of infected individuals is denoted by \(\sigma\). The parameter \(\alpha\) is the ratio at which the infected individuals contribute virus in terms of the environmental reservoirs, while \(\eta\) is the removal rate of the virus.

Existence analysis and positivity

We show the mathematical and biological feasibility of the epidemic problem that is under consideration. Since the sum of the humans population is define as \(N(t)=R_h(t)+I_h(t)+S_h(t)\), which implies that \(0\le N(t)\le \frac{\Pi }{d}\). Now, from the last equation of the system (1), we can re-write

$$\begin{aligned} \frac{dW(t)}{dt}+\eta W(t)=\alpha I_h(t), \end{aligned}$$

implies that

$$\begin{aligned} 0<W(t)\le \frac{\Pi \alpha }{d}\left\{ \frac{N(0)-1}{\eta -d}\left( e^{-dt}-e^{-\eta t}\right) +\frac{1}{\eta }\left( 1-e^{-\eta t}\right) \right\} +N(0)e^{-\eta t}. \end{aligned}$$

Clearly, W(t) is bounded; thus, all the model solutions are bounded as well as positively invariant, and the feasible region is given as

$$\begin{aligned} \Theta _1 = \bigg \{(S_h,I_h,R_h,W)\in \mathbb {R}_{+}^{4}:N\le \frac{\Pi }{d}\quad \text {and}\quad W\le \frac{\Pi \alpha }{d\eta }\bigg \}. \end{aligned}$$
(3)

In addition, let us assume that \(t_{+}>0\) and B represents Banach space, then

$$\begin{aligned} B = x^{1}(0,t_{+})\times x^{1}(0,t_{+})\times x^{1}(0,t_{+})\times x^{1}(0,t_{+}). \end{aligned}$$
(4)

We define the norm on B as \(\Vert \psi \Vert =\sum _{j=1}^{4}\Vert \psi _j\Vert ,~\text {where}~\psi =(\psi _1,\psi _2,\psi _3,\psi _4)\in B\) and suppose that the positive cone of \(x^{1}(0,t_{+})\) is \(B_{+}\), then from Eq.(4), \(B_{+}=x^{1}_{+}(0,t_{+})\times x^{1}_{+}(0,t_{+})\times x^{1}_{+}(0,t_{+})\times x^{1}_{+}(0,t_{+})\). Thus the state space of the proposed system (1), takes the following form

$$\begin{aligned} \Theta _2 = \bigg \{(S_h,I_h,R_h,W)\in B_{+}: 0\le N \le \frac{\Pi }{d}\quad \text {and}\quad W\le \frac{\Pi \alpha }{d\eta }\bigg \}. \end{aligned}$$
(5)

Let L is the linear operator and \(\phi =(S_h,I_h,R_h,W)\), then \(L\phi =(L_1,L_2,L_3,L_4)^{T}\), where

$$\begin{aligned} \begin{aligned} L_1&=\bigg (-\frac{dS_h}{dt}-dS_h,0,0,0,0\bigg ),~L_2=\bigg (0,-\frac{dI_h}{dt}-(\sigma +d+d_1)I_h,0,0\bigg ),\\ L_3&=\bigg (0,\sigma I_h,\frac{dR_h}{dt}-dR_h,0\bigg ),~L_4=\bigg (0, \alpha I_h,0,\frac{dW}{dt}-\eta W\bigg ), \end{aligned} \end{aligned}$$
(6)

and D(L) is the domain, then \(D(L)=\bigg \{\phi \in B:~S_h,I_h,R_h,W\in LC[0,t_{+}),\phi (0)=\big (S_h(0),I_h(0),R_h(0),W(0)\big )\bigg \}\). In which \(LC[0,t_{+})\) is the set of absolutely continues functions on \([0,t_{+})\). Now we define the nonlinear operator M, such that \(M:B\rightarrow B\) by

$$\begin{aligned} M\phi =\left( \begin{array}{c} \Pi -\big (\beta _1I_h+\beta _2W\big )S_h\\ \big (\beta _1I_h+\beta _2W\big )S_h\\ 0 \\ 0 \end{array} \right) . \end{aligned}$$
(7)

Let \(w(t)=(S_h(.,t),I_h(.,t),R_h(.,t),W(.,t))\), then the proposed system (1) can be written as

$$\begin{aligned} \frac{dw(t)}{dt}=L(w(t))+M(w(t)),~w(0)\in B, \end{aligned}$$
(8)

where \(w(0)=(s_h(0),I_h(0),R_h(0),W(0))^{T}\). Now to prove the existence of solution of the above system (8), we follow21,22, thus we have the following results.

Theorem 2.1

For each \(w(0)\in B_{+}\), there exist a maximal interval of existence \([0, t_0)\) and a unique continues mild solution \(w(t,w_0)\in B_{+}\), \(t\in [0, t_0)\) for Eq.(8), such that

$$\begin{aligned} w(t)=w(0)e^{Lt}+\int _0^{t}e^{L(t-\tau )}M(w(\tau ))d\tau . \end{aligned}$$
(9)

Proposition 2.1

Let us assume the initial conditions (2) of the proposed model, then the solutions of the proposed model are positive for all \(t_{+}\).

Proof

Let \(I\subset [0,t_{+})\), then the solution of the first equation of system (1) looks like

$$\begin{aligned} \begin{aligned} S_h(t)&=S_h(0)\exp \bigg \{-dt-\int _0^t\big (\beta _1I_h(x)+\beta _2W(x)\big )dx\bigg \}\\&\quad +\Pi \exp \bigg \{-dt-\int _0^t\big (\beta _1I_h(x)+\beta _2W(x)\big )dx\bigg \}\\&\quad \times \int _0^t\bigg \{dt+\int _0^x\big (\beta _1I_h(y)+\beta _2W(y)\big )dy\bigg \}dx>0. \end{aligned} \end{aligned}$$
(10)

Obviously the L.H.S of Eq. (10) is positive i.e., \(S_h(t)>0\). Similarly, it can be shown that the solution of the second, third, and fourth equations of the model (1) are non-negative. \(\square\)

Stability analysis

We investigated that the proposed epidemic problem is mathematically and biologically feasible. To discuss the temporal dynamics, we perform stability of the model with the aid of linear stability theory. We calculate the model equilibria and find the threshold quantity (basic reproductive number). We also derive the stability conditions to characterize whether the disease dies out or persists.

Threshold quantity and dynamical analysis of the model equilibria

To find the stability conditions, we study the qualitative behavior of the proposed epidemic problem. We first calculate the disease-free equilibrium and then the threshold quantity. For this, we set all the states equations of the model equating to zero except \(S_h=S_{h0}\), which gives \(X_0=(S_{h0},0,0,0,0)\) is the disease-free state of the considered model, where \(S_{h0}=\frac{\Pi }{d}\).

In addition, the threshold quantity represents the average number of secondary infections produced by a single infective whenever introduced into a susceptible population. This will identify what happens if an infected agent enters into the population of susceptible. One of the examples of the 2019-nCoV is given by Fig. 2, which suggests that the coronavirus spreads faster than the estimation of WHO. This study was organized between Jan 2020 and Feb 2020, and the estimate ranged from 1.4-6.49 with 3.28 of average and 2.79 of a median23. In this case, the estimated threshold value is high, and it is essential to control the pandemic by reducing the weight of the reproductive number. So here, we are going to figure out the threshold quantity (\(\mathcal {R}_0\)) of the model that is under consideration following the methodology of24. We then discuss the detailed sensitivity of the threshold parameter to find the relation of all parameters to the disease propagation. This type of analysis will be beneficial to optimize the value of the threshold quantity. Based on the next-generation matrix method, we calculate matrix F, and V is given by

Figure 2
figure 2

Estimated values of the threshold quantity \((\mathcal {R}_0\)) for the 2019-nCoV virus in China.

Full size image
$$\begin{aligned} F=\left( \begin{array}{cc} \beta _1S_{h0} & \beta _2S_{h0} \\ 0 & 0 \\ \end{array} \right) ,~V=\left( \begin{array}{cc} \sigma +d+d_1 & 0 \\ -\alpha & \eta \\ \end{array} \right) . \end{aligned}$$
(11)

The threshold quantity is the dominant eigenvalue of the matrix (\(\bar{H}=FV^{-1}\)), given by

$$\begin{aligned} \mathcal {R}_0=\frac{\beta _1\Pi }{d(\sigma +d+d_1)}+\frac{\alpha \beta _2\Pi }{\eta d(\sigma +d+d_1)}. \end{aligned}$$
(12)

From the above Eq. (12), it is clear that the threshold quantity for the proposed problem consists of two parts, which describe that there are two different transmission routes, one from the infected individuals and the other from reservoirs to the susceptible populations.

Sensitivity analysis

To find the relation of every parameter with the threshold quantity to the disease transmission and then its severeness, we use the formula \(X_{\gamma }=\frac{\gamma }{\mathcal {R}_0}\frac{\partial \mathcal {R}_0}{\partial \gamma }\) by following25. Once we calculate the indices of the threshold quantity, we can analyze the sensitivity of each parameter to the disease transmission and its prevalence. Using the above formula, we obtain the following indices

$$\begin{aligned} \begin{aligned}&X_{\beta _1}=\frac{\eta \beta _1}{\eta \beta _1+\alpha \beta _2}>0,~X_{\beta _2}=X_{\alpha }=\frac{\alpha \beta _2}{\eta \beta _1+\alpha \beta _2}>0,\\&X_{\sigma }=-\frac{\sigma \beta _1\eta +\alpha \beta _2}{\eta \beta _1+\alpha \beta _2}<0,~X_{\eta }=-\frac{\alpha \beta _2}{\eta \beta _1+\alpha \beta _2}<0. \end{aligned} \end{aligned}$$
(13)

It is very much clear from the Eq. (13) that the relation between \(\mathcal {R}_0\) and parameters \(S_{1}=\{\beta _1,\beta _2,\sigma \}\) is direct. But it is indirect with the parameters given by \(S_{2}=\{\sigma ,\eta \}\), which certify that if the parameters value of \(S_1\) increases, the value of the threshold parameter will also increase significantly. In contrast, a decrease in the parameters value of \(S_2\) causes increases in threshold quantity. Now to control the spread of coronavirus disease, we focus on optimizing the value of the parameters that have a direct relation with threshold quantity and maximizing the value of those parameters with an inverse relation. More precisely, we take some feasible values of the parameters to find the percentage of the sensitivity indices and their relative impact on the threshold quantity numerically. Let us assume the biological feasible values of the parameters, i.e., \(\beta _1=0.05\), \(\beta _2=0.000001231\), \(\sigma =0.09871\), \(\eta =0.01\) and \(\alpha =0.000398\) to find the numerical value of the sensitivity indices and their relative percentage impact on the threshold quantity. Using the above numerical values, we calculate the following indices, i.e., \(X_{\beta _1}=0.9999990210\), \(X_{\beta _2}=X_{\alpha }=0.0000979079\), \(X_{\sigma }=-0.0987008824\), \(X_{\eta }=0.0000979079\). The biological interpretation of these sensitivity indices shows that the increase or decrease in the value of \(\beta _1\), say, for example, 10 percent, would increase or decrease the value of \(\mathcal {R}_0\) by 9.999990210 percent as shown Fig. 4. In addition, if we increase or decrease the value of \(\beta _2\) and \(\alpha\) by 10 percent, collectively affects the value of the threshold quantity by 0.001958158 percent, see Fig. 3b. Similarly, from the sensitivity indices of \(\sigma\) and \(\eta\), we observe that the influence of these parameters on the threshold quantity is negative, in which an increase leads to a decrease in threshold quantity. For example, an increase of 10 percent in the value of \(\sigma\) and \(\alpha\) would decrease the value of the threshold quantity by 0.9870088245 and 0.0009790790415 respectively, for instance, see Fig. 3a,b.

Figure 3
figure 3

The graphical results show the sensitivity analysis of the threshold quantity (\(\mathcal {R}_0\)) and its relative impact against the variations of various epidemic parameters \(\{\beta _1,\beta _2,\alpha ,\sigma ,\eta \}\). For this analysis, the numerical values of the parameters used are \(\beta _1= 0.05\) \(\sigma =0.09871\), \(\eta =0.01\), \(\beta _2=0.00000123\), \(\Pi =0.0033907997\), \(\alpha =0.000398\), \(d_1=0.0404720925\), \(d=0.0003567816\).

Full size image
Figure 4
figure 4

The graphical results show the sensitivity analysis of the threshold quantity (\(\mathcal {R}_0\)) and its relative impact against the variation of various epidemic parameters \(\beta _1\). We use the numerical values of other parameters \(\Pi =0.0033907997\), \(\alpha =0.000398\), \(d_1=0.0404720925\), \(d=0.0003567816\).

Full size image

Dynamical analysis

To discuss the model dynamics at the disease-free state using the theory of dynamical systems, we state the following results.

Theorem 3.1

If \(\mathcal {R}_0\) is less than unity, i.e. \(\mathcal {R}_0<1\), then the model is stable asymptotically at \(X_0=(S_{h0},0,0,0,0)\).

Proof

Since the considered model reveals that \(R_h(t)\) is not explicitly involved in all other equations, we assume only the three equations of the model, while investigating the dynamics of the proposed epidemic problem. Let \(A(X_0)\) is the Jacobian the model (1) at \(X_0\), then

$$\begin{aligned} A(X_0)=\left( \begin{array}{ccc} -d & -\frac{\beta _1 \Pi }{d} & -\frac{\beta _2 \Pi }{d} \\ 0 & \frac{\beta _1 \Pi }{d}-\sigma -d-d_1 & \frac{\beta _2 \Pi }{d} \\ 0 & \alpha & -\eta \end{array} \right) . \end{aligned}$$
(14)

Obviously, \(A(X_0)\) has an eigenvalue, \(\lambda _1=-d_1<0\). For the remaining, we take the following reduced matrix

$$\begin{aligned} A(X_0)=\left( \begin{array}{cc} \frac{\beta _1 \Pi }{d}-\sigma -d-d_1 & \frac{\beta _2 \Pi }{d} \\ \alpha & -\eta \end{array} \right) . \end{aligned}$$
(15)

The local dynamics of the epidemic problem depend on the nature of the eigenvalues of the Jacobian matrix. So, \(\lambda _1\) is negative while the eigenvalues of (15) will be negative if the Routh- Hurwitz criteria hold. We noted that the trace and determinant of matrix (15) are respectively negative and positive if the \(\mathcal {R}_0\) is less than unity, such that

$$\begin{aligned} trace(A(X_0)) =-(\sigma +d+d_1)(1-\mathcal {R}_0)-\bigg (\eta +\frac{\alpha \beta _2 \Pi }{d}\bigg ) \end{aligned}$$
(16)

and

$$\begin{aligned} det(A(X_0)) = \eta (\sigma +d+d_1)(1-\mathcal {R}_0). \end{aligned}$$
(17)

Equations (16)–(17) implies that Routh–Hurwitz criteria satisfies under the condition that \(\mathcal {R}_0<1\), so we conclude that the model (1) is stable asymptotically in local sense if \(\mathcal {R}_0<1\).

For the global dynamics around \(X_0\), we define a function \(H_1:R_{+}^4\rightarrow R\), such that

$$\begin{aligned} H_1(S_h,I_h,R_h,W)=a_1(S_h-S_{h0})+a_2I_h+a_3R_h+a_4W. \end{aligned}$$
(18)

In the above Eq. (18), \(a_i\), \(i=1,2,3,4\) are arbitrary positive constants. We differentiate \(H_1\) temporally and assume the values of the positive constants, i.e., \(a_1=a_2=\sigma +d+d_1\), \(b_3=\eta \beta _1S_{h}+\beta _2\alpha S_h\), then arrives at the following equation

$$\begin{aligned} \frac{dH_1(t)}{dt}=-(\sigma +d+d_1)(S_h-S_{h0})-(1-\mathcal {R}_0)\sigma I_h. \end{aligned}$$

Obviously \(\frac{dH_1(t)}{dt}\) is negative if \(S_h>S_{h0}\) and \(\mathcal {R}_0<1\). Moreover, \(\frac{dH_1(t)}{dt}=0\) if \(S_h=S_{h0}\) and \(I_h=R_h=W=0\). Therefore, the well-known LaSalle invariance principle or simply invariance principle26 implies that \(X_{0}\) is globally asymptotically stable. \(\square\)

To determine the dynamics of the model at the endemic equilibrium, we assume that \(X_{*}\) be the endemic equilibrium of the proposed model, which is obtained by setting \(S_h=S_h^{*}\), \(I_h=I_h^{*}\), \(R_h=R_h^{*}\), \(W=W^{*}\) at the steady state, then the corresponding disease endemic equilibrium leads to \(X_{*}=(S_h^{*}, I_h^{*}, R_h^{*}, W^{*})\), where the following system of equation defines the components of the endemic state:

$$\begin{aligned} \begin{aligned} S_h^{*}&=\frac{\eta \left( \sigma +d+d_1\right) }{\eta \beta _1+\alpha \beta _2},\quad R_h^{*}=\frac{\sigma }{d}I_{h}^{*},\\ I_h^{*}&=\frac{\eta d\left( \sigma +d+d_1\right) \left( \mathcal {R}_0-1\right) }{\sigma \beta _1\eta +\sigma \beta _2\alpha +d\beta _1\eta +d\beta _2\alpha +d_1\beta _1\eta +d_1\beta _2\alpha },\\ W^{*}&= \frac{\alpha \eta \left( \sigma +d+d_1\right) \left( \mathcal {R}_0-1\right) }{\sigma \beta _1\eta +\sigma \beta _2\alpha +d\beta _1\eta +d\beta _2\alpha +d_1\beta _1\eta +d_1\beta _2\alpha }. \end{aligned} \end{aligned}$$
(19)

It can be noted from the above Eq. (19) that the endemic state of the model exists only in the case whenever \(\mathcal {R}_0>1\). Thus we state the following result.

Lemma 3.2

The endemic state \(X_{*}=(S_h^{*}, I_h^{*}, R_h^{*}, W^{*})\), of the proposed problem (1) exists only in the case \(\mathcal {R}_0>1\).

To discuss the model temporal dynamics around the endemic state we prove the following theorem.

Theorem 3.3

If the conditions \(\mathcal {R}_0>1\) and \(d+d_1>\alpha\) holds then \(X_{*}=(S_{h}^{*},I_{h}^{*},R_{h}^{*},W^{*})\) is asymptotically stable.

Proof

We use the theory of linear stability analysis to discuss the dynamics of the model at \(X_{*}\). Let \(A(X_{*})\) be the Jacobian matrix of the model (1) around \(X_{*}\), then

$$\begin{aligned} A(X_{*})=\left( \begin{array}{ccc} -(\beta _1I_{h}^{*}+\beta _2W^{*}+d) & -\beta _1S_{h}^{*} & -\beta _2S_{h}^{*} \\ \beta _1I_{h}^{*}+\beta _2W^{*} & \beta _1S_{h}^{*}-(\sigma +d+d_1) & \beta _2S_{h}^{*} \\ 0 & \alpha & -\eta \end{array} \right) . \end{aligned}$$
(20)

To discuss the local dynamics, we find the characteristic polynomial of (20) as

$$\begin{aligned} X(\lambda )=\lambda ^3+h_1\lambda ^2+h_2\lambda +h_3, \end{aligned}$$
(21)

where

$$\begin{aligned} \begin{aligned} h_1&=\frac{(\eta +d)(\eta \beta _1+\alpha \beta _2)+(\eta d \beta _1+\alpha \beta _2)(\mathcal {R}_0-1)+\alpha \beta _2(\sigma +d+d_1)}{\eta \beta _1+\alpha \beta _2},\\ h_2&=\frac{(\eta d^2\beta _1+\eta d\beta _1d_1+\eta ^2d\beta _1+\alpha \beta _2d+\alpha \beta _2d_1+\alpha \beta _2+\sigma \eta d\beta _1+\alpha \sigma \beta 2)(\mathcal {R}_0-1)}{\eta \beta _1+\alpha \beta _2}\\&\quad +\frac{\alpha \beta _2d(\sigma +d+d_1)+\eta d(\eta \beta _1+\alpha \beta _2)}{\eta \beta _1+\alpha \beta _2},\\ h_3&=\frac{(\eta ^2d^2\beta _1+\eta ^2d\beta _1d_1+\alpha \beta _2d+\alpha \beta _2d_1+\eta ^2d\beta _1\sigma +\alpha \beta _2\sigma )(\mathcal {R}_0-1)}{\eta \beta _1+\alpha \beta _2}. \end{aligned} \end{aligned}$$
(22)

According to linear stability analysis, the proposed model around endemic equilibrium (19) is stable if all eigenvalues of the Jacobian matrix (20) are negative. One may follow the Routh- Hurwitz criteria to discuss the properties of the Jacobian matrix (20). Roots (eigenvalues) of the above Eq. (21) are negative, if Routh-Hurwitz criteria \((H_3)\) i.e. \(h_1>0\), \(h_3>0\) and \(h_1h_2-h_3>0\) holds. It can be noted from Eq. (22) that \(h_1\) and \(h_3\) are positive under the condition that \(\mathcal {R}_0>1\), while \(h_1h_2-h_3>0\) is calculated and given by the following expression

$$\begin{aligned} \begin{aligned} h_1h_2-h_3&=((\eta +d)(\eta \beta _1+\alpha \beta _2)+(\eta d\beta _1+\alpha \beta _2)(\mathcal {R}_0-1)\\&\quad +\alpha \beta _2(\sigma +d+d_1))(\alpha \beta _2d(\sigma +d+d_1)+\eta d(\eta \beta _1+\alpha \beta _2))\\&\quad +(\eta d\beta _1(1-\eta )(d+d_1)+\eta ^2d\beta _1(1-\sigma )+\alpha \beta _2+\sigma \eta d\beta _1)(\mathcal {R}_0-1). \end{aligned} \end{aligned}$$
(23)

Clearly \(h_1h_2-h_3>0\), whenever \(\mathcal {R}_0>1\), so \(H_3\) satisfied if \(\mathcal {R}_0>1\). Thus, model (1) is asymptotically stable at \(X_{*}\) if \(\mathcal {R}_0>1\).

For the global dynamics, we assume a function \(H_2:R_{+}^4\rightarrow R\) by the following equation

$$\begin{aligned} H_2(S_h,I_h,W)=\frac{1}{2}\big \{(S_h-S_h^{*})+(I_h-I_h^{*})+(R_h-R_h^{*})+(W-W^{*})\big \}^2. \end{aligned}$$
(24)

We differentiate \(H_2\) and making use of some simple algebra, we get

$$\begin{aligned} \begin{aligned} \frac{dH_2(t)}{dt}=-&\big \{(S_h-S_h^{*})+(I_h-I_h^{*})+(R_h-R_h^{*})+(W-W^{*})\big \}\\&\big \{d(S_h-S_0)+(d+d_1-\alpha )I_h-dR_h-\eta W\big \}. \end{aligned} \end{aligned}$$

It could be noted from the above equation that \(\frac{dH_2(t)}{dt}\le 0\) if and only if \(S_h>S_{h}^{*}\), \(d+d1>\alpha\), and \(S_h=S_{h}^{*}\), \(I_h=I_h^{*}\), \(R_h=R_h^{*}\), \(W=W^{*}\). Thus we conclude that according to the invariance principle the disease endemic state \(X_{*}\) is globally asymptotically stable, if \(\mathcal {R}_0>a\) and \(d+d_1>\alpha\). \(\square\)

Optimal control

Keeping in mind the model dynamics and sensitivity of the threshold parameter, we analyse a control mechanism with the aid of optimal control theory. It can be noted that the maximum sensitivity index parameter in the epidemic parameters is \(\beta _1\), whose value is 0.9999990210. This indicates that the increase of 10 percent in \(\beta _1\) causes the increase of 9.999990210 percent in the value of \(\mathcal {R}_0\). Therefore, the very first step to control the spread needs much attention to minimize the value of this parameter. With the help of the control variable, isolation of infected and non-infected (\(u_1(t)\)), we can reduce and optimize the spreading in the community. Similarly, the index of \(\beta _2\) is 0.000097079, which affects the threshold quantity by 0.00097079 percent. To minimize the value of this parameter, we use the control variable, (\(u_2(t)\)), which physically represents personal protection, e.g., wearing of mass, washing hands, etc. The index of \(\alpha\) is 0.000097079, which variates the threshold quantity by 0.00097079 percent. We must reduce this by taking the \(u_4(t)\) control variable. Moreover, the collective variation of \(\sigma\) and \(\eta\) affect the threshold quantity by 0.9870088245 and 0.0009790790415. Here, we need to increase the values of these parameters because the increase would decrease the threshold quantity value. For this, we assume \(u_3(t)\), a control measure representing infected individuals’ treatment. Thus we are now in a position to take the set of control variables symbolized by u(t) and define as \(u(t)=\{u_1(t),u_2(t),u_3(t),u_4(t)\}\), where \(u_1(t)\) is used for the isolation of infected individuals from the community, while \(u_2(t)\) represents the personal protection including wearing of masks, washing hands, etc. Similarly, the control variable \(u_3(t)\) represents the treatment of infected individuals, as it is clear that Corona has no proper treatment, so this control means that treatment which strengthens the immunity system and \(u_4(t)\) is the control variable which is used for minimizing the value of the parameter \(\alpha\). Placing the control variables, we arrive at the following control problem:

$$\begin{aligned} J(u,x)=\min \int _0^T\bigg (\omega _1I_{h}(t)+\omega _2W(t)+\frac{1}{2}(\omega _3^2u_1^2(t)+\omega _4^2u_2^2(t)+\omega _5^2u_3^2(t)+\omega _6^2u_4^2(t))\bigg )dt, \end{aligned}$$
(25)

subject to the system, which is the extended version of Eq. (1) as defined by

$$\begin{aligned} \begin{aligned} \frac{dS_{h}(t)}{dt}&=\Pi -\beta _1S_h(t)I_h(t)(1-u_1(t))-\beta _2S_h(t)W(t)(1-u_2(t))-dS_h(t),\\ \frac{dI_{h}(t)}{dt}&=\beta _1I_h(t)S_h(t)(1-u_1(t)+\beta _2W(t)S_h(t)(1-u_2(t)\\&\quad -(\sigma +d+d_1+u_3(t))I_h(t),\\ \frac{dR_{h}(t)}{dt}&=(\sigma +u_3(t)) I_h(t)-dR_h(t),\\ \frac{dW(t)}{dt}&=\alpha I_h(t)-(\eta +u_4(t))W(t), \end{aligned} \end{aligned}$$
(26)

with

$$\begin{aligned} S_h(0)>0,\quad I_h(0),~R_h(0),~W(0)\ge 0. \end{aligned}$$
(27)

In the above Eq. (25), \(x=(S_h,I_h,R_h,W)\) and \(\omega _1\), \(\omega _2\), \(\omega _3\), \(\omega _4\), \(\omega _5\) and \(\omega _6\) are the weight constants in which \(\omega _1\) and \(\omega _2\) represent the relative cost of the infected population and reservoir, while \(\omega _3\), \(\omega _4\), \(\omega _5\) and \(\omega _6\) measure the associated cost of \(u_1(t)\), \(u_2(t)\), \(u_3(t)\) and \(u_4(t)\), respectively. The goal of our control problem is very clear from the objective functional (25): to control the disease by maximizing the recovered and reducing the infected, and the reservoir by taking into account the cost of controls. We obtain the cost function at \((u_1^{*},u_2^{*},u_3^{*},u_4^{*})\), such that

$$\begin{aligned} J(u_1^{*},u_2^{*},u_3^{*},u_4^{*})=\min \{J(u_1,u_2,u_3,u_4), u_i\in U,~\text {for}~i=1\ldots 4\} \end{aligned}$$
(28)

subject to the control system (26), where U in the above Eq. (28) is called control set and defined by the following equation

$$\begin{aligned} U:=\{(u_1,u_2,u_3,u_4)|u_i~\text {is Lebesgue measurable on}~[0,T],0\le u_i(t)\le 1, 1\ldots 4\}. \end{aligned}$$
(29)

In addition, first, we prove that such controls exist. Following the methodology given in27, demonstrates the existence of a solution for the proposed control system whenever the controls are bounded, and Lebesgue’s measurable along with the initial conditions are non-negative. Re-writing the control problem as

$$\begin{aligned} \frac{d\Phi }{dt}=A\Phi +B(\Phi ), \end{aligned}$$
(30)

In the above system (30), \(\Phi =(S_h,I_h,R_h,W)^T\), while \(B(\Phi )\) and A contain the non-linear and bounded co-efficient, and linearize matrix, respectively, such that

$$\begin{aligned} \begin{aligned} A&=\left( \begin{array}{cccc} -d & 0 & 0 & 0 \\ 0 & -(\sigma +d+d_1+u_3(t)) & 0 & 0 \\ 0 & (u_3(t)+\sigma ) & -d & 0 \\ 0 & 0 & \alpha & -(\eta +u_4(t)) \end{array} \right) ,\\ B(\Phi )&=\left( \begin{array}{cccc} \Pi -\beta _1 I_h(t)S_h(t)(1-u_1(t))-\beta _2 S_h(t)W(t)(1-u_1(t)) \\ \beta _1 S_h(t)I_h(t)(1-u_1(t))+\beta _2 S_h(t)W(t)(1-u_1(t))\\ 0 \\ 0 \end{array} \right) . \end{aligned} \end{aligned}$$

Setting \(L(\Phi )=A\Phi +F(\Phi )\), satisfies that

$$\begin{aligned} \begin{aligned} |F(\Phi _1)-F(\Phi _2)|&\le n_1|S_{h1}-S_{h2}|+n_2|I_{h1}-I_{h2}|+n_3|R_{h1}-R_{h2}|+n_4|W_{1}-W_{2}|,\\&\le N(|S_{h1}-S_{h2}|+|I_{h1}-I_{h2}|+|R_{h1}-R_{h2}|+|W_1-W_2|), \end{aligned} \end{aligned}$$

where \(N=\max (n_1,n_2,n_3,n_4)\). Also, we have

$$\begin{aligned} |L(\Phi _1)-L(\Phi _2)| \le M|\Phi _1-\Phi _2|, \end{aligned}$$

where \(M=\max (N,||A||)<\infty\), so L is uniformly Lipschitz continuous and \(S_h(t)\), \(I_h(t)\), \(R_h(t)\) and W(t) having non-negative values which implies the existence of solution to the system (26).

Thus, we describe the existence analysis with the aid of the following theorem regarding.

Theorem 4.1

There exists an optimal solution \(u^{*}=(u_1^{*},\ldots ,u_4^{*})\in U\) to system (25)–(29).

Proof

Obviously, the model states and control functions have non-negative values while set U is closed and convex. In addition, the control system is bounded implying compactness. Also, the integrand in Eq.(25) is convex, therefore it investigates the existence of optimal controls. \(\square\)

Optimality condition

We find the optimal solution to the control problem (25)–(26). We define the Hamiltonian and Lagrangian by assuming that if x is the state variable and u is a control function, then

$$\begin{aligned} L(x,u)=\omega _1I_h(t)+\omega _2W(t)+\frac{1}{2}\big (\omega _3^2u_1^2(t)+\omega _4^2u_2^2(t)+\omega _5^2u_3^2(t)+\omega _6^2u_4^2(t)\big ), \end{aligned}$$
(31)

and

$$\begin{aligned} H(x,u,\lambda )=L(x,u)+\lambda f(x,u), \end{aligned}$$
(32)

where \(\lambda =(\lambda _1,\ldots ,\lambda _4)\) and \(f=(f_1,\ldots ,f_4)\), and

$$\begin{aligned} \begin{aligned} f_1&=\Pi -\beta _1S_h(t)I_h(t)(1-u_1(t))-\beta _2S_h(t)W(t)(1-u_2(t))-dS_h(t),\\ f_2&=\beta _1I_h(t)S_h(t)(1-u_1(t)+\beta _2W(t)S_h(t)(1-u_2(t)\\&-(\sigma +d+d_1+u_3(t))I_h(t),\\ f_3&=(\sigma +u_3(t)) I_h(t)-dR_h(t),\\ f_4&=\alpha I_h(t)-(\eta +u_4(t))W(t). \end{aligned} \end{aligned}$$

We now use the well-known Pontryagin Maximum Principle28 to find the optimal solution. We assume that if \((x^{*},u^{*})\) is an optimal solution, then there exists a nontrivial function \(\lambda\), satisfying

$$\begin{aligned} \frac{dx^{*}(t)}{dt}=\frac{\partial H}{\partial \lambda },~ 0=\frac{\partial H}{\partial u},~\frac{d\lambda (t)}{dt}=\frac{\partial H}{\partial x}, \end{aligned}$$
(33)

and

$$\begin{aligned} H(x^{*},u^{*},\lambda )=\max _{u\in [0,1]\times [0,1]}H, \end{aligned}$$
(34)

which is the maximality condition and \(\lambda (T)=0\), the transversal condition holds. Using Eq. (33) to find the adjoint system and optimal value for controls with the aid of the result given below.

Theorem 4.2

We assume that \(S^{*}_h(t)\), \(I^{*}_h(t)\), \(R^{*}_h(t)\), \(W^{*}(t)\) are the optimal states and \((u_1^{*},u_2^{*},u_3^{*},u_4^{*})\) are optimal control functions for the control problem., the adjoint system looks like

$$\begin{aligned} \begin{aligned} \frac{d\lambda _1(t)}{dt}&=-\big (\lambda _2(t)-\lambda _1(t)\big )\big (\beta _1 I_h(t)(1-u_1(t))+\beta _2 W(t)(1-u_2(t))\big )+\lambda _1(t)d,\\ \frac{d\lambda _2(t)}{dt}&=-\omega _1-\big (\lambda _2(t)-\lambda _1(t)\big )\big (\beta _1 S_h(t)(1-u_1(t))\big )\\&\quad +\big (\lambda _2(t)-\lambda _3(t)\big )\big (\sigma +u_3(t)\big )+\lambda _2(t)(d+d_1)-\alpha \lambda _4(t),\\ \frac{d\lambda _3(t)}{dt}&=d\lambda _4(t),\\ \frac{d\lambda _4(t)}{dt}&=-\omega _2-\big (\lambda _2(t)-\lambda _1(t)\big )\big (\beta _2 S_h(t)(1-u_2(t))\big )+\lambda _4(t)\big (\eta +u_4(t)\big ), \end{aligned} \end{aligned}$$
(35)

with terminal conditions \(\lambda (T)=0\). Furthermore, the set of optimal control functions is defined by \(\bigg \{u_1^{*}=\max \{\min \{l_1,0\},1\},u_2^{*}=\max \{\min \{l_2,0\},1\},u_3^{*}=\max \{\min \{l_3,0\},1\},u_4^{*}=\max \{\min \{l_4,0\},1\}\bigg \}\), where

$$\begin{aligned} \begin{aligned}&l_1=\frac{1}{\omega _3^2}\big (\lambda _2(t)-\lambda _1(t)\big )\beta _1 I_h(t)S_h(t) ,~l_2=\frac{1}{\omega _4^2}\big (\lambda _2(t)-\lambda _1(t)\big )\beta _2W(t)S_h(t),\\&l_3=\frac{1}{\omega _5^2}\big (\lambda _2(t)-\lambda _3(t)\big )I_h(t),~l_4=\frac{1}{\omega _6^2}\lambda _4(t)W(t). \end{aligned} \end{aligned}$$
(36)

Proof

The adjoint variables (35) can be derived from the direct use of the Eq. (33), and the transversal conditions can be obtained from \(\lambda (T)=0\). In addition, to find the optimal value of the controls, we set \(\frac{\partial H}{\partial u}=0\), and solve for \(u_1\), \(u_1\), \(u_1\) and \(u_1\). \(\square\)

Numerical experiments

The verification of our analytical results will be analyzed with the help of some computational findings. The large-scale simulations are based on both quantitative and qualitative analysis. We take some of the parameters from the literature while some are adjusted based on the analytical results. We use purely numerical methods, i.e., Runge-Kutta method (RKM) of the 4th order to simulate the model based on numerical data along with time interval of \(0-500\) units. More precisely, we take the parameters value i.e., \(\beta _1=0.00005\), \(\sigma =0.09871\), \(\eta =0.01\), \(\beta _2=0.0000000123\), \(\Pi =0.093907997\), \(\alpha =0.0398\), \(d_1=0.00404720925\) and \(d=0.009567816\) to verify the analytical findings of the model at the disease-free state. We calculate \(X_0\) and \(\mathcal {R}_0\) for the above parameter’s value as (9.814987767, 0, 0, 0) and 0.004373289660. Thus, we perform the numerical experiments for the model (1) using the above parametric values, and obtain the results as depicted by Fig. 5, which ensures the verification of the result stated by Theorem 3.1. Perturbing the initial sizes of the compartmental population along with the parametric ratios and interval of time, the solution goes to the disease-free equilibrium irrespective of its initial sizes, which ensures that the model is stable at \(X_0\). Further, the theoretical interpretation states that in case of \(\mathcal {R}_0<1\), each solution curve of the susceptible population approximately taking 300 to 400 days to reach the equilibrium value 9.814987767 as shown in Fig. 5a. Similarly, the dynamics of infected and recovered population in the case when \(\mathcal {R}_0<1\) describe that the solution curves of \(I_h\) and \(R_h\) respectively take approximately 50 and 400 days to reach its equilibrium position as shown in Figs. 5b,c. Thus, the biological interpretation states that the elimination of the 2019-nCoV virus from the community will take more time almost 12 to 13 units if \(\mathcal {R}_0<1\). So it is essential to optimize and keep the value \(\mathcal {R}_0\) less than one.

Next, we assume another set of parameters value, i.e., \(\beta _1=0.0005\), \(\sigma =0.009871\), \(\eta =0.01\), \(\beta _2=0.0000123\), \(\Pi =0.93907997\), \(\alpha =0.398\), \(d_1=0.00404720925\) and \(d=0.009567816\) to study the dynamics of the considered problem at the endemic state (19). We calculate the associated threshold quantity and endemic equilibrium \(X_{*}\), defined as: \(\mathcal {R}_0=4.135362579\) and \(X_{*}=(23.73428588,30.31567421,31.27631426,1261.065047)\). Here, if \(\mathcal {R}_0>1\), we take the same initial sizes of the compartmental populations. We observed that the susceptible populations decrease or increase during the initial 100 days of the infection, however after that, there will be no effect, which guarantees the stability of the susceptible individuals as shown in Fig. 6a and therefore approaches to its equilibrium position 23.73428588. The dynamics of the infected population reveal that initially, the ratio of this population increases suddenly during the initial period of the infection, i.e., up to 50 days, while later decreases day by day and attains its endemic position, i.e., 30.31567421 as shown in the Fig. 6b. From this, it could be noted that the disease will reach its endemic position. The 30.31567421 ratio of the infected population will persist in the community if the proper control measure is not implemented. Moreover, we simulate the problem to study the recovered population and reservoir dynamics as shown in Fig. 6c,d. The dynamics of the recovered population describe that in the initial days between 0 to 450 days, the ratio of recovered individuals increases or decreases while then becoming stabilized and attaining its equilibrium position 31.27631426. To control the spread, it should be characterized to increase this ratio. Finally, the ratio of reservoir increases in the initial period of the infection, i.e., between 0 to 400 days according to the dynamics of the proposed problem, while then becoming stable and reaching to its equilibrium 1261.065047. Special attention will be required to sustain the reservoir’s position to zero to control COVID-19.

In order to discuss the impact of control analysis, we use the well-known purely numerical method Runge-Kutta method (RKM) of order 4th. We solve the optimality system, i.e., numerically, we discretize the model (1), the proposed control problem (26), and the adjoint system (35) along with the initial conditions (2)–(4), boundary conditions \(\lambda (T)=0\) and characterization of the control problem. Moreover, the value of the parameter used for this purpose are \(\beta _1=0.85\), \(\sigma =0.04\), \(\eta =0.25\), \(\beta _2=0.46\), \(\Pi =0.0221\), \(\alpha =0.398\), \(d_1=0.08\) and \(d=0.0693\), while the values of weight constants are assumed to be \(\omega _1=0.6\), \(\omega _2=0.9\), \(\omega _3=10.44\), \(\omega _4=0.2\), \(\omega _5=10.90\) and \(\omega _6=10.90\). More precisely, the control and without a control system will be solved by the forward RKM method of order 4th and then the adjoint system by backward RKM method of order 4th with the help of terminal conditions and the solution of system (1). The numerical experiments support the suggested control mechanism as shown in Figs. 7 and 8, while in the case if we ignore the control mechanism and don’t follow the suitable control measure in the form of precautions and treatment etc., it would significantly increase and grow up the spreading of the disease throughout the World. On the other hand, the results describe that the infection will be easily eliminated once the control mechanism imposes in the true spirit.

Cost effective optimal control analysis

We discuss the most cost-effective optimal control measure among the single and combined implementation of the control measures in order to optimally determine the spread of SARS-CoV-2 at the lowest possible cost. We use incremental cost effectiveness ratio (ICER) to explore cost-effectiveness analysis29,30. To avoid dissipation of available limited resources, ICER is calculated to compare any two competing measures for controlling the spread of disease or related problems as given by the ICER formula

$$\begin{aligned} ICER=\frac{\text {Change in total costs between control measures}}{\text {Change in total number of infection averted by control measures}}. \end{aligned}$$

The total cost for each of the single implementation and combined effort of the control measures is calculated from the objective functional (25), while the infection averted is obtained by calculating the difference between infectious individuals without and with control measures. Let \(C_i\) for \(i=1,2,3,4\) and \(C_{14}\) respectively represent the single implementation of each control measure \(u_i (t)\), and combined effort of the four control measures. Table 1 summarizes the ICER for each and combination of the control variables \(u_i (t)\) in an increasing order of the total infection averted.

Using the ICER using the above formula for each \(C_i\) \((i=1,2,3,4)\) and \(C_{14}\) shown in Table 1 are calculated, respectively,

Table 1 ICER in the order of COVID-19 cases averted by control measures.
Full size table

Comparing \(C_i\)’s, it is seen that \({ ICER(C_4)}\) is less than any other \({ ICER(C_i)}\) for \(i=1,2,3\). This means that \(C_i\) for \(i=1,2,3\) are more costly and less effective than \(C_4\). In other words, \(C_1\) dominates \(C_i\) for \(i=1,2,3\). Thus, single implementation of management control is removed from the list. As a consequence, \(C_4\) and \(C_{14}\) are assessed in Table 2 using similar procedure.

It is revealed in Table 2 that \(C_4\) is dominated by \(C_{14}\) since \({ ICER(C_4)}\) is greater than \({ ICER(C_{14})}\). This implies that \(C_{14}\) is less costly and more effective than \(C_4\). Hence, single implementation of preventive measure is excluded from the list. This shows that combined effort of the four control measures is the most cost effective intervention capable of diminishing the burden of the novel coronavirus optimally in the host population.

Table 2 Comparison between \(C_4\) and \(C_{14}\).
Full size table
Figure 5
figure 5

The dynamics of the model in the case whenever \(\mathcal {R}_0<1\) and \(X_0=(9.814987767,0,0,0)\). For this, we use the value of the parameters: \(\beta _1=0.00005\), \(\sigma =0.09871\), \(\eta =0.01\), \(\beta _2=0.0000000123\), \(\Pi =0.093907997\), \(\alpha =0.0398\), \(d_1=0.00404720925\) and \(d=0.009567816\).

Full size image
Figure 6
figure 6

The dynamics of the model compartmental population in the case whenever \(\mathcal {R}_0>1\) and \(X_{*}=(23.73428588,30.31567421,31.27631426,1261.065047)\). For this, we use the value of the parameters are \(\beta _1=0.0005\), \(\sigma =0.009871\), \(\eta =0.01\), \(\beta _2=0.0000123\), \(\Pi =0.93907997\), \(\alpha =0.398\), \(d_1=0.00404720925\) and \(d=0.009567816\).

Full size image
Figure 7
figure 7

The graphical visualization of the control problem with and without controls.

Full size image
Figure 8
figure 8

The plots visualizes the dynamics of time dependent controls measures.

Full size image

Conclusions

In the current work, we studied the pandemic trend of the 2019-nCoV virus transmission. We developed an epidemiological model with control mechanism strategies to eradicate/control the virus. We studied the fundamental axioms of the model to show the mathematical and biological feasibility. Once the model is well-possed, we find the threshold parameter (\(\mathcal {R}_0\)) and discuss the detailed sensitivity, which clarifies the role of every epidemic parameter involved in the combination of the proposed model. It would be noted that the key parameters that strengthen the spread of the disease are \(\beta _1\), which is approximately 99.99 percent, affecting the threshold quantity’s value. Moreover, the stability analysis is also performed to determine the conditions for stabilizing the disease’s exponential spread. Based on stability analysis, the various perspective of the results is investigated numerically, which states that one of the necessary steps is to optimize the threshold quantity value and keep less than unity. Similarly, if the disease persists, then the long run of the model suggests that it attains its endemic equilibrium, which is very high and not valid for the community. For this, special attention is needed to apply some control strategies to keep this minimum. Once we observed all these situations, we designed a control mechanism based on sensitivity analysis and the dynamics of the problem, which describe that the 2019-nCoV virus infection will be easily eliminated once the control program is applied truly. All the theoretical works are justified with the aid of numerical experiments. This would be quite easy for the readers to study the dynamics of the problem.

In future, we will modify the proposed model by replacing the mass action law with standard incidence to study its significance on the novel corona virus disease dynamics.