IntroductionTonsillitis, the inflammation of the palatine tonsils, is a highly prevalent upper respiratory tract infection among children and adolescents worldwide that imposes a substantial socioeconomic burden1,2,3 and transmitted via respiratory droplets with a typical incubation period of two to five days4, the condition presents clinically with a sore throat, fever, and swollen cervical lymph nodes5. Its etiology can be viral or bacterial, with the most significant bacterial pathogen being Streptococcus pyogenic (Group A Streptococcus), which accounts for an estimated 15–30% of related pharyngitis cases in children and 5–15% in adults1,2. The disease manifests in several forms: acute, with symptoms lasting up to two weeks1,2 recurrent, defined by multiple episodes a year that cause considerable morbidity by impacting school and work attendance6; and chronic, where persistent symptoms can diminish quality of life and create a public health concern7. As no vaccine is available, prevention relies on non-pharmaceutical interventions like good hygiene4. Management of established infections depends on the cause, with antibiotics being standard for bacterial cases, though effectiveness is contingent on local sensitivity patterns2,3 and in cases of frequent, severe recurrence, a tonsillectomy may be recommended to reduce long-term morbidity6.Mathematical modeling has become an indispensable tool for understanding the complex dynamics of infectious diseases and evaluating public health interventions8,9,10,11. Researchers successfully apply these quantitative frameworks to link biological transmission mechanisms with population-level infection dynamics12, thereby guiding epidemic control strategies13. Similar to the present study on tonsillitis, many of these models employ classical ordinary differential equations (ODEs) as a robust framework to describe the flow of individuals between health states14,15, a method that has been effectively applied to diverse diseases like tuberculosis16, Hepatitis B17,18, and pneumonia19. However, the field has also advanced to incorporate greater complexity. For instance, some models have explored intricate epidemic properties such as bi-modal20 and bi-susceptibility dynamics21 in COVID-19, or have included structural complexities like treatment and delay factors22,23. Furthermore, to capture memory effects in disease progression, some researchers have enhanced their models with fractional-order derivatives, including Atangana-Baleanu24,25 and Caputo-Fabrizio operators26 a level of complexity the current tonsillitis study acknowledges but approaches with a classical ODE framework. Despite this extensive body of work on a wide range of pathogens27, the modeling of tonsillitis has been largely overlooked, highlighting the novelty of the present study in adapting this well-established ODE-based framework to address a critical gap in infectious disease epidemiology.A prominent and powerful technique used in conjunction with these models is optimal control theory, rooted in Pontryagin’s Maximum Principle28,29, which is used to determine the most cost-effective allocation of intervention resources. This approach has been widely applied to design optimal strategies for managing major infectious global diseases, including HIV co-infections30, Ebola31, measles32, rubella33, and saw numerous applications during the COVID-19 pandemic25,34. A central and recurring conclusion across these studies, which strongly aligns with the findings of the tonsillitis model, is the superior efficacy of integrated or combined control strategies over single-intervention approaches, with investigations consistently demonstrating that a multi-pronged approach targeting different aspects of the transmission cycle yields the most significant reduction in disease burden20,24,28,30,32,35 While the present study employs the robust fourth-order Runge-Kutta method for its numerical simulations, other works in the reference list highlight a focus on enhancing the computational execution of these models. Researchers have dedicated efforts to developing and implementing computationally efficient numerical approaches specifically for complex optimal control problems in the context of the COVID-19 pandemic35,36,37 and Rubella33. Ultimately, despite the varied pathogens and increasing methodological sophistication, the broad body of work validates the use of optimal control to design effective, multi-faceted public health policies18,25,38, reinforcing the central conclusion of this tonsillitis study: that an integrated assault on the disease is unequivocally the most effective strategy.This study pioneers the application of compartmental modeling to the transmission dynamics of tonsillitis, addressing a critical gap in a field traditionally dominated by clinical and bacteriological research1,6. Our primary objective is to construct a robust quantitative framework that moves beyond descriptive studies2,3 to predict community-level spread and identify optimal intervention pathways. To achieve this, we introduce a novel, well-posed six-compartment model and conduct a multi-layered analytical investigation. We first establish the fundamental epidemic dynamics by deriving the basic reproduction number and rigorously proving the local and global stability of the system’s equilibria. A key contribution is our bifurcation analysis using center manifold theory, which confirms a forward bifurcation at the model basic reproduction number is one, thereby providing a clear, unambiguous target for eradication efforts. Furthermore, we perform a comprehensive sensitivity analysis to identify the most influential parameters such as the transmission and chronic treatment rates pinpointing the most effective levers for public health intervention. Building on this analytical foundation, the central contribution of this work is the formulation and analysis of an optimal control problem. By leveraging Pontryagin’s Maximum Principle with three time-dependent controls proactive prevention \(\:(w_1)\), acute treatment \(\:(w_2),\) and chronic therapy \(\:(w_3)\) our findings reveal a crucial insight: an integrated strategy that simultaneously suppresses new transmissions while clearing individuals from persistent infectious reservoirs is unequivocally superior to any single-intervention approach. Ultimately, this research provides public health stakeholders with a novel, evidence-based tool, validated by numerical simulations, to transition from reactive treatment towards proactive, integrated policies capable of mitigating the significant socioeconomic burden of tonsillitis.The format of the remaining part of this study is as follows. In “Compartmental model formulation”, we introduce the model’s formulation, including the six-compartment structure, the underlying assumptions, the flow diagram, and the governing system of ordinary differential equations. “Qualitative analysis of the model” is dedicated to the qualitative analysis of the model, where we prove the positivity and boundedness of solutions, derive the basic reproduction number, and establish the local and global stability of the disease-free equilibrium. A sensitivity analysis is conducted in what constitutes “Methods of sensitivity analysis and discussions”, identifying the most influential parameters on the basic reproduction number and visualizing their impact on disease transmission. Following this, “Optimal control problem analysis ” details the theoretical framework of the optimal control problem, applying Pontryagin’s Maximum Principle to characterize the intervention strategies. “Numerical simulations and discussions” presents the numerical simulations, which first validate the model’s convergence towards its stable states and then provide a comparative evaluation of the different single, dual, and triple control scenarios to determine their effectiveness. Finally, “Discussions, conclusions and future directions” offers the study’s conclusions, summarizing the key findings on the most effective intervention strategies and outlining potential directions for future research.Compartmental model formulationIn this study, we need to formulate a model with optimal control theory to capture the transmission dynamics of tonsillitis within a community. The formulated model is built upon a system of ordinary differential equations (ODEs), which is a standard and robust approach for tracking the flow of individuals between different health states over time. The total human population denoted by \(\:M\left(t\right)\) at time\(\:\:t\ge\:0\) is partitioned into five distinct, mutually exclusive compartments based on an individual’s disease status: Susceptible \(\:S\left(t\right),\) Acutely Infected \(\:A\left(t\right),\:\)Chronically Infected \(\:C\left(t\right),\) Treated \(\:T\left(t\right),\) and Recovered \(\:R\left(t\right)\) such that$$\:M\left(t\right)=\:S\left(t\right)+A\left(t\right)+C\left(t\right)+T\left(t\right)+R\left(t\right).$$(1)The flow of individuals between these compartments is governed by a set of carefully defined parameters representing biological and social processes. New individuals enter the population through a constant recruitment rate, \(\:\varLambda\:\) and are initially placed in the susceptible \(\:S\left(t\right)\) compartment. Susceptible individuals become infected through contact with infectious individuals, specifically those in the Acute \(\:A\left(t\right)\:\)and Chronic \(\:C\left(t\right)\:\)compartments. This transmission process is captured by the standard incidence force of infection, \(\:{\lambda\:}_{T}\), defined by:$$\:{\lambda\:}_{T}\left(t\right)=\frac{\beta\:}{M}\left(\:C\left(t\right)+\vartheta\:A\left(t\right)\right)$$(2)which is dependent on the transmission rate \(\:\beta\:\) and the relative infectiousness of acute cases\(\:\:\vartheta\:.\:\) Upon infection, a fraction \(\:(1-q)\) of individuals enters the Acute \(\:A\left(t\right)\) stage, characterized by a short-term, highly inflammatory condition, while the remaining fraction \(\:\left(q\right)\) progresses directly to the chronic \(\:C\left(t\right)\) stage. Individuals in the acute \(\:A\left(t\right)\) compartment can transition to the chronic \(\:C\left(t\right)\) compartment at a rate\(\:\:\delta\:\). Treatment is modeled as a pathway originating from the chronic C(t) compartment, where individuals move to the Treated T(t) compartment at a rate \(\:\gamma\:\), representing the initiation of medical intervention for persistent infections. Following treatment, individuals move from the treated \(\:T\left(t\right)\) to the Recovered \(\:R\left(t\right)\) compartment at a rate\(\:\:\sigma\:.\) A key feature of the model is the inclusion of waning immunity, where recovered \(\:R\left(t\right)\) individuals can lose their immunity at a rate \(\:\omega\:\) and return to the susceptible S(t) compartment, allowing for reinfection. Finally, all compartments are subject to a natural death rate \(\:\left(d\right),\) and the chronically infected \(\:C\left(t\right)\) compartment faces an additional disease-induced death rate\(\:{\:d}_{T}\). This comprehensive structure allows the model to simulate the complex interplay between infection, progression, treatment, and immunity that characterizes tonsillitis dynamics at the population level.The formulation of the tonsillitis transmission model is based on the following set of basic assumptions, which define the biological and demographic framework for the system of differential equations: The population is considered variable, with new individuals entering the susceptible class \(\:S\left(t\right)\) at a constant rate \(\:\varLambda\:\), representing births and immigration, all individuals, regardless of their disease state, are subject to a natural death rate \(\:d\) only individuals in the chronic infection compartment \(\:C\left(t\right)\) are assumed to experience an additional disease-induced mortality rate \(\:{d}_{T}\) this reflects the higher risk of complications and long-term health decline associated with persistent tonsillitis compared to the acute phase, the model assumes a homogeneously mixed population, meaning that every individual has an equal probability of coming into contact with any other individual, factors like age, geography, and social structure are not explicitly considered, both acutely infected \(\:A\left(t\right)\) and chronically infected \(\:C\left(t\right)\:\)individuals are capable of transmitting the disease to susceptible individuals, the model employs a standard incidence rate, which is more realistic than mass action in large populations where the number of contacts per person does not necessarily increase with population size, acutely infected individuals \(\:A\left(t\right)\) may have a different level of infectiousness compared to chronically infected individuals \(\:C\left(t\right),\) which is accounted for by the modification parameter \(\:\vartheta\:\:\)a susceptible individual transitions into one of two infectious states: a fraction \(\:(1-q)\) enters the acute stage \(\:A\left(t\right),\) while the remaining fraction \(\:q\) progresses directly to the chronic stage \(\:C\left(t\right),\) individuals in the acute compartment A(t) do not recover directly but instead progress to the chronic compartment \(\:C\left(t\right)\) at a rate \(\:\delta\:\), treatment is assumed to be sought primarily by individuals with chronic tonsillitis and hence individuals move from the chronic compartment \(\:C\left(t\right)\) to the treated compartment \(\:T\left(t\right)\) at a rate \(\:\gamma\:\) following successful treatment, they move from the treated \(\:T\left(t\right)\) to the recovered \(\:R\left(t\right)\) compartment at a rate \(\:\sigma\:\) immunity acquired after recovery is not permanent, recovered individuals R(t) lose their immunity at a rate \(\:\omega\:\) and return to the susceptible compartment \(\:S\left(t\right)\) making them eligible for reinfection.Based on the descriptions and model assumptions stated above, the population flow diagram that illustrates the dynamics of tonsillitis infectious disease transmission is provided in Fig. 1 below.Fig. 1Population flow diagram of the tonsillitis infection spreading dynamics where \(\:{\varvec{\lambda\:}}_{\varvec{T}}\) is described in the Eq. (2).Full size imageAccording to the population flow diagram described by Fig. 1 above, the mathematical model for the tonsillitis infectious disease transmission dynamics is represented by the following system of ordinary differential equations.$$\:\frac{dS}{dt}={\Lambda\:}+{\upomega\:}\text{R}-\left({\lambda\:}_{T}+d\right)S,$$$$\:\frac{dA}{dt}=(1-q){\lambda\:}_{T}S-(\delta\:+d)A,$$$$\:\frac{dC}{dt}=q{\lambda\:}_{T}S+\delta\:A-(\gamma\:+d+{d}_{T})C,$$(3)$$\:\frac{dT}{dt}=\gamma\:C-(\sigma\:+\text{d})T,$$$$\:\frac{dR}{dt}={\upsigma\:}\text{T}-(\omega\:+d)R,$$with initial data given by \(\:S\left(0\right)\ge\:0\),\(\:\:A\left(0\right)\ge\:0\),\(\:\:C\left(0\right)\ge\:0\), \(\:T\left(0\right)\ge\:0\), and \(\:R\left(0\right)\ge\:0\).Qualitative analysis of the modelBasic properties of the modelPositivity, Boundedness, Existence and Uniqueness of the Model (3) solutions: For compartmental epidemiological model to be both biologically and mathematically well posed we have to prove that the model solutions have to be non-negative bounded, exist and unique in the invariant region represented by:$$\:{\Omega\:}=\left\{\begin{array}{c}\left(S,\:A,C,T,R\right)\in\:{\mathbb{R}}_{+}^{5},N\le\:\frac{{\Lambda\:}}{d}\end{array}\right\}.$$(4)Theorem 1Given the initial conditions illustrated in Eq. (3), then the solutions \(\:\:S\left(t\right),\:A\left(t\right),C\left(t\right),T\left(t\right),,\:\) and\(\:\:R\left(t\right)\) of the tonsillitis infectious disease transmission dynamics model (3) are positive for all time \(\:\:t>0\).ProofThe theorem put out by Picard and Lindelöf and shown in the reference14 is the most often used and basic technique for proving that solutions to the initial value issue will exist and be unique when the functions involved are Lipschitz continuous locally. The functions specified on the right-hand side of model (3) are \(\:{\text{C}}^{1}\) on the topology of \(\:{\mathbb{R}}_{+}^{5}\), where the Picard–Lindelöf criterion is valid. This is the model we developed and presented in (3). Therefore, there are only three solutions for the tonsillitis model. As long as \(\:\text{y}\in\:{\left[0,{\infty\:}\right]}^{5}\:\)and\(\:\:\text{y}{\:}_{\text{i}}=\:0\), we can confirm that the tonsillitis model (3) functions represented on the right by \(\:{\text{h}}_{\text{i}}\left(S,\:A,C,T,R\right)\) are non-negative using the same methods as the research14 in claim A.1. This proved the essential assertion and is a basic characteristic in mathematical epidemiology models, where negative for the population state variables is illogical.Theorem 2The solutions on the model (3) are bounded in the region described by Eq. (4).ProofGiven the equation’s initial conditions, let \(\:\left(S,\:A,C,T,R\right)\in\:{\mathbb{R}}_{+}^{5}\) be arbitrary positive solutions of the dynamical system (3). When \(\:t\:\to\:\:\infty\:,\) use Birkhoff’s and Rota’s theorem39 to find the result provided by \(\:0\le\:M\left(t\right)\le\:\frac{{\Lambda\:}}{d}\). The sum of all differential equations shown in Eq. (3) becomes \(\:\frac{dM}{dt}\le\:{\Lambda\:}-dM\), where M is specified in (1). Thus, the region indicated by Eq. (4) is occupied by each of the positive solutions of the model (3).Invariance regionThe dynamics of the tonsillitis model described by system (3) are analyzed within a specific region of the state space that is both biologically and mathematically meaningful. This epidemiological feasible region, denoted by \(\:\varOmega\:\), is defined as:\(\:\:{\Omega\:}=\left\{\begin{array}{c}\left(S,\:A,C,T,R\right)\in\:{\mathbb{R}}_{+}^{5}:N\left(t\right)=S\left(t\right)\:+\:A\left(t\right)\:+\:C\left(t\right)\:+\:T\left(t\right)\:+\:R\left(t\right)\:\le\:\frac{{\Lambda\:}}{d}\end{array}\right\}.\)RemarkThe results of Theorem 1 (Positivity) and Theorem 2 (Boundedness) collectively establish that the region \(\:\varOmega\:\) is positively invariant under the flow of the system (3). This means that any solution trajectory that originates within\(\:\:\varOmega\:\), with non-negative initial conditions and a total population at or below the value \(\:\frac{{\Lambda\:}}{d}\), will remain within Ω for all future time \(\:t\:>\:0.\) Therefore, the model is well-posed in this region, ensuring that the population dynamics are confined to a biologically realistic domain where all state variables remain non-negative and bounded40.Theorem 3(existence and uniqueness of solutions).For any non-negative initial condition represented by \(\:X_0\:=\:\left(S\right(0),\:A(0),\:C(0),\:T(0),\:R(0\left)\right)\:\in\:\:{\Omega\:}=\left\{\begin{array}{c}\left(S,\:A,C,T,R\right)\in\:{\mathbb{R}}_{+}^{5},N\le\:\frac{{\Lambda\:}}{d}\end{array}\right\}\), there exists a unique solution for the tonsillitis model (3) given by solution \(\:X\left(t\right)\:=\:\left(S\right(t),\:A(t),\:C(t),\:T(t),\:R(t\left)\right)\) for all \(\:t\:\ge\:\:0.\).ProofLet the right-hand side of the system of Eq. (3) be represented by a vector function \(\:F\left(X\right),\) where \(\:X\:=\:(S,\:A,\:C,\:T,\:R)^{\rm T}\). The system can be written as the initial value problem \(\:\frac{dX}{dt}=F\left(X\right)\) with \(\:X\left(0\right)\:=\:X_0\). We will use the Picard-Lindelöf theorem (also known as the Cauchy-Lipschitz theorem), which guarantees the existence and uniqueness of a local solution if the function \(\:F\left(X\right)\) is locally Lipschitz continuous with respect to the state vector \(\:X\) in the domain \(\:\varOmega\:\). The vector function \(\:F\left(X\right)\:\)is given by the right hand side equations of the model (3). A sufficient condition for \(\:F\left(X\right)\) to be locally Lipschitz continuous is that it is continuously differentiable (of class \(\:C^{1}\)). We can verify this by examining the partial derivatives of each component of\(\:\:\:F\) with respect to each state variable (i.e., the Jacobian matrix of \(\:F\)). The components of \(\:F\left(X\right)\) are composed of sums, products, and quotients of the state variables. These operations result in functions that are continuously differentiable everywhere except where the denominator \(\:M\) is zero. Within our defined feasible region \(\:\varOmega\:\), the total population \(\:M\left(t\right)\) is bounded below. Since the recruitment rate \(\:\varLambda\:\) is positive, the total population \(\:M\left(t\right)\) cannot reach zero. The analysis of \(\:\frac{dM}{dt}=\varLambda\:\:-\:dM\:-{d}_{T}\:C\) shows that if \(\:M\) were to approach zero, \(\:\frac{dM}{dt}\) would be approximately \(\:\varLambda\:\:>\:0\), ensuring the population remains positive. Therefore, \(\:M\left(t\right)\) is strictly positive for all \(\:t\:>\:0\) within \(\:\varOmega\:\). Since \(\:M\:>\:0\), all partial derivatives \(\:\frac{\partial\:{F}_{i}}{\partial\:{X}_{j}}\) exist and are continuous functions of the state variables throughout the interior of the feasible region \(\:\varOmega\:\). This confirms that the vector field \(\:F\left(X\right)\) is continuously differentiable, and thus locally Lipschitz continuous, in Ω.Therefore, by the Picard-Lindelöf theorem, for any initial condition \(\:X_0\:\in\:\:\varOmega\:\), a unique local solution \(\:X\left(t\right)\) exists. Furthermore, as we have already shown in Theorem 2 that all solutions are bounded and remain within the compact set\(\:\:\varOmega\:\), these local solutions do not blow up in finite time. Therefore, the unique solution can be extended for all \(\:t\:\ge\:\:0\). This completes the proof that the model system (3) is well-posed.Disease-free equilibrium and local stabilityBy setting the right-hand side of the dynamical system (3) to zero so that \(\:A=C=T=R=0\), the tonsillitis disease transmission dynamics model (3) disease-free equilibrium point is found and following a few computation and simplification processes, it is provided by$$\:{E}_{T}^{0}=\left({S}^{0},{A}^{0},{C}^{0},{T}^{0},{R}^{0}\right)=\left(\frac{{\Lambda\:}}{d},\text{0,0},\text{0,0}\right)$$.Basic reproduction numberDefinitionFor the tonsillitis infection model (3), the basic reproduction number, represented by \(\:{\mathcal{R}}_{T}\), is the average number of secondary cases that a typical infective individual would create in a fully susceptible community30,34.We can determine the basic reproduction number \(\:{\mathcal{R}}_{T}\) of our system (3) by applying the next-generation matrix approach41. This method emphasizes the importance of distinguishing new infections from all other transitions within the population to compute the basic reproduction number. The infected categories are \(\:A\) and \(\:C\). We can express system (3) as: \(\:\dot{x}\mathcal{\:}=\mathcal{\:}\mathcal{F}\left(x\right)-\mathcal{\:}\mathcal{V}\mathcal{\:}\left(x\right)\), \(\:\mathcal{\:}\mathcal{V}\mathcal{\:}=\mathcal{\:}{\mathcal{V}}^{-}-{\mathcal{V}}^{+}\), and \(\:x\:=\:(S,\:\:A,\:\:C,\:\:T,\:\:R\:)\). Here, \(\:\mathcal{F}\) represents the rate at which new infections appear in each category, (\(\:\mathcal{V}\)) is the rate of movement into each category via all other processes, and \(\:{\mathcal{V}}^{+}\:\)indicates the rate at which infectious individuals transition out of each category. Utilizing the differential equations from system (3), the corresponding matrices for the new infection components, \(\:\mathcal{F}\left(x\right)\) and the other transition components, \(\:\mathcal{V}\left(x\right)\), are provided as follows:$$\:\mathcal{F}\left(x\right)=\left(\begin{array}{c}(1-q)\frac{\beta\:}{M}\left(C+\vartheta\:A\right)S\\\:q\frac{\beta\:}{M}\left(C+\vartheta\:A\right)S\\\:0\\\:0\\\:0\end{array}\right),$$and$$\:\mathcal{V}\left(x\right)=\left(\begin{array}{c}\left(\delta\:+d\right)A\:\\\:\left(\gamma\:+d+{d}_{T}\right)C-\delta\:A\\\:\left(\sigma\:+d\right)T-\gamma\:C\\\:\left(\omega\:+d\right)R-\sigma\:T\\\:dS-\omega\:-{\Lambda\:}\end{array}\right).$$(5)Then, evaluating the partial derivatives of (5) and considering in mind that the model (3) has two infected classes, namely \(\:\:A\:\)and\(\:\:C\), we obtain$$\:F=\left[\begin{array}{cc}(1-q)\frac{\beta\:\vartheta\:S}{M}&\:(1-q)\frac{\beta\:S}{M}\\\:\frac{q\beta\:\vartheta\:S}{M}&\:\frac{q\beta\:S}{M}\end{array}\right].$$At the disease-free equilibrium point, we have \(\:S\approx\:M\). Thus$$\:F=\left[\begin{array}{cc}\left(1-q\right)\beta\:\vartheta\:&\:(1-q)\beta\:\\\:q\beta\:\vartheta\:&\:q\beta\:\end{array}\right].$$Similarly, the partial derivatives of (5) with respect to \(\:A\:\)and\(\:\:C\) at \(\:{E}_{0}\) gives$$\:V=\left[\begin{array}{cc}\delta\:+d&\:0\\\:-{\updelta\:}&\:\gamma\:+d+{d}_{T}\end{array}\right].$$To get \(\:{V}^{-1}\), we use the adjoint matrix method.$$\:{V}^{-1}=\frac{1}{\text{det}\left(V\right)}adj\left(V\right)\:\text{w}\text{h}\text{e}\text{r}\text{e}\:{det}\left(V\right)=\left|\begin{array}{cc}\delta\:+d&\:0\\\:-{\updelta\:}&\:\gamma\:+d+{d}_{T}\end{array}\right|=\left(\delta\:+d\right)\left(\gamma\:+d+{d}_{T}\right).$$Then we have determined.\(\:\Rightarrow\:{V}^{-1}=\left[\begin{array}{cc}\frac{1}{\left(\delta\:+d\right)}&\:0\\\:\frac{\delta\:}{\left(\delta\:+d\right)\left(\gamma\:+d+{d}_{T}\right)}&\:\frac{1}{\left(\gamma\:+d+{d}_{T}\right)}\end{array}\right]\) and$$\:F{V}^{-1}=\left[\begin{array}{cc}\left(1-q\right)\beta\:\vartheta\:&\:(1-q)\beta\:\\\:q\beta\:\vartheta\:&\:q\beta\:\end{array}\right]\left[\begin{array}{cc}\frac{1}{\left(\delta\:+d\right)}&\:0\\\:\frac{\delta\:}{\left(\delta\:+d\right)\left(\gamma\:+d+{d}_{T}\right)}&\:\frac{1}{\left(\gamma\:+d+{d}_{T}\right)}\end{array}\right],$$$$\:\Rightarrow\:F{V}^{-1}=\left[\begin{array}{cc}\frac{\left(1-q\right)\beta\:\vartheta\:}{\left(\delta\:+d\right)}+\frac{(1-q)\beta\:\delta\:}{\left(\delta\:+d\right)\left(\gamma\:+d+{d}_{T}\right)}&\:\frac{(1-q)\beta\:}{\left(\gamma\:+d+{d}_{T}\right)}\\\:\frac{q\beta\:\vartheta\:}{\left(\delta\:+d\right)}+\frac{q\beta\:\delta\:}{\left(\delta\:+d\right)\left(\gamma\:+d+{d}_{T}\right)}&\:\frac{q\beta\:}{\left(\gamma\:+d+{d}_{T}\right)}\end{array}\right].$$Then, we need to find the eigenvalues of the matrix \(\:F{V}^{-1}\) with the characteristic equation\(\:\:\left|F{V}^{-1}-\lambda\:I\right|=0\) as$$\:\Rightarrow\:\left|\begin{array}{cc}\frac{\left(1-q\right)\beta\:\vartheta\:}{\left(\delta\:+d\right)}+\frac{(1-q)\beta\:\delta\:}{\left(\delta\:+d\right)\left(\gamma\:+d+{d}_{T}\right)}-\lambda\:&\:\frac{(1-q)\beta\:}{\left(\gamma\:+d+{d}_{T}\right)}\\\:\frac{q\beta\:\vartheta\:}{\left(\delta\:+d\right)}+\frac{q\beta\:\delta\:}{\left(\delta\:+d\right)\left(\gamma\:+d+{d}_{T}\right)}&\:\frac{q\beta\:}{\left(\gamma\:+d+{d}_{T}\right)}-\lambda\:\end{array}\right|=0,$$$$\begin{aligned}&\:\Rightarrow\:\left(\frac{\left(1-q\right)\beta\:\vartheta\:}{\left(\delta\:+d\right)}+\frac{(1-q)\beta\:\delta\:}{\left(\delta\:+d\right)\left(\gamma\:+d+{d}_{T}\right)}-\lambda\:\right)\:\left(\frac{q\beta\:}{\left(\gamma\:+d+{d}_{T}\right)}-\lambda\:\right)\\ & \quad -\frac{(1-q)\beta\:}{\left(\gamma\:+d+{d}_{T}\right)}\left(\frac{q\beta\:\vartheta\:}{\left(\delta\:+d\right)}+\frac{q\beta\:\delta\:}{\left(\delta\:+d\right)\left(\gamma\:+d+{d}_{T}\right)}\right)=0, \end{aligned}$$$$\begin{aligned}&\:\Rightarrow\:{\lambda\:}^{2}-\left(\frac{\left(1-q\right)\beta\:\vartheta\:}{\left(\delta\:+d\right)}+\frac{(1-q)\beta\:\delta\:}{\left(\delta\:+d\right)\left(\gamma\:+d+{d}_{T}\right)}+\frac{q\beta\:}{\left(\gamma\:+d+{d}_{T}\right)}\right)\lambda\:+\frac{q\beta\:}{\left(\gamma\:+d+{d}_{T}\right)}\left(\frac{\left(1-q\right)\beta\:\vartheta\:}{\left(\delta\:+d\right)}+\frac{(1-q)\beta\:\delta\:}{\left(\delta\:+d\right)\left(\gamma\:+d+{d}_{T}\right)}\right)\\ & \quad -\frac{(1-q)\beta\:}{\left(\gamma\:+d+{d}_{T}\right)}\left(\frac{q\beta\:\vartheta\:}{\left(\delta\:+d\right)}+\frac{q\beta\:\delta\:}{\left(\delta\:+d\right)\left(\gamma\:+d+{d}_{T}\right)}\right)=0, \end{aligned}$$$$\begin{aligned}&\:\Rightarrow\:{\lambda\:}^{2}-\left(\frac{\left(1-q\right)\beta\:\vartheta\:}{\left(\delta\:+d\right)}+\frac{(1-q)\beta\:\delta\:}{\left(\delta\:+d\right)\left(\gamma\:+d+{d}_{T}\right)}+\frac{q\beta\:}{\left(\gamma\:+d+{d}_{T}\right)}\right)\lambda\\ & \quad \:+\frac{q\left(1-q\right)\beta\:}{\left(\gamma\:+d+{d}_{T}\right)}\left(\frac{\beta\:\vartheta\:}{\left(\delta\:+d\right)}+\frac{\beta\:\delta\:}{\left(\delta\:+d\right)\left(\gamma\:+d+{d}_{T}\right)}\right)-\frac{q(1-q)\beta\:}{\left(\gamma\:+d+{d}_{T}\right)}\left(\frac{\beta\:\vartheta\:}{\left(\delta\:+d\right)} +\frac{\beta\:\delta\:}{\left(\delta\:+d\right)\left(\gamma\:+d+{d}_{T}\right)}\right)=0,\end{aligned}$$$$\begin{aligned}&\:\Rightarrow\:{\lambda\:}^{2}-\left(\frac{\left(1-q\right)\beta\:\vartheta\:}{\left(\delta\:+d\right)}+\frac{\left(1-q\right)\beta\:\vartheta\:\delta\:}{\left(\delta\:+d\right)\left(\gamma\:+d+{d}_{T}\right)}+\frac{q\beta\:}{\left(\gamma\:+d+{d}_{T}\right)}\right)\lambda \\ & \quad\:+\frac{q\left(1-q\right)\beta\:}{\left(\gamma\:+d+{d}_{T}\right)}\:\left(\frac{\beta\:\vartheta\:}{\left(\delta\:+d\right)} +\frac{\beta\:\delta\:}{\left(\delta\:+d\right)\left(\gamma\:+d+{d}_{T}\right)}-\frac{\beta\:\vartheta\:}{\left(\delta\:+d\right)} -\frac{\beta\:\delta\:}{\left(\delta\:+d\right)\left(\gamma\:+d+{d}_{T}\right)}\right)=0,\end{aligned}$$$$\:\Rightarrow\:{\lambda\:}^{2}-\left(\frac{\left(1-q\right)\beta\:\vartheta\:}{\left(\delta\:+d\right)}+\frac{\left(1-q\right)\beta\:\vartheta\:\delta\:}{\left(\delta\:+d\right)\left(\gamma\:+d+{d}_{T}\right)}+\frac{q\beta\:}{\left(\gamma\:+d+{d}_{T}\right)}\right)\lambda\:=0,$$$$\:\Rightarrow\:\lambda\:\left(\lambda\:-\left(\frac{\left(1-q\right)\beta\:\vartheta\:}{\left(\delta\:+d\right)}+\frac{\left(1-q\right)\beta\:\vartheta\:\delta\:}{\left(\delta\:+d\right)\left(\gamma\:+d+{d}_{T}\right)}+\frac{q\beta\:}{\left(\gamma\:+d+{d}_{T}\right)}\right)\right)=0,$$\(\:\Rightarrow\:{\lambda\:}_{1}=0\) and \(\:{\lambda\:}_{2}=\frac{\left(1-q\right)\beta\:\vartheta\:}{\left(\delta\:+d\right)}+\frac{\left(1-q\right)\beta\:\vartheta\:\delta\:}{\left(\delta\:+d\right)\left(\gamma\:+d+{d}_{T}\right)}+\frac{q\beta\:}{\left(\gamma\:+d+{d}_{T}\right)}\) thus the spectral radius of the matrix \(\:F{V}^{-1}\) is given by \(\:{\mathcal{R}}_{T}={max}\left[{\lambda\:}_{1},\:{\lambda\:}_{2}\right]={\lambda\:}_{2}\). Therefore, based on the approach41 the basic reproduction number of the model is$$\:{\mathcal{R}}_{T}=\frac{\left(1-q\right)\beta\:\vartheta\:}{\delta\:+d}+\frac{\left(1-q\right)\beta\:\vartheta\:\delta\:}{\left(\delta\:+d\right)\left(\gamma\:+d+{d}_{T}\right)}+\frac{q\beta\:}{\gamma\:+d+{d}_{T}}.$$(6)RemarkThe basic reproduction number \(\:{\mathcal{R}}_{T}\) of the tonsillitis model (3) given by Eq. (6) is a crucial epidemiological threshold, can be calculated using several established methods. The method chosen for this study is the Next-Generation Matrix (NGM) method, a systematic and widely adopted approach for compartmental models, as formalized by van den Driessche & Watmough (2002)41. The key feature of this method is its clear conceptual separation of the model dynamics into two distinct processes: the rate of new infections (represented by the matrix \(\:F\)) and the rates of all other transitions between compartments, such as progression, recovery, and death (represented by the matrix \(\:V\)). The basic reproduction number is then elegantly defined as the spectral radius (the largest eigenvalue) of the “next-generation” matrix\(\:,\:F{V}^{-1},\) this matrix epidemiologically represents the expected number of new infections in the next generation of infected individuals, making the resulting \(\:{\mathcal{R}}_{T}\) value highly interpretable. Alternative method for determining the stability threshold also exists. For instance, the survival function method calculates \(\:{\mathcal{R}}_{T}\) by integrating the product of the survival probability of an infected individual and their infectivity over their entire infectious period. For this study, the Next-Generation Matrix method was deliberately chosen for several compelling reasons. First, its structured approach is ideally suited for models with multiple infectious compartments, such as our tonsillitis model with both Acute \(\:A\left(t\right)\) and Chronic \(\:C\left(t\right)\) stages. The NGM method systematically handles the complex interactions and transitions between these states. Second, it provides a direct, explicit, and biologically meaningful formula for \(\:{\mathcal{R}}_{T}\), which breaks down the overall transmission into contributions from each infectious pathway. This interpretability is invaluable for the subsequent sensitivity analysis, allowing us to pinpoint which parameters most significantly drive disease spread. Finally, its standardization in the field of mathematical epidemiology ensures that our results are robust and directly comparable to the vast body of literature on other infectious diseases.Theorem 4The disease-free equilibrium point \(\:{E}_{0}^{T}\) of the system of ordinary differential Eq. (3) is locally asymptotically stable if \(\:{\mathcal{R}}_{T}1\).ProofInitially at\(\:\:t=0\),\(\:\:S\left(0\right)\:>0\), \(\:A\left(0\right)\:>0\), \(\:C\left(0\right)>0\),\(\:\:T\left(0\right)\ge\:0\), and \(\:R\left(0\right)\ge\:0\). The Jacobian matrix associated with the model equations at the disease-free equilibrium point \(\:{E}_{0}^{T}\:=(\frac{{\Lambda\:}}{d},\:0,\:0,\:0,\:0)\) is given by:$$\:J\left({E}_{0}^{T}\:\right)=\left[\begin{array}{ccccc}-d&\:-\beta\:\vartheta\:&\:-\beta\:&\:0&\:\omega\:\\\:0&\:\left(1-q\right)\beta\:\vartheta\:-(\delta\:+d)&\:(1-q)\beta\:&\:0&\:0\\\:0&\:q\beta\:\vartheta\:+\delta\:&\:q\beta\:-\left(\gamma\:+d+{d}_{T}\right)&\:0&\:0\\\:0&\:0&\:\gamma\:&\:-({\upsigma\:}+\text{d})&\:0\\\:0&\:0&\:0&\:\sigma\:&\:-({\upomega\:}+\text{d})\end{array}\right].$$$$\:\Rightarrow\:J=\left[\begin{array}{ccccc}-d&\:-\beta\:\vartheta\:&\:-\beta\:&\:0&\:\omega\:\\\:0&\:{b}_{22}&\:{b}_{23}&\:0&\:0\\\:0&\:{b}_{32}&\:{b}_{33}&\:0&\:0\\\:0&\:0&\:\gamma\:&\:-({\upsigma\:}+\text{d})&\:0\\\:0&\:0&\:0&\:\sigma\:&\:-({\upomega\:}+\text{d})\end{array}\right]$$$$\text{w}\text{h}\text{e}\text{r}\text{e}\:{b}_{22}=\left(1-q\right)\beta\:\vartheta\:-(\delta\:+d),\:{b}_{23}=(1-q)\beta\:,\:{b}_{32}=q\beta\:\vartheta\:+\delta\:\:\text{a}\text{n}\text{d}\:{b}_{33}=\:q\beta\:-\left(\gamma\:+d+{d}_{T}\right).$$Let us focus on the disease sub-system and since the last two rows (\(\:T\) and \(\:R\) equations) form a lower triangular block with clearly negative diagonal entries \(\:-(\sigma\:+d)\) and \(\:-(\omega\:+d)\) their eigenvalues are negative. Hence, local stability depends only on the upper-left \(\:3\times\:3\) sub-matrix involving \(\:A\) an \(\:C\) given by$$\:{J}_{3}=\left[\begin{array}{ccc}-d&\:-\beta\:\vartheta\:&\:-\beta\:\\\:0&\:{b}_{22}&\:{b}_{23}\\\:0&\:{b}_{32}&\:{b}_{33}\end{array}\right].$$Then the characteristics polynomial of \(\:{J}_{3}\) is given by \(\:\left(-d-\lambda\:\right)\left({b}_{22}-\lambda\:\right)\left({b}_{33}-\lambda\:\right)=0\) and the second degree characteristics polynomial \(\:{\lambda\:}^{2}+{a}_{1}\lambda\:+{a}_{0}=0\) where \(\:{a}_{1}=-({b}_{22}+{b}_{33})\) and\(\:\:{a}_{0}={b}_{22}{b}_{33}-{b}_{23}{b}_{32}\). Now let us apply the Routh-Hurwitz stability criteria for the second-degree polynomial \(\:{\lambda\:}^{2}+{a}_{1}\lambda\:+{a}_{0}=0\) and the disease-free equilibrium point is locally asymptotically stable if and only if \(\:{a}_{1}>0\) and\(\:{\:a}_{0}>0.\).In order to verify the above result we substituted the model parameters into\(\:\:{a}_{1}\) and\(\:{\:a}_{0}\) as \(\:{a}_{1}=-(\left(1-q\right)\beta\:\vartheta\:-(\delta\:+d)+q\beta\:-\left(\gamma\:+d+{d}_{T}\right))\:\text{a}\text{n}\text{d}\) \({a}_{0}=(\left(1-q\right)\beta\:\vartheta\:-(\delta\:+d\left)\right)(q\beta\:-\left(\gamma\:+d+{d}_{T}\right))-\left(\right(1-q\left)\beta\:\right)(q\beta\:\vartheta\:+\delta\:\:).\)Then,\(\:\:{a}_{1}=\left(\delta\:+\gamma\:+d+{d}_{T}\right)-\beta\:[\left(1-q\right)\vartheta\:+q]\) and hence \(\:{a}_{1}>0\) whenever \(\:{\mathcal{R}}_{T}0\), whenever \(\:{\mathcal{R}}_{T}1\), and hence model (3) has unique endemic equilibrium point if and only if \(\:{\mathcal{R}}_{T}>1.\)Global stability of disease-free equilibrium (DFE)The global asymptotic stability (GAS) of the disease-free state variable of the model is investigated using the theorem that is developed by Castillo-Chavez et al. (2002)42. The model is re-written as follows:$$\:\frac{dX}{dt}=F(X,\:I)$$$$\:\frac{dI}{dt}=G(X,\:I),\:G\left(X,0\right)=0,$$(7)where the components of the column-vector \(\:X\in\:{\mathbb{R}}^{m}\) denote the uninfected individuals \(\:(S,\:T,\:R,\:P)\) and the components of \(\:I\in\:{\mathbb{R}}^{n}\) denote the infected individuals \(\:(A,\:C)\). \(\:{E}_{0}^{T}\:=({X}^{0},\:0)\), denotes the disease-free equilibrium of the system (3). The equilibrium point \(\:{E}_{0}^{T}\:=({X}^{0},\:\:0)\) is a globally asymptotically stable equilibrium for this system provided that \(\:{\mathcal{R}}_{T}\:0\) by the convexity of logarithm function and it is zero if and only if for\(\:\:X={X}^{*}\) for each compartment means for the endemic equilibrium point \(\:{E}^{*}\). Now let us differentiate \(\:V\) along the solutions of the dynamical system which yields the result:$$\:\frac{dV}{dt}=\sum\:_{X\in\:\left(\text{S},\:\text{A},\text{C},\text{T},\text{R}\right)}(1-\frac{X}{{X}^{*}})\frac{dX}{dt},$$Substitute the model equations and use the endemic equilibrium conditions:$$\:{\Lambda\:}+{\upomega\:}{\text{R}}^{*}=\left({\lambda\:}_{T}^{*}+d\right){S}^{*},$$$$\:\left(1-q\right){\lambda\:}_{T}^{*}{S}^{*}=(\delta\:+d){A}^{*},$$$$\:q{\lambda\:}_{T}^{*}{S}^{*}+\delta\:{A}^{*}=(\gamma\:+d+{d}_{T}){C}^{*},$$$$\:\gamma\:{C}^{*}=(\sigma\:+\text{d}){T}^{*},$$$$\:{\upsigma\:}{\text{T}}^{*}=(\omega\:+d){R}^{*}.$$Let \(\:{\lambda\:}_{T}=\frac{\beta\:}{M}(C+\vartheta\:A)\) is used to compute \(\:{\lambda\:}_{T}^{*}=\frac{\beta\:}{M}({C}^{*}+\vartheta\:{A}^{*})\) then the derivative simplifies to a sum where inter-compartmental transitions cancel each other. We choose all Lyapunov weights \(\:{c}_{S}={c}_{A}={c}_{C}={c}_{T}={c}_{R}=1\) to ensure this cancellation, guided by the infection flow:\(\:R\to\:S\): which is balanced if \(\:{c}_{S}\)\(\:={c}_{R}\), \(\:T\to\:\:R\): balanced if \(\:{c}_{T}={c}_{R}\), \(\:C\to\:T:\) balanced if \(\:{c}_{C}={c}_{T},\) \(\:A\to\:C:\) balanced if \(\:{c}_{A}={c}_{C}\). The result gives us:\(\:\:\frac{dV}{dt}\le\:{\lambda\:}_{T}S\left[\frac{{S}^{*}}{S}-\left(1-q\right)\frac{{A}^{*}}{A}-q\frac{{C}^{*}}{C}\right]+\)negative terms and using the arithmetic mean–geometric mean \(\:(AM\--GM)\) inequality and logarithmic convexity we have the following results: \(\:x-1-lnx\ge\:0\) and \(\:(1-\frac{a}{b})(a-b)\le\:0\). Hence, each term in \(\:\frac{dV}{dt}\) is non-positive, and the whole derivative satisfies: \(\:\frac{dV}{dt}\le\:0\), with equality if and only if \(\:S={S}^{*},A={A}^{*},C={C}^{*},T={T}^{*},R={R}^{*}\), i.e., only at the endemic equilibrium.Thus, by LaSalle’s Invariance Principle44, all tonsillitis model (3) solutions trajectories in the interior of \(\:\varOmega\:\) converge to \(\:{E}^{*}\) and hence the endemic equilibrium \(\:{E}^{*}\) is globally asymptotically stable when \(\:{\mathcal{R}}_{T}>1\). This completes the proof.Methods of sensitivity analysis and discussionsThe parameter values used in this study were estimated through a synthesis of epidemiological evidence from published literature on tonsillitis spreading dynamics and reasonable assumptions informed by clinical practice, addressing the absence of primary community data. Transmission rate \(\:(\beta\:\:=\:0.50/day)\) was calibrated to achieve a basic reproduction number (\(\:{\mathcal{R}}_{T}\)) of 1.5–2.0, consistent with streptococcal pharyngitis models4,8, reflecting moderate transmissibility through respiratory droplets. Disease progression parameters were grounded in clinical evidence: acute duration \(\:(\delta\:\:\approx\:\:14\:\)days\(\:)\)1,2, chronic progression (q = 20%)6,7, and treatment timelines (σ⁻¹ = 7 days for antibiotics; \(\:\gamma\:^{-1}\:=\:30\:\)days for chronic management)2,3,6. Mortality\(\:{\:(d}_{T}=\:10^{-6}/day\)) was minimized based on modern care outcomes1,2,4, while loss of immunity (\(\:\omega\:^{-1}\:=\:365\:\)days) aligned with recurrent infection patterns4. This parameterization strategy ensures the model’s epidemiological plausibility despite data limitations, providing a robust foundation for evaluating control strategies.Table 1 Model parameters descriptions and values used for numerical simulation results.Full size tableSensitivity analysisDefinitionThe normalized sensitivity index (elasticity) of \(\:{\mathcal{R}}_{T}\) with respect to a parameter \(\:p\) is defined as:\(\:{SI}_{p}^{{\mathcal{R}}_{T}}=\frac{\partial\:{\mathcal{R}}_{T}}{\partial\:\text{p}}\times\:\frac{p}{{\mathcal{R}}_{T}}\)45In this section, we used the model parameter values stated in Table 1 above and the normalizing forward sensitivity index formula to determine the most influential model parameters on the tonsillitis infection spreading dynamics. The results are revealed by Table 2; Fig. 3 below respectively.Table 2 Sensitivity indices of for \(\:{\mathcal{R}}_{\varvec{T}}\).Full size tableThe sensitivity analysis of \(\:{\mathcal{R}}_{T}\) using the specified parameter values reveals that the transmission rate \(\:(\beta\:=0.9)\) and relative infectiousness of acute cases \(\:(\vartheta\:=0.5)\) are the dominant parameters with the highest direct positive impact on tonsillitis transmission, exhibiting sensitivity indices of \(\:+1.000\) and \(\:+0.9999\) respectively. This indicates that a 1% increase in \(\:\beta\:\) would increase\(\:{\mathcal{\:}\mathcal{R}}_{T}\) by 1%, while a 1% increase in \(\:\vartheta\:\) would raise \(\:{\mathcal{R}}_{T}\) by 0.999%, confirming their roles as primary amplifiers of transmission. The fraction progressing to chronic infection \(\:(q=0.25)\) also demonstrates moderate negative sensitive impact, indicating that interventions reducing chronic progression could moderately suppress transmission. Notably, chronic treatment rate \(\:\left(\gamma\:\right)\) has a highly indirect impact suggesting limited leverage for control via these parameters. These results prioritize prevention strategies targeting \(\:\beta\:\:\:\)(hygiene campaigns) and \(\:\vartheta\:\) (early acute case isolation) as the most efficient pathways for transmission reduction, while highlighting that reducing chronic treatment rate \(\:\left(\gamma\:\right)\) offers complementary indirect benefits.Fig. 3Simulation graph for sensitivity indices.Full size imageSimilarly, the results obtained from Fig. 3 above investigate the sensitivity indices of the parameters in the tonsillitis model (3). The transmission rate of tonsillitis, denoted by \(\:\beta\:\), and the modification factor of the contagious rate, \(\:\vartheta\:,\) are identified as the most sensitive parameters with a positive relationship to the basic reproduction number \(\:{\mathcal{R}}_{T}\). On the other hand, the tonsillitis treatment rate \(\:\gamma\:\) is also sensitive parameters but exhibit a negative relationship with \(\:{\mathcal{R}}_{T}\). According to the epidemiology of tonsillitis, in order to reduce and control the spread of the disease within the community, emphasis should be placed on minimizing the transmission rate \(\:\beta\:\) and enhancing the treatment rate\(\:\:\gamma\:\).Impacts of sensitive parameters on \(\:{\mathcal{R}}_{\varvec{T}}\)Fig. 4Effect of \(\:\beta\:\) on \(\:{\mathcal{R}}_{T}\).Full size imageThe graphical results illustrated by Fig. 4 above show the direct and linear impact of the tonsillitis transmission rate \(\:\left(\beta\:\right)\) on the basic reproduction number (\(\:{\mathcal{R}}_{T}\)), providing a clear visual contrast between the calculated transmission potential and the critical epidemic threshold. The blue curve, representing the calculated \(\:{\mathcal{R}}_{T}\), steadily rises as\(\:\:\beta\:\) increases, demonstrating that even small changes in transmission behavior directly amplify the disease’s spread. In contrast, the red line at \(\:{\mathcal{R}}_{T}=\:1\) serves as the crucial epidemiological tipping point. Below the critical\(\:\:\beta\:\) value of approximately 0.35, the blue curve lies beneath the red line, signifying a state where \(\:{\mathcal{R}}_{T}\:1\), meaning the infection has the potential to become endemic, as each case leads to more than one new infection, ensuring sustained community transmission. Therefore, the intersection of these two lines pinpoints the exact public health target: interventions must keep the transmission rate \(\:\beta\:\) below this critical value to prevent an outbreak.Fig. 5Effect of \(\:\:\vartheta\:\)on \(\:{\mathcal{R}}_{T}\).Full size imageThis graphical results described by Fig. 5 above demonstrates the influence of tonsillitis acute-phase infectiousness \(\:\vartheta\:\:\)on the overall transmission dynamics, contrasting the calculated epidemic potential against the fixed critical threshold. The blue curve, representing the calculated \(\:{\mathcal{R}}_{T}\), starts from a baseline value driven solely by chronic transmission and increases linearly, showing that as acute cases become more contagious, they contribute more significantly to the spread. This rising potential is contrasting to the constant red line at \(\:{\mathcal{R}}_{T}=\:1\), which represents the tipping point for an epidemic. Below the critical \(\:\vartheta\:\) value of approximately 0.30, the blue curve is subordinate to the red threshold; here \(\:{\mathcal{R}}_{T}\:1\) territory. This transition indicates that the added infectivity from acute cases is now sufficient to fuel a self-sustaining outbreak, turning a declining disease into an endemic one. The intersection, therefore, pinpoints the precise level of acute-phase contagiousness that public health measures must stay below to maintain control.Fig. 6Effect of \(\:\gamma\:\:\)on \(\:{\mathcal{R}}_{T}\).Full size imageThe simulation curve described by Fig. 6 above mainly highlights the powerful, inverse relationship between the chronic tonsillitis infection treatment rate \(\:\gamma\:\) and the spread of tonsillitis, contrasting the calculated \(\:{\mathcal{R}}_{T}\) with the epidemic threshold. The blue curve, representing the calculated \(\:{\mathcal{R}}_{T}\), starts high when the treatment rate is low and decreases sharply as treatment becomes more rapid, demonstrating that effective management of chronic cases is a potent control measure. This is juxtaposed against the static red line at \(\:{\mathcal{R}}_{T}=\:1\), the tipping point for an epidemic. To the right of the critical γ value of approximately 0.03, the blue curve lies below the red line (\(\:{\mathcal{R}}_{T}\:1\)), allowing chronic cases to persist and sustain transmission, which would lead to an endemic state. Therefore, the intersection point precisely defines the minimum chronic treatment rate that public health programs must achieve to successfully control the disease.Optimal control problem analysisIn this section, since the sensitivity analysis detailed results discussed in Sect. 4 above verified that some of the model parameters such as the transmission rate and the treatment rate are sensitive to tonsillitis infection spreading dynamics we consider to formulate and analyze an optimal control problem. The control mechanisms to manage the tonsillitis infection can be divided into (i) preventive measures (while we cannot completely avert tonsillitis, we can minimize the risk by adhering to good hygiene practices, such as frequently washing our hands particularly before touching our nose or mouth and avoiding the sharing of food, drinks, or utensils with individuals who are ill) and (ii) treatment for affected individuals (pertaining to both acute and chronic tonsillitis), indicating that the management of tonsillitis varies based on its underlying cause. Although the symptoms of viral and bacterial tonsillitis can often overlap, their treatment strategies differ. In the following sections, we shall detail the specific control measures that can be applied for each of these two types of interventions. Our quantitative and qualitative research has shown that, as previously mentioned, implementing a range of control measures is vital for reducing the effects of tonsillitis infections. To determine the most effective control strategy for decreasing the incidence of tonsillitis across various socioeconomic groups, we intend to carry out an optimal control analysis. To reduce the costs associated with implementing the controls, we assume that the control parameter \(\:{w}_{1}\:\) represents preventive measures for managing tonsillitis infections. Likewise, to keep implementation costs low, we also assume that the control parameters \(\:{w}_{2}\:\)and \(\:{w}_{3}\:\)signify treatment measures for managing acute tonsillitis and chronic tonsillitis, respectively, by enhancing treatments for both conditions over time. We then construct a suitable optimal control function aimed at minimizing the cost of these interventions while adhering to the model presented in (3). In this section, we revisit and assess the optimal control model (7) to establish the most effective control strategy that decreases the number of individuals infected with tonsillitis, all while incurring the least total intervention costs. The goal is to determine the optimal values \(\:{W}^{*}=\left({w}_{1}^{*},\:{w}_{2}^{*},{w}_{3}^{*}\right)\) for the controls \(\:W=\left({w}_{1},\:{w}_{2},{w}_{3}\right)\), such that the corresponding state trajectories \(\:{S}^{*}\left(t\right),\:\:{A}^{*}\left(t\right),\:\:{C}^{*}\left(t\right),\:\:{T}^{*}\left(t\right),\:\:{R}^{*}\left(t\right))\) serve as solutions to the system (3) over the intervention time interval\(\:\left[0,\:{T}_{f}\right]\) with the initial conditions specified in (3) and minimize the objective functional. This aligns with the execution of effective treatment and prevention strategies for those infected with tonsillitis within a community, aiming to enhance recovery periods through appropriate treatment and counseling, ensuring that \(\:0\:\le\:{w}_{1},\:{w}_{2},{w}_{3}\le\:1\). The model system (3) with optimal control is now presented as follows.$$\:\frac{dS}{dt}={\Lambda\:}+\omega\:R-\left((1-{w}_{1}){\lambda\:}_{T}+d\right)S,$$$$\:\frac{dA}{dt}=(1-q)(1-{w}_{1}){\lambda\:}_{T}S-\left(\delta\:\right(1+{w}_{2})+d)A,$$$$\:\frac{dC}{dt}=q(1-{w}_{1}){\lambda\:}_{T}S+\delta\:(1+{w}_{2})A-\left(\gamma\:{(1+w}_{3}\right)+d+{d}_{T})C,$$$$\:\frac{dT}{dt}=\gamma\:(1+{w}_{3})C-(\sigma\:+\text{d})T,$$(8)$$\:\frac{dR}{dt}={\upsigma\:}\text{T}-(\omega\:+d)R,$$with initial data given by \(\:S\left(0\right)\ge\:0\),\(\:\:A\left(0\right)\ge\:0\),\(\:\:C\left(0\right)\ge\:0\), \(\:T\left(0\right)\ge\:0\), and \(\:R\left(0\right)\ge\:0\).For this, minimizing the objective functional is our optimum control issue which is given by$$\:J({w}_{1},\:{w}_{2},{w}_{3})\:={\int\:}_{0}^{{T}_{f}}\left({{\Theta\:}}_{1}A+{{\Theta\:}}_{2}C+\frac{1}{2}\sum\:_{i=1}^{3}{\text{Y}}_{i}{w}_{i}^{2}\right)dt,$$(9)where\(\:\:I\left({E}_{p}^{*},\:{w}_{1},\:{w}_{2},{w}_{3}\right)\:={{\Theta\:}}_{1}A+{{\Theta\:}}_{2}C+\frac{{\text{Y}}_{1}}{2}{w}_{1}^{2}+\frac{{\text{Y}}_{2}}{2}{w}_{2}^{2}+\frac{{\text{Y}}_{3}}{2}{w}_{3}^{2}\), the current cost at time t, is assessed, with \(\:{T}_{f}\) representing the final time. The coefficients \(\:{{\Theta\:}}_{1},\:\:\)and \(\:{{\Theta\:}}_{2}\:\)are positive weight constants, while \(\:\frac{{\text{Y}}_{1}}{2},\:\frac{{\text{Y}}_{2}}{2}\)and \(\:\:\frac{{\text{Y}}_{3}}{2}\:\) measure the relative costs associated with the controls \(\:{w}_{1},\:{w}_{2},\)and\(\:{\:w}_{3}\), respectively, and ensure the integrand’s units are balanced. In the cost functional, the term\(\:{{\Theta\:}}_{1}A\) pertains to the expenses associated with the acute tonsillitis class, while \(\:{{\Theta\:}}_{2}C\) pertains to the expenses linked to the chronic tonsillitis class in a community context. Additionally, the cost functional J represents the total expenses stemming from tonsillitis infections and the strategies for their control. The permissible control functions are defined by$$\:{{\Omega\:}}_{w}=\left\{{w}_{1}\left(t\right),\:{w}_{2}\left(t\right),{w}_{3}\left(t\right)\in\:{L}^{3}:0\le\:{w}_{1},\:{w}_{2},{w}_{3}\le\:1,\:t\in\:\left[0,{T}_{f}\:\right]\right\}.$$(10)More specifically, we look for the best control pair represented by:$$\:J\left({w}_{1}^{*},\:{w}_{2}^{*},{w}_{3}^{*}\right)=\underset{{{\Omega\:}}_{w}}{\text{min}}J\left({w}_{1},\:{w}_{2},{w}_{3}\right).$$(11)According to Pontryagin’s Maximum Principle28,46,47,48, if \(\:{W}^{*}\left(t\right)\in\:{{\Omega\:}}_{w}\:\) is optimal for dynamical system (9) with the initial value (9) and (10) with fixed final time\(\:\:{T}_{f}\), then there exists a non-trivial absolutely continuous mapping\(\:\:{\Phi\:}:\left[0,\:{T}_{f}\right]\to\:{\mathbb{R}}^{6}\),\(\:\:{\Phi\:}\:=\left({{\Phi\:}}_{1}\left(t\right),\:\:{{\Phi\:}}_{2}\left(t\right),\:\:{{\Phi\:}}_{3}\left(t\right),\:\:{{\Phi\:}}_{4}\left(t\right),\:\:{{\Phi\:}}_{5}\left(t\right))\right)\) called the adjoint vector, such that.(1)The Hamiltonian function is defined as$$\:H\:={{\Theta\:}}_{1}A+{{\Theta\:}}_{2}C+\frac{1}{2}\sum\:_{i=1}^{3}{\text{Y}}_{i}{w}_{i}^{2}+\sum\:_{i=1}^{5}{{\Phi\:}}_{i}{\text{{\rm\:Y}}}_{i}.$$(12)where \(\:{\text{{\rm\:Y}}}_{i}\:\) stands for the right-hands side of model (8) which is the \(\:{i}{\rm th}\) state variable equation.(2)The control system$$\:\frac{dS}{dt}=\frac{\partial\:H}{\partial\:{{\Phi\:}}_{1}},\:\:\frac{dA}{dt}=\frac{\partial\:H}{\partial\:{{\Phi\:}}_{2}},\:\:\:\frac{dC}{dt}=\frac{\partial\:H}{\partial\:{{\Phi\:}}_{3}},\:\:\frac{dT}{dt}=\frac{\partial\:H}{\partial\:{{\Phi\:}}_{4}},\:\:\frac{dR}{dt}=\frac{\partial\:H}{\partial\:{{\Phi\:}}_{5}}\: .$$(13)(3)The adjoint system$$\:\frac{d{{\Phi\:}}_{1}}{dt}=-\frac{\partial\:H}{\partial\:S},\:\:\frac{d{{\Phi\:}}_{2}}{dt}=-\frac{\partial\:H}{\partial\:A},\:\:\:\frac{d{{\Phi\:}}_{3}}{dt}=-\frac{\partial\:H}{\partial\:C},\:\:\frac{d{{\Phi\:}}_{4}}{dt}=-\frac{\partial\:H}{\partial\:T},\:\:\frac{d{{\Phi\:}}_{5}}{dt}=-\frac{\partial\:H}{\partial\:R}\:.$$(14)(4)The optimality condition$$\:H\left({E}_{p}^{*},w,{{\Phi\:}}^{*}\right)=\underset{w\in\:{{\Omega\:}}_{w}}{\text{min}}\left({E}_{p}^{*},{w}^{*},{{\Phi\:}}^{*}\right)$$(15)holds for almost all \(\:t\in\:[0,\:{T}_{f}].\).(5)Moreover, the transversality condition$$\:{{\Phi\:}}_{i}\left({T}_{f}\right)=0,\:i=1,\:2,\:3,\dots\:,5,$$(16)also holds. We discuss the characterization of optimal controls and adjoint variables in the next result.Theorem 6Let \(\:{W}^{*}=\:\left({w}_{1}^{*},\:{w}_{2}^{*},{w}_{3}^{*}\right)\) be the optimal control value for the proposed control strategies \(\:W=\:\left({w}_{1},\:{w}_{2},{w}_{3}\right)\) and \(\:\left({S}^{*}\left(t\right),{A}^{*}\left(t\right),{C}^{*}\left(t\right),{T}^{*}\left(t\right),{R}^{*}\left(t\right))\right)\) be the associated unique optimal solutions of the optimal control problem (8) with initial condition and objective functional (9) with fixed final time \(\:{T}_{f}\) (10). Then there exists adjoint function \(\:{{\Phi\:}}_{i}^{*}\left(t\right),\:i\:=\:1,\:...\:,5\) satisfying the following canonical equations:$$\:\frac{d{{\Phi\:}}_{1}}{dt}={{\Phi\:}}_{1}\left(\left(1-{w}_{1}\right){\lambda\:}_{T}^{*}+d\right)-{{\Phi\:}}_{2}\left(1-q\right)\left(1-{w}_{1}\right){\lambda\:}_{T}^{*}-{{\Phi\:}}_{3}q\left(1-{w}_{1}\right){\lambda\:}_{T}^{*},$$$$\:\:\frac{d{{\Phi\:}}_{2}}{dt}=-{{\Theta\:}}_{1}-{{\Phi\:}}_{2}\left(\delta\:{(1+w}_{2}\right)+d)+{{\Phi\:}}_{3}\delta\:\left(1+{w}_{2}\right),$$$$\:\:\frac{d{{\Phi\:}}_{3}}{dt}=-{{\Theta\:}}_{2}-{{\Phi\:}}_{3}\left(\gamma\:\left(1+{w}_{3}\right)+d+{d}_{T}\right)+{{\Phi\:}}_{4}\gamma\:\left(1+{w}_{3}\right),$$$$\:\frac{d{{\Phi\:}}_{4}}{dt}={{\Phi\:}}_{4}\left(\sigma\:+d\right)-{{\Phi\:}}_{5}\sigma\:,$$$$\:\:\frac{d{{\Psi\:}}_{5}}{dt}={{\Phi\:}}_{1}\omega\:-{{\Phi\:}}_{5}(\omega\:+d),$$$$\:{{\Phi\:}}_{i}^{*}\left({T}_{f}\right)=0,\:i=1,\:2,\:\dots\:,5.\:$$(17)The control functions \(\:{w}_{1}^{*}\left(t\right),\:\:{w}_{2}\left(t\right)\:\)and\(\:\:{w}_{3}^{*}\left(t\right)\) are obtained by setting \(\:\frac{\partial\:H}{\partial\:{w}_{i}}=0,\:i=\text{1,2},3\) i.e.$$\:\frac{\partial\:H}{\partial\:{w}_{1}}={Y}_{1}{w}_{1}+{\lambda\:}_{T}S\left({{\Phi\:}}_{1}-{{\Phi\:}}_{2}\right)+q{\lambda\:}_{T}S\left({{\Phi\:}}_{2}-{{\Phi\:}}_{3}\right)=0$$$$\:{w}_{1}^{\circ\:}=\frac{\left[{{\Phi\:}}_{2}\right(1-q)+{{\Phi\:}}_{3}q-{{\Phi\:}}_{1}]{\lambda\:}_{T}S}{{Y}_{1}}$$$$\:\frac{\partial\:H}{\partial\:{w}_{2}}={Y}_{2}{w}_{2}-\delta\:A\left({{\Phi\:}}_{3}-{{\Phi\:}}_{2}\right)=0$$$$\:{w}_{2}^{\circ\:}=\frac{\delta\:A\left({{\Phi\:}}_{3}-{{\Phi\:}}_{2}\right)}{{Y}_{2}}$$$$\:\frac{\partial\:H}{\partial\:{w}_{3}}=\:{Y}_{3}{w}_{3}-\gamma\:C\left({{\Phi\:}}_{4}-{{\Phi\:}}_{3}\right)=0$$\(\:{w}_{3}^{\circ\:}=\frac{\gamma\:C\left({{{\Phi\:}}_{4}-{\Phi\:}}_{3}\right)}{{Y}_{3}}\), then the control solution denoted by$$\:{w}_{i}^{*}=\left\{\begin{array}{c}0,\:\:if\:{w}_{i}