1 Introduction

Mathematical models are convenient to recognize the nature of an infection when it arrives a community and to explore under which circumstances it will be continued or wiped out from that community. Presently, COVID-19 is of main alarm to governments, researchers and to all of the mankind. This is due to the high percentage of the infection range and the very high number of deaths that happened. The very first case of infectious disease COVID-19 is reported in Wuhan city of China. Moreover, it spread almost everywhere in the world. On January 2020, the WHO approved its outburst as a public health disaster of international concern. On April 2020, 47, 249 people had died and the total of 936,237 people are tested positive for COVID-19 [1]. COVID-19 ranges due to the interaction of individuals with infected individual when they sneeze or cough. COVID-19 is a respiratory disease with slight to temperate signs like fever, dry cough and tiredness. In severe cases, it causes the difficulty in breathing [1]. Also, few people having mild symptoms of COVID-19 disease may recuperate themselves if they evade to contact with infected cases and maintain good sanitization. COVID-19 is a key health warning to individuals with a previous medical history and also to individuals who are 60 years or above in age (elderly population). This was described by Li et al. [2] who considered the average age of 425 people infected with COVID-19 in city Wuhan, China, and almost half of the patients were above the age of 60. Governments of every country in the world are taking numerous defensive measures to control the spread of COVID-19. Till date, there is no effective vaccine present in the world to fight against the COVID-19 virus. Hence, stopping the spread is the only way to fight against this virus. Defensive actions include sanitizing hands regularly, keeping the distance of 1 meter at least from a person coughing or sneezing, and maintaining social distance. Numerous mathematical models have been established so far to report various challenges in foreseeing the outburst of COVID-19 disease. The author in [3] has used the elementary SIR-model to discover the real dimensions of the epidemic. Peng et al. [4] analyzed the situation of COVID-19 in China by expressing the SEIR dynamical system. Furthermore, he have foretold that the situation will be under control at the start of April. Sun et al. [5] argued numerous parts of COVID-19 condition in China which assistances recognize the casualty rate and spread rate of COVID-19 and control the epidemic transmission. In the situation of COVID-19 prevalent, experience to disease plays a dynamic role in the transmission of the COVID-19. The author in [6] has studied numerous models and determined that with the basic reproduction number being 2 and bearing in mind the 14-day contagious period, if an infected individual stays for 9 or more hours with others, he could infect others. If the contact time is 18 h, the model endorses total shield with more than 70 usefulness to the attendees of the communal gathering. The authors in [7] have discussed the importance of travel restriction to control the COVID-19 disease. In addition, they explained that these travel restriction delayed the progression of disease in Wuhan city of China. A related consequence was also revealed by Kucharski et al. [8]. They exposed that with the journey limitation in Wuhan, the daily reproduction number decayed from 2.35 to 1.05. Furthermore, there is correspondingly a model calculated by Tang et al. [9] and Tang et al. [10] wherein they separated the subpopulation into isolated and unisolated classes to know the spread threat of the epidemic.

Khan et al. [11] have studied the dynamics of COVID-19 with quarantined and isolation by considering the number of real statistical cases reported in China. The authors in [12] have studied the dynamics of a mathematical model by using classical Caputo fractional derivative. Oud et al. [13] have studied a fractional order mathematical model for COVID-19 dynamics with isolation, quarantine, and environmental viral load. For the study of some interesting models related to COVID-19 pandemic, we refer a reader to [14,15,16,17]. From the evaluations so far, we have noticed that isolation plays a huge role in controlling the disease spread. Here, we have considered a mathematically compartmental model for the persons already anguish from hypertension or diabetes who are at a greater chance of getting infected with COVID-19 (see [1]). Moreover, the authors in [1] have considered a model by means of Z-control applied to the isolated class to get the essential exposure state in command to make the system free of chaos.

Table 1 Definitions of parameters and their respective values
Full size table
Fig. 1
figure 1

Flow diagram representing the shifting of any individual from one class to other class

Full size image

Limited number of basic models related to the Z-control contains wide-ranging model by Samanta [18] and prey-predator model by Alzahrani et al. [19]. With three conjointly exclusive classes, namely, susceptible class, tested for hypertension or diabetes (S); exposed or unprotected class (E); and quarantine class (Q). Values of parameters considered in the preparation of our mathematical system are provided in Table 1. In our model, the susceptible class is considered as a population tested for hypertension or suffering from diabetes. Here, the birth rate for new individuals is represented by B,  and all the population in their particular sections undergo death at the uniform rate \(\mu \). Let us represent the saturated rate by \(\frac{\beta _2 SE}{a+S}\), where a is the constant of saturation and \(\beta _2\) is the infection force as shown in Fig. 1. With the help of Fig. 1, we get the following system of nonlinear differential equations [1]:

$$\begin{aligned} \begin{aligned} \frac{\mathrm{d}S}{\mathrm{d}t}&=BS-\beta _{1}SE-\frac{\beta _{2}SE}{a+S}-\frac{\beta _{4}SQ}{b+S}-\mu S^{2}, \\ \frac{\mathrm{d}E}{\mathrm{d}t}&=\beta _{1}SE+\frac{\beta _{2}SE}{a+S}-\frac{\beta _{3}EQ}{d+E}-\mu E,\\ \frac{\mathrm{d}Q}{\mathrm{d}t}&=\frac{c_{1}\beta _{3}QE}{d+E}+\frac{c_{2}\beta _{4}SQ}{b+S}-\mu Q, \end{aligned} \end{aligned}$$
(1)

where, \(S\ge 0, \) \(E\ge 0,\) \(Q\ge 0.\)

In [20], Singh et al. show that the discretization by using Euler’s scheme is not appropriate for every continuous-time dynamical. Furthermore, it violates accuracy of the numerical method and bifurcations occur for larger values of step size used in Euler’s scheme. In order to eliminate this lack, a different discretization technique can be applied as follows: Supposing that the ordinary evolution rates in each of the compartments vary at the fixed pause of time. Then by using the technique of piecewise constant arguments(see [21]) for differential equations, system (1) can be written as:

$$\begin{aligned} \begin{aligned} \frac{1}{S(t)}\frac{\mathrm{d}S(t)}{\mathrm{d}t}&=B-\beta _{1}E[(t)]-\frac{\beta _{2}E[(t)]}{a+S[(t)]}-\frac{\beta _{4}Q[(t)]}{b+S[(t)]}-\mu S[(t)], \\ \frac{1}{E(t)}\frac{\mathrm{d}E(t)}{\mathrm{d}t}&=\beta _{1}S[(t)]+\frac{\beta _{2}S[(t)]}{a+S[(t)]}-\frac{\beta _{3}Q[(t)]}{d+E[(t)]}-\mu ,\\ \frac{1}{Q(t)}\frac{\mathrm{d}Q(t)}{\mathrm{d}t}&=\frac{c_{1}\beta _{3}E[(t)]}{d+E[(t)]} +\frac{c_{2}\beta _{4}S[(t)]}{b+S[(t)]}-\mu , \end{aligned} \end{aligned}$$
(2)

where \(0<t<\infty \) and [t] represents the integer part of t. Furthermore, on an interval \([n, n + 1)\) with \(n\in {\mathbb {W}}\) one can get the next system by integrating system (2) for \(t\in [n, n + 1)\), \(n\in {\mathbb {W}}.\)

$$\begin{aligned} \begin{aligned} S(t)&=S_{n}e^{[{B-\beta _{1}E_{n}-\frac{\beta _{2}E_{n}}{a+S_{n}}-\frac{\beta _{4}Q_{n}}{b+S_{n}}-\mu S_{n}}](t-n)}, \\ E(t)&=E_{n}e^{[{\beta _{1}S_{n}+\frac{\beta _{2}S_{n}}{a+S_{n}}-\frac{\beta _{3}Q_{n}}{d+E_{n}}-\mu }](t-n)},\\ Q(t)&=Q_{n}e^{[{\frac{c_{1}\beta _{3}E_{n}}{d+E_{n}}+\frac{c_{2}\beta _{4}S_{n}}{b+S_{n}}-\mu }](t-n)}. \end{aligned} \end{aligned}$$
(3)

by letting \(t \longrightarrow n + 1\) , we get the following discrete-time mathematical model from system (3):

$$\begin{aligned} \begin{aligned} S_{n+1}&=S_{n}e^{B-\beta _{1}E_{n}-\frac{\beta _{2}E_{n}}{a+S_{n}}-\frac{\beta _{4}Q_{n}}{b+S_{n}}-\mu S_{n}}, \\ E_{n+1}&=E_{n}e^{\beta _{1}S_{n}+\frac{\beta _{2}S_{n}}{a+S_{n}}-\frac{\beta _{3}Q_{n}}{d+E_{n}}-\mu },\\ Q_{n+1}&=Q_{n}e^{\frac{c_{1}\beta _{3}E_{n}}{d+E_{n}}+\frac{c_{2}\beta _{4}S_{n}}{b+S_{n}}-\mu }. \end{aligned} \end{aligned}$$
(4)

The aim of our study in this article is to explain the boundedness of every solution of the system (4). To discuss the local stability of system (4) about each of its equilibrium point. Moreover, the existence of Neimark–Sacker bifurcation and chaos for one and only equilibrium point of system (4) is scrutinized. In concern to control the Neimark–Sacker bifurcation and, a modified hybrid control technique is implemented [22] in Sect. 6. In final section, some numerical examples are provided.

2 Boundedness of system (4)

In this section, the boundedness of every positive solution \((S_{n},E_{n},Q_{n})\) is proved. For this purpose, the next lemma is presented.

Lemma 2.1

[23] Assume that \(\lambda _{n}\) fulfills \(\lambda _{n+1}\le \lambda _{n} \mathrm{exp}(a(1-b\lambda _{n}))\) for every \(n\in [n_{1}, \infty )\) with \(\lambda _{0}>0,\) where \(a,b>0.\) Then,

$$\begin{aligned} \lim _{n \longrightarrow \infty }\mathrm{sup}{\lambda _{n}}\le \frac{1}{ab}\mathrm{exp}(a-1). \end{aligned}$$

Lemma 2.2

Every solution \((S_{n},E_{n},Q_{n})\) of system (4) is bounded uniformly if for a finite \(N\in \mathfrak {R},\) we have

$$\begin{aligned} \lim _{n \longrightarrow \infty }\mathrm{sup}{E_{n}}\le N. \end{aligned}$$

Proof

Assume that \(S_{0}>0,\) \(E_{0}>0\) and \(Q_{0}>0\) then each solution \((S_{n},E_{n},Q_{n})\) of system (4) satisfies \(S_{n}>0,\) \(E_{n}>0\) and \(Q_{n}>0\) for each \(n\ge 0.\) Firstly, by taking into account the positivity of solutions of (4) and from first equation of system (4) it can be seen that

$$\begin{aligned} S_{n+1}= & {} S_{n}\mathrm{exp}\left( {B-\beta _{1}E_{n}-\frac{\beta _{2}E_{n}}{a+S_{n}} -\frac{\beta _{4}Q_{n}}{b+S_{n}}-\mu S_{n}}\right) \\\le & {} S_{n}\mathrm{exp}\left( {B-\beta _{1}E_{n}-\frac{\beta _{2}E_{n}}{a+S_{n}}-\mu S_{n}}\right) \\\le & {} S_{n}\mathrm{exp}({B-\mu S_{n}})\\= & {} S_{n}\mathrm{exp}\left( {B(1-\frac{\mu }{B}S_{n})}\right) . \end{aligned}$$

Next, by applying Lemma 2.3 we get

$$\begin{aligned} \lim _{n \longrightarrow \infty }\mathrm{sup}{ S_{n}}\le \frac{1}{\mu }\mathrm{exp}(B-1)=\digamma (\mathrm{say}). \end{aligned}$$

Now, from third equation of system (4) it can be seen that

$$\begin{aligned} Q_{n+1}= & {} Q_{n}\mathrm{exp}\left( {\frac{c_{1}\beta _{3}E_{n}}{d+E_{n}} +\frac{c_{2}\beta _{4}S_{n}}{b+S_{n}}-\mu }\right) \\\le & {} Q_{n}\mathrm{exp}\left( {\frac{c_{1}\beta _{3}N}{d+N}+\frac{c_{2}\beta _{4} ({B-\beta _{1}N-\frac{\beta _{2}N}{a+\digamma }-\frac{\beta _{4}Q_{n}}{b+\digamma }})}{\mu b}}\right) \\\le & {} Q_{n}\mathrm{exp}\left( {\frac{c_{1}\beta _{3}N}{d+N}+\frac{c_{2}\beta _{4}B}{\mu b} -\frac{c_{2}\beta _{4}^{2}B}{\mu b(b+\digamma )}Q_{n}}\right) \\= & {} Q_{n}\mathrm{exp}\left( \frac{(\mu N b c_{1} \beta _{3}+c_{2}\beta _{4} B(d+N))}{\mu b(d+N)}\left( 1-\frac{c_{2}\beta _{4}^{2} B(d+N)Q_{n}}{(b+\digamma )(\mu N b c_{1} \beta _{3}+c_{2}\beta _{4} B(d+N))}\right) \right) . \end{aligned}$$

Hence, by applying Lemma 2.3 we get

$$\begin{aligned} \lim _{n \longrightarrow \infty }\mathrm{sup}{Q_{n}}\le \frac{1}{\frac{c_{2}\beta _{4}^{2}B}{\mu b(b+\digamma )}}\mathrm{exp}\left( \frac{(\mu N b c_{1} \beta _{3}+c_{2}\beta _{4} B(d+N))}{\mu b(d+N)}-1\right) . \end{aligned}$$

\(\square \)

By using the method of mathematical induction, one can prove the next result.

Lemma 2.3

Assume that \(0<S_{0}<\frac{1}{\mu }\mathrm{exp}(B-1),\) \(0<E_{0}<N \) and

$$\begin{aligned} 0<Q_{0}<\frac{1}{\frac{c_{2}\beta _{4}^{2}B}{\mu b(b+\digamma )}}\mathrm{exp}\left( \frac{(\mu N b c_{1} \beta _{3}+c_{2}\beta _{4} B(d+N))}{\mu b(d+N)}-1\right) , \end{aligned}$$

then the set \( [0,\frac{1}{\mu }\mathrm{exp}(B-1)]\times [0,N]\times [0,\frac{1}{\frac{c_{2}\beta _{4}^{2}B}{\mu b(b+\digamma )}}\mathrm{exp}\left( \frac{(\mu N b c_{1} \beta _{3}+c_{2}\beta _{4} B(d+N))}{\mu b(d+N)}-1\right) ]\) remains invariant for every solution \((S_{n},E_{n},Q_{n})\) of the system (4).

3 Existence of equilibrium points and local stability of system (4)

In this section, we contemplate probable equilibrium points (S, E, Q) of system (4), which can be acquired by considering the next system:

$$\begin{aligned} \begin{aligned} S^*&=S^*e^{B-\beta _{1}E^*-\frac{\beta _{2}E^*}{a+S^*}-\frac{\beta _{4}Q^*}{b+S^*}-\mu S^*}, \\ E^*&=E^*e^{\beta _{1}S^*+\frac{\beta _{2}S^*}{a+S^*}-\frac{\beta _{3}Q^*}{d+E^*}-\mu },\\ Q^*&=Q^*e^{\frac{c_{1}\beta _{3}E^*}{d+E^*}+\frac{c_{2}\beta _{4}S^*}{b+S^*}-\mu }. \end{aligned} \end{aligned}$$
(5)

On solving system (5), one can obtain six equilibrium points (0, 0, 0), \((\frac{B}{\mu },0,0),\) \((0,\frac{d \mu }{ c_1 \beta _3-\mu },\)\(\frac{d \mu c_1}{\mu -c_1 \beta _3}),\) \((\frac{\mu b}{(c_{2}\beta _{4}-\mu )}, 0,\frac{bc_{2}(B c_{2}\beta _{4}-b\mu ^{2}-B \mu )}{(c_{2}\beta _{4}-\mu )^2}),\) \(({\bar{S}}, {\bar{E}}, 0),\) and the unique positive equilibrium point \((S^*,E^*,Q^*)\).

Remark 3.1

The equilibrium point \((0,\frac{d \mu }{ c_1 \beta _3-\mu },\frac{d \mu c_1}{\mu -c_1 \beta _3})\) ceased to exist as one of the components from \(\frac{d \mu }{ c_1 \beta _3-\mu }\) or \(\frac{d \mu c_1}{\mu -c_1 \beta _3}\) remains negative for every \(c_1, \beta _3, \mu >0.\)

Let

$$\begin{aligned} F_J=\begin{array}{ccc} \left[ \begin{array}{ccc} j_{11}&{}j_{12}&{}j_{13}\\ j_{21}&{} j_{22}&{}j_{23}\\ j_{31}&{}j_{32}&{}j_{33} \end{array} \right] \end{array} \end{aligned}$$

be the variational matrix evaluated at \((S^*,E^*,Q^*)\). Then, characteristic polynomial \({\mathbb {H}}(\omega )\) of matrix \(F_J\) is:

$$\begin{aligned} {\mathbb {H}}(\omega )=\omega ^3-A_1\omega ^2+A_2\omega -A_3, \end{aligned}$$
(6)

where

$$\begin{aligned} A_1= & {} (j_{11}+j_{22}+j_{33}), \\ A_2= & {} J_{11}+J_{22}+J_{33}, \end{aligned}$$

and

$$\begin{aligned} A_3=det(F_J). \end{aligned}$$

Where \(J_{11}\), \(J_{22}\) and \(J_{33}\) are minor determinants of Jacobian matrix \(F_J.\) Firstly, we explore the stability analysis of the trivial fixed point (0, 0, 0). The Jacobian matrix \(F_{J_{1}}\) about equilibrium point (0, 0, 0), is given by;

$$\begin{aligned} F_{J_{1}}=\left( \begin{array}{ccc} e^B &{} 0 &{} 0 \\ 0 &{} e^{-\mu } &{} 0 \\ 0 &{} 0 &{} e^{-\mu } \end{array} \right) . \end{aligned}$$

In addition, \(F_{J_{1}}\) has three eigenvalues, namely \(\tau _1=e^B,\) \(\tau _2=e^{-\mu }\) and \(\tau _3=e^{-\mu }\), such that \(|\tau _1|>1\) and \(|\tau _2|=|\tau _3|<1\) remains true for all parametric values. Hence, we conclude the following proposition about the local stability of (4) about (0, 0, 0).

Proposition 3.1

Let (0, 0, 0) be a equilibrium point of (4), then (0, 0, 0) remains unstable for every \(B,\mu >0\).

Fig. 2
figure 2

Stable region for \((\frac{B}{\mu },0,0)\) for \(a=2, b=10, \beta _1=0.99, d=0.4, \beta _2=0.01, \beta _3=0.80, \beta _4=0.60, c_1=1, c_2=2\) , \(0<B<1 ,0<\mu <2\) and \(S_0=0.995, E_0=0.38855, Q_0=0.3455\)

Full size image
Fig. 3
figure 3

Stable and unstable regions for \(\left( \frac{b \mu }{c_2 \beta _4-\mu },0,\frac{b c_2 \left( B c_2 \beta _4-\mu (B+b \mu )\right) }{\left( \mu -c_2 \beta _4\right) {}^2}\right) \) for \(a=2, b=10, \beta _1=0.99, d=0.4, \beta _2=0.01, \beta _3=0.80, \beta _4=0.60, c_1=1, c_2=2\) , \(0<B<5 ,0<\mu <1\) and \(S_0=0.995, E_0=0.38855, Q_0=0.3455\)

Full size image

From Proposition 3.1, one can observe that the extinction equilibrium point is mathematically unstable and biologically it is not possible to hold whenever any one of the three classes from (4) exists. Next, our goal is to explore the local stability of system (4) about \((\frac{B}{\mu },0,0)\). The matrix of variation \(F_{J_{2}}\) evaluated about \((\frac{B}{\mu },0,0)\) can be calculated as:

$$\begin{aligned} F_{J_{2}}=\left( \begin{array}{ccc} 1-B &{} -\frac{B \beta _1}{\mu }-\frac{B \beta _2}{B+a \mu } &{} -\frac{B \beta _4}{B+b \mu } \\ 0 &{} e^{-\mu +\frac{B \beta _1}{\mu }+\frac{B \beta _2}{B+a \mu }} &{} 0 \\ 0 &{} 0 &{} e^{-\mu -\frac{B c_2 \beta _4}{B+b \mu }} \end{array} \right) . \end{aligned}$$

Let us assume that \({\mathbb {H}}(\tau )\) is characteristic polynomial of Jacobian matrix \(F_{J_{2}}\). Then, characteristic roots of \({\mathbb {H}}(\tau )=0\) are given by \(\tau _1=1-B,\) \(\tau _2=\frac{1}{e^{\mu -\frac{B \beta _1}{\mu }-\frac{B \beta _2}{B+a \mu }}}\) and \(\tau _3=\frac{1}{e^{\mu +\frac{B c_2 \beta _4}{B+b \mu }}}.\) Hence, we have the following proposition related to the local stability of (4) about \((\frac{B}{\mu },0,0)\) (Fig. 2, 3).

Proposition 3.2

Let \((\frac{B}{\mu },0,0)\) be a equilibrium point of (4) then;

  1. (i)

    \((\frac{B}{\mu },0,0)\) is a stable equilibrium point if \(|\tau _1|,|\tau _2|,|\tau _3|<1\) if and only if

    $$\begin{aligned} 0<B<2 \ and \ \mu ^2(B+a\mu )>B(\beta _{1}(B+a\mu )+\mu \beta _{2}). \end{aligned}$$
  2. (ii)

    \((\frac{B}{\mu },0,0)\) is a unstable equilibrium point if

    $$\begin{aligned} B>2. \end{aligned}$$

Consider the equilibrium point \(\left( \frac{b \mu }{c_2 \beta _4-\mu },0,\frac{b c_2 \left( B c_2 \beta _4-\mu (B+b \mu )\right) }{\left( \mu -c_2 \beta _4\right) {}^2}\right) \) of system (4) and suppose \(F_{J_{3}}\) be the matrix of variation for system (4) about \(\left( \frac{b \mu }{c_2 \beta _4-\mu },0,\frac{b c_2 \left( B c_2 \beta _4-\mu (B+b \mu )\right) }{\left( \mu -c_2 \beta _4\right) {}^2}\right) \) . Then, \(F_{J_{3}}\) has the following mathematical form:

$$\begin{aligned} F_{J_{3}}=\left( \begin{array}{ccc} 1+\frac{\mu (B+b \mu )}{c_2 \beta _4}+\frac{2 b \mu ^2}{\mu -c_2 \beta _4} &{} b \mu \left( \frac{\beta _1}{\mu -c_2 \beta _4}+\frac{\beta _2}{a \mu -b \mu -a c_2 \beta _4}\right) &{} -\frac{\mu }{c_2} \\ 0 &{} e^{-\mu +\frac{b \mu \beta _1}{-\mu +c_2 \beta _4}+\frac{b \mu \beta _2}{(-a+b) \mu +a c_2 \beta _4}+\frac{b c_2 \beta _3 \left( \mu (B+b \mu )-B c_2 \beta _4\right) }{d \left( \mu -c_2 \beta _4\right) {}^2}} &{} 0 \\ \frac{e^{-2 \mu } \left( \mu (B+b \mu )-B c_2 \beta _4\right) }{\beta _4} &{} -\frac{b e^{-2 \mu } c_1 c_2 \beta _3 \left( \mu (B+b \mu )-B c_2 \beta _4\right) }{d \left( \mu -c_2 \beta _4\right) {}^2} &{} e^{-2 \mu } \end{array} \right) . \end{aligned}$$

Assume that \({\mathbb {H}}(\psi )\) is characteristic polynomial of Jacobian matrix \(F_{J_{3}}\). Then, characteristic roots of \({\mathbb {H}}(\psi )=0\) are given by

$$\begin{aligned} \psi _1= & {} \frac{1}{e^{\mu -\frac{b \mu \beta _1}{c_2 \beta _4-\mu }-\frac{b \mu \beta _2}{(b-a) \mu +a c_2 \beta _4}-\frac{b c_2 \beta _3 \left( \mu (B+b \mu )-B c_2 \beta _4\right) }{d \left( \mu -c_2 \beta _4\right) {}^2}}}, \\ \psi _2,\psi _3= & {} \frac{\left( a_{11}\pm \sqrt{a_{12}+a_{13}}\right) }{e^{2 \mu }2 c_2 \beta _4 \left( \mu -c_2 \beta _4\right) } , \end{aligned}$$

where

$$\begin{aligned} \left\{ \begin{array}{ll} a_{11}\!=\!e^{2 \mu } \mu ^2 (B\!+\!b \mu )\!+\!\mu \left( 1\!+\!e^{2 \mu } (1\!-\!B\!+\!b \mu )\right) c_2 \beta _4\!-\!\left( 1\!+\!e^{2 \mu }\right) c_2^2 \beta _4^2, \\ a_{12}\!=\!\left( e^{2 \mu } \mu ^2 (B\!+\!b \mu )\!+\!c_2 \beta _4 \left( \mu \!+\!e^{2 \mu } \mu (1\!-\!B\!+\!b \mu )\!-\!\left( 1\!+\!e^{2 \mu }\right) c_2 \beta _4\right) \right) {}^2, \\ \ \ \ and\\ a_{13}\!=\!4 e^{2 \mu } c_2 \beta _4 \left( c_2 \beta _4\!-\!\mu \right) \left( \mu ^2 (1\!+\!\mu ) (B+b \mu )\!-\! c_2 \beta _4 \left( \mu (B\!-\!1\!+\!2 B \mu \!+\!b (\mu \!-\!1) \mu )\!+\!(B \mu \!-\!1) c_2 \beta _4\right) \right) . \end{array} \right. \end{aligned}$$
(7)

Proposition 3.3

Let \(\left( \frac{b \mu }{c_2 \beta _4-\mu },0,\frac{b c_2 \left( B c_2 \beta _4-\mu (B+b \mu )\right) }{\left( \mu -c_2 \beta _4\right) {}^2}\right) \) be a equilibrium point of (4) and \(\psi _1,\) \(\psi _2\) and \(\psi _3\) are characteristic roots for \(F_{J_{3}}\), then \(\left( \frac{b \mu }{c_2 \beta _4-\mu },0,\frac{b c_2 \left( B c_2 \beta _4-\mu (B+b \mu )\right) }{\left( \mu -c_2 \beta _4\right) {}^2}\right) \) remains stable if

$$\begin{aligned} \mu >\frac{b \mu \beta _1}{c_2 \beta _4-\mu }+\frac{b \mu \beta _2}{(b-a) \mu +a c_2 \beta _4}+\frac{b c_2 \beta _3 \left( \mu (B+b \mu )-B c_2 \beta _4\right) }{d \left( \mu -c_2 \beta _4\right) {}^2} \end{aligned}$$

and

$$\begin{aligned} |\left( a_{11}\pm \sqrt{a_{12}+a_{13}}\right) |< e^{2 \mu }2 c_2 \beta _4 \left( \mu -c_2 \beta _4\right) \end{aligned}$$

with \(\mu >c_2 \beta _4\) and \(a_{11}, a_{12}, a_{13}\) are given in (7).

It is prominent that the equilibrium points of mathematical system (4), that is, the solution of system (5) cannot be unique, however, for biological aims; we are concerned to the positive results of (5). Hence, we do not care about exactly how many results are there of system (5). From the system (5), we get

$$\begin{aligned} E^*=-d\left( 1+\frac{c_{1}\beta _{3}(b+S^*)}{\mu (b+S^*)+c_{2}\beta _{4}+S^*-c_{1} \beta _{3}(b+S^*)}\right) \end{aligned}$$

and

$$\begin{aligned} Q^*=\frac{-(d+E^*)}{\beta _{3}}\left( \mu -\left( \beta _{1}+\frac{\beta _{2}}{(a+S^*)} \right) S^*\right) , \end{aligned}$$

where \(S^*\) is one of the roots of cubic polynomial:

$$\begin{aligned} g(t)=A_{11}t^{3}+B_{11}t^{2}+C_{11}t+D_{11}, \end{aligned}$$
(8)

where

$$\begin{aligned} \left\{ \begin{array}{ll} A_{11}=\mu (c_{1}\beta _{3}+c_{2}\beta _{4}-\mu ), \\ B_{11}=(B-(a+b)\mu )(\mu -c_{1}\beta _{3})-(B-a\mu )c_{2}\beta _{4}+d\beta _{1} (\mu +(c_{1}-c_{2})\beta _{4}), \\ C_{11}=((a+b)B-ab\mu )(\mu -c_{1}\beta _{3})-(d\mu c_{1}+aB c_{2})\beta _{4} +(\mu +(c_{1}-c_{2})\beta _{4})d\beta _{2}\\ +((a+b)\mu +a(c_{1}-c_{2})\beta _{4})d\beta _{1},\\ D_{11}=(ab(B+d\beta _{1})+db\beta _{3})\mu -ac_{1}(bB\beta _{3}+d\mu \beta _{4}). \end{array} \right. \end{aligned}$$
(9)

We are looking for the unique positive equilibrium point of system (5). For this, we have the following Descartes rule of signs (see [24]).

Lemma 3.1

[24] Let \(f_{1}(x)=b_{n}x^{n}+b_{n-1}x^{n-1}+b_{n-2}x^{n-2}+.....+b_{1}x+b_{0}\) be a polynomial function with real coefficients. Then, the number of positive real roots of \(f_{1}\) is either the same as the number of sign changes for \(f_{1}(x)\) or less than by a positive even integer. Moreover, if f(x) has only one variation in sign, then \(f_{1}\) has exactly one positive real root.

Using Lemma 3.1, we have the following result for existence of unique positive real root of polynomial g(t) given in (8).

Lemma 3.2

Polynomial g(t) in (8) has unique positive real root if one of the following conditions hold:

\(i. \ \ A_{1}<0, B_{11}<0, C_{11}<0, D_{11}>0.\)

\(ii. \ \ A_{11}<0, B_{11}>0, C_{11}>0, D_{11}>0.\)

Lemma 3.3

Let \(S^*\) is one of the roots of cubic polynomial (8). Then, under the conditions of lemma 3.2, system (1) has one and only positive fixed point if the following conditions are satisfied:

$$\begin{aligned} c_{1}\beta _{3}+c_{2}\beta _{4}<\mu <(c_{2}-c_{1})\beta _{4} \end{aligned}$$

and

$$\begin{aligned} B<\mu . \end{aligned}$$

Proof

Positive equilibrium point of (4) fulfills algebraic equations same to its continuous counterpart (1). In case of system (1), the evidence is specified in [1].\(\square \)

4 Local stability analysis of system (4) about its endemic equilibrium

We consider the Jacobian matrix \(F_{J_{4}}\) of system (4) about unique positive equilibrium point \((S^* ,E^* ,Q^* ).\)

$$\begin{aligned} F_{J_{4}}=\left( \begin{array}{ccc} 1+S^* \left( \frac{\beta _2 E^*}{\left( a+S^*\right) ^2-\mu }+\frac{\beta _4 Q^*}{\left( b+S^*\right) ^2}\right) &{} -S^* \left( \beta _1+\frac{\beta _2}{a+S^*}\right) &{} -\frac{\beta _4 S^*}{b+S^*} \\ \frac{E^* \left( a \beta _2+\beta _1 \left( a+S^*\right) ^2\right) }{\left( a+S^*\right) ^2} &{} \frac{\left( d+E^*\right) ^2+\beta _3 E^* Q^*}{\left( d+E^*\right) ^2} &{} -\frac{\beta _3 E^*}{d+E^*} \\ -\frac{b c_2 \beta _4 Q^*}{\left( b+S^*\right) ^2} &{} \frac{d c_1 \beta _3 Q^*}{\left( d+E^*\right) ^2} &{} 1 \end{array} \right) . \end{aligned}$$
(10)

Moreover, assume that (14) be the characteristic equation of (10). Then, in direction to study the stability analysis of one and only positive equilibrium point, we have the next theorem, which provides us freely confirmable necessary and sufficient situations for all the roots of real polynomial of degree 3 to have magnitude smaller than one (see, Theorem 5 from [25]).

Theorem 4.1

Assume the next third-degree polynomial equation

$$\begin{aligned} {\mathbb {H}}(\omega )=\omega ^3-A_1\omega ^2+A_2\omega -A_3, \end{aligned}$$
(11)

where \(A_1, A_2, A_3 \in \mathfrak {R}.\) Then, necessary and sufficient conditions that all roots of (11) lie inside an open disk are given as follows:

$$\begin{aligned} \left\{ \begin{array}{ll} |A_1+A_3|<1+A_2, \\ |A_1-3A_3|<3-A_2, \\ {A_3}^{2}+A_2-A_1 A_3<1. \end{array} \right. \end{aligned}$$
(12)

Theorem 4.2

The unique positive equilibrium point \((S^* ,E^* ,Q^* )\) of system (4) is locally asymptotically stable if the following conditions are satisfied:

$$\begin{aligned} \left\{ \begin{array}{ll} |A_1+A_3|<1+A_2, \\ |A_1-3A_3|<3-A_2, \\ {A_3}^{2}+A_2-A_1 A_3<1. \end{array} \right. \end{aligned}$$

where

$$\begin{aligned} \left\{ \begin{array}{ll} A_1=3+\frac{\beta _3 E^* Q^*}{\tau ^2}+\left( -\mu +\frac{\beta _2 E^*}{\rho ^2}+\frac{\beta _4 Q^*}{\sigma ^2}\right) S^*, \\ A_2=3+\frac{d c_1 \beta _3^2 E^* Q^*}{\tau ^3}+S^* \left( \frac{\left( \rho ^3 \beta _1^2+\rho \left( 2+(a+\rho ) \beta _1\right) \beta _2+a \beta _2^2\right) E^*}{\rho ^3}-2 \mu +\beta _4 Q^* \left( \frac{\sigma -b c_2 \beta _4}{\sigma ^3}+\frac{1}{\left( \sigma ^*\right) ^2}\right) \right) \\ +\frac{\beta _3 E^* Q^* \left( \rho ^2 \beta _4 Q^* S^*+\left( 2 \rho ^2+\left( -\mu \rho ^2+\beta _2 E^*\right) S^*\right) \left( \sigma ^*\right) ^2\right) }{\rho ^2 \tau ^2 \left( \sigma ^*\right) ^2}, \\ A_3=\frac{1}{\rho ^3 \sigma ^3 \tau ^3}\left( \rho ^3 \sigma ^3 \left( \tau ^3+\beta _3 \left( \tau +d c_1 \beta _3\right) E^* Q^*\right) +\left( \sigma ^3 \tau ^3 \left( \rho ^3 \beta _1^2+\rho \left( 1+(a+\rho ) \beta _1\right) \beta _2+a \beta _2^2\right) E^*-\mu \rho ^3 \right) \right) \\ -\frac{\rho }{\rho ^3 \sigma ^3 \tau ^3} \left( b \rho ^2 \tau ^3 c_2 \beta _4^2+\sigma ^3 \beta _3 \left( \tau +d c_1 \beta _3\right) E^* \left( \mu \rho ^2-\beta _2 E^*\right) \right) \\ +\frac{1}{\rho ^3 \sigma ^3 \tau ^3}\left( \sigma \tau \beta _4 \left( -\rho ^2 \tau ^2+\left( b \rho \tau c_2 \left( \rho \beta _1+\beta _2\right) +d \sigma c_1 \left( \rho ^2 \beta _1+a \beta _2\right) \right) \beta _3 E^*\right) Q^*\right) \\ +\frac{\rho ^3 \beta _3 \beta _4 }{\rho ^3 \sigma ^3 \tau ^3}\left( d \sigma c_1 \beta _3+\tau \left( \sigma -b c_2 \beta _4\right) \right) E^* \left( Q^*\right) ^2 S^*,\\ a+S^*=\rho ,\\ b+S^*=\sigma ,\\ \ \ and\\ d+E^*=\tau .\\ \end{array} \right. \end{aligned}$$
(13)

Proof

The Jacobian matrix for system (4) evaluated at unique positive fixed point \((S^* ,E^* ,Q^* )\) is given by

$$\begin{aligned} F_{J_{4}}=\left( \begin{array}{ccc} 1+S^* \left( \frac{\beta _2 E^*}{\left( a+S^*\right) ^2-\mu }+\frac{\beta _4 Q^*}{\left( b+S^*\right) ^2}\right) &{} -S^* \left( \beta _1+\frac{\beta _2}{a+S^*}\right) &{} -\frac{\beta _4 S^*}{b+S^*} \\ \frac{E^* \left( a \beta _2+\beta _1 \left( a+S^*\right) ^2\right) }{\left( a+S^*\right) ^2} &{} \frac{\left( d+E^*\right) ^2+\beta _3 E^* Q^*}{\left( d+E^*\right) ^2} &{} -\frac{\beta _3 E^*}{d+E^*} \\ -\frac{b c_2 \beta _4 Q^*}{\left( b+S^*\right) ^2} &{} \frac{d c_1 \beta _3 Q^*}{\left( d+E^*\right) ^2} &{} 1 \end{array} \right) . \end{aligned}$$

Then, characteristic polynomial \({\mathbb {H}}(\omega )\) from matrix \(F_{J_{4}}\) is given by:

$$\begin{aligned} {\mathbb {H}}(\omega )=\omega ^3-A_1\omega ^2+A_2\omega -A_3, \end{aligned}$$
(14)

where \(A_1, A_2\) and \( A_3\) are given in (13). Now, by using Theorem 4.1, the one and only positive equilibrium point \((S^* ,E^* ,Q^* )\) is locally asymptotically stable if the following inequalities are satisfied:

$$\begin{aligned} \left\{ \begin{array}{ll} |A_1+A_3|<1+A_2, \\ |A_1-3A_3|<3-A_2, \\ {A_3}^{2}+A_2-A_1 A_3<1. \end{array} \right. \end{aligned}$$

\(\square \)

5 Period-doubling bifurcation

In this part of article, we examine that quarantined free fixed point of system (4) experience period-doubling bifurcation. For this bifurcation concept, center manifold theorem is applied after the application of normal forms to show the existence and direction of such kind of bifurcation. Freshly, period-doubling bifurcation associated with discrete-time models has been explored by many authors [26,27,28,29,30]

Under the suppositions that \(B>\mu \) and \(Q_{n}\rightarrow 0, \forall n, \) we study the following set:

$$\begin{aligned} \Phi _{1}=\{a,B,d,\beta _1,\beta _2 \in \mathfrak {R}^{+}:\mu =f({\overline{S}}, {\overline{E}}, a, B, d, \beta _1, \beta _2)\}, \end{aligned}$$

Then, the quarantined free fixed point \(({\overline{S}},{\overline{E}},0)\) of system (4) experiences period-doubling bifurcation such that \(\mu \) is taken as bifurcation parameter, and it varies in a slight neighborhood of \({\hat{\mu }},\) which is given as

$$\begin{aligned} {\hat{\mu }}=f({\overline{S}}, {\overline{E}}, a, B, d, \beta _1, \beta _2). \end{aligned}$$

In addition, system (4) is characterized evenly with the next two-dimensional map:

$$\begin{aligned} \left( \begin{array}{c} S \\ E \\ \end{array} \right) \rightarrow \left( \begin{array}{c} Se^{B-\beta _{1}E-\frac{\beta _{2}E}{a+S}-\mu S} \\ Ee^{\beta _{1}S+\frac{\beta _{2}S}{a+S}-\mu }\\ \end{array} \right) . \end{aligned}$$
(15)

To discuss and analyze the period-doubling bifurcation for steady state \(({\overline{S}}, {\overline{E}}, 0)\) of (15), we suppose that \(a, B, d, \beta _1, \beta _2, {\hat{\mu }} \in \Phi _{1}\). Then, it follows that

$$\begin{aligned} \left( \begin{array}{c} S \\ E \\ \end{array} \right) \rightarrow \left( \begin{array}{c} Se^{B-\beta _{1}E-\frac{\beta _{2}E}{a+S}-{\hat{\mu }} S} \\ Ee^{\beta _{1}S+\frac{\beta _{2}S}{a+S}-{\hat{\mu }}}\\ \end{array} \right) . \end{aligned}$$
(16)

Taking \({\bar{\mu }}\) as small parameter for bifurcation, then the perturbation of mapping (15) can be described by the next map:

$$\begin{aligned} \left( \begin{array}{c} S \\ E \\ \end{array} \right) \rightarrow \left( \begin{array}{c} Se^{B-\beta _{1}E-\frac{\beta _{2}E}{a+S}-({\hat{\mu }}+{\bar{\mu }}) S} \\ Ee^{\beta _{1}S+\frac{\beta _{2}S}{a+S}-({\hat{\mu }}+{\bar{\mu }})}\\ \end{array} \right) \end{aligned}$$
(17)

where \(|{\bar{\mu }}|\ll 1\), is a small parameter for perturbation.

Taking \(x=S-{\overline{S}}\) and \(y=E-{\overline{E}}\). Then, from (17)we obtained the next map whose fixed point is at (0, 0) ;

$$\begin{aligned} \left( \begin{array}{c} x \\ y \\ \end{array} \right) \rightarrow \left( \begin{array}{cc} \theta _{11} &{} \theta _{12} \\ \phi _{21}&{} \phi _{22}\\ \end{array} \right) \left( \begin{array}{c} x \\ y \\ \end{array} \right) +\left( \begin{array}{c} f_1(x,y, {\bar{\mu }}) \\ f_2(x,y,{\bar{\mu }})\\ \end{array} \right) , \end{aligned}$$
(18)

where

$$\begin{aligned} f_1(x,y,{\bar{\mu }})= & {} \theta _{13}x^2+\theta _{14}xy+\theta _{15}x{\bar{\mu }}+ \theta _{16}y^2+\theta _{17}{\bar{\mu }}y+\theta _{18}{\bar{\mu }}^2+\theta _{19}x^3+\theta _{20}x^2y\\&+\theta _{21}{\bar{\mu }}x^2+\theta _{22}xy^2+\theta _{23}xy{\bar{\mu }}+\theta _{24}x{\bar{\mu }}^2 +\theta _{25}y{\bar{\mu }}^2+\theta _{26}{\bar{\mu }}^3+O\left( (|x|+|y|+|{\bar{\mu }}|)^4\right) ,\\ f_2(x,y, {\bar{\eta }})= & {} \phi _{13}x^2+\phi _{14}xy+\phi _{15}y{\bar{\mu }}+ \phi _{16}{\bar{\mu }}^2+\phi _{17}x^3+\phi _{18}x^2y+\phi _{19}{\bar{\mu }}x^2+\phi _{20}{\bar{\mu }}xy\\&+\phi _{21}{\bar{\mu }}^2x+\phi _{22}y{\bar{\mu }}^2+\phi _{23}{\bar{\mu }}^3+O\left( (|x|+|y|+|{\bar{\mu }}|)^4\right) ,\\ \end{aligned}$$

with

$$\begin{aligned} \left\{ \begin{array}{ll} \theta _{11}={\overline{S}} \left( \frac{{\overline{E}} \beta _2}{(a+{\overline{S}})^2}-{\bar{\mu }} \right) ,\ \theta _{12}=-{\overline{S}} \left( \beta _1+\frac{\beta _2}{a+{\overline{S}}}\right) ,\ \theta _{13}=\frac{1}{2} \left( {\bar{\mu }} ({\overline{S}} {\bar{\mu }}-2)+\frac{{\overline{E}} \beta _2 \left( 2 (a+{\overline{S}}) (a-{\overline{S}} (a+{\overline{S}}) {\bar{\mu }})+{\overline{E}} {\overline{S}} \beta _2\right) }{(a+{\overline{S}})^4}\right) ,\\ \theta _{14}= \left( -\beta _1-\frac{\beta _2}{a+{\overline{S}}}\right) +\frac{{\overline{S}} \beta _2}{(a+{\overline{S}})^2}+{\overline{S}} \left( \frac{\beta _2{\overline{E}}}{(a+{\overline{S}})^2}-{\bar{\mu }}\right) \left( -\beta _1-\frac{\beta _2}{a+{\overline{S}}}\right) ,\ \theta _{15}={\overline{S}} \left( {\overline{S}} {\bar{\mu }} -2-\frac{{\overline{E}} {\overline{S}} \beta _2}{(a+{\overline{S}})^2}\right) ,\\ \theta _{16}=\frac{{\overline{S}} \left( (a+{\overline{S}}) \beta _1+\beta _2\right) {}^2}{2 (a+{\overline{S}})^2},\ \theta _{17}=\frac{{\overline{S}}^2 \left( (a+{\overline{S}}) \beta _1+\beta _2\right) }{a+{\overline{S}}}, \ \theta _{18}=\frac{{\overline{S}}^3}{2},\\ \theta _{19}=\frac{1}{6} \left( {\bar{\mu }} ^2 (3-{\overline{S}} {\bar{\mu }})+\frac{{\overline{E}} \beta _2 \left( 3 (a+{\overline{S}})^2 \left( -2 a-2 a (a+{\overline{S}}) {\bar{\mu }} +{\overline{S}} (a+{\overline{S}})^2 {\bar{\mu }}^2\right) +{\overline{E}} \beta _2 \left( -3 (a+{\overline{S}}) (-a+{\overline{S}}+{\overline{S}} (a+{\overline{S}}) {\bar{\mu }})+{\overline{E}} {\overline{S}} \beta _2\right) \right) }{(a+{\overline{S}})^6}\right) ,\\ \theta _{20}=\frac{\beta _2}{(a+{\overline{S}})^2}+\left( \frac{\beta _2{\overline{E}}}{(a+{\overline{S}})^2}-{\bar{\mu }}\right) \left( -\beta _1-\frac{\beta _2}{a+{\overline{S}}}\right) -\frac{{\overline{S}}\beta _2}{(a+{\overline{S}})^3}-\frac{S\beta _2E \left( -\beta _1-\frac{\beta _2}{a+S}\right) }{(a+{\overline{S}})^3}+\frac{S \left( \frac{\beta _2{\overline{E}}}{(a+{\overline{S}})^2}-{\bar{\mu }} \right) \beta _2}{(a+{\overline{S}})^2}\\ +\frac{1}{2} {\overline{S}} \left( \frac{\beta _2{\overline{E}}}{(a+{\overline{S}})^2}-{\bar{\mu }} \right) {}^2 \left( -\beta _1-\frac{\beta _2}{a+S}\right) ,\ \theta _{21}=-1+\frac{E S^2 \beta _2}{(a+S)^3}+2 S \left( {\bar{\mu }}-\frac{e \beta _2}{(a+{\overline{S}})^2}\right) -\frac{1}{2} {\overline{S}}^2 \left( {\bar{\mu }}-\frac{{\overline{E}} \beta _2}{(a+{\overline{S}})^2}\right) {}^2 ,\\ \theta _{22}=-\frac{\left( (a+{\overline{S}}) \beta _1+\beta _2\right) \left( (a+{\overline{S}})^3 (-1+S {\bar{\mu }}) \beta _1+(a+{\overline{S}}) \left( -a+{\overline{S}}+{\overline{S}} (a+{\overline{S}}) {\bar{\mu }} -{\overline{E}} {\overline{S}} \beta _1\right) \beta _2-{\overline{E}} {\overline{S}} \beta _2^2\right) }{2 (a+{\overline{S}})^4} ,\\ \theta _{23}=2 {\overline{S}}\left( \beta _1+\frac{\beta _2}{a+{\overline{S}}}\right) -\frac{{\overline{S}}^2}{(a+{\overline{S}})^2}+{\overline{S}}^2 \left( \frac{\beta _2{\overline{E}}}{(a+{\overline{S}})^2}-{\bar{\mu }} \right) \left( \beta _1+\frac{\beta _2}{a+{\overline{S}}}\right) ,\ \theta _{24}=\frac{1}{2} {\overline{S}}^2 \left( 3-{\overline{S}} {\bar{\mu }} +\frac{{\overline{E}} {\overline{S}} \beta _2}{(a+{\overline{S}})^2}\right) , \\ \theta _{25}=-\frac{{\overline{S}}^3 \left( (a+{\overline{S}}) \beta _1+\beta _2\right) }{2 (a+{\overline{S}})}, \ \theta _{26}=-\frac{{\overline{S}}^4}{6}, \\ \phi _{11}={\overline{E}} \left( \beta _1+\frac{a \beta _2}{(a+{\overline{S}})^2}\right) ,\ \phi _{12}=1, \ \phi _{13}= \frac{{\overline{E}} \left( (a+{\overline{S}})^4 \beta _1^2+2 a (a+S) \left( -1+(a+S) \beta _1\right) \beta _2+a^2 \beta _2^2\right) }{2 (a+S)^4},\ \phi _{14}=\beta _1+\frac{a \beta _2}{(a+{\overline{S}})^2},\\ \phi _{15}=-1,\ \phi _{16}=\frac{{\overline{E}}}{2}, \ \phi _{17}=\frac{{\overline{E}} \left( -6 {\overline{S}} (a+{\overline{S}})^2 \beta _2+6 (a+{\overline{S}})^3 \beta _2-6 a (a+{\overline{S}}) \beta _2 \left( (a+{\overline{S}})^2 \beta _1+a \beta _2\right) +\left( (a+{\overline{S}})^2 \beta _1+a \beta _2\right) {}^3\right) }{6 (a+{\overline{S}})^6},\\ \phi _{18}=\frac{-2 a (a+{\overline{S}}) \beta _2+\left( (a+{\overline{S}})^2 \beta _1+a \beta _2\right) {}^2}{2 (a+{\overline{S}})^4},\ \phi _{19}=-\frac{{\overline{E}} \left( (a+{\overline{S}})^4 \beta _1^2+2 a (a+{\overline{S}}) \left( (a+{\overline{S}}) \beta _1-1\right) \beta _2+a^2 \beta _2^2\right) }{2 (a+{\overline{S}})^4}, \\ \phi _{20}=-\beta _1-\frac{a \beta _2}{(a+{\overline{S}})^2},\ \phi _{21}= \frac{1}{2} {\overline{E}} \left( \beta _1+\frac{a \beta _2}{(a+{\overline{S}})^2}\right) ,\ \phi _{22}=\frac{1}{2},\ \phi _{23}=-\frac{{\overline{E}}}{6}. \end{array} \right. \end{aligned}$$

Assume that \(\xi _{1}, \xi _{2}\) are eigenvalues for system (15), then we have the following translation;

$$\begin{aligned} \left( \begin{array}{c} x \\ y \\ \end{array} \right) =M\left( \begin{array}{c} u \\ v \\ \end{array} \right) , \end{aligned}$$
(19)

where

$$\begin{aligned} M=\left[ \begin{array}{cc} -{\overline{S}} \left( \beta _1+\frac{\beta _2}{a+{\overline{S}}}\right) &{} -{\overline{S}} \left( \beta _1+\frac{\beta _2}{a+{\overline{S}}}\right) \\ {\overline{S}} {\bar{\mu }} -\frac{{\overline{E}} {\overline{S}} \beta _2}{(a+{\overline{S}})^2}-1&{} \xi _{2}-{\overline{S}} \left( \frac{{\overline{E}} \beta _2}{(a+{\overline{S}})^2}-{\bar{\mu }} \right) \\ \end{array} \right] \end{aligned}$$

be a nonsingular matrix. By applying transformation (19), the map (18) can be written as:

$$\begin{aligned} \left( \begin{array}{c} u \\ v \\ \end{array} \right) \rightarrow \left( \begin{array}{cc} -1 &{} 0 \\ 0&{} \xi _{2}\\ \end{array} \right) \left( \begin{array}{c} u \\ v \\ \end{array} \right) +\left( \begin{array}{c} f(u,v,{\bar{\mu }}) \\ g(u,v,{\bar{\mu }})\\ \end{array} \right) , \end{aligned}$$
(20)

where

$$\begin{aligned} f(u,v,{\bar{\mu }})= & {} \left( {\frac{ \left( \xi _{{2}}-\theta _{{11}} \right) \theta _{{26} }}{\theta _{{12}} \left( \xi _{{2}}+1 \right) }}-{\frac{\phi _{{23}}}{ \xi _{{2}}+1}} \right) {{\bar{\mu }}}^{3}+ \left( {\frac{ \left( \xi _{{2}}- \theta _{{11}} \right) \theta _{{24}}}{\theta _{{12}} \left( \xi _{{2}}+1 \right) }}-{\frac{\phi _{{21}}}{\xi _{{2}}+1}} \right) {{\bar{\mu }}}^{2}x\\&+\left( {\frac{ \left( \xi _{{2}}-\theta _{{11}} \right) \theta _{{25} }}{\theta _{{12}} \left( \xi _{{2}}+1 \right) }}-{\frac{\phi _{{22}}}{ \xi _{{2}}+1}} \right) {{\bar{\mu }}}^{2}y+ \left( {\frac{ \left( \xi _{{2}}- \theta _{{11}} \right) \theta _{{18}}}{\theta _{{12}} \left( \xi _{{2}}+1 \right) }}-{\frac{\phi _{{16}}}{\xi _{{2}}+1}} \right) {{\bar{\mu }}}^{2}\\&+\left( {\frac{ \left( \xi _{{2}}-\theta _{{11}} \right) \theta _{{21} }}{\theta _{{12}} \left( \xi _{{2}}+1 \right) }}-{\frac{\phi _{{19}}}{ \xi _{{2}}+1}} \right) {\bar{\mu }}\,{x}^{2}+ \left( {\frac{ \left( \xi _{{2} }-\theta _{{11}} \right) \theta _{{23}}}{\theta _{{12}} \left( \xi _{{2}}+ 1 \right) }}-{\frac{\phi _{{20}}}{\xi _{{2}}+1}} \right) {\bar{\mu }}\,xy\\&+{\frac{ \left( \xi _{{2}}-\theta _{{11}} \right) \theta _{{15}}{\bar{\mu }}\,x}{ \theta _{{12}} \left( \xi _{{2}}+1 \right) }}+ \left( {\frac{ \left( \xi _{{2}}-\theta _{{11}} \right) \theta _{{17}}}{\theta _{{12}} \left( \xi _{{2}}+1 \right) }}-{\frac{\phi _{{15}}}{\xi _{{2}}+1}} \right) y {\bar{\mu }}\\&+ \left( {\frac{ \left( \xi _{{2}}-\theta _{{11}} \right) \theta _ {{19}}}{\theta _{{12}} \left( \xi _{{2}}+1 \right) }}-{\frac{\phi _{{17} }}{\xi _{{2}}+1}} \right) {x}^{3}+ \left( {\frac{ \left( \xi _{{2}}- \theta _{{11}} \right) \theta _{{20}}}{\theta _{{12}} \left( \xi _{{2}}+1 \right) }}-{\frac{\phi _{{18}}}{\xi _{{2}}+1}} \right) {x}^{2}y\\&+\left( {\frac{ \left( \xi _{{2}}-\theta _{{11}} \right) \theta _{{13} }}{\theta _{{12}} \left( \xi _{{2}}+1 \right) }}-{\frac{\phi _{{13}}}{ \xi _{{2}}+1}} \right) {x}^{2}+{\frac{ \left( \xi _{{2}}-\theta _{{11}} \right) \theta _{{22}}x{y}^{2}}{\theta _{{12}} \left( \xi _{{2}}+1 \right) }}\\&+ \left( {\frac{ \left( \xi _{{2}}-\theta _{{11}} \right) \theta _{{14}}}{\theta _{{12}} \left( \xi _{{2}}+1 \right) }}-{\frac{ \phi _{{14}}}{\xi _{{2}}+1}} \right) xy+{\frac{ \left( \xi _{{2}}- \theta _{{11}} \right) \theta _{{16}}{y}^{2}}{\theta _{{12}} \left( \xi _{ {2}}+1 \right) }}\\&+O\left( (|u|+|v|+|{\bar{\mu }} |)^4\right) , \\ g(u,v,{\bar{\mu }})= & {} \left( {\frac{ \left( 1+\theta _{{11}} \right) \theta _{{26}}}{\theta _ {{12}} \left( \xi _{{2}}+1 \right) }}+{\frac{\phi _{{23}}}{\xi _{{2}}+1} } \right) {{\bar{\mu }}}^{3}+ \left( {\frac{ \left( 1+\theta _{{11}} \right) \theta _{{24}}}{\theta _{{12}} \left( \xi _{{2}}+1 \right) }}+{\frac{ \phi _{{21}}}{\xi _{{2}}+1}} \right) {{\bar{\mu }}}^{2}x\\&+ \left( {\frac{ \left( 1+\theta _{{11}} \right) \theta _{{25}}}{\theta _{{12}} \left( \xi _{{2}}+1 \right) }}+{\frac{\phi _{{22}}}{\xi _{{2}}+1}} \right) { {\bar{\mu }}}^{2}y+ \left( {\frac{ \left( 1+\theta _{{11}} \right) \theta _{{18 }}}{\theta _{{12}} \left( \xi _{{2}}+1 \right) }}+{\frac{\phi _{{16}}}{ \xi _{{2}}+1}} \right) {{\bar{\mu }}}^{2}\\&+ \left( {\frac{ \left( 1+\theta _{{11 }} \right) \theta _{{21}}}{\theta _{{12}} \left( \xi _{{2}}+1 \right) }}+ {\frac{\phi _{{19}}}{\xi _{{2}}+1}} \right) {\bar{\mu }}\,{x}^{2}+ \left( { \frac{ \left( 1+\theta _{{11}} \right) \theta _{{23}}}{\theta _{{12}} \left( \xi _{{2}}+1 \right) }}+{\frac{\phi _{{20}}}{\xi _{{2}}+1}} \right) {\bar{\mu }}\,xy\\&+{\frac{ \left( 1+\theta _{{11}} \right) \theta _{{15} }{\bar{\mu }}\,x}{\theta _{{12}} \left( \xi _{{2}}+1 \right) }}+ \left( {\frac{ \left( 1+\theta _{{11}} \right) \theta _{{17}}}{\theta _{{12}} \left( \xi _{{2}}+1 \right) }}+{\frac{\phi _{{15}}}{\xi _{{2}}+1}} \right) y {\bar{\mu }}+ \left( {\frac{ \left( 1+\theta _{{11}} \right) \theta _{{19}}}{ \theta _{{12}} \left( \xi _{{2}}+1 \right) }}+{\frac{\phi _{{17}}}{\xi _{ {2}}+1}} \right) {x}^{3}\\&+ \left( {\frac{ \left( 1+\theta _{{11}} \right) \theta _{{20}}}{\theta _{{12}} \left( \xi _{{2}}+1 \right) }}+{ \frac{\phi _{{18}}}{\xi _{{2}}+1}} \right) {x}^{2}y+ \left( {\frac{ \left( 1+\theta _{{11}} \right) \theta _{{13}}}{\theta _{{12}} \left( \xi _{{2}}+1 \right) }}+{\frac{\phi _{{13}}}{\xi _{{2}}+1}} \right) {x}^ {2}+{\frac{ \left( 1+\theta _{{11}} \right) \theta _{{22}}x{y}^{2}}{ \theta _{{12}} \left( \xi _{{2}}+1 \right) }}\\&+ \left( {\frac{ \left( 1+ \theta _{{11}} \right) \theta _{{14}}}{\theta _{{12}} \left( \xi _{{2}}+1 \right) }}+{\frac{\phi _{{14}}}{\xi _{{2}}+1}} \right) xy+{\frac{ \left( 1+\theta _{{11}} \right) \theta _{{16}}{y}^{2}}{\theta _{{12}} \left( \xi _{{2}}+1 \right) }}+O\left( (|u|+|v|+|{\bar{\mu }}|)^4\right) , \end{aligned}$$

where,

$$\begin{aligned} x=\theta _{12}(u+v),\ \ \ \ \ \ \ \ \ y=-(1+\theta _{11})u+(\xi _{2}-\theta _{11}) v. \end{aligned}$$

Assume that \(\mho ^c(0,0,0)\) be the center manifold of (20) intended at (0, 0) in a smallest neighborhood of \({\bar{\eta }}=0.\) Then, \(\mho ^c(0,0,0)\) can be estimated as follows:

$$\begin{aligned} \mho ^c(0,0,0)=\left\{ (u,v,{\bar{\mu }})\in {\mathbb {R}}^3:v=M_{11}u^2+M_{12}u{\bar{\mu }}+M_{13}{\bar{\mu }}^2 +O\left( (|{\bar{\mu }}|+|u|)^3\right) \right\} , \end{aligned}$$

where

$$\begin{aligned} M_{11}= & {} \frac{1}{1-\xi _{{2}}} \left( \left( {\frac{ \left( 1+ \theta _{{11}} \right) \theta _{{13}}}{\theta _{{12}} \left( \xi _{{2}}+1 \right) }}+{\frac{\phi _{{13}}}{\xi _{{2}}+1}} \right) {\theta _{{12}}} ^{2}- \left( {\frac{ \left( 1+\theta _{{11}} \right) \theta _{{14}}}{ \theta _{{12}} \left( \xi _{{2}}+1 \right) }}+{\frac{\phi _{{14}}}{\xi _{ {2}}+1}} \right) \theta _{{12}} \left( 1+\theta _{{11}} \right) \right) \\&+\frac{1}{1-\xi _{{2}}}{ \frac{ \left( 1+\theta _{{11}} \right) \theta _{{16}} \left( 1+\theta _ {{11}} \right) ^{2}}{\theta _{{12}} \left( \xi _{{2}}+1 \right) }}, \\ M_{12}= & {} {\frac{1}{1-\xi _{{2}}} \left( {\frac{ \left( 1+\theta _{{11} } \right) \theta _{{15}}}{\xi _{{2}}+1}}- \left( {\frac{ \left( 1+ \theta _{{11}} \right) \theta _{{17}}}{\theta _{{12}} \left( \xi _{{2}}+1 \right) }}+{\frac{\phi _{{15}}}{\xi _{{2}}+1}} \right) \left( 1+ \theta _{{11}} \right) \right) } , \\ M_{13}= & {} {\frac{1}{1-\xi _{{2}}} \left( {\frac{ \left( 1+\theta _{{11} } \right) \theta _{{18}}}{\theta _{{12}} \left( \xi _{{2}}+1 \right) }}+{ \frac{\phi _{{16}}}{\xi _{{2}}+1}} \right) } . \end{aligned}$$

Now, the map restricted to set \(\mho ^c(0,0,0)\) is described as follows:

$$\begin{aligned} G: u\rightarrow -u+k_{11}u^2+k_{12}u{\bar{\mu }}+k_{13}u^2{\bar{\mu }}+k_{14}u{\bar{\mu }}^2 +k_{15}u^3+O\left( (|u|+|{\bar{\mu }}|)^4\right) , \end{aligned}$$

where

$$\begin{aligned} k_{11}= & {} \left( {\frac{ \left( \xi _{{2}}-\theta _{{11}} \right) \theta _{{13}}}{\theta _{{12}} \left( \xi _{{2}}+1 \right) }}-{\frac{ \phi _{{13}}}{\xi _{{2}}+1}} \right) {\theta _{{12}}}^{2}- \left( { \frac{ \left( \xi _{{2}}-\theta _{{11}} \right) \theta _{{14}}}{\theta _ {{12}} \left( \xi _{{2}}+1 \right) }}-{\frac{\phi _{{14}}}{\xi _{{2}}+1} } \right) \theta _{{12}} \left( 1+\theta _{{11}} \right) \\&+{\frac{ \left( \xi _{{2}}-\theta _{{11}} \right) \theta _{{16}} \left( 1+ \theta _{{11}} \right) ^{2}}{\theta _{{12}} \left( \xi _{{2}}+1 \right) } }, \\ k_{12}= & {} {\frac{ \left( \xi _{{2}}-\theta _{{11}} \right) \theta _{{15} }}{\xi _{{2}}+1}}- \left( {\frac{ \left( \xi _{{2}}-\theta _{{11}} \right) \theta _{{17}}}{\theta _{{12}} \left( \xi _{{2}}+1 \right) }}-{ \frac{\phi _{{15}}}{\xi _{{2}}+1}} \right) \left( 1+\theta _{{11}} \right) , \\ k_{13}= & {} \left( {\frac{ \left( \xi _{{2}}-\theta _{{11}} \right) \theta _{{21}}}{\theta _{{12}} \left( \xi _{{2}}+1 \right) }}-{\frac{ \phi _{{19}}}{\xi _{{2}}+1}} \right) {\theta _{{12}}}^{2}- \left( { \frac{ \left( \xi _{{2}}-\theta _{{11}} \right) \theta _{{23}}}{\theta _ {{12}} \left( \xi _{{2}}+1 \right) }}-{\frac{\phi _{{20}}}{\xi _{{2}}+1} } \right) \theta _{{12}} \left( 1+\theta _{{11}} \right) \\&+{\frac{ \left( \xi _{{2}}-\theta _{{11}} \right) \theta _{{15}}M_{{1}}}{\xi _{{2 }}+1}}+ \left( {\frac{ \left( \xi _{{2}}-\theta _{{11}} \right) \theta _{{17}}}{\theta _{{12}} \left( \xi _{{2}}+1 \right) }}-{\frac{ \phi _{{15}}}{\xi _{{2}}+1}} \right) \left( \xi _{{2}}-\theta _{{11}} \right) M_{{1}}\\&+2\, \left( {\frac{ \left( \xi _{{2}}-\theta _{{11}} \right) \theta _{{13}}}{\theta _{{12}} \left( \xi _{{2}}+1 \right) }}-{ \frac{\phi _{{13}}}{\xi _{{2}}+1}} \right) {\theta _{{12}}}^{2}M_{{2}}+ \left( {\frac{ \left( \xi _{{2}}-\theta _{{11}} \right) \theta _{{14} }}{\theta _{{12}} \left( \xi _{{2}}+1 \right) }}-{\frac{\phi _{{14}}}{ \xi _{{2}}+1}} \right) \theta _{{12}} \left( \xi _{{2}}-\theta _{{11}} \right) M_{{2}}\\&- \left( {\frac{ \left( \xi _{{2}}-\theta _{{11}} \right) \theta _{{14}}}{\theta _{{12}} \left( \xi _{{2}}+1 \right) }}-{ \frac{\phi _{{14}}}{\xi _{{2}}+1}} \right) \theta _{{12}}M_{{2}} \left( 1+\theta _{{11}} \right) -2\,{\frac{ \left( \xi _{{2}}-\theta _{{11}} \right) \theta _{{16}} \left( 1+\theta _{{11}} \right) \left( \xi _{{2 }}-\theta _{{11}} \right) M_{{2}}}{\theta _{{12}} \left( \xi _{{2}}+1 \right) }} \\ k_{14}= & {} \left( {\frac{ \left( \xi _{{2}}-\theta _{{11}} \right) \theta _{{24}}}{\theta _{{12}} \left( \xi _{{2}}+1 \right) }}-{\frac{ \phi _{{21}}}{\xi _{{2}}+1}} \right) \theta _{{12}}- \left( {\frac{ \left( \xi _{{2}}-\theta _{{11}} \right) \theta _{{25}}}{\theta _{{12}} \left( \xi _{{2}}+1 \right) }}-{\frac{\phi _{{22}}}{\xi _{{2}}+1}} \right) \left( 1+\theta _{{11}} \right) \\&+{\frac{ \left( \xi _{{2}}-\theta _{{11}} \right) \theta _{{15}}M_{{2}}}{\xi _{{2}}+1}}+ \left( { \frac{ \left( \xi _{{2}}-\theta _{{11}} \right) \theta _{{17}}}{\theta _ {{12}} \left( \xi _{{2}}+1 \right) }}-{\frac{\phi _{{15}}}{\xi _{{2}}+1} } \right) \left( \xi _{{2}}-\theta _{{11}} \right) M_{{2}}\\&+2\, \left( {\frac{ \left( \xi _{{2}}-\theta _{{11}} \right) \theta _{{13}}}{\theta _{{12}} \left( \xi _{{2}}+1 \right) }}-{\frac{\phi _{{13}}}{\xi _{{2}}+1 }} \right) {\theta _{{12}}}^{2}M_{{3}}-2\,{ \frac{ \left( \xi _{{2}}-\theta _{{11}} \right) \theta _{{16}} \left( 1+\theta _{{11}} \right) \left( \xi _{{2}}-\theta _{{11}} \right) M_{{3} }}{\theta _{{12}} \left( \xi _{{2}}+1 \right) }}\\&+ \left( {\frac{ \left( \xi _{{ 2}}-\theta _{{11}} \right) \theta _{{14}}}{\theta _{{12}} \left( \xi _{{2} }+1 \right) }}-{\frac{\phi _{{14}}}{\xi _{{2}}+1}} \right) \theta _{{12} } \left( \xi _{{2}}-\theta _{{11}} \right) M_{{3}}- \left( {\frac{\left( \xi _{{2}}-\theta _{{11}} \right) \theta _{{14}}}{\theta _{{12}} \left( \xi _{{2}}+1 \right) }}-{\frac{\phi _{{14}}}{\xi _{{2}}+1}} \right) \theta _{{12}}M_{{3}} \left( 1+\theta _{{11}} \right) , \\ k_{15}= & {} {\frac{ \left( \xi _{{2}}-\theta _{{11}} \right) \theta _{{26} }}{\theta _{{12}} \left( \xi _{{2}}+1 \right) }}-{\frac{\phi _{{23}}}{ \xi _{{2}}+1}}+{\frac{ \left( \xi _{{2}}-\theta _{{11}} \right) \theta _ {{15}}M_{{3}}}{\xi _{{2}}+1}}+ \left( {\frac{ \left( \xi _{{2}}- \theta _{{11}} \right) \theta _{{17}}}{\theta _{{12}} \left( \xi _{{2}}+1 \right) }}-{\frac{\phi _{{15}}}{\xi _{{2}}+1}} \right) \left( \xi _{{2 }}-\theta _{{11}} \right) M_{{3}}. \end{aligned}$$

Next, we have the following real numbers:

$$\begin{aligned} L_{11}= & {} \left( \frac{\partial ^2f_{1}}{\partial u \partial {\bar{\mu }}}+\frac{1}{2}\frac{\partial G}{\partial {\bar{\mu }}}\frac{\partial ^2G}{\partial u^2}\right) _{(0,0)} =\frac{\left( \left( 1+\theta _{11}\right) \theta _{17}-\theta _{12} \theta _{15}\right) \left( \theta _{11}-\xi _2\right) +\left( 1+\theta _{11}\right) \theta _{12} \phi _{15}}{\theta _{12} \left( 1+\xi _2\right) } \ne 0,\\ L_{12}= & {} \left( \left( \frac{1}{2}\frac{\partial ^2G}{\partial u^2}\right) ^2+\frac{1}{6}\frac{\partial ^3G}{\partial u^3}\right) _{(0,0)}={k_{11}}^2+{k_{15}}\ne 0. \end{aligned}$$

Hence, by aforementioned study we have the next conclusive theorem related to the existence of period-doubling bifurcation of system (4) about \(({\overline{S}}, {\overline{E}}, 0).\)

Theorem 5.1

If \(L_{12}\ne 0\), then system (4) undergoes period-doubling bifurcation at the isolation free equilibrium \(({\overline{S}}, {\overline{E}}, 0)\) when parameter \(\mu \) varies in small neighborhood of \({\bar{\mu }}.\) Furthermore, if \(L_{12}>0\), then the period-two orbits that bifurcate from \(({\overline{S}}, {\overline{E}}, 0)\) are stable, and if \(L_{12}<0\), then these orbits are unstable.

6 Neimark–Sacker bifurcation

This section is related to the bifurcation analysis of the system (4) about \((S^* ,E^* ,Q^* ).\) Where all conditions for existence and positivity of \((S^* ,E^* ,Q^* )\) are given in Lemmas 2.3 and 3.3. Here, we will discuss the Neimark–Sacker bifurcation experienced by system (4) about \((S^* ,E^* ,Q^* )\) under some mathematical conditions. Bifurcation is the mathematical phenomena produced in any system due to creation of very small change in stability of system. Mathematically, bifurcation arises whenever parameters are varied in very least neighborhood of equilibrium point. Moreover, for further study of bifurcation theory and to understand this surprising behavior of a discrete-time mathematical system one can see [31,32,33,34]. We deliberate the Neimark–Sacker Bifurcation for positive equilibrium point \((S^* ,E^* ,Q^* )\) of system (4) by using an obvious criterion for Neimark–Sacker Bifurcation and compelling \(\mu \) as a bifurcation parameter. Due to appearance of Neimark–Sacker Bifurcation, closed invariant circles are formed. Equally, one can find some lonely orbits of periodic performance alongside with paths that cover the invariant circle tightly. The bifurcation can be supercritical or subcritical causing in a stable or unstable closed invariant curve, correspondingly. In command to study the Neimark–Sacker bifurcation in system (4), we have the next obvious standard of Hopf bifurcation [35]. By means of this standard, one can catch the presence of Neimark–Sacker bifurcation deprived of finding the eigenvalues.

Lemma 6.1

(see [35]) Consider an m-dimensional discrete dynamical system \(Y_{K+1}=g_{\alpha }(Y_{K})\), where \(\alpha \in \mathfrak {R}\) is bifurcation parameter. Let \(Y^*\) be a fixed point of \(g_{\alpha }\) and the characteristic polynomial for Jacobian matrix \(J(Y^*)=(s_{ij})_{m\times m}\) of m-dimensional map \(g_{\alpha }\) is given by

$$\begin{aligned} {\mathbb {H}}_{\alpha }(\kappa )=\kappa ^{m}+b_1\kappa ^{m-1}+b_2\kappa ^{m-2} +........+b_{m-1}\kappa +b_{m} \end{aligned}$$
(21)

where \(b_{i}=b_{i}(\alpha ,s),\) \(i=1,2,3,....,m\) and s is control parameter or another parameter to be determined. Let \({\Delta }_{0}^{\pm }(\alpha , s)=1,{\Delta }_{1}^{\pm }(\alpha , s),{\Delta }_{2}^{\pm }(\alpha , s),.....,{\Delta }_{m}^{\pm }(\alpha , s).\) be the sequence of determinants defined by \({\Delta }_{i}^{\pm }(\alpha , s)=det(N_{1}\pm N_{2}), i=1,2,3,....,m\) where

$$\begin{aligned} N_{1}= & {} \left( \begin{array}{ccccc} 1 &{} b_{1} &{} b_{2} &{} .... &{} b_{i-1} \\ 0 &{} 1 &{} b_{1} &{} .... &{} b_{i-2} \\ 0 &{} 0 &{} 1 &{} .... &{} b_{i-3} \\ .... &{} .... &{} .... &{} .... &{} .... \\ 0 &{} 0 &{} 0 &{} .... &{} 1 \\ \end{array} \right) \end{aligned}$$
(22)
$$\begin{aligned} N_{2}= & {} \left( \begin{array}{ccccc} b_{m-i+1} &{} b_{m-i+2}&{}.... &{} b_{m-1} &{} b_{m} \\ b_{m-i+2} &{} b_{m-i+3} &{} .... &{} b_{m} &{} 0 \\ b_{m-i+3} &{} b_{m-i+4} &{} .... &{} 0 &{} 0 \\ .... &{} .... &{} .... &{} .... &{} .... \\ b_{m} &{} 0 &{} 0 &{} .... &{} 0 \\ \end{array} \right) . \end{aligned}$$
(23)

Moreover, the following conditions are satisfied:

\(C_{1}\):

Eigenvalue assignment: \({\Delta }_{m-1}^{-}(\alpha _{0},s)=0,{\Delta }_{m-1}^{+}(\alpha _{0},s)>0, {\mathbb {H}}_{\alpha _{0}}(1)>0\), \((-1)^m{\mathbb {H}}_{\alpha _{0}}(-1)>0\), \({\Delta }_{i}^{\pm }(\alpha _{0},s)>0,i=m-3,m-5,m-7,....,1\) or \(i=m-3,m-5,m-7,....,2\) when m is even or odd, respectively.

\(C_{2}\):

Transversality condition: \(\left[ \frac{d}{d \alpha }({\Delta }_{m-1}^{-}(\alpha ,s))\right] _{\alpha =\alpha _{0}}\ne 0\)

\(C_{3}\):

Nonresonance condition: \(cos(\frac{2 \pi }{n})\ne \psi ,\) or resonance condition \(cos(\frac{2 \pi }{n})= \psi ,\) where \(n=3,4,5,....,\) and \(\psi =\frac{-1+0.5 {\mathbb {H}}_{\alpha _{0}}(1){\Delta }_{m-3}^{-}(\alpha _{0},s)}{{\Delta }_{m-2}^{+} (\alpha _{0},s)}.\) Then, Neimark–Sacker bifurcation exits at \(\alpha _{0}.\)

The following result shows that system (4) undergoes Neimark–Sacker bifurcation if we take \(\alpha \) as bifurcation parameter.

Theorem 6.1

The unique positive equilibrium point of system (4) undergoes Neimark–Sacker bifurcation if the following conditions hold:

$$\begin{aligned} \left\{ \begin{array}{ll} 1-A_{2}+A_{3}(A_{1}-A_{3})=0, \\ 1+A_{2}-A_{3}(A_{1}+A_{3})>0, \\ 1+A_{1}+A_{2}+A_{3}>0,\\ 1-A_{1}+A_{2}-A_{3}>0. \end{array} \right. \end{aligned}$$

where, \(A_{1},A_{2}\) and \(A_{3}\) are given in (13).

Proof

According to Lemma 6.1, for \(m=3\), we have in (14) the characteristic polynomial of system (4) evaluated at its unique positive equilibrium, then we obtain the following equalities and inequalities:

$$\begin{aligned} \left\{ \begin{array}{ll} {\Delta }_{2}^{-}(\mu )=1-A_{2}+A_{3}(A_{1}-A_{3})=0, \\ {\Delta }_{2}^{+}(\mu )=1+A_{2}-A_{3}(A_{1}+A_{3})>0,\\ {\mathbb {H}}_{\mu }(1)=1+A_{1}+A_{2}+A_{3}>0,\\ (-1)^3{\mathbb {H}}_{\mu }(-1)=1-A_{1}+A_{2}-A_{3}>0. \end{array} \right. \end{aligned}$$

\(\square \)

Remark 6.1

The arithmetical rule for \(\mu \) at which Neimark–Sacker bifurcation arises can be established by finding the solutions of the equation \({\Delta }_{2}^{-}(\mu )=0\).

7 Modified hybrid control of Bifurcation and Chaos

Control of bifurcation and chaos in mathematical models is considered as a key element for population models mainly when these models are associated with biological interactions and breeding of different species. Generally discrete-time mathematical systems are more complex to analyze as compared to continuous one. It is necessary that the population does not experience any irregular situation. Hence, to avoid these irregularities a chaos controlling technique must be implemented. In this part of manuscript, we study a feedback control strategy with parameter perturbation to move unstable and irregular trajectories toward the stable trajectories. The most useful and well-known method in the field of chaos is given by Ott et al. [36] to control period-doubling bifurcation, which is known as OGY method. Latter on, numerous strategic control methods are developed (see [22, 37]). Here, we present a modified hybrid control method to control the Neimark–Sacker bifurcation and chaos. Furthermore, this mathematical method is well applicable to every discrete-time system experiencing the period-doubling bifurcation and chaos. Originally, hybrid method was proposed by Liu et al. [27]. Moreover, it was developed to control the period-doubling bifurcation (see [26,27,28,29,30]). Here, we have used the modified hybrid control technique [22] to control Neimark–Sacker bifurcation. Moreover, this technique is much better then old techniques of control. Consider the following n-dimensional discrete dynamical system:

$$\begin{aligned} Z_{n+1}=g(Z_n,\ell ) \end{aligned}$$
(24)

with \(Z_n \in \mathfrak {R}^n\), \(n \in Z.\) Suppose \(\ell \in \mathfrak {R}\) is parameter for which system (24) experiences the bifurcation. The purpose of proposing the modified technique for controlling the bifurcation is to regain the maximum range of stable region in (24) by reduction in length of unstable region. Hence, we present the following generalized hybrid control technique by applying state feedback along with parameter perturbation;

$$\begin{aligned} Z_{n+k}=\theta ^3 g^{(\hbar )}(Z_n, \ell )+ (1-\theta ^3)Z_n \end{aligned}$$
(25)

where \(\hbar >0\) is in Z and \(0<{\theta }<1\) is parameter for controlling the bifurcation appearing in (25). In addition, \( g^{(\hbar )}\) is kth iterative value of g(.). By application of (25) on system (4), we get the following system;

$$\begin{aligned} \begin{aligned} S_{n+1}&=\theta ^3\left( S_{n}e^{B-\beta _{1}E_{n}-\frac{\beta _{2}E_{n}}{a+S_{n}}-\frac{\beta _{4}Q_{n}}{b+S_{n}}-\mu S_{n}}\right) +(1-\theta ^3)S_{n}, \\ E_{n+1}&=\theta ^3\left( E_{n}e^{\beta _{1}S_{n}+\frac{\beta _{2}S_{n}}{a+S_{n}}-\frac{\beta _{3}Q_{n}}{d+E_{n}}-\mu }\right) +(1-\theta ^3)E_{n},\\ Q_{n+1}&=\theta ^3\left( Q_{n}e^{\frac{c_{1}\beta _{3}E_{n}}{d+E_{n}}+\frac{c_{2}\beta _{4}S_{n}}{b+S_{n}}-\mu }\right) +(1-\theta ^3)Q_{n} \end{aligned} \end{aligned}$$
(26)

Furthermore, the system (26) and system (4) have same constant solutions. Additionally, the Jacobian matrix of (26) about \((S^* ,E^* ,Q^* )\) is given as follows:

$$\begin{aligned} \left( \begin{array}{ccc} 1+ \theta ^3 S^* \left( \mu -\frac{\beta _2 e^*}{\left( a+S^*\right) ^2}-\frac{\beta _4 Q^*}{\left( b+S^*\right) ^2}\right) &{} -\theta ^3 S^* \left( \beta _1+\frac{\beta _2}{a+S^*}\right) &{} -\frac{ \theta ^3 \beta _4 S^*}{b+S^*} \\ \frac{\theta ^3 E^* \left( a \beta _2+\beta _1 \left( a+S^*\right) ^2\right) }{\left( a+S^*\right) ^2} &{} 1-\theta ^3+\frac{\theta ^3 \left( \left( d+E^*\right) ^2+\beta _3 E^* Q^*\right) }{\left( d+E^*\right) ^2} &{} -\frac{\theta ^3 \beta _3 E^*}{d+E^*} \\ -\frac{b \theta ^3 c_2 \beta _4 Q^*}{\left( b+S^*\right) ^2} &{} \frac{d \theta ^3 c_1 \beta _3 Q^*}{\left( d+e^*\right) ^2} &{} 1 \end{array} \right) . \end{aligned}$$
(27)

The one and only positive equilibrium point \((S^* ,E^* ,Q^* )\) of the controlled system (26) is locally asymptotically stable, if all solutions of the characteristic polynomial of (27) lie inside \(D_{1}.\) Where \(D_{1}\) is an open unit disk.

8 Numerical simulation

In this part of article, the numerical analysis for the dynamics of (4) is provide. Moreover, this study is the direct verification of our theoretical analysis and analytic results which we proved in previous sections.

Example 8.1

Assume that \(a=2, b=10, \beta _1=0.99, d=0.4, \beta _2=0.01, \beta _3=0.80, \beta _4=0.60, c_1=1, c_2=2\) and \(B>0.\) Then, the discrete-time mathematical system (4) takes the following form:

$$\begin{aligned} \begin{aligned} S_{n+1}&=S_{n}e^{B-0.99E_{n}-\frac{0.01E_{n}}{2+S_{n}}-\frac{0.60Q_{n}}{10+S_{n}}-\mu S_{n}}, \\ E_{n+1}&=E_{n}e^{0.99S_{n}+\frac{0.01S_{n}}{2+S_{n}}-\frac{0.80Q_{n}}{0.4+E_{n}}-\mu },\\ Q_{n+1}&=Q_{n}e^{\frac{0.80E_{n}}{0.4+E_{n}}+\frac{2(0.60S)_{n}}{10+S_{n}}-\mu }, \end{aligned} \end{aligned}$$
(28)

where \(\mu >0\) and \(S_0=0.91059, E_0=0.38854, Q_0=0.3549 \) are initial conditions. In this case, the graphical behavior of both population variables is shown in (Fig. 4). In (Fig. 4c), maximum Lyapunov exponents for the existence of bifurcation are given. In (Fig. 5), some phase portraits are given for variations of \(\mu >0\). Hence, it can be easily seen that there exists the Neimark–Sacker bifurcation for large range of bifurcation parameter \(\mu .\) For aforementioned values of parameters, one can obtain the Jacobian matrix \(F_{J_{3}}\) as follows:

$$\begin{aligned} F_{J_{3}}=\left( \begin{array}{ccc} 1+0.91059 (0.0021294 -\mu ) &{} -0.904613 &{} -0.0500756 \\ 0.286336 &{} 1.1728 &{} -0.335248 \\ -0.0357759 &{} 0.239551 &{} 1 \end{array} \right) . \end{aligned}$$

Moreover, the characteristic equation \({\mathbb {H}}(\psi )=0\) for \(F_{J_{3}}\) has the following coefficients:

$$\begin{aligned} A_{1}= & {} 3.17946 -0.91059\,\mu , \\ A_{2}= & {} 3.68735-1.97853\,\mu , \\ A_{3}= & {} 1.58237-1.1377\,\mu . \end{aligned}$$

Moreover, the equilibrium point \((S^*, E^*, Q^*)\) will be locally asymptotically stable for

$$\begin{aligned} 0.37858611389735736<\mu <1.0675487017977925. \end{aligned}$$

In addition, the one and only positive equilibrium point \((S^*, E^*, Q^*)\) experiences the Neimark–Sacker bifurcation whenever

$$\begin{aligned} 1.0675487017977952<\mu <2.3456458284152606. \end{aligned}$$
Fig. 4
figure 4

Existence of Neimark–Sacker bifurcation in system (4) for \(a=2, b=10, \beta _1=0.99, d=0.4, \beta _2=0.01, \beta _3=0.80, \beta _4=0.60, c_1=1, c_2=2\) , \(B>0 ,\mu >0\) and \(S_0=0.91059, E_0=0.38854, Q_0=0.3549\)

Full size image
Fig. 5
figure 5

Phase portraits of system (4) for \(a=2, b=10, \beta _1=0.99, d=0.4, \beta _2=0.01, \beta _3=0.80, \beta _4=0.60, c_1=1, c_2=2\) , \(B>0 ,\mu >0\) and \(S_0=0.91059, E_0=0.38854, Q_0=0.3549\)

Full size image

Example 8.2

Assume that \(a=2, b=10, \beta _1=0.99, d=0.4, \beta _2=0.01,\) \(B>\mu \) and \(Q_{n}\rightarrow 0, \forall n.\) Then, the discrete-time mathematical system (4) takes the following form:

$$\begin{aligned} \begin{aligned} S_{n+1}&=S_{n}e^{B-0.99E_{n}-\frac{0.01E_{n}}{2+S_{n}}-\mu S_{n}}, \\ E_{n+1}&=E_{n}e^{0.99S_{n}+\frac{0.01S_{n}}{2+S_{n}}-\mu } \end{aligned} \end{aligned}$$
(29)

where \(S_0=1.4391059, E_0=1.138854\) are initial conditions. In this case, the graphical behavior of population variables \(S_{n}\) and \(E_{n}\) is, respectively, shown in Fig. 6a and b. It is clearly seen that in the absence of isolated compartment the system (4) experiences the chaos for higher values of death parameter. In addition, the isolation-free equilibrium point \(({\overline{S}}, {\overline{E}}, 0)\) experiences the period-doubling bifurcation whenever we have

$$\begin{aligned} 1.26997017977952<\mu <1.9887756458284152. \end{aligned}$$
Fig. 6
figure 6

Existence of period-doubling bifurcation in system (4) for \(a=2, b=10, \beta _1=0.99, d=0.4, \beta _2=0.01,\) \(B>\mu \) and \(Q_{n}\rightarrow 0, \forall n\) and \(S_0=1.4391059, E_0=1.138854\)

Full size image

Example 8.3

Assume that \(a=2, b=10, \beta _1=0.99, d=0.4, \beta _2=0.01, \beta _3=0.80, \beta _4=0.60, c_1=1, c_2=2,\) \(B>0\) and \(\mu >0\) with initial conditions \(S_0=1.35, E_0=0.31, Q_0=0.75\) in system (4). Then, in this case the graphical behavior of susceptible population is shown in (Fig. 7). Furthermore, it can be seen that the system (4) undergoes period-doubling bifurcation for higher values of \(\mu \) (see Fig. 7a). Additionally, the existence of chaos for susceptible population can be seen from (Fig. 7b).

Fig. 7
figure 7

Period-doubling bifurcation and chaos in system (4) for \(a=2, b=10, \beta _1=0.99, d=0.4, \beta _2=0.01, \beta _3=0.80, \beta _4=0.60, c_1=1, c_2=2\) , \(B>0 ,\mu >0\) and \(S_0=1.35, E_0=0.31, Q_0=0.75\)

Full size image

Example 8.4

Assume that \(a=2, b=10, \beta _1=0.99, d=0.4, \beta _2=0.01, \beta _3=0.80, \beta _4=0.60, c_1=1, c_2=2\) and \(B>0.\) Then, the discrete-time mathematical system (26) takes the following form:

$$\begin{aligned} \begin{aligned} S_{n+1}&=\theta ^{3}\left( S_{n}e^{B-0.99E_{n}-\frac{0.01E_{n}}{2+S_{n}}-\frac{0.60Q_{n}}{10+S_{n}}-\mu S_{n}}\right) +(1-\theta ^{3})S_{n}, \\ E_{n+1}&=\theta ^{3}\left( E_{n}e^{0.99S_{n}+\frac{0.01S_{n}}{2+S_{n}}-\frac{0.80Q_{n}}{0.4+E_{n}}-\mu }\right) +(1-\theta ^{3})E_{n},\\ Q_{n+1}&=\theta ^{3}\left( Q_{n}e^{\frac{0.80E_{n}}{0.4+E_{n}}+\frac{2(0.60S)_{n}}{10+S_{n}}-\mu }\right) +(1-\theta ^{3})Q_{n}, \end{aligned} \end{aligned}$$
(30)

where \(\mu >0\) ,\(S_0=0.91059, E_0=0.38854, Q_0=0.3549 \) and \(\theta \in [0,1].\) In this case, the graphical behavior of both population variables is shown in (Fig. 8). Hence, it can be easily seen that the Neimark–Sacker bifurcation has been controlled for large range of control parameter \(\theta .\)

Fig. 8
figure 8

Control diagrams for Neimark–Sacker bifurcation in system (26) for \(a=2, b=10, \beta _1=0.99, d=0.4, \beta _2=0.01, \beta _3=0.80, \beta _4=0.60, c_1=1, c_2=2\) , \(B>0 ,\mu >0\) and \(S_0=0.91059, E_0=0.38854, Q_0=0.3549\)

Full size image

Example 8.5

Assume that \(a=2, b=10, \beta _1=0.99, d=0.4, \beta _2=0.01, \beta _3=0.80, \beta _4=0.60, c_1=1, c_2=2\) and \(B>0\) in system (26). Then, for \(\mu >0\) with initial conditions \(S_0=1.35, E_0=0.31,\) and \(Q_0=0.75\) the graphical behavior of susceptible population is shown in (Fig. 9). Hence, it can be easily seen that the period-doubling bifurcation is controlled for large range of control parameter \(\theta .\) Additionally, from (Fig. 9b) it can be seen that the chaos in system (26) has been controlled effectively large range of control parameter \(\theta .\)

Fig. 9
figure 9

Control diagrams for period-doubling bifurcation and chaos in system (26) for \(a=2, b=10, \beta _1=0.99, d=0.4, \beta _2=0.01, \beta _3=0.80, \beta _4=0.60, c_1=1, c_2=2\) , \(B>0 ,\mu >0\) and \(S_0=1.35, E_0=0.31, Q_0=0.75\)

Full size image

Example 8.6

Assume that \(a=2, b=10, \beta _1=0.99, d=0.4, \beta _2=0.01, \beta _3=0.80, \beta _4=0.60, c_1=1, c_2=2\) and \(B>0.\) Then, the continuous-time mathematical system (1) takes the following form:

$$\begin{aligned} \begin{aligned} \frac{\mathrm{d}S}{\mathrm{d}t}&=BS-0.99SE-\frac{0.01SE}{2+S}-\frac{0.60SQ}{10+S}-\mu S^{2}, \\ \frac{\mathrm{d}E}{\mathrm{d}t}&=0.99SE+\frac{0.01SE}{2+S}-\frac{0.80EQ}{0.4+E}-\mu E,\\ \frac{\mathrm{d}Q}{\mathrm{d}t}&=\frac{1(0.80)QE}{0.4+E}+\frac{2(0.60)SQ}{10+S}-\mu Q. \end{aligned} \end{aligned}$$
(31)

From (Fig. 11) it can be seen that the system (31) is stable whenever the death rate \(\mu \le 1.\) Furthermore, in (Fig. 10) the stable behavior of each population variable, namely, S, E and Q is represented in (Fig. 10a–c), respectively. In (Fig. 10), some stable phase portraits are given. Moreover, the system (31) will remain unstable whenever the death parameter \(\mu \) is taken as \(\mu > 1.\)

Fig. 10
figure 10

Two-dimensional plots and phase portraits for system (31) for \(a=2, b=10, \beta _1=0.99, d=0.4, \beta _2=0.01, \beta _3=0.80, \beta _4=0.60, c_1=1, c_2=2\) , \(B>0 ,\mu >0\) and \(S_0=1.35, E_0=0.31, Q_0=0.75\)

Full size image
Fig. 11
figure 11

Three-dimensional phase portraits for system (31) for \(a=2, b=10, \beta _1=0.99, d=0.4, \beta _2=0.01, \beta _3=0.80, \beta _4=0.60, c_1=1, c_2=2\) , \(B>0 ,\mu >0\) and \(S_0=1.35, E_0=0.31, Q_0=0.75\)

Full size image

9 Concluding remarks

We study the qualitative behavior of a discrete-time mathematical model for the population suffering from hypertension or diabetes exposed to COVID-19. The continuous-time counterpart of our model was modeled and analyzed in [1] with Z-control. In our study, we discretize the model by using piecewise constant arguments and provides a better graphical and theoretical analysis of model. In addition, in this paper, we considered the reported cases from the early stages of pandemic to the number of cases in the start of year 2020 in country India [1]. Firstly, by assuming the condition that exposed population(E) will remain finite, the boundedness of every positive solution of system (4) is discussed and an explicit lemma for the boundedness of every positive solution of system (4) is provided in Sect. 2. It is shown that there exist six equilibria for system (4). The stability of system (4) is discussed about each of its equilibrium point. Additionally, we have evaluated some mathematical results concerned to existence of one and only positive equilibrium point and some conditions are developed for local asymptotic stability of positive equilibrium point. It is shown that in the absence of quarantined compartment, the system (4) undergoes period-doubling bifurcation and chaos (see Sect. 9). In order to show the complexity in system (4), the existence of Neimark–Sacker bifurcation for one and only positive equilibrium point is proved mathematically. Through numerical study, we show that system (4) undergoes Neimark–Sacker for wide range of bifurcation parameter \(\mu .\) Moreover, it is shown that for discrete-time mathematical system (4) and its continuous counterpart (1), if we take death parameter \(\mu \le 1\), then both systems are stable and for \(\mu >1\) these systems are unstable. A stability comparison of continuous-time mathematical (1) system with its discrete counterpart (4) for death parameter \(\mu \) is given in Examples 8.1 and 8.6. It is easy to see that our discrete-time mathematical system undergoes chaos when an individual exposed to COVID-19 does not quarantined himself, that is, for \(Q_{n}\rightarrow 0, \forall n,\) the system (4) must be chaotic but for system (1) if \(Q(t)\rightarrow 0\), then system (1) reduces to a two-dimensional continuous-time system, in which chaos is ceased to exist (see [38, 39]). Thus, our discrete-time mathematical system (4) will remain bounded whenever E is finite and (4) become unbounded when E is unbounded. In addition, an individual which is exposed to COVID-19 must have to quarantine himself such that the system (4) remains stable.