Neimark–Sacker bifurcation of a chemotherapy treatment of glioblastoma multiform (GBM)

Abstract

In this paper, we propose a system of differential equations with piecewise constant arguments to describe the growth of GBM under chemotherapeutic treatment and the interaction among the glial cells, the cancer cells, and the chemotherapeutic agents. In this work, the cancer cells are considered as two populations: the sensitive cancer cells and the resistant cancer cells. The sensitive tumor cells produce a population that is known as the resistant cell population, where this population has more resistance to the drug treatment than the sensitive tumor cell population. We analyze at first the local and global stability of the positive equilibrium point by considering the Schur–Cohn criteria and constructing a suitable Lyapunov function, respectively. Moreover, we use the center manifold theorem and bifurcation theory to show that the model undergoes Neimark–Sacker bifurcation. To investigate the case for the extinction of the tumor population, we consider the Allee threshold at time t. Simulation results support the theoretical study.

1 Introduction

Cancer is one of the greatest killers in the world, and the control of tumor growth requires special attention [1]. The most common malignant intrinsic primary tumors of the adult human brain are the gliomas [2]. This name is coined because cancer attacks the glial cells, leaving the neurons intact. Glial cells are responsible for delivering nutrients to neurons, providing structural support to them, and controling the biochemical compositions of the fluid surrounding the neurons [3, 4]. Their death would certainly affect the neurons.

To construct a mathematical model about GBM and the treatment process,we have to consider the destroying effect of the tumor population as well as the dosage of the therapy. In addition, to model biological phenomena, we have to consider the overlapping and non-overlapping generation cases of populations, such as:

• (a) For an overlapping generation of a single species, it is preferred to use a model of differential equations.

• (b) If there is a non-overlapping generation of a single species, then it is convenient to construct a model with a difference equation [5,6,7].

However, for both time situations, continuous and discrete, there are some population dynamics in the ecosystem, which combine the properties of both differential and difference equations.

In this study, there are two events to realize; the first event is the continuity of tumor growth where the resting and converting time of the tumor population occur in discrete times. The second event is the treatment process, which also occurs in discrete times, while the interaction between glial cells and tumor cells is in continuous-time.

For both biological events, it is suitable to construct a model with piecewise constant arguments. Some of these studies can be considered in [8,9,10,11,12,13,14], while other studies in mathematical modeling can be seen in [15,16,17,18,19,20,21,22,23,24].

In this paper, we construct a model as a system of differential equations with piecewise constant arguments, where we use the continuous-time to explain the tumor growth and the interaction between the glial cells and the tumor cells. However, to explain the mutation from the sensitive cells to the resistant cells and the treatment process, we consider the discrete-time.

In our model, glioblastoma multiform produces after a specific time another tumor population which is more resistant to the treatment than the sensitive tumor cells. Both cancer cells attack the glial cells. A cancerous cell never returns to normal, resulting in invasion and destruction of surrounding healthy tissue by cancer cells [25,26,27,28]. The cancer cells only attack the glial cells. The chemotherapeutic agent behaves like a predator acting on all cells. Thus, the model shows interaction among glial cells, cancer cells (sensitive-resistant), and the chemotherapeutic agent.

This biological phenomenon is as follows:

$$ \textstyle\begin{cases} \frac{dG}{dt} = r_{1} G ( t ) ( K_{1} - \alpha _{1} G ( t ) - \alpha _{2} G ( [\!\![t ]\!\!]) ) - \mu _{1} G ( t ) S ( [\!\![t ]\!\!]) - \mu _{3} G ( t ) R ( [\!\![t ]\!\!]) \\ \hphantom{\frac{dG}{dt} =}{}- \frac{P _{1} G ( t ) C ( [\!\![t ]\!\!])}{K_{1} + G ( [\!\![t ]\!\!])}, \\ \frac{{dS}}{{dt}} ={pS} ( t ) + r_{2} S ( t ) ( K_{2} - \beta _{1} S ( t ) - \beta _{2} S ( [\!\![t ]\!\!]) ) - \mu _{2} S ( t ) G ( [\!\![t ]\!\!]) -\rho S ( t ) R ( [\!\![t ]\!\!]) \\ \hphantom{\frac{{dS}}{{dt}} =}{}- \frac{P_{2} S ( t ) C ( [\!\![t ]\!\!])}{K _{2} + S ( [\!\![t ]\!\!])}, \\ \frac{{dR}}{{dt}} = r_{3} R ( t ) ( K_{3} - \gamma _{1} R ( t ) - \gamma _{2} R ( [\!\![t ]\!\!]) ) - \mu _{4} R ( t ) G ( [\!\![t ]\!\!]) +\rho S ( [\!\![t ]\!\!]) R ( t ) \\ \hphantom{\frac{{dR}}{{dt}} =}{}- \frac{P _{3} R ( t ) C ( [\!\![t ]\!\!])}{K_{3} + R ( [\!\![t ]\!\!])}, \\ \frac{{dC}}{{dt}} =\sigma - \omega _{1} C ( t ) - \omega _{2} C ( [\!\![t ]\!\!]), \end{cases} $$

(1.1)

where the parameters denote positive numbers and $[\!\![t ]\!\!]$ is the integer part of $t\in [0,\infty)$. $G(t)$ represents the glial cells concentration (in kg/m³), $S(t)$ represents the sensitive cancerous cells concentration (in kg/m³), $R(t)$ is the resistant cancerous cells concentration (in kg/m³), and $C ( t )$ is the concentration of the chemotherapeutic agent (in mg/m²). The first three equations of system (1.1) shows a logistic growth of each population and the interaction between the glial and cancer cells. The second and third equation of this system involves a converting rate from sensitive cells to resistant cells, while the chemotherapeutic agent effects all three classes $G(t)$, $S(t)$ and $R(t)$. The last equation of system (1.1) describes the dynamics of the chemotherapeutic agent, presenting an exponential decay in concentration. Figure 1 shows the schematic representation of the interaction between the glial cells and the tumor cells, the effet of the chemotherapeutic agent to each class, and the mutation of sensitive tumor cells to resistant cells.

Figure 1

Schematic representation of the interaction between the glial cells and the tumor cells, the chemotherapy agent, and the mutation of the sensitive tumor cells to the resistant cells

First of all, we want to find the positive equilibrium point $\varLambda _{1} =( \overline{G}, \overline{S}, \overline{R}, \overline{C})$ of system (1.1) to work on the scenario of interaction. Since (1.1) is a system of differential equations with piecewise constant arguments, the solution is obtained in each subinterval of the form $[n,n+1)$, $n \in \mathbb{N}$. For $t\in [n,n+1)$, we have $[\!\![t ]\!\!]=n$, and system (1.1) becomes

$$ \textstyle\begin{cases} \frac{dG}{dt} =G(t) ( \xi _{1} - \alpha _{1} r_{1} G ( n ) ), \\ \frac{dS}{dt} =S(t) ( \xi _{2} - \beta _{1} r_{2} S ( n ) ), \\ \frac{dR}{dt} =R(t) ( \xi _{3} - \gamma _{1} r_{3} R ( n ) ), \\ \frac{dC}{dt} =\sigma - \omega _{1} C ( t ) - \omega _{2} C(n), \end{cases} $$

(1.2)

where

$$\begin{aligned}& \xi _{1} (n)= r_{1} K_{1} - \alpha _{2} r _{1} G ( n ) - \mu _{1} S ( n ) - \mu _{3} R ( n ) - \frac{P_{1} C(n)}{K_{1} +G(n)}, \\& \xi _{2} (n)=p+ r_{2} K_{2} - \beta _{2} r _{2} S ( n ) - \mu _{2} G ( n ) -\rho R ( n ) - \frac{P_{2} C(n)}{K_{2} +S(n)}, \\& \xi _{3} (n)= r_{3} K_{3} - \gamma _{2} r _{3} R ( n ) - \mu _{4} G ( n ) +\rho S ( n ) - \frac{P_{3} C ( n )}{K_{3} +R ( n )}. \end{aligned}$$

Solving system (1.2), we get, for $t\rightarrow n+1$,

$$ \textstyle\begin{cases} G ( n+1 ) = \frac{G ( n ) \xi _{1} (n)}{ ( \xi _{1} (n)- \alpha _{1} r_{1} G(n) ) e^{- \xi _{1} (n)} + \alpha _{1} r_{1} G ( n )}, \\ S ( n+1 ) = \frac{S ( n ) \xi _{2} (n)}{ ( \xi _{2} (n)- \beta _{1} r_{2} S ( n ) ) e^{- \xi _{2} (n)} + \beta _{1} r_{2} S ( n )}, \\ R ( n+1 ) = \frac{R ( n ) \xi _{3} (n)}{ ( \xi _{3} (n)- \gamma _{1} r_{3} R ( n ) ) e^{- \xi _{3} (n)} + \gamma _{1} r_{3} R ( n )}, \\ C ( n+1 ) = \frac{\sigma ( 1- e^{- \omega _{1}} )}{ \omega _{1}} + \frac{ ( ( \omega _{1} + \omega _{2} ) e ^{- \omega _{1}} - \omega _{2} ) C(n)}{\omega _{1}}. \end{cases} $$

(1.3)

Let $( \overline{G}, \overline{S}, \overline{R}, \overline{C} )$ be an equilibrium point of (1.3). Then, if we assume that $\overline{\xi _{i}} \neq 0$, $i=1,2,3$, then the fixed points G̅, S̅, and R̅ satisfy

$$ \overline{G} ( \overline{\xi _{1}} - \alpha _{1} r_{1} \overline{G} ) =0, \qquad \overline{S} ( \overline{\xi _{2}} - \beta _{1} r_{2} \overline{ S} ) =0, \qquad \overline{R} ( \overline{\xi _{3}} - \gamma _{1} r_{3} \overline{ R} ) =0, \quad \mbox{and}\quad \overline{C} = \frac{\sigma }{\omega _{1} + \omega _{2}}. $$

Henceforth, we consider $\varLambda _{1} = ( \overline{G}, \overline{S}, \overline{R}, \overline{C} )$ of (1.3) to analyze the interaction among the cells.

This paper organized as follows: In Sect. 2, we analyze the boundedness character of the solutions of (1.3) and obtain conditions for a single semi-cycle behavior. Moreover, we prove that the growth of cancer during the treatment process does not show a period two behavior. In Sect. 3, we discuss the local and global stability of the positive equilibrium point of (1.3) by applying the Schur–Cohn criteria and considering a suitable Lyapunov function. By using the center manifold theorem and bifurcation theory, we show in Sect. 4 that the model undergoes Neimark–Sacker bifurcation. In Sect. 5, we discuss the case for a low-density size of the tumor cell population. This scenario is designed to represent the early detection case, where we focused on the strong Allee effect study to consider the threshold effect. Finally, in Sect. 6 we provide our remarks and results using simulation.

2 Boundedness character and the periodic behavior

In this section, we analyze the boundedness character of the system and obtain conditions for possible periodic behavior.

Theorem 2.1

Assume that $\{ G ( n ),S ( n ), R ( n ), C(n) \} _{n=0}^{ \infty }$ is a positive solution of system (1.3). Then the following statements are true.

• (i) If

  $$ \textstyle\begin{cases} r_{1} K_{1} - \alpha _{2} r_{1} G ( n ) - \mu _{1} S ( n ) - \mu _{3} R ( n ) - \frac{P _{1} C(n)}{K_{1} +G(n)} > \alpha _{1} r_{1} G ( n ), \\ p+ r_{2} K_{2} - \beta _{2} r_{2} S ( n ) - \mu _{2} G ( n ) -\rho R ( n ) - \frac{P _{2} C(n)}{K_{2} +S(n)} > \beta _{1} r_{2} S ( n ), \\ r_{3} K_{3} - \gamma _{2} r_{3} R ( n ) - \mu _{4} G ( n ) +\rho S ( n ) - \frac{P_{3} C ( n )}{K_{3} +R ( n )} > \gamma _{1} r_{3} R ( n ), \\ \frac{\sigma }{\omega _{1} + \omega _{2}} >C ( n ), \end{cases} $$

  (2.1)

  then the solution of system (1.3) increases monotonically.

• (ii) If

  $$ \textstyle\begin{cases} \alpha _{1} r_{1} G ( n ) > r _{1} K_{1} - \alpha _{2} r_{1} G ( n ) - \mu _{1} S ( n ) - \mu _{3} R ( n ) - \frac{P_{1} C(n)}{K _{1} +G(n)} >0, \\ p+ r_{2} K_{2} - \beta _{2} r_{2} S ( n ) - \mu _{2} G ( n ) -\rho R ( n ) - \frac{P _{2} C(n)}{K_{2} +S(n)} > \xi _{2} >0, \\ r_{3} K_{3} - \gamma _{2} r_{3} R ( n ) - \mu _{4} G ( n ) +\rho S ( n ) - \frac{P_{3} C ( n )}{K_{3} +R ( n )} > \xi _{3} >0, \\ C ( n ) > \frac{\sigma }{\omega _{1} + \omega _{2}}, \end{cases} $$

  (2.2)

  then the solution of system (1.3) decreases monotonically.

Proof

(i) Assume that $\{ G ( n ),S ( n ), R ( n ),C(n) \} _{n=0}^{\infty }$ is a positive solution of system (1.3). Then we have

$$ \textstyle\begin{cases} \frac{G ( n+1 )}{G ( n )} = \frac{\xi _{1} (n)}{ ( \xi _{1} (n)- \alpha _{1} r_{1} G(n) ) e^{- \xi _{1} (n)} + \alpha _{1} r_{1} G ( n )} >1, \\ \frac{S ( n+1 )}{S ( n )} = \frac{\xi _{2} (n)}{ ( \xi _{2} (n)- \beta _{1} r_{2} S ( n ) ) e^{- \xi _{2} (n)} + \beta _{1} r_{2} S ( n )} >1, \\ \frac{R ( n+1 )}{R ( n )} = \frac{\xi _{3} (n)}{ ( \xi _{3} (n)- \gamma _{1} r_{3} R ( n ) ) e^{- \xi _{3} (n)} + \gamma _{1} r_{3} R ( n )} >1, \\ \frac{C ( n+1 )}{C(n)} = \frac{\sigma }{\omega _{1} C(n)} ( 1- e^{- \omega _{1}} ) + \frac{ ( \omega _{1} + \omega _{2} ) e^{- \omega _{1}} - \omega _{2}}{\omega _{1}} >1 \end{cases} $$

(2.3)

if (2.1) holds.

(ii) The proof is similar to Theorem 2.1(i) and is omitted. □

Theorem 2.2

Assume that $\{ G ( n ),S ( n ), R ( n ), C(n) \} _{n=0}^{ \infty }$ is a positive solution of system (1.3). Then the following statements are true.

• (i) Let Theorem 2.1(i) hold. Then

  $$ \begin{gathered} G ( n ) < \frac{K_{1}}{\alpha _{1} ( 1- e ^{- r_{1} K_{1}} )},\qquad S ( n ) < \frac{ K_{2}}{\beta _{1} ( 1- e^{- r_{2} K _{2}} )}, \\ R ( n ) < \frac{K_{3}}{\gamma _{1} ( 1- e^{- r_{3} K_{3}} )},\quad \textit{and}\quad C ( n ) < \frac{\sigma ( 1+ e^{- \omega _{1}} )}{\omega _{1}}. \end{gathered} $$

• (ii) Let Theorem 2.1(ii) hold. If $\xi _{i} ( n ) > \ln 2$ (for $i=1, 2, 3$), then system (1.3) is upper bounded with

  $$\begin{aligned}& G ( n ) < G ( 0 ) \exp \Biggl( \alpha _{1} r_{1} \sum _{i=0}^{n-1} G ( i ) \Biggr), \\& S ( n ) < S ( 0 ) \exp \Biggl( \beta _{1} r_{2} \sum _{i=0}^{n-1} S ( i ) \Biggr), \\& R ( n ) < R ( 0 ) \exp \Biggl( \gamma _{1} r_{3} \sum _{i=0}^{n-1} R ( i ) \Biggr), \end{aligned}$$

  and

  $$ C ( n ) < \biggl( \frac{\omega _{1} + \omega _{2}}{\omega _{1}} \biggr)^{n} C ( 0 ) + \frac{\sigma }{\omega _{1}} \sum_{i=1}^{n} \biggl( \frac{\omega _{1} + \omega _{2}}{\omega _{1}} \biggr) ^{i-1}. $$

Proof

(i) Assume that the conditions in Theorem 2.1(i) hold. Since

$$ r_{1} K_{1} > r_{1} K_{1} - \alpha _{2} r_{1} G ( n ) - \mu _{1} S ( n ) - \mu _{3} R ( n ) - \frac{P_{1} C(n)}{K_{1} +G(n)} > \alpha _{1} r _{1} G ( n ), $$

we obtain

$$\begin{aligned} G ( n+1 ) &= \frac{\xi _{1} ( n ) G ( n )}{ ( \xi _{1} ( n ) - \alpha _{1} r _{1} G ( n ) ) e^{- \xi _{1} ( n )} + \alpha _{1} r_{1} G ( n )} \\ &< \frac{r_{1} K_{1} G ( n )}{ ( - \alpha _{1} r_{1} G ( n ) ) e ^{- \xi _{1}} + \alpha _{1} r_{1} G ( n )} \\ &< \frac{r_{1} K_{1} G(n)}{ ( - \alpha _{1} r_{1} G(n) ) e^{- r_{1} K_{1}} + \alpha _{1} r_{1} G(n)} < \frac{K_{1}}{\alpha _{1} ( 1- e^{- r_{1} K_{1}} )}. \end{aligned}$$

(2.4)

Similarly to (2.4), we have

$$ S ( n+1 ) < \frac{K_{2}}{\beta _{1} ( 1- e ^{- r_{2} K_{2}} )}\quad \text{and}\quad R ( n+1 ) < \frac{K_{3}}{\gamma _{1} ( 1- e^{- r_{3} K_{3}} )}. $$

Moreover, from (2.1), we get

$$ C ( n+1 ) < \frac{\sigma ( 1+ e^{- \omega _{1}} )}{ \omega _{1}}, $$

which completes the proof of part (i).

(ii) Assume that Theorem 2.1(ii) holds. Then we obtain

$$\begin{aligned} G ( n+1 ) &= \frac{\xi _{1} ( n ) G ( n )}{ ( \xi _{1} ( n ) - \alpha _{1} r _{1} G ( n ) ) e^{- \xi _{1} ( n )} + \alpha _{1} r_{1} G ( n )} \\ &< \frac{\alpha _{1} r_{1} G ( n )}{ \alpha _{1} r_{1} ( 1- e^{- \xi _{1} ( n )} )} = \frac{G(n) e^{\xi _{1} ( n )}}{ e^{\xi _{1} ( n )} -1} < \frac{G(n) e ^{\alpha _{1} r_{1} G ( n )}}{ e^{\xi _{1} ( n )} -1} < G ( n ) e^{\alpha _{1} r_{1} G ( n )}. \end{aligned}$$

For $n=1, 2, \dots$, we get

$$ G ( n ) < G ( 0 ) \exp \Biggl( \alpha _{1} r_{1} \sum _{i=0}^{n-1} G ( i ) \Biggr). $$

(2.5)

Similarly, we obtain

$$ S ( n ) < S ( 0 ) \exp \Biggl( \beta _{1} r_{2} \sum _{i=0}^{n-1} S ( i ) \Biggr)\quad \text{and} \quad R ( n ) < R ( 0 ) \exp \Biggl( \gamma _{1} r _{3} \sum _{i=0}^{n-1} R ( i ) \Biggr). $$

At last, from

$$ C ( n+1 ) = \frac{\sigma ( 1- e^{- \omega _{1}} ) + ( ( \omega _{1} + \omega _{2} ) e^{- \omega _{1}} - \omega _{2} ) C(n)}{\omega _{1}} < \frac{\sigma + ( \omega _{1} + \omega _{2} ) C(n)}{\omega _{1}}. $$

Since $\frac{\omega _{1} + \omega _{2}}{\omega _{1}} >1$, we have for $n=1, 2, 3,\dots$ the upper bound

$$ C ( n ) < \biggl( \frac{\omega _{1} + \omega _{2}}{\omega _{1}} \biggr)^{n} C ( 0 ) + \frac{\sigma }{\omega _{1}} \sum_{i=1}^{n} \biggl( \frac{\omega _{1} + \omega _{2}}{\omega _{1}} \biggr) ^{i-1}, $$

(2.6)

which completes the proof of part (ii). □

Theorem 2.3

Let $\{ G ( n ),S ( n ), R ( n ), C(n) \} _{n=0}^{\infty }$ be a positive solution of system (1.3), which consists of a single semi-cycle. Furthermore, assume that Theorem 2.1 holds. If

$$ \begin{aligned} &\xi _{1} ( n-1 ) < \ln \biggl( \frac{ \overline{G}}{G ( n-1 )} \biggr),\qquad \xi _{2} ( n-1 ) < \ln \biggl( \frac{\overline{ S}}{S ( n-1 )} \biggr)\quad \textit{and} \\ &\xi _{3} ( n-1 ) < \ln \biggl( \frac{\overline{ R}}{R ( n-1 )} \biggr) \end{aligned} $$

(2.7)

for all $n\geq 1$, then $\{ G ( n ),S ( n ), R ( n ), C(n) \} _{n=0}^{\infty }$ converges monotonically to the positive equilibrium point $\varLambda _{1} = ( \overline{G}, \overline{S}, \overline{R}, \overline{C} )$.

Proof

Suppose that $0< G ( n-1 ) < \overline{G}$ for all $n\geq 1$. In this case, we have

$$ G ( n-1 ) < G ( n ) = \frac{G ( n-1 ) \xi _{1} (n-1)}{ ( \xi _{1} (n-1)- \alpha _{1} r_{1} G(n-1) ) e^{- \xi _{1} (n-1)} + \alpha _{1} r_{1} G ( n-1 )}, $$

(2.8)

where

$$ \bigl( 1- e^{- \xi _{1} (n-1)} \bigr) \cdot \bigl( \xi _{1} (n-1)- \alpha _{1} r_{1} G( n-1) \bigr) >0 $$

holds for the conditions in (2.1). Moreover, from

$$ G ( n ) = \frac{G ( n-1 ) \xi _{1} (n-1)}{ ( \xi _{1} (n-1)- \alpha _{1} r_{1} G(n-1) ) e^{- \xi _{1} (n-1)} + \alpha _{1} r_{1} G ( n-1 )} < \overline{G}, $$

we can write

$$ G ( n-1 ) \bigl( \xi _{1} ( n-1 ) - \alpha _{1} r_{1} \overline{G} \bigr) < \overline{G} \bigl( \xi _{1} ( n-1 ) - \alpha _{1} r_{1} G ( n-1 ) \bigr) e^{- \xi _{1} ( n-1 )}, $$

(2.9)

which holds for

$$ \xi _{1} ( n-1 ) - \alpha _{1} r_{1} \overline{ G} < \xi _{1} ( n-1 ) - \alpha _{1} r _{1} G ( n-1 ) $$

and

$$ G ( n-1 ) < \overline{G} \cdot e^{- \xi _{1} ( n-1 )} \quad \Rightarrow \quad \xi _{1} ( n-1 ) < \ln \biggl( \frac{\overline{ G}}{G ( n-1 )} \biggr). $$

(2.10)

Thus we obtain

$$ 0< G ( n-1 ) < G ( n ) < \overline{G}\quad \text{for all } n\geq 1. $$

Similarly, we can prove that

$$ 0< S ( n-1 ) < S ( n ) < \overline{S}\quad \text{for all } n \geq 1 $$

and

$$ 0< R ( n-1 ) < R ( n ) < \overline{R}\quad \text{for all } n\geq 1. $$

Let us consider $0< N ( n-1 ) < \overline{ N}$ for all $n\geq 1$. Thus, we have

$$ N ( n-1 ) < N ( n ) = N ( n-1 ) e^{\theta - \frac{P_{4} C ( n-1 )}{K_{4} +N ( n-1 )}}. $$

(2.11)

Finally, assume that $0< C ( n-1 ) < \overline{ C}$ for all $n\geq 1$. In this case, we have

$$ C ( n-1 ) < C ( n ) = \frac{\sigma - \omega _{2} C ( n-1 )}{\omega _{1}} + \biggl( \frac{ ( \omega _{1} + \omega _{2} ) C ( n-1 ) -\sigma }{ \omega _{1}} \biggr) e^{- \omega _{1}}, $$

(2.12)

where we obtain

$$ \bigl( ( \omega _{1} + \omega _{2} ) C ( n-1 ) -\sigma \bigr) \bigl( 1- e^{- \omega _{1}} \bigr) < 0. $$

(2.13)

From (2.1), since $\frac{\sigma }{\omega _{1} + \omega _{2}} >C ( n )$, the above inequality (2.13) holds. The proof is completed. □

Theorem 2.4

Let $\{ G ( n ),S ( n ), R ( n ), C(n) \} _{n=0}^{\infty }$ be a positive solution of system (1.3). If

$$ \begin{aligned} &G ( n ) < \sqrt{\frac{P_{1} C ( n )}{ ( \alpha _{1} + \alpha _{2} ) r_{1}}} - K_{1},\qquad S ( n ) < \sqrt{ \frac{P_{2} C ( n )}{ ( \beta _{1} + \beta _{2} ) r_{2}}} - K_{2},\quad \textit{and} \\ & R ( n ) < \sqrt{ \frac{P _{3} C ( n )}{ ( \gamma _{1} + \gamma _{2} ) r_{3}}} - K_{3}, \end{aligned} $$

(2.14)

then system (1.3) has no positive solutions of prime period two.

Proof

Let

$$ \dots , \varpi , \psi , \varpi , \psi , \dots $$

(2.15)

be period two solutions to the $( G ( n ) )_{n=0}^{\infty }$ class such that $\varpi \neq \psi$. Then we have

$$ \begin{aligned}[b] \varpi &= \biggl({\psi \biggl( r_{1} K_{1} - \alpha _{2} r_{1} \psi - \mu _{1} S ( n ) - \mu _{3} R ( n ) - \frac{P_{1} C(n)}{K _{1} +\psi } \biggr)} \biggr) \\ &\quad {} \Big/ \biggl( \biggl( r_{1} K_{1} - \alpha _{2} r_{1} \psi - \mu _{1} S ( n ) - \mu _{3} R ( n ) - \frac{P _{1} C(n)}{K_{1} +\psi } - \alpha _{1} r_{1} \psi \biggr) \\ &\quad {} \cdot e^{- ( r_{1} K_{1} - \alpha _{2} r_{1} \psi - \mu _{1} S ( n ) - \mu _{3} R ( n ) - \frac{P_{1} C(n)}{K_{1} +\psi } )} + \alpha _{1} r_{1} \psi \biggr) \end{aligned} $$

(2.16)

and

$$ \begin{aligned}[b] \psi &= \biggl({\varpi \biggl( r_{1} K_{1} - \alpha _{2} r_{1} \varpi - \mu _{1} S ( n ) - \mu _{3} R ( n ) - \frac{P_{1} C(n)}{K _{1} +\varpi } \biggr)} \biggr) \\ &\quad {} \Big/ \biggl( \biggl( r_{1} K_{1} - \alpha _{2} r_{1} \varpi - \mu _{1} S ( n ) - \mu _{3} R ( n ) - \frac{P_{1} C(n)}{K_{1} +\varpi } - \alpha _{1} r_{1} \varpi \biggr) \\ &\quad {}\cdot e^{- ( r_{1} K_{1} - \alpha _{2} r_{1} \varpi - \mu _{1} S ( n ) - \mu _{3} R ( n ) - \frac{P_{1} C(n)}{K_{1} +\varpi } )} + \alpha _{1} r _{1} \varpi \biggr). \end{aligned} $$

(2.17)

Note that

$$ \begin{gathered}[b] \psi \biggl( r_{1} K_{1} - \alpha _{2} r _{1} \psi - \mu _{1} S ( n ) - \mu _{3} R ( n ) - \frac{P_{1} C(n)}{K_{1} + \psi } \biggr) \\ \qquad {}-\varpi \biggl( r_{1} K_{1} - \alpha _{2} r_{1} \psi - \mu _{1} S ( n ) - \mu _{3} R ( n )- \frac{P _{1} C(n)}{K_{1} +\psi } - \alpha _{1} r_{1} \psi \biggr) \\ \qquad {}\cdot e^{- ( r_{1} K_{1} - \alpha _{2} r_{1} \psi - \mu _{1} S ( n ) - \mu _{3} R ( n ) - \frac{P_{1} C(n)}{K_{1} +\psi } )} \\ \quad = \alpha _{1} r_{1} \varpi \psi \\ \quad = \varpi \biggl( r_{1} K_{1} - \alpha _{2} r_{1} \varpi - \mu _{1} S ( n ) - \mu _{3} R ( n ) - \frac{P_{1} C(n)}{K _{1} +\varpi } \biggr) \\ \qquad {}-\psi \biggl( r_{1} K_{1} - \alpha _{2} r_{1} \varpi - \mu _{1} S ( n ) - \mu _{3} R ( n ) - \frac{P_{1} C ( n )}{K_{1} +\varpi } - \alpha _{1} r_{1} \varpi \biggr) \\ \qquad {} \cdot e^{- ( r_{1} K_{1} - \alpha _{2} r_{1} \varpi - \mu _{1} S ( n ) - \mu _{3} R ( n ) - \frac{P_{1} C ( n )}{K_{1} +\varpi } )}. \end{gathered} $$

(2.18)

Let us set

$$\begin{aligned} F ( z ) &=z \biggl( r_{1} K_{1} - \alpha _{2} r_{1} z - \mu _{1} S ( n ) - \mu _{3} R ( n ) - \frac{P_{1} C ( n )}{K_{1} +z} - \alpha _{1} r_{1} z \biggr) \\ &\quad {}\cdot \bigl( 1- e^{- ( r_{1} K _{1} - \alpha _{2} r_{1} z - \mu _{1} S ( n ) - \mu _{3} R ( n ) - \frac{P_{1} C ( n )}{K_{1} +z} )} \bigr). \end{aligned}$$

It is easy to see that for $\varLambda _{1}$ we have $F ( \overline{G} ) =0$. Since

$$ \begin{aligned}[b] \frac{dF ( z )}{dz} &= \biggl( \bigl( \xi _{1} ( z ) - \alpha _{1} r_{1} z \bigr) + \biggl( \frac{P_{1} C ( n ) ( K_{1} +z-1 )}{ ( K_{1} +z ) ^{2}} - \alpha _{2} r_{1} z- \alpha _{1} r_{1} z \biggr) \biggr) \cdot \bigl( 1- e^{- \xi _{1} ( z )} \bigr) \\ &\quad {}+ \biggl( \frac{P_{1} C ( n )}{ ( K_{1} +z )^{2}} - \alpha _{2} r_{1} \biggr) e ^{- \xi _{1} ( z )} \cdot z \bigl( \xi _{1} ( z ) - \alpha _{1} r_{1} z \bigr) \\ &>0 \end{aligned} $$

(2.19)

if

$$ G ( n ) =z< \sqrt{\frac{P_{1} C ( n )}{ ( \alpha _{1} + \alpha _{2} ) r_{1}}} - K_{1} $$

(2.20)

since $1- e^{- \xi _{1} ( z )} >0$ and $\xi _{1} ( z ) - \alpha _{1} r_{1} z>0$. Similarly, we obtain

$$ S ( n ) < \sqrt{\frac{P_{2} C ( n )}{ ( \beta _{1} + \beta _{2} ) r_{2}}} - K_{2}\quad \text{and}\quad R ( n ) < \sqrt{\frac{P _{3} C ( n )}{ ( \gamma _{1} + \gamma _{2} ) r_{3}}} - K_{3}. $$

(2.21)

Finally, from

$$ \textstyle\begin{cases} \nu = \frac{\sigma ( 1- e^{- \omega _{1}} )}{\omega _{1}} + ( \frac{ ( \omega _{1} + \omega _{2} ) e^{- \omega _{1}} - \omega _{2}}{\omega _{1}} ) \hat{\nu }, \\ \hat{\nu } = \frac{\sigma ( 1- e^{- \omega _{1}} )}{\omega _{1}} + ( \frac{ ( \omega _{1} + \omega _{2} ) e^{- \omega _{1}} - \omega _{2}}{\omega _{1}} ) \nu, \end{cases} $$

(2.22)

where $\nu \neq \hat{\nu }$, we get

$$ ( \nu - \hat{\nu } ) \biggl( \frac{ ( \omega _{1} - \omega _{2} ) + ( \omega _{1} + \omega _{2} ) e^{- \omega _{1}}}{\omega _{1}} \biggr) =0, $$

(2.23)

which holds for $\nu = \hat{\nu }$. This contradicts the assumption. The proof is completed. □

3 Local and global stability analysis

In this section, we use the Schur–Cohn criterion to show the local asymptotic stability of the positive equilibrium point.

The Jacobian matrix for the positive equilibrium point $\varLambda _{1} =( \overline{G}, \overline{S}, \overline{R}, \overline{C})$ is

$$ J ( \varLambda _{1} ) = \begin{pmatrix} a_{11} & a_{12}&a_{13} & a_{14}\\ a_{21} & a_{22}&a_{23} & a_{24}\\ a_{31} & a_{32}&a_{33} & a_{34}\\ 0 & 0& 0 & a_{44} \end{pmatrix}, $$

(3.1)

where

$$\begin{aligned}& a_{11} = \frac{ ( P_{1} \overline{C} - \alpha _{2} r_{1} ( K_{1} + \overline{G} )^{2} ) + ( ( \alpha _{1} + \alpha _{2} ) r_{1} ( K_{1} + \overline{G} )^{2} - P_{1} \overline{C} ) e ^{- \overline{\xi _{1}}}}{\alpha _{1} r_{1} ( K_{1} + \overline{G} )^{2}}, \\& a_{12} = \frac{\mu _{1} ( e ^{- \overline{\xi _{1}}} -1 )}{\alpha _{1} r_{1}}, \qquad a_{13} = \frac{\mu _{3} ( e^{- \overline{\xi _{1}}} -1 )}{ \alpha _{1} r_{1}}, \\& a_{14} =- \frac{P_{1} ( 1- e^{- \overline{ \xi _{1}}} )}{\alpha _{1} r_{1} ( K_{1} + \overline{G} )}, \qquad a_{21} = \frac{\mu _{2} ( e^{- \overline{\xi _{2}}} -1 )}{ \beta _{1} r_{2}}, \\& a_{22} = \frac{ ( P_{2} \overline{C} - \beta _{2} r_{2} ( K_{2} + \overline{S} )^{2} ) + ( ( \beta _{1} + \beta _{2} ) r_{2} ( K_{2} + \overline{S} )^{2} - P_{2} \overline{C} ) e ^{- \overline{\xi _{2}}}}{\beta _{1} r_{2} ( K_{2} + \overline{S} )^{2}}, \\& a_{23} = \frac{\rho ( e^{- \overline{\xi _{2}}} -1 )}{\beta _{1} r_{2}}, \qquad a_{24} =- \frac{P_{2} ( 1- e^{- \overline{\xi _{2}}} )}{ \beta _{1} r_{2} ( K_{2} + \overline{S} )}, \\& a_{31} = \frac{\mu _{4} ( e^{- \overline{\xi _{3}}} -1 )}{ \gamma _{1} r_{3}},\qquad a_{32} = \frac{\rho ( 1- e ^{- \overline{\xi _{3}}} )}{\gamma _{1} r_{3}}, \\& a_{33} = \frac{ ( P_{3} \overline{C} - \gamma _{2} r_{3} ( K_{3} + \overline{R} )^{2} ) + ( ( \gamma _{1} + \gamma _{2} ) r_{3} ( K_{3} + \overline{R} )^{2} - P_{3} \overline{C} ) e ^{- \overline{\xi _{3}}}}{\gamma _{1} r_{3} ( K_{3} + \overline{R} )^{2}}, \\& a_{34} =- \frac{P_{3} ( 1- e^{- \overline{\xi _{3}}} )}{ \gamma _{1} r_{3} ( K_{3} + \overline{R} )}, \end{aligned}$$

and

$$a_{44} = \frac{- \omega _{2} + ( \omega _{1} + \omega _{2} ) e^{- \omega _{1}}}{\omega _{1}}. $$

The characteristic equation around the positive equilibrium point is given by

$$ \lambda _{1} = \frac{- \omega _{2} + ( \omega _{1} + \omega _{2} ) e^{- \omega _{1}}}{\omega _{1}} $$

(3.2)

and

$$ \lambda ^{3} + a_{2} \lambda ^{2} + a_{1} \lambda + a_{0} =0, $$

(3.3)

where

$$\begin{aligned}& a_{2} =- a_{11} - a_{22} - a_{33}, \\& a_{1} = a_{11} a_{22} + a_{11} a_{33} + a_{22} a_{33} - a _{13} a_{31} - a_{23} a_{32} - a_{12} a_{21}, \\& a_{0} = a_{13} a_{31} a_{22} + a_{23} a_{32} a_{11} + a_{12} a_{21} a_{33} - a_{11} a_{22} a_{33} - a_{21} a_{32} a_{13} - a_{31} a _{12} a_{23}. \end{aligned}$$

From (3.2) $\lambda _{1}$ is always inside the unit disk since $\omega _{1} > 0$.

Theorem 3.1

([29])

The characteristic polynomial

$$ P ( \lambda ) = \lambda ^{3} + a_{2} \lambda ^{2} + a_{1} \lambda + a_{0} $$

(3.4)

has all its roots inside the unit disk $( \vert \lambda \vert <1 )$ if and only if

• (i) $P ( 1 ) =1+ a_{2} + a_{1} + a_{0} >0$ and $( -1 )^{3} P ( -1 ) =1- a_{2} + a_{1} - a_{0} >0$;

• (ii) $D_{2}^{+} =1+ a_{1} - a_{0} ^{2} - a_{0} a_{2} >0$;

• (iii) $D_{2}^{-} =1- a_{1} - a_{0} ^{2} + a_{0} a_{2} >0$.

Theorem 3.2

Let $\varLambda _{1} =( \overline{G}, \overline{S}, \overline{R}, \overline{C})$ be a positive equilibrium point of system (1.3). Assume that

$$ \begin{gathered} \frac{\mu _{1} \mu _{2}}{\mu _{3} \mu _{4}} > \frac{ ( P_{2} \overline{C} - ( \beta _{1} + \beta _{2} ) r_{2} ( K_{2} + \overline{S} )^{2} ) ( K_{3} + \overline{R} )^{2}}{ ( P_{3} \overline{C} - ( \gamma _{1} + \gamma _{2} ) r_{3} ( K_{3} + \overline{R} )^{2} ) ( K_{2} + \overline{S} ) ^{2}}, \\ \frac{\mu _{1}}{\mu _{2}} < \frac{\rho ^{2} ( P_{1} \overline{C} - \alpha _{2} r_{1} ( K_{1} + \overline{G} )^{2} ) ( K_{2} + \overline{S} )^{2}}{ ( P_{2} \overline{C} - \beta _{2} r_{2} ( K_{2} + \overline{S} )^{2} ) ( K_{1} + \overline{G} ) ^{2}},\end{gathered} $$

(3.5)

where $\overline{C} > \frac{ ( \alpha _{1} + \alpha _{2} ) r_{1} ( K_{1} + \overline{G} )^{2}}{P_{1}} > \frac{ ( \gamma _{1} + \gamma _{2} ) r_{3} ( K_{3} + \overline{R} )^{2}}{P_{3}} > \frac{ ( \beta _{1} + \beta _{2} ) r_{2} ( K_{2} + \overline{S} ) ^{2}}{P_{2}}$. If

$$ r_{1} > \frac{1}{\alpha _{1}} \biggl( \frac{\mu _{3} \mu _{4}}{\gamma _{1} r_{3}} + \frac{\mu _{1} \mu _{2}}{\beta _{1} r_{2}} \biggr) > \frac{ ( P_{3} \overline{C} - \gamma _{2} r_{3} ( K_{3} + \overline{R} )^{2} ) ( P_{1} \overline{C} - ( \alpha _{1} + \alpha _{2} ) r_{1} ( K _{1} + \overline{G} )^{2} )}{\alpha _{1} \gamma _{1} r_{3} ( K_{1} + \overline{G} )^{2} ( K _{3} + \overline{R} )^{2}}, $$

(3.6)

and the conditions

$$ \overline{\xi _{2}} < \ln \biggl( \frac{ ( P_{2} \overline{C} - ( \beta _{1} + \beta _{2} ) r_{2} ( K_{2} + \overline{S} )^{2} ) - \mu _{2} \mu _{4} \rho ( K _{2} + \overline{S} )^{2}}{ ( P_{2} \overline{C} - \beta _{2} r_{2} ( K_{2} + \overline{S} )^{2} ) - \mu _{2} \mu _{4} \rho ( K_{2} + \overline{S} )^{2}} \biggr) $$

and

$$ \overline{\xi _{1}} > \ln \biggl( \frac{\rho ( P_{1} \overline{C} - ( \alpha _{1} + \alpha _{2} ) r _{1} ( K_{1} + \overline{G} )^{2} ) - \mu _{1} \mu _{4} ( K_{1} + \overline{G} )^{2}}{\rho ( P_{1} \overline{C} - \alpha _{2} r_{1} ( K_{1} + \overline{G} )^{2} ) - \mu _{1} \mu _{4} ( K_{1} + \overline{G} ) ^{2}} \biggr) $$

hold, then the positive equilibrium point is locally asymptotically stable.

Proof

We want to consider at first Theorem 3.1(i), where we have

$$ \begin{aligned}[b] &1+ ( - a_{11} - a_{22} - a_{33} ) + ( a_{11} a_{22} + a_{11} a_{33} + a_{22} a_{33} - a_{13} a_{31} - a_{23} a_{32} - a _{12} a_{21} ) \\ &\qquad {}+ ( a_{13} a_{31} a_{22} + a_{23} a_{32} a _{11} + a_{12} a_{21} a_{33} - a_{11} a_{22} a_{33} - a_{21} a_{32} a _{13} - a_{31} a_{12} a_{23} ) \\ &\quad >0 \end{aligned} $$

(3.7)

and

$$ \begin{aligned}[b] &1- ( - a_{11} - a_{22} - a_{33} ) + ( a_{11} a_{22} + a_{11} a_{33} + a_{22} a_{33} - a_{13} a_{31} - a_{23} a_{32} - a _{12} a_{21} ) \\ &\qquad {} - ( a_{13} a_{31} a_{22} + a_{23} a_{32} a _{11} + a_{12} a_{21} a_{33} - a_{11} a_{22} a_{33} - a_{21} a_{32} a _{13} - a_{31} a_{12} a_{23} ) \\ &\quad >0. \end{aligned} $$

(3.8)

Considering (3.7) and (3.8) together, we obtain

$$ 1+ a_{11} a_{22} + a_{11} a_{33} + a_{22} a_{33} > a_{13} a_{31} + a _{23} a_{32} + a_{12} a_{21}. $$

(3.9)

Inequality (3.9) can be written in the following form:

$$\begin{aligned}& \begin{gathered} 1+ \frac{ ( P_{1} \overline{C} - \alpha _{2} r_{1} ( K_{1} + \overline{G} )^{2} ) - ( P_{1} \overline{C} - ( \alpha _{1} + \alpha _{2} ) r _{1} ( K_{1} + \overline{G} )^{2} ) e ^{- \overline{\xi _{1}}}}{\alpha _{1} r_{1} ( K_{1} + \overline{G} )^{2}} \\ \qquad {}\cdot \frac{ ( P_{2} \overline{C} - \beta _{2} r_{2} ( K_{2} + \overline{S} ) ^{2} ) - ( P_{2} \overline{C} - ( \beta _{1} + \beta _{2} ) r_{2} ( K_{2} + \overline{S} ) ^{2} ) e^{- \overline{\xi _{2}}}}{\beta _{1} r_{2} ( K_{2} + \overline{S} )^{2}} \end{gathered} \\& \qquad {}+ \frac{ ( P_{1} \overline{C} - \alpha _{2} r_{1} ( K_{1} + \overline{G} )^{2} ) - ( P_{1} \overline{C} - ( \alpha _{1} + \alpha _{2} ) r _{1} ( K_{1} + \overline{G} )^{2} ) e ^{- \overline{\xi _{1}}}}{\alpha _{1} r_{1} ( K_{1} + \overline{G} )^{2}} \\& \qquad {}\cdot \frac{ ( P_{3} \overline{C} - \gamma _{2} r_{3} ( K_{3} + \overline{R} )^{2} ) - ( P_{3} \overline{C} - ( \gamma _{1} + \gamma _{2} ) r_{3} ( K_{3} + \overline{R} )^{2} ) e^{- \overline{\xi _{3}}}}{ \gamma _{1} r_{3} ( K_{3} + \overline{R} )^{2}} \\& \qquad {}+ \frac{ ( P_{2} \overline{C} - \beta _{2} r_{2} ( K_{2} + \overline{S} )^{2} ) - ( P_{2} \overline{C} - ( \beta _{1} + \beta _{2} ) r_{2} ( K_{2} + \overline{S} )^{2} ) e^{- \overline{ \xi _{2}}}}{\beta _{1} r_{2} ( K_{2} + \overline{S} )^{2}} \\& \qquad {}\cdot \frac{ ( P_{3} \overline{C} - \gamma _{2} r_{3} ( K_{3} + \overline{R} )^{2} ) - ( P_{3} \overline{C} - ( \gamma _{1} + \gamma _{2} ) r_{3} ( K_{3} + \overline{R} )^{2} ) e^{- \overline{\xi _{3}}}}{\gamma _{1} r_{3} ( K _{3} + \overline{R} )^{2}} \\& \qquad {}+ \frac{\rho ^{2} ( 1- e^{- \overline{\xi _{2}}} ) ( 1- e^{- \overline{ \xi _{3}}} )}{\beta _{1} r_{2} \gamma _{1} r_{3}} \\& \quad > \frac{ \mu _{3} \mu _{4} ( 1- e^{- \overline{\xi _{1}}} ) ( 1- e^{- \overline{\xi _{3}}} )}{\alpha _{1} r _{1} \gamma _{1} r_{3}} + \frac{\mu _{1} \mu _{2} ( 1- e^{- \overline{\xi _{1}}} ) ( 1- e^{- \overline{ \xi _{2}}} )}{\alpha _{1} r_{1} \beta _{1} r_{2}}, \end{aligned}$$

(3.10)

which is reordered such as

$$\begin{aligned}& \biggl( 1+ \frac{ ( P_{1} \overline{C} - \alpha _{2} r _{1} ( K_{1} + \overline{G} )^{2} ) ( P_{2} \overline{C} - \beta _{2} r_{2} ( K_{2} + \overline{S} ) ^{2} )}{\alpha _{1} \beta _{1} r_{1} r_{2} ( K_{1} + \overline{G} )^{2} ( K_{2} + \overline{S} )^{2}} \\& \quad {}+ \frac{ ( P_{1} \overline{C} - \alpha _{2} r _{1} ( K_{1} + \overline{G} )^{2} ) ( P_{3} \overline{C} - \gamma _{2} r_{3} ( K_{3} + \overline{R} )^{2} )}{\alpha _{1} \gamma _{1} r_{1} r_{3} ( K_{1} + \overline{G} )^{2} ( K_{3} + \overline{R} )^{2}} \\& \quad {} + \frac{ ( P_{2} \overline{C} - \beta _{2} r_{2} ( K_{2} + \overline{S} )^{2} ) ( P_{3} \overline{C} - \gamma _{2} r_{3} ( K_{3} + \overline{R} )^{2} )}{\beta _{1} \gamma _{1} r_{2} r_{3} ( K_{2} + \overline{S} )^{2} ( K _{3} + \overline{R} )^{2}} - \frac{\mu _{3} \mu _{4}}{\alpha _{1} \gamma _{1} r_{1} r_{3}} - \frac{\mu _{1} \mu _{2}}{ \alpha _{1} \beta _{1} r_{1} r_{2}} \biggr) \\& \quad {}+ \biggl( - \frac{ ( P_{2} \overline{C} - \beta _{2} r _{2} ( K_{2} + \overline{S} )^{2} ) ( P_{1} \overline{C} - ( \alpha _{1} + \alpha _{2} ) r _{1} ( K_{1} + \overline{G} )^{2} )}{\alpha _{1} \beta _{1} r_{1} r_{2} ( K_{1} + \overline{G} )^{2} ( K_{2} + \overline{S} )^{2}} \\& \quad {}- \frac{ ( P_{3} \overline{C} - \gamma _{2} r_{3} ( K_{3} + \overline{R} )^{2} ) ( P_{1} \overline{C} - ( \alpha _{1} + \alpha _{2} ) r_{1} ( K_{1} + \overline{G} )^{2} )}{ \alpha _{1} \gamma _{1} r_{1} r_{3} ( K_{1} + \overline{G} )^{2} ( K_{3} + \overline{R} )^{2}} + \frac{\mu _{3} \mu _{4}}{\alpha _{1} \gamma _{1} r_{1} r_{3}} + \frac{\mu _{1} \mu _{2}}{\alpha _{1} \beta _{1} r_{1} r_{2}} \biggr) e^{- \overline{\xi _{1}}} \\& \quad {}+ \biggl( - \frac{ ( P_{1} \overline{C} - \alpha _{2} r _{1} ( K_{1} + \overline{G} )^{2} ) ( P_{2} \overline{C} - ( \beta _{1} + \beta _{2} ) r_{2} ( K_{2} + \overline{S} )^{2} )}{\alpha _{1} \beta _{1} r_{1} r_{2} ( K_{1} + \overline{G} )^{2} ( K _{2} + \overline{S} )^{2}} \\& \quad {}- \frac{ ( P_{3} \overline{C} - \gamma _{2} r_{3} ( K_{3} + \overline{R} )^{2} ) ( P_{2} \overline{C} - ( \beta _{1} + \beta _{2} ) r_{2} ( K_{2} + \overline{S} )^{2} )}{ \beta _{1} \gamma _{1} r_{2} r_{3} ( K_{2} + \overline{S} )^{2} ( K_{3} + \overline{R} )^{2}} + \frac{\mu _{1} \mu _{2}}{\alpha _{1} \beta _{1} r_{1} r_{2}} \biggr) e^{- \overline{\xi _{2}}} \\& \quad {}+ \biggl( - \frac{ ( P_{1} \overline{C} - \alpha _{2} r _{1} ( K_{1} + \overline{G} )^{2} ) ( P_{3} \overline{C} - ( \gamma _{1} + \gamma _{2} ) r _{3} ( K_{3} + \overline{R} )^{2} )}{\alpha _{1} \gamma _{1} r_{1} r_{3} ( K_{1} + \overline{G} ) ^{2} ( K_{3} + \overline{R} )^{2}} \\& \quad {}- \frac{ ( P_{2} \overline{C} - \beta _{2} r_{2} ( K_{2} + \overline{S} ) ^{2} ) ( P_{3} \overline{C} - ( \gamma _{1} + \gamma _{2} ) r_{3} ( K_{3} + \overline{R} ) ^{2} )}{\beta _{1} \gamma _{1} r_{2} r_{3} ( K _{2} + \overline{S} )^{2} ( K_{3} + \overline{R} ) ^{2}} + \frac{\mu _{3} \mu _{4}}{\alpha _{1} \gamma _{1} r_{1} r_{3}} \biggr) e^{- \overline{\xi _{3}}} \\& \quad {}+ \biggl( \frac{ ( P_{1} \overline{C} - ( \alpha _{1} + \alpha _{2} ) r_{1} ( K_{1} + \overline{G} ) ^{2} ) ( P_{2} \overline{C} - ( \beta _{1} + \beta _{2} ) r_{2} ( K_{2} + \overline{S} ) ^{2} )}{\alpha _{1} \beta _{1} r_{1} r_{2} ( K_{1} + \overline{G} )^{2} ( K_{2} + \overline{S} )^{2}} - \frac{\mu _{1} \mu _{2}}{\alpha _{1} \beta _{1} r_{1} r_{2}} \biggr) e^{- \overline{\xi _{2}}} e^{- \overline{\xi _{1}}} \\& \quad {}+ \biggl( \frac{ ( P_{1} \overline{C} - ( \alpha _{1} + \alpha _{2} ) r_{1} ( K_{1} + \overline{G} ) ^{2} ) ( P_{3} \overline{C} - ( \gamma _{1} + \gamma _{2} ) r_{3} ( K_{3} + \overline{R} ) ^{2} )}{\alpha _{1} \gamma _{1} r_{1} r_{3} ( K _{1} + \overline{G} )^{2} ( K_{3} + \overline{R} ) ^{2}} - \frac{\mu _{3} \mu _{4}}{\alpha _{1} \gamma _{1} r_{1} r _{3}} \biggr) e^{- \overline{\xi _{3}}} e^{- \overline{ \xi _{1}}} \\& \quad {}+ \biggl( \frac{ ( P_{2} \overline{C} - ( \beta _{1} + \beta _{2} ) r_{2} ( K_{2} + \overline{S} ) ^{2} ) ( P_{3} \overline{C} - ( \gamma _{1} + \gamma _{2} ) r_{3} ( K_{3} + \overline{R} ) ^{2} )}{\beta _{1} \gamma _{1} r_{2} r_{3} ( K _{2} + \overline{S} )^{2} ( K_{3} + \overline{R} ) ^{2}} \biggr) e^{- \overline{\xi _{3}}} e^{- \overline{ \xi _{2}}} >0. \end{aligned}$$

(3.11)

Since

$$ \overline{C} > \frac{\alpha _{2} r_{1} ( K_{1} + \overline{G} )^{2}}{P_{1}} > \frac{\gamma _{2} r_{3} ( K_{3} + \overline{R} )^{2}}{P_{3}} > \frac{\beta _{2} r_{2} ( K_{2} + \overline{S} )^{2}}{P_{2}}, $$

(3.12)

we have

$$ \biggl( \frac{ ( P_{2} \overline{C} - ( \beta _{1} + \beta _{2} ) r_{2} ( K_{2} + \overline{S} ) ^{2} ) ( P_{3} \overline{C} - ( \gamma _{1} + \gamma _{2} ) r_{3} ( K_{3} + \overline{R} ) ^{2} )}{\beta _{1} \gamma _{1} r_{2} r_{3} ( K _{2} + \overline{S} )^{2} ( K_{3} + \overline{R} ) ^{2}} \biggr) e^{- \overline{\xi _{3}}} e^{- \overline{ \xi _{2}}} >0. $$

Furthermore, from (3.11), we obtain

$$ \frac{ ( P_{1} \overline{C} - ( \alpha _{1} + \alpha _{2} ) r_{1} ( K_{1} + \overline{G} )^{2} ) ( P_{3} \overline{C} - ( \gamma _{1} + \gamma _{2} ) r_{3} ( K_{3} + \overline{R} )^{2} )}{ \alpha _{1} \gamma _{1} r_{1} r_{3} ( K_{1} + \overline{G} )^{2} ( K_{3} + \overline{R} )^{2}} - \frac{\mu _{3} \mu _{4}}{\alpha _{1} \gamma _{1} r_{1} r_{3}} < 0, $$

if

$$ \mu _{3} \mu _{4} > \frac{ ( P_{1} \overline{C} - ( \alpha _{1} + \alpha _{2} ) r_{1} ( K_{1} + \overline{G} )^{2} ) ( P_{3} \overline{C} - ( \gamma _{1} + \gamma _{2} ) r_{3} ( K_{3} + \overline{R} )^{2} )}{ ( K_{1} + \overline{G} ) ^{2} ( K_{3} + \overline{R} )^{2}}, $$

(3.13)

while from

$$ \frac{ ( P_{1} \overline{C} - ( \alpha _{1} + \alpha _{2} ) r_{1} ( K_{1} + \overline{G} )^{2} ) ( P_{2} \overline{C} - ( \beta _{1} + \beta _{2} ) r_{2} ( K_{2} + \overline{S} )^{2} )}{ \alpha _{1} \beta _{1} r_{1} r_{2} ( K_{1} + \overline{G} ) ^{2} ( K_{2} + \overline{S} )^{2}} - \frac{\mu _{1} \mu _{2}}{\alpha _{1} \beta _{1} r_{1} r_{2}} < 0 $$

we get

$$ \mu _{1} \mu _{2} > \frac{ ( P_{1} \overline{C} - ( \alpha _{1} + \alpha _{2} ) r_{1} ( K_{1} + \overline{G} )^{2} ) ( P_{2} \overline{C} - ( \beta _{1} + \beta _{2} ) r_{2} ( K_{2} + \overline{S} )^{2} )}{ ( K_{1} + \overline{G} ) ^{2} ( K_{2} + \overline{S} )^{2}}. $$

(3.14)

Considering (3.13) and (3.14) together, we have

$$ \frac{\mu _{1} \mu _{2}}{\mu _{3} \mu _{4}} > \frac{ ( P_{2} \overline{C} - ( \beta _{1} + \beta _{2} ) r_{2} ( K_{2} + \overline{S} )^{2} ) ( K_{3} + \overline{R} )^{2}}{ ( P_{3} \overline{C} - ( \gamma _{1} + \gamma _{2} ) r_{3} ( K_{3} + \overline{R} )^{2} ) ( K_{2} + \overline{S} ) ^{2}}. $$

(3.15)

Moreover, from

$$ \frac{\mu _{3} \mu _{4}}{\alpha _{1} \gamma _{1} r_{1} r_{3}} + \frac{ \mu _{1} \mu _{2}}{\alpha _{1} \beta _{1} r_{1} r_{2}} < 1, $$

we have

$$ r_{1} > \frac{1}{\alpha _{1}} \biggl( \frac{\mu _{3} \mu _{4}}{\gamma _{1} r_{3}} + \frac{\mu _{1} \mu _{2}}{\beta _{1} r_{2}} \biggr). $$

(3.16)

Since $\overline{\xi _{i}} >0$ ($i=1, 2, 3$), it is obvious that

$$\begin{aligned}& \begin{gathered} \frac{ ( P_{1} \overline{C} - \alpha _{2} r_{1} ( K _{1} + \overline{G} )^{2} ) ( P_{2} \overline{C} - \beta _{2} r_{2} ( K_{2} + \overline{S} )^{2} )}{ \alpha _{1} \beta _{1} r_{1} r_{2} ( K_{1} + \overline{G} ) ^{2} ( K_{2} + \overline{S} )^{2}} \\ \quad {}- \frac{ ( P_{1} \overline{C} - \alpha _{2} r_{1} ( K_{1} + \overline{G} )^{2} ) ( P_{2} \overline{C} - ( \beta _{1} + \beta _{2} ) r_{2} ( K_{2} + \overline{S} ) ^{2} )}{\alpha _{1} \beta _{1} r_{1} r_{2} ( K_{1} + \overline{G} )^{2} ( K_{2} + \overline{S} )^{2}} e^{- \overline{\xi _{2}}} >0, \end{gathered} \\& \frac{ ( P_{1} \overline{C} - \alpha _{2} r_{1} ( K _{1} + \overline{G} )^{2} ) ( P_{3} \overline{C} - \gamma _{2} r_{3} ( K_{3} + \overline{R} )^{2} )}{ \alpha _{1} \gamma _{1} r_{1} r_{3} ( K_{1} + \overline{G} )^{2} ( K_{3} + \overline{R} )^{2}} \\& \quad {}- \frac{ ( P_{3} \overline{C} - \gamma _{2} r_{3} ( K _{3} + \overline{R} )^{2} ) ( P_{1} \overline{C} - ( \alpha _{1} + \alpha _{2} ) r_{1} ( K_{1} + \overline{G} )^{2} )}{\alpha _{1} \gamma _{1} r_{1} r_{3} ( K_{1} + \overline{G} )^{2} ( K _{3} + \overline{R} )^{2}} e^{- \overline{\xi _{1}}} >0, \end{aligned}$$

and

$$ \begin{gathered} \frac{ ( P_{2} \overline{C} - \beta _{2} r_{2} ( K _{2} + \overline{S} )^{2} ) ( P_{3} \overline{C} - \gamma _{2} r_{3} ( K_{3} + \overline{R} )^{2} )}{ \beta _{1} \gamma _{1} r_{2} r_{3} ( K_{2} + \overline{S} )^{2} ( K_{3} + \overline{R} )^{2}} \\ \quad {}- \frac{ ( P_{2} \overline{C} - \beta _{2} r_{2} ( K _{2} + \overline{S} )^{2} ) ( P_{3} \overline{C} - ( \gamma _{1} + \gamma _{2} ) r_{3} ( K_{3} + \overline{R} )^{2} )}{\beta _{1} \gamma _{1} r_{2} r_{3} ( K_{2} + \overline{S} )^{2} ( K _{3} + \overline{R} )^{2}} e^{- \overline{\xi _{3}}}>0. \end{gathered} $$

Considering (3.11), we have

$$ \begin{gathered}[b] \biggl( \frac{\mu _{1} \mu _{2}}{\alpha _{1} \beta _{1} r_{1} r_{2}} - \frac{ ( P_{1} \overline{C} - \alpha _{2} r_{1} ( K_{1} + \overline{G} )^{2} ) ( P_{2} \overline{C} - ( \beta _{1} + \beta _{2} ) r_{2} ( K_{2} + \overline{S} )^{2} )}{\alpha _{1} \beta _{1} r_{1} r_{2} ( K_{1} + \overline{G} )^{2} ( K_{2} + \overline{S} )^{2}} \biggr) e^{- \overline{\xi _{2}}} \\ \quad > \biggl( \frac{\mu _{1} \mu _{2}}{\alpha _{1} \beta _{1} r_{1} r_{2}} - \frac{ ( P_{1} \overline{C} - ( \alpha _{1} + \alpha _{2} ) r_{1} ( K_{1} + \overline{G} )^{2} ) ( P_{2} \overline{C} - ( \beta _{1} + \beta _{2} ) r_{2} ( K_{2} + \overline{S} )^{2} )}{ \alpha _{1} \beta _{1} r_{1} r_{2} ( K_{1} + \overline{G} ) ^{2} ( K_{2} + \overline{S} )^{2}} \biggr) e ^{- \overline{\xi _{2}}} e^{- \overline{\xi _{1}}}, \end{gathered} $$

(3.17)

which holds for (3.14). Similarly, from (3.13), we have

$$ \begin{gathered}[b] \biggl( \frac{\mu _{3} \mu _{4}}{\alpha _{1} \gamma _{1} r_{1} r _{3}} - \frac{ ( P_{2} \overline{C} - \beta _{2} r_{2} ( K_{2} + \overline{S} )^{2} ) ( P_{3} \overline{C} - ( \gamma _{1} + \gamma _{2} ) r _{3} ( K_{3} + \overline{R} )^{2} )}{\beta _{1} \gamma _{1} r_{2} r_{3} ( K_{2} + \overline{S} ) ^{2} ( K_{3} + \overline{R} )^{2}} \biggr) e ^{- \overline{\xi _{3}}} \\ \quad > \biggl( \frac{\mu _{3} \mu _{4}}{\alpha _{1} \gamma _{1} r_{1} r_{3}} - \frac{ ( P_{1} \overline{C} - ( \alpha _{1} + \alpha _{2} ) r_{1} ( K_{1} + \overline{G} )^{2} ) ( P_{3} \overline{C} - ( \gamma _{1} + \gamma _{2} ) r_{3} ( K_{3} + \overline{R} )^{2} )}{\alpha _{1} \gamma _{1} r_{1} r_{3} ( K_{1} + \overline{G} )^{2} ( K _{3} + \overline{R} )^{2}} \biggr) e^{- \overline{\xi _{3}}} e^{- \overline{\xi _{1}}}. \end{gathered} $$

(3.18)

Finally, from

$$ r_{1} > \frac{1}{\alpha _{1}} \biggl( \frac{\mu _{3} \mu _{4}}{\gamma _{1} r_{3}} + \frac{\mu _{1} \mu _{2}}{\beta _{1} r_{2}} \biggr) > \frac{ ( P_{3} \overline{C} - \gamma _{2} r_{3} ( K_{3} + \overline{R} )^{2} ) ( P_{1} \overline{C} - ( \alpha _{1} + \alpha _{2} ) r_{1} ( K _{1} + \overline{G} )^{2} )}{\alpha _{1} \gamma _{1} r_{3} ( K_{1} + \overline{G} )^{2} ( K _{3} + \overline{R} )^{2}}. $$

(3.19)

This completes part (i) of the theorem.

From (ii) and (iii), we will consider the case for

$$ D_{2}^{+} =1+ a_{1} - a_{0}^{2} - a_{0} a_{2} >0\quad \text{and}\quad D_{2}^{-} =1- a_{1} - a_{0}^{2} + a_{0} a_{2} >0. $$

From

$$ \begin{gathered} a_{0} ( a_{0} + a_{2} ) < 1+ a_{1} \quad \text{and}\quad 2+ a_{0} ( a _{2} - a_{0} ) >1+ a_{1} \\ \quad \Rightarrow\quad 2+ a_{0} ( a_{2} - a_{0} ) > a_{0} ( a_{0} + a_{2} ) \quad \Rightarrow \quad a_{0} < 1. \end{gathered} $$

In this case, we have to show that $a_{0} <1$, that is,

$$ a_{13} a_{31} a_{22} + a_{23} a_{32} a_{11} + a_{12} a_{21} a_{33} < 1+ a_{11} a_{22} a_{33} + a_{21} a_{32} a_{13} + a_{31} a_{12} a_{23}, $$

where we obtain

$$ \begin{gathered}[b] \mu _{2} \mu _{4} \rho ( K_{2} + \overline{S} )^{2} \bigl( 1- e^{- \overline{\xi _{2}}} \bigr) \\ \quad > \bigl( P_{2} \overline{C} - \beta _{2} r_{2} ( K_{2} + \overline{S} ) ^{2} \bigr) - \bigl( P_{2} \overline{C} - ( \beta _{1} + \beta _{2} ) r_{2} ( K_{2} + \overline{S} ) ^{2} \bigr) e^{- \overline{\xi _{2}}} \end{gathered} $$

(3.20)

and

$$\begin{aligned}& \mu _{1} \mu _{4} ( K_{1} + \overline{G} )^{2} \bigl( 1- e^{- \overline{\xi _{1}}} \bigr) \\& \quad < \rho \bigl( \bigl( P_{1} \overline{C} - \alpha _{2} r_{1} ( K_{1} + \overline{G} )^{2} \bigr) - \bigl( P_{1} \overline{C} - ( \alpha _{1} + \alpha _{2} ) r _{1} ( K_{1} + \overline{G} )^{2} \bigr) e ^{- \overline{\xi _{1}}} \bigr), \end{aligned}$$

(3.21)

which holds for

$$\begin{aligned}& \bigl( P_{2} \overline{C} - \beta _{2} r_{2} ( K_{2} + \overline{S} )^{2} \bigr) - \mu _{2} \mu _{4} \rho ( K _{2} + \overline{S} )^{2} \\& \quad < \bigl( \bigl( P_{2} \overline{C} - ( \beta _{1} + \beta _{2} ) r_{2} ( K_{2} + \overline{S} )^{2} \bigr) - \mu _{2} \mu _{4} \rho ( K_{2} + \overline{S} ) ^{2} \bigr) e^{- \overline{\xi _{2}}} \end{aligned}$$

and

$$ \begin{gathered} \rho \bigl( P_{1} \overline{C} - \alpha _{2} r_{1} ( K _{1} + \overline{G} )^{2} \bigr) - \mu _{1} \mu _{4} ( K _{1} + \overline{G} )^{2} \\ \quad > \bigl( \rho \bigl( P_{1} \overline{C} - ( \alpha _{1} + \alpha _{2} ) r _{1} ( K_{1} + \overline{G} )^{2} \bigr) - \mu _{1} \mu _{4} ( K_{1} + \overline{G} )^{2} \bigr) e ^{- \overline{\xi _{1}}}, \end{gathered} $$

if

$$ \frac{\mu _{1}}{\mu _{2}} < \frac{\rho ^{2} ( P_{1} \overline{C} - \alpha _{2} r_{1} ( K_{1} + \overline{G} )^{2} ) ( K_{2} + \overline{S} )^{2}}{ ( P_{2} \overline{C} - \beta _{2} r_{2} ( K_{2} + \overline{S} ) ^{2} ) ( K_{1} + \overline{G} )^{2}}, $$

and

$$ \overline{\xi _{2}} < \ln \biggl( \frac{ ( P_{2} \overline{C} - ( \beta _{1} + \beta _{2} ) r_{2} ( K_{2} + \overline{S} )^{2} ) - \mu _{2} \mu _{4} \rho ( K _{2} + \overline{S} )^{2}}{ ( P_{2} \overline{C} - \beta _{2} r_{2} ( K_{2} + \overline{S} )^{2} ) - \mu _{2} \mu _{4} \rho ( K_{2} + \overline{S} )^{2}} \biggr), $$

and

$$ \overline{\xi _{1}} > \ln \biggl( \frac{\rho ( P_{1} \overline{C} - ( \alpha _{1} + \alpha _{2} ) r _{1} ( K_{1} + \overline{G} )^{2} ) - \mu _{1} \mu _{4} ( K_{1} + \overline{G} )^{2}}{\rho ( P_{1} \overline{C} - \alpha _{2} r_{1} ( K_{1} + \overline{G} )^{2} ) - \mu _{1} \mu _{4} ( K_{1} + \overline{G} ) ^{2}} \biggr). $$

This completes the proof. □

Theorem 3.3

Let $\varLambda _{1} =( \overline{G}, \overline{S}, \overline{R}, \overline{C})$ be the positive equilibrium point of system (1.3) and assume that the conditions in Theorem 2.1(i) and Theorem 3.2 hold. Then the positive equilibrium point of system (1.3) is globally asymptotically stable if

$$\begin{aligned}& 0< \xi _{1} ( n ) < \ln \biggl( \frac{2 \overline{G} -G(n)}{G(n)} \biggr)\quad \textit{and}\quad G ( n ) < \overline{G}, \\& 0< \xi _{2} ( n ) < \ln \biggl( \frac{2 \overline{ S} - S (n)}{S (n)} \biggr)\quad \textit{and}\quad S(n) < \overline{S}, \end{aligned}$$

and

$$ 0< \xi _{3} ( n ) < \ln \biggl( \frac{2 \overline{ R} - R (n)}{R (n)} \biggr)\quad \textit{and}\quad R ( n ) < \overline{R}. $$

Proof

Let $\varLambda _{1} =( \overline{G}, \overline{S}, \overline{R}, \overline{N}, \overline{C})$ be the positive equilibrium point of system (1.3), and let us consider a Lyapunov function $V ( n )$ defined by

$$ V ( n ) = \bigl( X(n)- \overline{X} \bigr)^{2}, \quad n=0, 1, 2, \dots , $$

(3.22)

where $X ( n ) = ( G ( n ), S ( n ), R ( n ),C ( n ) )$, and $\varLambda _{1} =( \overline{G}, \overline{S}, \overline{R}, \overline{C} )$.

The change along the solutions of the system is

$$\begin{aligned} \Delta V ( n ) &=V ( n+1 ) -V ( n ) = \bigl( X ( n+1 ) - \overline{X} \bigr)^{2} - \bigl( X ( n ) - \overline{X} \bigr)^{2} \\ &= \bigl( X ( n+1 ) -X ( n ) \bigr) \bigl( X ( n+1 ) +X ( n ) -2 \overline{X} \bigr). \end{aligned}$$

(3.23)

From the first equation of system (1.3), we have

$$\begin{aligned} \Delta V_{1} ( n ) &= \bigl( G ( n+1 ) -G ( n ) \bigr) \bigl( G ( n+1 ) +G ( n ) -2 \overline{G} \bigr) \\ &= \biggl( \frac{G ( n ) \xi _{1} (n)}{ ( \xi _{1} (n)- \alpha _{1} r_{1} G ( n ) ) \cdot e^{- \xi _{1} (n)} + \alpha _{1} r_{1} G ( n )} -G ( n ) \biggr) \\ &\quad {}\times \biggl( \frac{G ( n ) \xi _{1} (n)}{ ( \xi _{1} (n)- \alpha _{1} r _{1} G ( n ) ) \cdot e^{- \xi _{1} (n)} + \alpha _{1} r_{1} G ( n )} +G ( n ) -2 \overline{G} \biggr) \\ &= \frac{1}{ ( T ( n ) )^{2}} \cdot \bigl( G ( n ) \xi _{1} (n)-T ( n ) \cdot G ( n ) \bigr) \cdot \bigl( G ( n ) \xi _{1} (n)+ \bigl( G ( n ) -2 \overline{G} \bigr) T ( n ) \bigr), \end{aligned}$$

where $T ( n ) = ( \xi _{1} (n)- \alpha _{1} r_{1} G ( n ) ) \cdot e^{- \xi _{1} (n)} + \alpha _{1} r_{1} G ( n )$. In this case, if Theorem 2.1(i) holds, then

$$ G ( n ) \xi _{1} (n)-T ( n ) \cdot G ( n ) =G ( n ) \bigl( \xi _{1} (n)- \alpha _{1} r_{1} G ( n ) \bigr) \bigl( 1- e^{- \xi _{1} (n)} \bigr) >0. $$

(3.24)

So, we have to obtain

$$ G ( n ) \xi _{1} (n)+ \bigl( G ( n ) -2 \overline{G} \bigr) T ( n ) < 0 $$

(3.25)

to get $\Delta V_{1} ( n ) <0$, which holds for

$$ 0< \xi _{1} (n)< \ln \biggl( \frac{2 \overline{G} -G(n)}{G(n)} \biggr) \quad \text{and}\quad G ( n ) < \overline{G}. $$

(3.26)

This implies that $\lim_{n\rightarrow \infty } G ( n ) = \overline{G}$. Similarly, we can obtain the conditions

$$ 0< \xi _{2} ( n ) < \ln \biggl( \frac{2 \overline{ S} - S (n)}{S (n)} \biggr)\quad \text{and}\quad S(n) < \overline{S} $$

(3.27)

for $\Delta V_{2} ( n ) <0$ and

$$ 0< \xi _{3} ( n ) < \ln \biggl( \frac{2 \overline{ R} - R (n)}{R (n)} \biggr)\quad \text{and}\quad R ( n ) < \overline{R} $$

(3.28)

for $\Delta V_{3} ( n ) <0$. This completes the proof. □

Example 1

We use a system of differential equations with piecewise constant arguments to explain the growth of GBM and interactions among the glial cells, the cancer cells, and the chemotherapeutic agents. The below demonstrated graphs are as follows: the blue graph shows the sensitive tumor cells, the red represents the resistant cells, the yellow represents the glial cells, and the green shows the treatment agents.

The initial conditions are $G ( 0 ) =0.99$, $C ( 0 ) =0.01$, $S ( 0 ) =0.3$, and $R ( 0 ) =0.035$.

The values of the parameters used in this example are from Table 1. In Fig. 2 and Fig. 3 the chemotherapy agent rate for infusion is $\sigma =0.5$. For $( r_{1}, \frac{G(n+1)}{G(n)} )$, $( r_{2}, \frac{S(n+1)}{S(n)} )$, $( r_{3}, \frac{R(n+1)}{R(n)} )$, with an unchanged dosage of chemotherapy, we noticed that after a specific time, the cancer cells cover the glial cells, and the interaction between the cancer cells and glial cells goes to an undesired case. Figure 3 shows the population growth against the fitness of the cells $( G(n), \frac{G(n+1)}{G(n)} )$, $( S(n), \frac{S(n+1)}{S(n)} )$, $( R(n), \frac{R(n+1)}{R(n)} )$, and $( C(n), \frac{C(n+1)}{C(n)} )$.

Figure 2

Bifurcation diagram of (1.3), where $\sigma =0.5$

Figure 3

Dynamical behavior of the population growth due to the fitness of each class

Table 1 Description of the parameters according to the literature [14]

Figure 4 and Fig. 5 show the impact of increasing the dosage of the chemotherapeutic treatment to $\sigma =50$. We realize that we could minimize the population of the resistant cells while still there is a strong interaction between the glial cells and the tumor cells.

Figure 4

Bifurcation diagram of (1.3), where $\sigma =0.5$

Figure 5

Dynamical behavior of the population growth due to the fitness of each class

4 Analysis of Neimark–Sacker bifurcation

For Neimark–Sacker bifurcation, one has to show that

$$ D_{2}^{+} =1+ a_{1} - a_{0}^{2} - a_{0} a_{2} >0\quad \text{and}\quad D_{2}^{-} =1- a_{1} - a_{0}^{2} + a_{0} a_{2} =0, $$

(4.1)

while the conditions in Theorem 3.2(i) hold. From (4.1), we show at first that

$$ D_{2}^{+} =1+ a_{1} - a_{0} ( a_{0} + a_{2} )\quad \text{and}\quad D _{2}^{-} =1- a_{1} - a_{0} ( a_{0} - a_{2} ), $$

(4.2)

where $a_{0} =0$ and $a_{1} =1$. From Theorem 3.2(i), it is proven that

$$ 1+ a_{11} a_{22} + a_{11} a_{33} + a_{22} a_{33} > a_{13} a_{31} + a _{23} a_{32} + a_{12} a_{21}, $$

and since $a_{11} + a_{22} + a_{33}>1$, we have also

$$ a_{11} + a_{22} + a_{33} + a_{11} a_{22} a_{33} + a_{21} a_{32} a_{13} + a_{31} a_{12} a_{23} - a_{13} a_{31} a_{22} - a_{23} a_{32} a_{11} - a_{12} a_{21} a_{33} >0. $$

In this case, we will consider the case, where

$$ a_{13} a_{31} a_{22} + a_{23} a_{32} a_{11} + a_{12} a_{21} a_{33} = a _{11} a_{22} a_{33} + a_{21} a_{32} a_{13} + a_{31} a_{12} a_{23}. $$

(4.3)

From (4.3), we have

$$\begin{aligned}& a_{23} a_{32} a_{11} = a_{31} a_{12} a_{23} \\& \quad \Rightarrow \quad \bigl( \bigl( \mu _{1} \mu _{4} + ( \alpha _{1} + \alpha _{2} ) r_{1} \bigr) ( K_{1} + \overline{G} )^{2} - \rho P_{1} \overline{C} \bigr) e^{- \overline{\xi _{1}}} \\& \hphantom{\quad \Rightarrow \quad }\quad = ( \mu _{1} \mu _{4} + \alpha _{2} r_{1} ) ( K_{1} + \overline{G} )^{2} - \rho P_{1} \overline{C}, \end{aligned}$$

where we obtain

$$ \begin{gathered}[b] \overline{\xi _{1}} = \ln \biggl( \frac{ ( \mu _{1} \mu _{4} + ( \alpha _{1} + \alpha _{2} ) r_{1} ) ( K_{1} + \overline{G} )^{2} - \rho P_{1} \overline{C}}{ ( \mu _{1} \mu _{4} + \alpha _{2} r_{1} ) ( K _{1} + \overline{G} )^{2} - \rho P_{1} \overline{C}} \biggr) \\ \quad \text{for } \overline{C} < \frac{ ( \alpha _{2} r _{1} + \mu _{1} \mu _{4} ) ( K_{1} + \overline{G} ) ^{2}}{\rho P_{1}}.\end{gathered} $$

(4.4)

Furthermore, from

$$\begin{aligned}& a_{13} a_{31} a_{22} = a_{11} a_{22} a_{33} \\& \quad \Rightarrow\quad \mu _{3} \mu _{4} \bigl( 1- e^{- \overline{\xi _{1}}} \bigr) \bigl( 1- e^{- \overline{\xi _{3}}} \bigr) ( K_{1} + \overline{G} )^{2} ( K_{3} + \overline{R} )^{2} \\& \hphantom{\quad \Rightarrow\quad}\quad > \bigl( \bigl( P_{1} \overline{C} - \alpha _{2} r_{1} ( K_{1} + \overline{G} )^{2} \bigr) - \bigl( P_{1} \overline{C} - ( \alpha _{1} + \alpha _{2} ) r_{1} ( K_{1} + \overline{G} ) ^{2} \bigr) e^{- \overline{\xi _{1}}} \bigr) \\& \hphantom{\quad \Rightarrow\quad}\qquad {}\times \bigl( \bigl( P_{3} \overline{C} - \gamma _{2} r_{3} ( K_{3} + \overline{R} )^{2} \bigr) - \bigl( P_{3} \overline{C} - ( \gamma _{1} + \gamma _{2} ) r_{3} ( K_{3} + \overline{R} )^{2} \bigr) e^{- \overline{\xi _{3}}} \bigr), \end{aligned}$$

we obtain

$$ \overline{\xi _{1}} = \ln \biggl( \frac{ ( ( \alpha _{1} + \alpha _{2} ) r_{1} + \mu _{3} ) ( K_{1} + \overline{G} )^{2} - P_{1} \overline{C}}{ ( \alpha _{2} r_{1} + \mu _{3} ) ( K_{1} + \overline{G} ) ^{2} - P_{1} \overline{C}} \biggr)\quad \text{for } \overline{C} < \frac{ ( \alpha _{2} r_{1} + \mu _{3} ) ( K_{1} + \overline{G} )^{2}}{P_{1}} $$

(4.5)

and

$$ \overline{\xi _{3}} = \ln \biggl( \frac{ ( ( \gamma _{1} + \gamma _{2} ) r_{3} + \mu _{4} ) ( K_{3} + \overline{R} )^{2} - P_{3} \overline{C}}{ ( \gamma _{2} r_{3} + \mu _{4} ) ( K_{3} + \overline{R} ) ^{2} - P_{3} \overline{C}} \biggr)\quad \text{for } \overline{C} < \frac{ ( \gamma _{2} r_{3} + \mu _{4} ) ( K_{3} + \overline{R} )^{2}}{P_{3}}. $$

(4.6)

Considering (4.4) and (4.5) together, we get $\mu _{3} = \mu _{1} \mu _{4}$ and $\rho <1$. Finally, from

$$ a_{12} a_{21} a_{33} = a_{11} a_{22} a_{33} \quad \Rightarrow\quad a_{12} a_{21} = a_{11} a_{22}, $$

where we get

$$ \overline{\xi _{1}} = \ln \biggl( \frac{ ( ( \alpha _{1} + \alpha _{2} ) r_{1} + \mu _{1} ) ( K_{1} + \overline{G} )^{2} - P_{1} \overline{C}}{ ( \alpha _{2} r_{1} + \mu _{1} ) ( K_{1} + \overline{G} ) ^{2} - P_{1} \overline{C}} \biggr)\quad \text{for } \overline{C} < \frac{ ( \alpha _{2} r_{1} + \mu _{1} ) ( K_{1} + \overline{G} )^{2}}{P_{1}} $$

(4.7)

and

$$ \overline{\xi _{2}} = \ln \biggl( \frac{ ( ( \beta _{1} + \beta _{2} ) r_{2} + \mu _{2} ) ( K_{2} + \overline{S} )^{2} - P_{2} \overline{C}}{ ( \beta _{2} r_{2} + \mu _{2} ) ( K_{1} + \overline{S} ) ^{2} - P_{2} \overline{C}} \biggr)\quad \text{for } \overline{C} < \frac{ ( \beta _{2} r_{2} + \mu _{2} ) ( K_{1} + \overline{S} )^{2}}{P_{2}}, $$

(4.8)

where $\mu _{3} = \mu _{2} = \mu _{1} \mu _{4}$. This completes the proof that $a_{0} =0$. To consider $a_{1} =1$, we have to prove only that

$$ a_{22} a_{33} =1+ a_{23} a_{32}. $$

(4.9)

From (4.9), we obtain

$$\begin{aligned}& \bigl( \bigl( P_{2} \overline{C} - \beta _{2} r_{2} ( K_{2} + \overline{S} )^{2} \bigr) \bigl( P_{3} \overline{C} - \gamma _{2} r_{3} ( K_{3} + \overline{R} )^{2} \bigr) - \bigl( \rho ^{2} - \beta _{1} r_{2} \gamma _{1} r _{3} \bigr) \bigr) \\& \quad {}- \bigl( - \bigl( P_{2} \overline{C} - \beta _{2} r_{2} ( K_{2} + \overline{S} )^{2} \bigr) \bigl( ( \gamma _{1} + \gamma _{2} ) r_{3} ( K_{3} + \overline{R} )^{2} - P_{3} \overline{C} \bigr) + \rho ^{2} ( K_{2} + \overline{S} )^{2} ( K_{3} + \overline{R} )^{2} \bigr) e^{- \overline{\xi _{3}}} \\& \quad {}- \bigl( - \bigl( P_{3} \overline{C} - \gamma _{2} r_{3} ( K_{3} + \overline{R} )^{2} \bigr) \bigl( ( \beta _{1} + \beta _{2} ) r_{2} ( K_{2} + \overline{S} )^{2} -P_{2} \overline{C} \bigr) + \rho ^{2} ( K_{2} + \overline{S} )^{2} ( K_{3} + \overline{R} )^{2} \bigr) e^{- \overline{\xi _{2}}} \\& \quad {}+ \bigl( \bigl( P_{2} \overline{C} - ( \beta _{1} + \beta _{2} ) r_{2} ( K_{2} + \overline{S} )^{2} \bigr) \bigl( P_{3} \overline{C} - ( \gamma _{1} + \gamma _{2} ) r_{3} ( K_{3} + \overline{R} )^{2} \bigr) \\& \quad {} + \rho ^{2} ( K_{2} + \overline{S} )^{2} ( K_{3} + \overline{R} )^{2} \bigr) e^{- \overline{\xi _{2}}} e^{- \overline{\xi _{3}}} =0, \end{aligned}$$

where

$$ \begin{aligned}[b] \overline{\xi _{3}} &= \ln \bigl( \bigl( \bigl( ( \beta _{1} + \beta _{2} ) r_{2} ( K_{2} + \overline{S} ) ^{2} - P_{2} \overline{C} \bigr) \bigl( ( \gamma _{1} + \gamma _{2} ) r_{3} ( K_{3} + \overline{R} ) ^{2} - P_{3} \overline{C} \bigr) \\ &\quad {} + \rho ^{2} ( K_{2} + \overline{S} )^{2} ( K_{3} + \overline{R} )^{2}\bigr) \\ &\quad {}\big/\bigl( \rho ^{2} ( K_{2} + \overline{S} )^{2} ( K_{3} + \overline{R} )^{2} \\ &\quad {} - \bigl( ( \beta _{1} + \beta _{2} ) r_{2} ( K_{2} + \overline{S} )^{2} - P _{2} \overline{C} \bigr)\bigl( \gamma _{2} r_{3} ( K _{3} + \overline{R} )^{2} - P_{3} \overline{C} \bigr)\bigr) \bigr) \end{aligned} $$

(4.10)

and for $\rho = \sqrt{\beta _{1} r_{2} \gamma _{1} r_{3}}$, we get

$$ \overline{\xi _{3}} = \ln \biggl( \frac{\rho ^{2} ( K_{2} + \overline{S} )^{2} ( K_{3} + \overline{R} )^{2} - ( P_{2} \overline{C} - \beta _{2} r_{2} ( K_{2} + \overline{S} )^{2} ) ( ( \gamma _{1} + \gamma _{2} ) r_{3} ( K_{3} + \overline{R} ) ^{2} - P_{3} \overline{C} )}{ ( P_{2} \overline{C} - \beta _{2} r_{2} ( K_{2} + \overline{S} )^{2} ) ( P_{3} \overline{C} - \gamma _{2} r_{3} ( K_{3} + \overline{R} )^{2} )} \biggr). $$

(4.11)

From (4.6), (4.10), and (4.11), we have

$$ \overline{C} = \frac{r_{3} ( \gamma _{1} - \gamma _{2} ) ( K_{3} + \overline{R} )^{2}}{P_{3}}, \qquad 2 \gamma _{2} = \gamma _{1},\quad \text{and}\quad \frac{r_{2}}{ r_{3}} = \frac{P_{2} ( \gamma _{1} - \gamma _{2} ) ( K_{3} + \overline{R} )^{2}}{P_{3} ( \beta _{1} + \beta _{2} ) ( K_{2} + \overline{S} )^{2}}, $$

which completes the proof.

The theorem is in this case as follows.

Theorem 4.1

Let $\varLambda _{1} =( \overline{G}, \overline{S}, \overline{R}, \overline{C})$ be a positive equilibrium point of system (1.3) and assume that the conditions in Theorem 3.2(i) hold. Furthermore, let $\overline{C} = \frac{r_{3} ( \gamma _{1} - \gamma _{2} ) ( K_{3} + \overline{R} )^{2}}{P _{3}}$, $2 \gamma _{2} = \gamma _{1}$, $\rho = \sqrt{\beta _{1} r_{2} \gamma _{1} r_{3}}$, and $\frac{r_{2}}{r _{3}} = \frac{P_{2} ( \gamma _{1} - \gamma _{2} ) ( K _{3} + \overline{R} )^{2}}{P_{3} ( \beta _{1} + \beta _{2} ) ( K_{2} + \overline{S} )^{2}}$. If

$$\begin{aligned}& \overline{\xi _{1}} = \ln \biggl( \frac{ ( \mu _{1} \mu _{4} + ( \alpha _{1} + \alpha _{2} ) r_{1} ) ( K_{1} + \overline{G} )^{2} - \rho P_{1} \overline{C}}{ ( \mu _{1} \mu _{4} + \alpha _{2} r_{1} ) ( K _{1} + \overline{G} )^{2} - \rho P_{1} \overline{C}} \biggr), \end{aligned}$$

(4.12)

$$\begin{aligned}& \overline{\xi _{2}} = \ln \biggl( \frac{ ( ( \beta _{1} + \beta _{2} ) r_{2} + \mu _{2} ) ( K_{1} + \overline{S} )^{2} - P_{2} \overline{C}}{ ( \beta _{2} r_{2} + \mu _{2} ) ( K_{1} + \overline{S} ) ^{2} - P_{2} \overline{C}} \biggr), \end{aligned}$$

(4.13)

and

$$ \overline{\xi _{3}} = \ln \biggl( \frac{ ( ( \gamma _{1} + \gamma _{2} ) r_{3} + \mu _{4} ) ( K_{1} + \overline{R} )^{2} - P_{3} \overline{C}}{ ( \gamma _{2} r_{3} + \mu _{4} ) ( K_{1} + \overline{R} ) ^{2} - P_{3} \overline{C}} \biggr), $$

(4.14)

where $\mu _{3} = \mu _{2} = \mu _{1} \mu _{4}$ and $\rho <1$. In this case, Eq. (1.3) undergoes Neimark–Sacker bifurcation.

5 Stability analysis for a model with Allee effect

In 1931, Allee [30] introduced a significant point in population dynamics, related to the low-density size of the species. He showed that the Allee effect occurs when the population growth rate reduces at low population size.

It is known that the logistic model assumes that the per capita growth rate declines monotonically when the density increases. However, for a low-density model, the logistic equation is incapable of representing the real biological meaning of the suggested phenomena. Logistic equations with Allee effect are capable of representing the real meaning of the events, where it can be seen that for low-density size the per capita growth rate gives a humped curve up to a maximum intermediate density and then declines again.

In literature, Allee effect can be divided into two main types based on whether or not the population exhibits a critical population size or density: (i) strong Allee effect and (ii) weak Allee effect.

A population exhibiting a strong Allee effect will have a critical population size below the population growth rate which becomes negative: when the population density hits a number below this threshold, it will go to extinction without any further aid. A population exhibiting a weak Allee effect will possess a reduced per capita growth rate at lower population density or size [31,32,33,34,35].

Let (1.1) be as follows:

$$ \textstyle\begin{cases} \frac{dG}{dt} = r_{1} G ( t ) ( K_{1} - \alpha _{1} G ( t ) - \alpha _{2} G ( [\!\![t ]\!\!]) ) - \mu _{1} G ( t ) S ( [\!\![t ]\!\!]) - \mu _{3} G ( t ) R ( [\!\![t ]\!\!]) \\ \hphantom{\frac{dG}{dt} =}{}- \frac{P _{1} G ( t ) C ( [\!\![t ]\!\!])}{K_{1} + G ( [\!\![t ]\!\!])}, \\ \frac{{dS}}{{dt}} = ( {pS} ( t ) + r_{2} S ( t ) ( K_{2} - \beta _{1} S ( t ) - \beta _{2} S ( [\!\![t ]\!\!]) ) - \mu _{2} S ( t ) G ( [\!\![t ]\!\!]) -\rho S ( t ) R ( [\!\![t ]\!\!]) \\ \hphantom{\frac{{dS}}{{dt}} =}{}- \frac{P_{2} S ( t ) C ( [\!\![t ]\!\!])}{K _{2} + S ( [\!\![t ]\!\!])} ) ( \frac{S ( [\!\![t ]\!\!])}{K_{2}^{*}} -1 ), \\ \frac{{dR}}{{dt}} = ( r_{3} R ( t ) ( K_{3} - \gamma _{1} R ( t ) - \gamma _{2} R ( [\!\![t ]\!\!]) ) - \mu _{4} R ( t ) G ( [\!\![t ]\!\!]) +\rho S ( [\!\![t ]\!\!]) R ( t ) \\ \hphantom{\frac{{dR}}{{dt}} = }{}- \frac{P _{3} R ( t ) C ( [\!\![t ]\!\!])}{K_{3} + R ( [\!\![t ]\!\!])} ) ( \frac{R ( [\!\![t ]\!\!])}{K _{3}^{*}} -1 ), \\ \frac{{dC}}{{dt}} =\sigma - \omega _{1} C ( t ) - \omega _{2} C ( [\!\![t ]\!\!]), \end{cases} $$

(5.1)

where $K_{2}^{*}$ and $K_{3}^{*}$ are the thresholds of both tumor populations $S(t)$ and $R ( t )$. Below the threshold level, the growth rate decreases and the tumor population goes to extinction. This situation describes a strong Allee effect. In the following, we consider the conditions for a strong Allee effect.

The discrete solution of Eq. (5.1) on the interval $t \in [ n,n+1 )$ can be obtained as

$$ \textstyle\begin{cases} G ( n+1 ) = \frac{G ( n ) \xi _{1} (n)}{ ( \xi _{1} (n)- \alpha _{1} r_{1} G(n) ) e^{- \xi _{1} (n)} + \alpha _{1} r_{1} G ( n )}, \\ S ( n+1 ) = \frac{S ( n ) \xi _{2} (n)}{ ( \xi _{2} (n)- \beta _{1} r_{2} S ( n ) ) e^{- ( \frac{ \overline{S}}{K_{2}^{*}} -1 ) \xi _{2} (n)} + \beta _{1} r_{2} S ( n )}, \\ R ( n+1 ) = \frac{R ( n ) \xi _{3} (n)}{ ( \xi _{3} (n)- \gamma _{1} r_{3} R ( n ) ) e^{- ( \frac{\overline{R}}{K_{3}^{*}} -1 ) \xi _{3} (n)} + \gamma _{1} r_{3} R ( n )}, \\ C ( n+1 ) = \frac{\sigma ( 1- e^{- \omega _{1}} )}{ \omega _{1}} + \frac{ ( ( \omega _{1} + \omega _{2} ) e ^{- \omega _{1}} - \omega _{2} ) C(n)}{\omega _{1}}. \end{cases} $$

(5.2)

The following theorems are given without proves since they are similar to Theorem 3.2 and Theorem 3.3.

Theorem 5.1

Let $\varLambda _{1} =( \overline{G}, \overline{S}, \overline{R}, \overline{C})$ be a positive equilibrium point of system (5.1). Assume that

$$ \begin{gathered} \frac{\mu _{1} \mu _{2}}{\mu _{3} \mu _{4}} > \frac{ ( P_{2} \overline{C} - ( \beta _{1} + \beta _{2} ) r_{2} ( K_{2} + \overline{S} )^{2} ) ( K_{3} + \overline{R} )^{2}}{ ( P_{3} \overline{C} - ( \gamma _{1} + \gamma _{2} ) r_{3} ( K_{3} + \overline{R} )^{2} ) ( K_{2} + \overline{S} ) ^{2}},\\ \frac{\mu _{1}}{\mu _{2}} < \frac{\rho ^{2} ( P_{1} \overline{C} - \alpha _{2} r_{1} ( K_{1} + \overline{G} )^{2} ) ( K_{2} + \overline{S} )^{2}}{ ( P_{2} \overline{C} - \beta _{2} r_{2} ( K_{2} + \overline{S} )^{2} ) ( K_{1} + \overline{G} ) ^{2}}, \end{gathered} $$

(5.3)

where $\overline{C} > \frac{ ( \alpha _{1} + \alpha _{2} ) r_{1} ( K_{1} + \overline{G} )^{2}}{P_{1}} > \frac{ ( \gamma _{1} + \gamma _{2} ) r_{3} ( K_{3} + \overline{R} )^{2}}{P_{3}} > \frac{ ( \beta _{1} + \beta _{2} ) r_{2} ( K_{2} + \overline{S} ) ^{2}}{P_{2}}$. If

$$ r_{1} > \frac{1}{\alpha _{1}} \biggl( \frac{\mu _{3} \mu _{4}}{\gamma _{1} r_{3}} + \frac{\mu _{1} \mu _{2}}{\beta _{1} r_{2}} \biggr) > \frac{ ( P_{3} \overline{C} - \gamma _{2} r_{3} ( K_{3} + \overline{R} )^{2} ) ( P_{1} \overline{C} - ( \alpha _{1} + \alpha _{2} ) r_{1} ( K _{1} + \overline{G} )^{2} )}{\alpha _{1} \gamma _{1} r_{3} ( K_{1} + \overline{G} )^{2} ( K _{3} + \overline{R} )^{2}}, $$

(5.4)

and the conditions

$$ \overline{\xi _{2}} < \ln \biggl( \frac{ ( P_{2} \overline{C} - ( \beta _{1} + \beta _{2} ) r_{2} ( K_{2} + \overline{S} )^{2} ) - \mu _{2} \mu _{4} \rho ( K _{2} + \overline{S} )^{2}}{ ( P_{2} \overline{C} - \beta _{2} r_{2} ( K_{2} + \overline{S} )^{2} ) - \mu _{2} \mu _{4} \rho ( K_{2} + \overline{S} )^{2}} \biggr) ^{\frac{K_{2}^{*}}{\overline{S} - K_{2}^{*}}} $$

and

$$ \overline{\xi _{1}} > \ln \biggl( \frac{\rho ( P_{1} \overline{C} - ( \alpha _{1} + \alpha _{2} ) r _{1} ( K_{1} + \overline{G} )^{2} ) - \mu _{1} \mu _{4} ( K_{1} + \overline{G} )^{2}}{\rho ( P_{1} \overline{C} - \alpha _{2} r_{1} ( K_{1} + \overline{G} )^{2} ) - \mu _{1} \mu _{4} ( K_{1} + \overline{G} ) ^{2}} \biggr)^{\frac{K_{3}^{*}}{\overline{R} - K_{3}^{*}}} $$

hold, then the positive equilibrium point is locally asymptotically stable.

Theorem 5.2

Let $\varLambda _{1} =( \overline{G}, \overline{S}, \overline{R}, \overline{C})$ be the positive equilibrium point of system (5.1) and assume that the conditions in Theorem 2.1(i) and Theorem 5.1 hold. Then the positive equilibrium point of system (5.1) is globally asymptotically stable if

$$\begin{aligned}& 0< \xi _{1} ( n ) < \ln \biggl( \frac{2 \overline{G} -G(n)}{G(n)} \biggr)\quad \textit{and}\quad G ( n ) < \overline{G}, \end{aligned}$$

(5.5)

$$\begin{aligned}& 0< \xi _{2} ( n ) < \ln \biggl( \frac{2 \overline{ S} - S (n)}{S (n)} \biggr)^{\frac{K_{2} ^{*}}{\overline{S} - K_{2}^{*}}}\quad \textit{and}\quad S( n) < \overline{S}, \end{aligned}$$

(5.6)

and

$$ 0< \xi _{3} ( n ) < \ln \biggl( \frac{2 \overline{ R} - R (n)}{R (n)} \biggr)^{\frac{K_{3} ^{*}}{\overline{R} - K_{3}^{*}}}\quad \textit{and}\quad R ( n ) < \overline{R}. $$

(5.7)

Example 2

For early detection of tumor cells, we incorporate the Allee threshold to get results for the extinction case. The threshold for both cancer cells is given as $K_{2}^{*} = K_{3}^{*} =0.5\mbox{ mm}^{3}$, which corresponds to a tumor population in day 17 (see [36]). Applying the same values as given in Example 1 with $\sigma =0.5$, we noticed that for an early detection scenario, the dosage is enough to keep the sensitive cancer cells in a specific density so that they do not convert to a population of resistant tumor cells. Moreover, the glial cells are insignificantly affected by the dosage of the treatment process.

It is known that the Allee effect plays an essential role in the stability analysis of equilibrium points of population dynamics. Generally, an Allee effect has a stabilizing effect on population dynamics and shows realistic results for low population densities, as demonstrated in Fig. 6.

Figure 6

Bifurcation diagram of (1.3), where $\sigma =0.5$

6 Conclusion

In this paper, we construct a model as a system of differential equations with piecewise constant arguments that was established in [14] as a system of differential equations. We consider the overlapping and non-overlapping events in the model and incorporate both continuous- and discrete-time to reach to a more realistic interpretation for the cancer growth and the treatment process. Regarding the resting time and continuous growth of the tumor density, and emphasizing the state that the dosage of the drug changes in discrete-time according to the findings of the treatment, we designed a hybrid (continuous and discrete) time model.

First of all, we obtained the positive equilibrium point of system (1.3) to consider the case, where the interaction between the cells is active and the treatment process started. We showed in Sect. 2 that the system has positive solutions and obtained a permanent interval under specific conditions. It is seen that the solutions have single semi-cycle behavior. However, it is also proven that the solutions of the system do not have a period two behavior. For the obtained conditions, it is shown that system (1.3) has only monotonic increasing or decreasing behaviors.

In Sect. 3, we analyzed the local and global stability of the system by applying the Schur–Cohn criteria and using a suitable Lyapunov function, respectively. By considering the data of Table 1 and the theorems in Sect. 3, we show in Example 1 that, for a chemotherapy agent rate of $\sigma =0.5$, the interaction between the tumor cells and glial cells still exists. However, increasing this rate to $\sigma =50$ could only decrease the effectiveness of the resistant tumor cell population, while still a strong interaction between the sensitive tumor cells and glial cells continues. Thus, it was recognized that for late detection of GBM, the treatment and the interaction would be painful and with challenges. In Sect. 4, we proved that the system undergoes Neimark–Sacker bifurcation (or Hopf bifurcation for map) when the bifurcation parameter exceeds a critical value.

To consider the tumor for a low-density rate, we incorporate to the sensitive and resistant tumor formulations Allee function at time t since they are formulated as logistic equations. Analyzing the problem for an extinction case, we apply the strong Allee effect and add the threshold effect to the system. It was shown in Example 2 that for early detection of GBM in day 17, a chemotherapy agent rate of $\sigma =0.5$ could be enough to keep the sensitive tumor cells under a specific density so that it does not produce the resistance cell population.

The future direction of this paper is to establish the model as a fractional order differential equation system, and for this study some essential works such as [15,16,17,18,19,20,21,22,23,24] will be considered to highlight the theory of the study.