Nonlinear dynamics of a SIRI model incorporating the impact of information and saturated treatment with optimal control.

In this article, we propose and analyze an infectious disease model with reinfection and investigate disease dynamics by incorporating saturated treatment and information effect. In the model, we consider the case where an individual's immunity deteriorates and they become infected again after recovering. According to our findings, multiple steady states and backward bifurcation may occur as a result of treatment saturation. Further, if treatment is available for all, the disease will be eradicated provided R 0 < 1 ; however, because limited medical resources caused saturation in treatment, the disease may persist even if R 0 < 1 . The global stability of the unique endemic steady state is established using a geometric approach. We also establish certain conditions on the transmission rate for the occurrence of periodic oscillations in the model system. Among nonlinear dynamics, we show supercritical Hopf bifurcation, bi-stability, backward Hopf bifurcation, and double Hopf bifurcation. To illustrate and validate our theoretical results, we present numerical examples. We found that when disease information coverage is high, infection cases fall considerably, and the disease persists when the reinfection rate is high. We then extend our model by incorporating two time-dependent controls, namely inhibitory interventions and treatment. Using Pontryagin's maximum principle, we prove the existence of optimal control paths and find the optimal pair of controls. According to our numerical simulations, the second control is less effective than the first. Furthermore, while implementing a single intervention at a time may be effective, combining both interventions is most effective in reducing disease burden and cost.

Introduction

Infectious disease models had been to the core of scientific studies in the last century, and researchers have used them for various diseases, such as tuberculosis (TB) [1, 2], human immunodeficiency virus (HIV) [3, 4], influenza [5, 6], malaria [7-9], and the recent one COVID-19 [10, 11]. Infectious diseases outbreaks have always heavily impacted the social and economic well-being of the affected countries and have proved to be havoc, especially for third world countries. Therefore, the cessation and control of such diseases are one of the top priorities worldwide specifically for the countries with poor medical/health care or those having limited resources. In this scenario, mathematical models are extremely helpful for the proper understanding and exploration of such diseases’ dynamics, their prevention, and control. Of late, there is a growing interest among researchers for diseases which show a peculiar property of reinfection. There are many infectious diseases in the community for which individuals become reinfected even after fully recovering and taking treatment. The study of reinfection epidemic models is useful in controlling the prevalence of these diseases. One of the widely used mathematical models is SIR [12] which is used for the disease analysis in terms of its propagation and evolution by dividing the total population into three compartments: susceptible (S), infective (I), and recovered (R). As there is possibility of reinfection which indicates future probability for the recovered individuals still being vulnerable to reinfection, SIRI models are immediate possible extensions to the SIR models.

A closed demographic SIRI model for herpes viral infections considering usual bilinear incidence was first developed and analyzed by Tudor [13] in 1990. Moreira and Wang [14] extended the work of Tudor [13] and considered a SIRI model with the nonlinear incidence rate. Van Den Driessche and Zou [15] proposed a SIRI modeling approach for the relapsing phenomenon in infectious diseases using an integro-differential equation. In [16], Castillo et al. considered the SIRI model for the dynamics of tobacco and alcohol use like a contagious disease with relapse. Leon [17] established the global stability conditions of two epidemiological models with relapse having bilinear and standard incidence rates. Naz et al. [18] presented a SIRI and a tuberculosis model with demographic progress and then analyzed the exact solutions and first integrals for these models. Martins and Pinto [19] considered a SIRI epidemic model to analyze the influence of awareness programs on an individual’s decision for vaccination. They introduced the evolutionary vaccination strategy for the model and established its bi-stability. Recently, Buonomo [11] studied a SIRI model incorporating behavioral changes toward information-based vaccination. Many other SIRI models have been reported in the literature considering various aspects; see for example [20-33] and the references therein.

Medical treatment, isolation, quarantine, vaccination, etc., are the commonly used disease control methods. It is important to note that treating all infected individuals or vaccinating all susceptible individuals is a challenging task, especially for the countries with limited medical resources or awareness. Thus, when a disease outbreak occurs in a region, not all infected individuals receive medical treatment when hospitalization reaches its capacity. A saturated treatment function in the mathematical modeling has been used in the literature to represent this situation. In 2004, Wang and Ruan [34] examined an epidemic model with a constant treatment rate of infective and studied the effect of medical resources on disease progression. In this case, they showed that model system possesses homoclinic bifurcation, subcritical Hopf bifurcation, and saddle-node bifurcation. Wang [35] modified the constant treatment rate and considered the non-smooth treatment rate according to available infective individuals. The author observed backward bifurcations in the system due to insufficient treatment capacity and found that system possesses bistable endemic equilibria. The author also showed that is only necessary and not sufficient to eradicate the disease. Recently, there have been few more studies appeared in the literature involving different types of saturated treatment function in the epidemic models; for details, one can see [36-43] and references therein. Most of these studies exhibit rich nonlinear dynamics of the system due to saturation in treatment.

In infectious disease modeling, the incidence rate plays a significant role to ensure that the model gives the appropriate dynamics of disease progression, evolution, and spread. In classical epidemiological models, the standard and bilinear incidence rates have been frequently used. Recently, many researchers are showing their interest in other forms of nonlinear incidence rates and the nonlinear dynamics in epidemiological models having multiple endemic steady states. Capasso and Serio [44] generalized the Kermack–McKendrick deterministic model [12, 45] and incorporated a saturated incidence rate which represents the inhibition effect due to the behavioral changes of susceptible in the presence of the infective in the population. For large I,  (saturation level), and for small I, . Liu et al. in [46] proposed a more general nonlinear incidence rate and later on many authors used this which can be found in [47-52]. Xiao and Ruan [53] studied an epidemic model with non-monotonic incidence rate and established the stability results for the model considered. The non-monotonicity of incidence rate function g(I) corresponds to the possible psychological impacts, that is, g(I) is increasing and decreasing for small and large I,  respectively. One can understand that the force of infection is reduced as I increases because for a larger number of infected individuals present in the population, the healthy individual may reduce the interactions per unit time. This may seem reasonable as in low density infective infection rate may be higher but once population has information about the disease prevalence, it will self-induce behavioral changes and try to reduce contacts which will eventually reduce the infection rate. In this work, we have chosen this non-monotonic incidence rate.

As we mentioned earlier, diseases not only cause mortality and morbidity but also have their economic consequences. Gupta et al. [54] studied the economic impact of severe acute respiratory syndrome (SARS) in Toronto city, and similarly in [55, 56] the authors analyzed the impact of COVID-19 on the tourism industry. Moreover, there is associated financial cost involved in implementing any intervention to control the spread or eradicate a disease. Using one’s resources in optimal way is of utmost importance for policymakers, and optimal control theory may provide suitable answers to these questions. The choice and suitability of control strategies may vary from one disease to another depending on behavior and availability of interventions. For example, Kumar and Srivastava [57] used treatment and vaccination as controls in an SVIR model; they observed that among all strategies, the simultaneous use of vaccination and treatment is most effective and less expensive. In [58-60], the authors studied optimal control for different TB models with exogenous reinfection. Hattaf et al. [58] observed that the prevention of exogenous reinfection is much effective in reducing disease prevalence, and the authors in [59, 60] studied how the control measures should be implemented to reduce the number of actively infected individuals while controlling the cost of implementing the intervention. For our study, we consider isolation with precautions and medical treatment as two interventions and study their impact on the dynamics of the disease. We also find the optimal control profile of both the interventions with respect to predefined cost which should be minimal along with disease load.

The layout of the article is as follows. In Sect. 2, we propose a SIRI model with treatment where the incidence rate is non-monotonic function. The stability analysis of disease-free and unique endemic steady states is performed. Further, we obtain the parametric conditions for the existence of transcritical and backward bifurcation. We also observe that system undergoes Hopf bifurcation and exhibits the existence of periodic oscillation under certain conditions. In Sect. 3, we provide numerical examples to analyze and further illustrate the significance of the theoretical results obtained in Sect. 2. In Sect. 4, we define our optimal control problem by introducing two time-dependent control variables to the original model. The existence and characterization of the optimal controls are derived with respect to a cost functional, and the results are discussed with the help of numerical simulations. Finally, we present the concluding remarks and discussion.

Mathematical model

We partition the total population (N) into susceptible (S), infective (I), and recovered (R) subclasses. For our model, we consider the non-monotonic incidence rate f(S, I) defined aswhich is a generalized version of the standard and bilinear incidence rate incorporating the inhibitory effects of healthy individuals due to the presence of the disease [53, 61]. Treatment is an important aspect to control or eradicate infectious diseases. Hence, we employ h(I) as the treatment rate function in our model, which is defined as

Parameters and b are saturation constants, i.e., indicates the information-induced sensitivity of individuals to the level of infection, b denotes the extent of infective being delayed for treatment, and a is the treatment rate.

We propose the following SIRI model with an inhibitory effect on transmission and saturated treatment. We also consider reinfection of recovered individual which follows mass action. The schematic diagram of the model system is shown in Fig. 1. The mathematical representation of the model is given by the following nonlinear ordinary differential equations:with initial conditions Here, is an inflow rate in susceptible population, is a natural mortality rate, is the per individual natural recovery rate of the infective individual, is the transmission rate of infection, and is the rate of reinfection for recovered individuals. Here, is a reduction parameter for the actual infection rate because the recovered individual will have reduced susceptibility due to acquired immunity ( denotes full immunity). All the considered parameters are nonnegative.

Fig. 1

Schematic flow diagram of the disease model with treatment and information

Positivity and boundedness

From the model system (1), we render the following:Note that all the rates are nonnegative on the boundary plane of the nonnegative cone of As the vector fields’ directions are inward on the boundary planes, so if we start with any interior part of this cone then the solutions will remain within this cone for all the future time. This ensures the positivity of all the populations in the system.

Theorem 1

The set is positively invariant and global attractive set of the model system (1).

Proof

To prove this, we note that the total population On differentiating, we get By solving this rate equation, we obtain the solution N(t) asIf then It means that all solutions with initial conditions in remain in for all future times and hence the set is positively invariant for the model system (1). If then from (2), we observe that Hence, we conclude that the set is globally attractive for the model system (1).

Remark 1

Existence of global attractor implies existence of absorbing set.

Existence of disease-free steady state and the basic reproduction number

Equating the right-hand side of the model system (1) to zero, we obtain the disease-free steady state which always exists.

Now, using next-generation matrix method [62], the basic reproduction number for the model system (1) is given as

Sensitivity analysis of

As the qualitative behavior of the system is -dependent, listing the involved parameters in order of their impact on the value of is worthwhile. This investigation will help in finding out the most effective parameters pertaining to the disease control. We define the parameter spaceWe can easily verify that Now, following Chitnis et al. [63], we give the following definition.

Definition 1

[63] The normalized forward sensitivity index of with respect to a parameter say is defined as

This definition gives an assessment of the normalized change in when one parameter changes its value while the other parameters remain still. The analytical formulations for each parameter’s sensitivity index are intricate and seldom ever helpful in revealing any information. Positive indexing of a parameter represents that is increasing with respect to that parameter, whereas a negative indexing indicates that is decreasing. Using the expression of and Eq. (3), we obtainThe sensitivity indices are obtained by taking the initial values mentioned in Table 1. The index indicates that decreases by 12.462% with an increase of 10% in the -value. From Fig. 2 and Table 1, we conclude that , , and are very influential parameters, of which and have maximum positive correlation with whereas has maximum negative correlation with However, apart from and a also have considerable negative correlation with . The sensitivity indices from Table 1 show that the epidemic threshold quantity is more sensitive to four key parameters and a. Therefore, we provide three contour plots for the variation of with arguments as and as shown in Fig. 3a–c, respectively.

Table 1

Initial value and sensitivity indices of the parameters associated with

Parameters	Initial value	Sensitivity indices
β	0.1	+1
Λ	0.99	+1
μ	0.098	−1.2462
γ	0.1	−0.2513
a	0.2	−0.5025

Fig. 2

Sensitivity plot for with respect to the associated parameters

Fig. 3

Contour plot of as a function of transmission rate and the parameters and a, in (a), (b), and (c), respectively

Existence of endemic steady states

In this section, we find the endemic steady state of the model system (1). Here, and is a positive real root of the following equation:where and

The coefficients A and B are always negative and depending on for a positive set of parameters. Also, we observe that the coefficients C and D may be positive or negative based on parameter values. Thus, f(I) may have more than one positive real roots. Further, from Eq. (4), for we obtain and and hence the continuity property ensures the existence of a positive such that Thus there exists at least one positive endemic steady state for Various possibilities of endemic steady states are given in Table 2, using Descartes’ rule of the sign.

Table 2

Maximum number of positive real roots

Coefficients	R_0<1	R_0>1
A	----	----
B	----	----
C	-++-	-++-
D	-+-+	-+-+
E	----	++++
Maximum number of positive real roots	0   2   2   2	1   1   3   1

Theorem 2

When or or the model system (1) has a unique endemic steady state for

When or or the model system (1) has maximum two endemic steady states and for

When the model system (1) has maximum three endemic steady states and for

Local stability of steady states

Theorem 3

The disease-free steady state of the model system (1) is locally asymptotically stable when and is unstable when .

The Jacobian matrix of the model system (1) at the disease-free steady state is given byThe characteristic equation of is given asClearly, the eigenvalues are and The eigenvalue is negative (positive) for Hence, when the disease-free steady state is locally asymptotically stable and is unstable for

Theorem 4

When or or and the model system (1) consists of a unique endemic steady state and it is locally asymptotically stable whenever and where and are defined in the proof.

The Jacobian matrix of the model system (1) at the endemic steady state is given byThe characteristic equation of is given asHere, and By the Routh–Hurwitz criterion, all the roots of (5) have negative real parts provided and Hence, the unique endemic steady state is locally asymptotically stable if it exists and and

Global stability of disease-free steady state

Now, we establish the global stability of the disease-free steady state using the method given by Castillo-Chavez et al. [64].

Theorem 5

The disease-free steady state is globally asymptotically stable when and

Consider and , similar to as they are defined in [64]. Hence, the model system (1) can be written asThe disease-free steady state of the system (1) is with Clearly, is globally asymptotically stable for subsystem as when . Further, we havewhere and We observe that when Further, is a Metzler matrix in as for . Hence, the theorem follows by [64].

Remark 2

From the expression of we observe that may not be nonnegative when and So, in that case, the disease-free steady state may not be globally stable; there may exist multiple endemic steady states (via backward bifurcation) for .

Global stability of endemic steady state

Lemma 1

The model system (1) is uniformly persistent, that is, there exists a positive constant c,  such that

It is clear from Theorem 3 that is unstable when Using the uniform persistence result [65], the instability of the disease-free steady state assures the uniform persistence of the system when Hence, the model system (1) is uniformly persistent.

Theorem 6

The unique endemic steady state of the model system (1), provided it exists, is globally asymptotically stable when and where is defined in the proof.

Using the geometric approach as well as notations as used in [66], we show the global stability of the unique endemic steady state We consider the following subsystem of the model system (1) asThe Jacobian matrix corresponding to the subsystem (6) is given asThe second additive compound matrix [67, 68] of J isAssuming the functionwe getDefineso that From the second equation of model system (6), we can easily findBy using the above expression in we getNow, following [66], we havewhere denotes the Lozinskii measure with respect to the norm and are matrix norms with respect to the vector norm. So, we haveUsing Lemma 1, Eq. (7) becomesChoose so that we haveIntegrating Eq. (8), we getwhich givesor,as I(t) is bounded and Hence, if Thus, the system (6) is globally asymptotically stable for that is, as Now, using the first equation of the model system (1) in the limiting form, we havethat is,After solving the above equation, we get as Hence, the theorem follows.

Direction of bifurcation at

Theorem 7

If backward bifurcation occurs at for the model system (1), i.e., multiple endemic steady states exist for

If transcritical bifurcation occurs at i.e., the disease-free steady state changes its stability from stable to unstable as exceeds unity and a unique endemic steady state appears and is stable for

We prove the theorem by using the method of Castillo-Chavez and Song [69]. We consider as the bifurcation parameter and by we get The model system (1) can be rewritten as:The Jacobian matrix at is given byThis Jacobian matrix has two negative eigenvalues and one simple zero eigenvalue at . The right and left eigenvectors corresponding to the zero eigenvalue are and respectively. Further,Now using Theorem 4.1 of [69], we note that if and hence the model system (1) undergoes backward bifurcation at Further, if then and transcritical bifurcation occurs at Hence, the theorem follows.

Impact of reinfection parameter on backward bifurcation

From the expression of defined in Theorem 7, it can be observed that is dependent on the parameter This particular parameter has a significant role to avoid the case of backward bifurcation. Note that is a decreasing function of Therefore, if we increase gradually, keeping all other parameters fixed as then the threshold value decreases strictly. This observation is shown in Fig. 4. Moreover, the model system (1) undergoes backward bifurcation when the pair lies on the upper portion of the curve . For the pair which lies on the lower side of the curve, the case of backward bifurcation does not arise and only forward bifurcation occurs. In Fig. 4, BB and TB stand for backward bifurcation and transcritical bifurcation, respectively.

Fig. 4

Plot of as a function of and dividing the region in plane where backward bifurcation occurs. BB and TB stand for backward bifurcation and transcritical bifurcation, respectively

Saddle-node (SN) bifurcation

Using Sotomayor’s theorem, we derive the transversality condition for saddle-node bifurcation taking as a bifurcation parameter for model system (1). We have seen that the occurrence and destruction of an endemic steady state depends on the parameter Let is the threshold value of for the occurrence of saddle-node bifurcation. Hence, the model system (1) has two steady states ( and ) for no steady state for and collision () of two steady states at . The point of collision is represented as Define where , and are defined asLet the Jacobian matrix has a simple eigenvalue with eigenvector , and that the transpose of the Jacobian matrix has an eigenvector to the eigenvalue Now differentiating f with respect to parameter we getAt bifurcation point, we getNow, we haveAlso, we findThen, the model system (1) experiences a saddle-node bifurcation at , if the following transversality conditions [70] are satisfied:

Existence of periodic solutions through Hopf bifurcation around an endemic steady state

Now, we further explore the nonlinear dynamics of the model system (1). Considering as a bifurcation parameter, we observe that the coefficients of Eq. (5) are a smooth function of . A simple Hopf bifurcation exists around the endemic steady state, if the following conditions hold: To show the existence of Hopf bifurcation, we use the criterion derived by Liu [71]. We obtain a pair of purely imaginary eigenvalues of J at whenever the coefficients of the characteristic Eq. (5) satisfy and

The Jacobian matrix J at the has a negative real eigenvalue and a pair of purely imaginary eigenvalues,

(transversality condition), where denotes the eigenvalue of and is the threshold bifurcation value of parameter

Now, for the second condition, we rewrite the characteristic Eq. (5) as is a pair of purely imaginary eigenvalues corresponding to Now, differentiating Eq. (13) with respect to we getAfter simplifying the above equation, we getSubstituting we getThus, if where We summarize the above discussion in the following theorem.

Theorem 8

Periodic oscillations appear in the model system (1) from the unique endemic steady state at provided and along with the transversality condition:

Numerical simulations

In this section, we analyze and explore our analytical results in the previous section using MATLAB. For this purpose, we select a hypothetical set of parameters that satisfy the stated conditions and exhibit the dynamical properties.

Example 1

In this example, we choose a set of parameters as so that For this set of parameters, it can be verified that and hence according to Theorem 2 there exists a unique endemic steady state for while no endemic steady state for We note that as defined in Theorem 7. Figure 5a depicts the existence of transcritical bifurcation as crosses unity in this range of It is observed that the endemic steady state is stable.

Fig. 5

a Existence of transcritical bifurcation b Dynamics of solution trajectories for different initial values at

Further, we fix the incidence rate (so that ) and plot the solution trajectories in Fig. 5b corresponding to different initial conditions (6, 0.9, 1.2),  (4, 0.7, 6),  (4, 0.2, 3),  (4, 0.04, 2),  and (6, 0.04, 4). We observe that if we start from any of these initial values, the solution trajectories converge to endemic steady state which implies that is stable.

Example 2

In this example, we choose a set of parameters as so that For this set of parameters, it can be verified that C changes its sign from positive to negative and and hence according to Theorem 2 there are two endemic steady states and for while a unique endemic steady state for We also note that as defined in Theorem 7. Figure 6a depicts the existence of backward bifurcation as crosses unity. It is observed that the endemic steady states are stable and are unstable for this range of

Fig. 6

a Existence of backward bifurcation b Phase portrait for and c Phase portrait for and

Further, we show the phase diagram of for two fixed values of from the above range, where in one case we have two endemic steady states and in other case we have unique steady state. First, with fixed we note that . In this case, the endemic steady state is locally stable and is unstable. Corresponding to the same the solution trajectories with initial conditions (560, 130, 25),  (120, 130, 50),  (60, 63, 25),  (350, 50, 50),  and (150, 50, 100) are plotted in Fig. 6b, which are approaching to either the disease-free steady state or the stable showing bi-stability. Further, the solution trajectories corresponding to (for which ) with initial conditions (150, 50, 100),  (500, 50, 80),  (120, 130, 50) (60, 63, 25),  and (500, 100, 80) are plotted in Fig. 6c which are approaching the unique stable endemic steady state

Example 3

In this example, we choose the set of parameters as so that For this set of parameters, it can be verified that and hence according to Theorem 2 there exist two endemic steady states for as shown in Fig. 7a. It is observed that are always unstable and changes its behavior from unstable to stable for this range of .

Fig. 7

a Existence of oscillatory solution on a branch of multiple endemic steady states for b Plot for the real part of eigenvalues corresponding to the endemic steady state where it changes its stability from unstable to stable for

At (), there exist two endemic steady states and where is unstable. The endemic steady state is unstable when and stable when At and for the characteristic equation is given bywhere and The eigenvalues of the characteristic equation at are which shows the existence of a pair of purely imaginary eigenvalues at at . We have made a plot for real part of complex eigenvalues of in the small neighborhood (0.012016224, 0.012056224) of and we see that system changes its behavior from unstable to stable. From Fig. 7b, we can see that at we have All the conditions of Theorem 8 are satisfied. Using the approach and formula given in [72], we find the first Lyapunov coefficient () for determining the direction of Hopf bifurcation. We obtain which implies that the Hopf bifurcation is subcritical and hence the unstable limit cycle bifurcating from . A periodic orbit around at is obtained which is shown in Fig. 7a inset. Further, a phase diagram of this periodic solution and the corresponding oscillatory solutions of the infective populations at are shown in Fig. 8a and b, respectively. Further, in Fig. 8c we have shown that for regains its stability.

Fig. 8

a Existence of periodic orbit around the endemic steady state at b Oscillations in infective population for corresponding to periodic orbit c Stability of infective population for

In this example, we also discuss the case of saddle-node bifurcation. We have seen that the model system has two endemic steady states at As decreases the value of parameter the two endemic steady states collide at and denoted as The Jacobian matrix has a simple eigenvalue and , are the eigenvector of J and respectively. Both the transversality conditions and are satisfied. Hence, the model system (1) experiences saddle-node bifurcation at We observe that for the considered set of parameters in Example 3, is unstable and is unstable and stable when and respectively.

Example 4

This example shows very informative and rich dynamics of the considered model. Here, we discuss the occurrence of multiple endemic steady states for , their stability behavior and the appearance of Hopf bifurcation. We choose the set of parameters as so that For this set of parameters, it can be verified that and hence according to Theorem 2 there exist either one or three endemic steady states for as shown in Fig. 9a. This figure is plotted on semi-log scale. We have unique endemic steady state for three endemic steady states and for and again unique endemic steady state for It is observed that and are always stable and unstable, respectively while changes its behavior from stable to unstable.

Fig. 9

a Existence of multiple endemic steady states for . The graph is drawn on semi-log scale. b Bifurcation diagram which shows the stability and instability as the parameter is varied c Existence of periodic orbit around the endemic steady state at

At (), the unique endemic steady state is stable when and unstable when For the characteristic equation at is given bywhere and The eigenvalues of the characteristic equation for are We have made a plot for the real part of complex eigenvalues of in the small neighborhood (0.12526624, 0.13117) of and we see that system changes its behavior from stable to unstable. From Fig. 11c, we can see that at and all the conditions of Theorem 8 are satisfied which shows the existence of Hopf bifurcation around and Theorem 8 ensures the existence of periodic solutions around the endemic steady state at Here, we have which implies that the Hopf bifurcation is supercritical and hence the stable limit cycle bifurcating from . We make a bifurcation diagram in -I-R plane in Fig. 9b in the neighborhood of to show the change in behavior and formation of closed orbits as crosses , and changes from stable to unstable. The periodic orbit around at is shown in Fig. 9c, and the corresponding oscillatory solutions of the infective populations at are shown in Fig. 11a. We also plot the corresponding two-dimensional bifurcation diagram of all three populations and R vs. in Fig. 10. Further, we consider within the region between purple vertical lines in Fig. 9a and draw the trajectories with different initial conditions, as shown in Fig. 11b. From this figure, it is clear that all trajectories converging to hence is stable and are unstable.

Fig. 11

a Oscillations in infective population at corresponding to periodic orbit b Dynamics of solution trajectories with different initial values at c Plot for the real part of eigenvalues corresponding to the endemic steady state where it changes its stability from stable to unstable for

Fig. 10

a Depiction of Hopf bifurcation diagram for the susceptible population as varies b Depiction of Hopf bifurcation diagram for the infective population as varies c Depiction of Hopf bifurcation diagram for the recovered population as varies

Hopf–Hopf bifurcation and stability switch

In the next example, we have shown the validation of Theorem 8 for a set of parameters where periodic oscillations occur around a unique endemic steady state and also stability switches occur in the model system.

Example 5

We choose the following set of parameters: so that For this set of parameters, it can be verified that and hence according to Theorem 2 there exists a unique endemic steady state for For , we get (basic reproduction number corresponding to ) and a unique endemic steady state At the corresponding characteristic equation becomeswhere and The eigenvalues of the characteristic equation at are From Fig. 12b, we observe that . We notice that loses its stability at Also, all the conditions of Theorem 8 are satisfied. At we have , so for the system undergoes a supercritical Hopf bifurcation and hence the stable limit cycle bifurcating from .

Fig. 12

a Plot for the real part of eigenvalues corresponding to the endemic steady state where it changes its stability from stable to unstable and then from unstable to stable for b Plot for the real part of eigenvalues corresponding to the endemic steady state where it changes its stability from stable to unstable and then from unstable to stable for

As we further increase bifurcation parameter we notice that the model system (1) regains its stability at At we have a unique endemic steady state and For and The eigenvalues of the characteristic equation at are From Fig. 12b, we observe that Also, all the conditions of Theorem 8 are satisfied. At we have so, for the system undergoes a supercritical Hopf bifurcation and hence the stable limit cycle bifurcating from . Hence, the Hopf bifurcation exists and the model system (1) exhibits periodic oscillations. The real part of eigenvalues corresponding to unique endemic steady state are plotted in Fig. 12a which shows the change in sign twice as varies. Figure 12b shows a plot for the real part of complex eigenvalues of in the small neighborhood (0.000962, 0.002442) of in which both and lies. Occurrence of periodic orbit around and are shown in Fig. 13b and c, respectively. Further, we consider and to show the stability of as shown in Figs. 13a and 14a, respectively. From the expression of we see that is directly proportional to so that we provide the bifurcation plot with instead of bifurcation parameter In Fig. 14b, we show the stability switching behavior of system with respect to the parameter and the corresponding two-dimensional bifurcation diagram of infective population versus is shown in Fig. 14c. Therefore, from the above discussion, we observe that at and the system goes through the Hopf bifurcation or double Hopf bifurcation leading to stability switch. Via first Hopf bifurcation, the endemic steady state loses its stability and periodic solutions bifurcate at while via second Hopf bifurcation, periodic oscillation vanishes and the system regains its stability at The case where the endemic steady state changes its behavior from unstable to stable with the increment in parameter value is named as backward Hopf bifurcation [73]. So, in this example, the occurrence of second Hopf bifurcation at is the case of backward Hopf bifurcation. Overall, the system shows stability switch via two Hopf bifurcations where in between two stable endemic steady states, the system possesses the oscillatory solutions which increase in size, reach the maximum, and then decrease in size before vanishing. The phenomenon is referred as endemic bubble and has been observed in a few other works and has been a point of attraction recently [39, 74].

Fig. 13

a Stability of infective population for b Periodic orbit around at (c) Periodic orbit around at

Fig. 14

a Stability of infective population for b Bifurcation diagram which shows the stability and instability as the parameter is varied. From stable steady state periodic solutions bifurcate and then ultimately it stabilizes again as increases c Depiction of double Hopf bifurcation diagram for the infective population as varies

Impact of information and reinfection on infective population

This section is dedicated to show the impact of information and reinfection on the force of infection for the considered model. We take three different values of and the other parameters are same as considered in Example 2. The solution trajectories corresponding to infective population are plotted in Fig. 15a, this shows for a large value of , the incidence will be reduced because of less interaction between infectives and susceptibles. In other words, the prevalence of the disease will reduce as the information toward disease increases.

Fig. 15

a Effect of different information levels on infective population b Effect of different reinfection levels on infective population

Now, to illustrate the significance of reinfection on the force of infection for the considered model, we take three different values of and the other parameters are same as taken in Example 2. The solution trajectories corresponding to infective population are plotted in Fig. 15b. This shows that the incidence will increase as reinfection increases.

From the above discussion, we conclude that the parameter is inversely proportional to the force of infection and the parameter is directly proportional to the force of infection.

Optimal control problem

This section is dedicated to the optimal control analysis of our model. By introducing two time-dependent controls and in model system (1), we get the following model representing the controlled dynamical system:with initial conditions Considering economic restrictions on healthcare system, the intervention policies have to be restricted [75]; therefore, both the controls are taken as bounded in [0, 1]:

Inhibitory interventions administered to susceptible population, which represents disease prevention efforts by implementing precautionary isolation for susceptible individuals. Increasing implies a decrease in and this ultimately reduces the force of infection f(S, I). If we take , then disease prevention is fully effective, whereas, if , we get the original version of the equation without control.

Treatment with medication administered to infective population. The coefficient with treatment rate function h(I), represents efforts of providing treatment and medication to infective individuals. The value represents full treatment efforts, while, denotes unavailability of treatment. An increase in leads to an increase in the rate of recovery due to treatment.

Cost designing

Before talking about the existence of optimal controls, let us first discuss different components of costs which we have to minimize. The objective cost functional with respect to system (14) is given byThis cost functional represents the total cost incurred, which is taken as the weighted sum of the cost components: the cost due to disease prevalence(we consider this to be directly proportional to the infective population, reflecting the social and economic burden of disease prevalence), the cost of implementing inhibitory interventions,and the cost of availing treatment is considered asThe cost associated with includes the cost of isolation, sanitization, and other precautionary measures taken to avoid infection, whereas the cost associated with includes the cost of hospitalization and treatment of active patients, professionals for their care, and the entire drug course. The second-order nonlinearity in both of these cost components represents a nonlinear increase in cost while covering a larger population with these controls. From now on, we will use and in place of and for convenience.

Our main goal is to control disease burden while keeping costs to a minimum, that is, we need to determine a pair of optimal controls () such thatwhere is the set of admissible controls. The positive coefficients in the integrand of (15), and are the cost-balancing weight constants. The choice of these weights is based on the importance and size of the respective cost components.

Existence and characterization of optimal controls

Here, we will obtain the necessary and sufficient conditions for the optimal control problem. First, we discuss the sufficient condition for the existence of optimal controls, referring to the conditions stated in Theorem 2.1 in [76], Theorem 4.1 in [77], and related corollaries. Then, we study the characterization of the optimal controls with the help of the celebrated Pontryagin’s maximum principle (PMP) [76] and necessary conditions are derived.

Existence of optimal controls

Since and are bounded, the state variables of the model system (14) are bounded in , as discussed in Section (2.1).

Theorem 9

(Sufficient Condition) Let be the right-hand side of the system (14) and be the integrand of the objective functional (15), then there exist a pair of optimal controls if the following conditions are satisfied:

is of class and there exist positive constants , , such that , , ,

The feasible solution set with controls in is non-empty and bounded,

,

The set is bounded, closed, and convex,

is convex in

Proof

Write vector aswhere and . Further, it is easy to see that all the three terms of are differentiable functions whose derivatives are continuous, i.e, class of functions and . We also havewhereandAs state variables S, I, R, and parameters are bounded, there exist positive constants , , such that , , and which means condition (1) holds.

Condition (1) ensures the existence of unique solution to the model system (14) corresponding to the pair of controls in , that is, condition (2) holds.

We can write asand hence condition (3) holds.

Condition (4) holds using the definitions of closed, bounded, and convex sets. For proving the convexity of , we have to establish the following inequalitywhere and , are two control vectors. We also haveandNow, using both of these we havethusHence, this completes the proof.

Characterization of optimal controls

For using the Pontryagin’s maximum principle (PMP), we write the Hamiltonian bywhere is known as the vector of adjoint variables which satisfieswith transversality conditions .

Theorem 10

(Necessary Condition) With existence of a pair of optimal controls and corresponding state variables , , that optimize the objective cost functional (15), there exist , and , satisfyingalong with transversality conditions, , , and . Then, the pair of optimal controls can be characterized by the following piecewise continuous functions:

It is proved by direct application of Pontryagin’s Maximum Principle for bounded controls. Simplifying the canonical equations (16) gives us the adjoint system (17) along with transversality conditions , , and . Then, by solving the optimality conditions, stated belowwith and in the admissible set , we obtain a pair of optimal controls as stated above (18). Hence, this completes the proof.

Numerical simulations and results

For the purpose of simulations and discussion of results, we solved our optimal control problem using the forward backward sweep method (FBSM). The solution process using the FBSM [76] includes solving the state system (14) with initial conditions in forward direction (time) by using the fourth-order Runge–Kutta (RK4) method and application of current state solution to solve the adjoint system (17) with transversality conditions in backward direction (time), again, using a RK4 method. Then, both the state and adjoint solutions are utilized for updating the value of controls following the characterization discussed above (18) and this process is repeated iteratively. The iterations are continued till the values of state variables, adjoint variables, and controls converge sufficiently [78].

For numerically simulating the optimal control problem, we choose the set of parameters: , , , , , , and with initial condition (4000, 800, 10). The weights for the objective functional are chosen as , and which depicts the different importance and efforts required for different controls. We are employing the controls for a time period of 130 days.

Only control, i.e., implementing only inhibitory control

We begin the discussion with application of single control . Figure 16a and b shows optimal control profile for and its effect on the infective class I(t) of the population respectively. Figure 16a shows that the control requires full efforts for almost whole time period (126 days), while, in Fig. 16b, it is clearly observed that applying only control helps us to reduce the disease burden (green curve) significantly as compared to the disease burden with no controls (red curve). If we observe the effect of applying , at the 130th day of execution, in total 540 cases are averted.

Fig. 16

a Optimal control profile for when only is applied b Plot of I(t) with only control and without controls

Only control, i.e., implementing treatment only

Figure 17a shows that employment of only requires full efforts for initial 60 days before it slowly reduces to zero. The effect of applying only is shown in Fig. 17b; the number of cases (green curve) declines suddenly at the start and then remains below the red curve (without controls). It shows that application of only is also fruitful in decreasing disease prevalence throughout the time period.

Fig. 17

a Optimal control profile for when only is applied b Plot of I(t) with only control and without controls

Both and controls, i.e., applying both inhibitory and treatment controls

Now, we observe the employment of both the controls and at the same time. The optimal control path for (when both controls are applied) is shown in Fig. 18a; here, employment of requires full efforts for 127 days. Similarly, Fig. 18b shows that optimal application of requires full efforts for 96 days before it reduces to zero. Figure 18c shows the effect of applying both the controls at the same time; infectives (green curve) are now reduced to a further extent. At the end of the assumed time period, we have a difference of total 619 cases between the model with controls (green curve) and without controls (red curve).

Fig. 18

a Optimal control profile for when both controls are applied b Optimal control profile for when both controls are applied c Plot of I(t) with both controls and without controls

Comparative study

Finally, here we compare all three control strategies: use of only , use of only , and simultaneous use of both , with the case of model without controls. In Fig. 19a, it is clear that the implementation of control (when both controls are used) requires slightly more efforts than that of the case when only is used. A similar depiction can also be seen for the implementation of control in Fig. 19b; this increase in the time when full efforts (controls at upper bound) required corresponds to the fact that we are using both controls at a time with same resources and economic conditions. In Fig. 19c, we can observe that simultaneous implementation of both controls (solid green curve) is the most effective strategy as compared to the single control implementation. The cost profiles shown in Fig. 19d also justify this observation; the green curve (cost when both controls are applied) shows that simultaneous use of both controls is most cost-effective among all strategies.

Fig. 19

a Optimal control profiles for when both controls are applied and when only is applied b Optimal control profiles for when both controls are applied and when only is applied c Plot of I(t) depicting the effect of various strategies d Plot of cost profiles for only , only , both and , and no controls

Conclusion

In this article, we propose and analyze an infectious disease model with reinfection and explore the dynamics of the disease by incorporating saturated treatment and effect of information. We have proposed a general SIRI model that incorporates a specific aspect of reinfection which is a typical case of many infectious diseases such as malaria, TB, HIV, influenza, and COVID-19. For our model, we have obtained the case of multiple steady states for as well as We have proved local stability results for all the obtained steady states. Moreover, we have applied a geometric approach to derive global stability of unique endemic steady state. We have seen that in Theorem 7, saturation in treatment causes backward bifurcation in the model system and also observed that disease persists even for We have provided a condition on in Theorem 8 for the occurrence of periodic oscillation in the system. We have discussed different kinds of bifurcations such as supercritical Hopf bifurcation, backward Hopf bifurcation, double Hopf bifurcation, and also bi-stability of steady states in a particular region. Further, we have shown the impact of information parameter and reinfection parameter on the force of infection numerically. Our model analysis depicts that information and treatment of disease plays an important role to understand the dynamics of infectious disease. We have constructed an optimal control problem for our basic SIRI model where inhibitory interventions and treatment with medication are assumed as controls since both of these have their own significance and easiness in employment. After formulation of objective cost functional and establishment of existence of the control paths, we have obtained the pair of optimal controls using Pontryagin’s maximum principle. Using the forward backward sweep method, we have solved our optimal control problem and discussed the outcomes numerically by assuming appropriate model parameters. We have observed that implementation of single control at a time may work well but simultaneous application of both the control interventions is most effective in minimizing the disease prevalence and also the cost of implementation.