Dynamic analysis and optimization of syphilis spread: Simulations, integrating treatment and public health interventions

Muhammad Farman; Khadija Jamil; Kottakkaran Sooppy Nisar; Yasir Nadeem Anjam; Muhammad Umer Saleem; Evren Hincal

doi:10.1515/nleng-2025-0100

Article Open Access

Dynamic analysis and optimization of syphilis spread: Simulations, integrating treatment and public health interventions

Muhammad Farman , Khadija Jamil , Kottakkaran Sooppy Nisar , Yasir Nadeem Anjam , Muhammad Umer Saleem and Evren Hincal

Published/Copyright: July 1, 2025

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information

From the journal Nonlinear Engineering Volume 14 Issue 1

Abstract

Throughout global history, the population has faced unprecedented challenges due to infectious spreads. Addressing the need to mitigate these infections requires well-directed and comprehensive efforts. Sexually transmitted infections, including syphilis, remain significant global health concerns affecting both developed and developing nations. Syphilis, resulting from the transmission of the Treponema pallidum bacterium through sexual contact, is estimated to affect around 12 million individuals annually worldwide. The objectives are achieved by launching the syphilis model, dividing the entire population into six compartments. Additionally, an ideal control plan is presented that integrates the most practicable medical measures to lower the quantity of afflicted persons encapsulate the dynamics of the prevailing degree of syphilis in a population. The devised model is validated by verifying essential features such as positivity, invariant region, and equilibrium points of the point for feasibility of solutions. This investigation focuses on examining the local stability of the syphilis model with a specific emphasis on considering limited observations a critical aspect of epidemic models. The reproductive number R 0 has been calculated to evaluate its impact across various sub-compartments playing a pivotal role in determining community-wide transmission rates. A sensitivity analysis of the models parameters has been performed. To gain numerical solutions, the advanced and well-established numerical technique nonstandard finite difference is employed to provide insights into the genuine behavior of the model. Additionally, to assist in achieving the fundamental aim of this research, an optimal control strategy is induced by considering control variables, namely, an educational awareness campaign and treatment protocols intended to reduce the prevalence of infected individuals. The purposes are attained by employing the Pontryagin maximum principle through mathematically and modeling.

Keywords: syphilis disease model; stability analysis; sensitivity analysis; nonstandard finite difference; optimal control strategy

1 Introduction

Syphilis, an infectious disease transmitted through sexual contact and caused by the bacterium Treponema pallidum, remains a significant public health challenge worldwide [1,2]. This ailment accounts for an estimated 12–13 million cases annually, leading to almost 200,000 deaths [3]. The prevalence of syphilis varies across different countries, with figures ranging from 3.5 to 30 per 100,000 individuals in nations such as Germany, the United Kingdom, the Netherlands, and Denmark [4]. Unfortunately, there is a lack of data on syphilis in developing nations, although Downing’s research [5] provides insights, particularly regarding Nigeria. Research by the WHO reveals that approximately 2 million expectant mothers contract syphilis infection annually, with transmission from mother to fetus occurring in 80% of these cases [2,6,7].

Syphilis, often referred to as the “great imitator” by Sir William Osler [8], is transmitted through direct contact with ulcers caused by the bacterium. If left untreated, the infection progresses through various phases: the initial phase, typically lasting 3–6 weeks; the subsequent phase, lasting 4–10 weeks; the dormant phase, with an average duration of 10 weeks to 2 years; and finally, the advanced (tertiary) phase, extending beyond 2 years [4,9]. During the subsequent stage, several symptoms emerge, including fever, skin rash, lymphadenopathy, and genital or perineal condyloma latum. In the dormant phase, the risk of damage to internal organs becomes a concern, potentially leading to significant health complications and even death. Clinical expressions decrease during this time, and the infection can only be detected through serologic examination. The final (tertiary) phase may manifest as gummatous malady, cardiovascular complications, or involvement of the central nervous system. The treatment for syphilis involves the use of antibiotics such as penicillin, which can be effective at any time while the infection is still active. In the initial phases, a single injection of penicillin may be sufficient for a cure. On the contrary, those whose syphilis does not go away after a year might need to take further medication. It is essential to remember that even those who have recovered from syphilis infection are susceptible to reinfection because they do not have a lifetime immunity [10,11]. Hence, the risk of reinfection and relapse remains present at any stage, even following the successful treatment of a prior infection [12,13].

Several deterministic mathematical models pertaining to syphilis have been created and examined. Garnett et al. introduced a model that considers different stages of syphilis but does not distinguish between early and late latency, incorporating therapy where treated individuals return to susceptibility or an immune class based on the stage of treatment [14]. Pourbohloul et al. formulated a model with 210 differential equations based on the system of syphilis dynamics. This model integrates various stages of syphilis and stratifies the population based on gender, sexual behavior, and age. Fenton et al. conducted a comprehensive review of existing mathematical models related to syphilis, encompassing studies until 2008. Their analysis provided a detailed summary of the collected information and its implications [15]. Epidemiological models play a crucial role in devising effective control strategies by providing a comprehensive understanding of the complex dynamics of infectious diseases [16,17]. These models help identify key factors essential for understanding the transmission and management of infections [18,19]. Mathematical models, with their fundamental assumptions and principles, aid in characterizing variables for different infectious diseases and assessing the effects of potential treatments [20]. In works referenced as [21,22], study the impact of different ways of treatments. These models included various categories within the population, accounting for factors such as age, gender, and sexual activity [23]. An ordinary differential equations (ODEs) model considering partial immunization and vaccination, assuming the availability of an effective vaccine. They demonstrated the existence of backward bifurcation for specific input variables [24].

Verma et al. [25] have employed optimal control in mathematical models to study the dynamics of COVID-19, considering lockdown effects and providing epidemiological insights. Verma [26] conducted another investigation that employed a nonlinear smoking model to analyze the interactions between smokers and quiters. First integrals obtained and symmetries in a generalized coupled Lane–Emden system were investigated in the study of Muatjetjeja and Khalique [27], improving the understanding of the system dynamics and applications. Symmetry groups for the system are identified [28], which facilitates the use of symmetry-based methods to simplify and solve nonlinear differential equations. Find invariant solutions and conservation rules by applying symmetry analysis [29]. Classify the system and establish conservation rules in the study of Muatjetjeja [30] by applying symmetry approaches, which offer crucial insights into the behavior and solutions of this nonlinear differential system. The work of Omame and Abbas [31] involved with optimal control, additionally, an optimal control strategy introduced by Omame et al. [32] addressed interactions between a different epidemic model. Various researchers have applied optimal control methods to model and forecast the effectiveness of health monitoring strategies, including studies on several diseases [33], and optimal studied for public health is studied in the study of Omame et al. [34]. Sensitivity analysis and optimality control of infectious disease transmission have also been explored by Nadeem Anjam et al. [35]. Stimulated by the implications of the above-documented issues, this study explores the prevalence of syphilis disease in a community. The article presents a proposed syphilis model classified into six compartments. First, we will evaluate the validity of the proposed model by establishing fundamental properties, namely, positivity, equilibrium points. Next, we develop the syphilis model using control variables provided through a qualitative optimal control concerning two control variables: an educational awareness campaign and treatment protocols intended to reduce the prevalence of infected individuals.

The work is structured as follows: Section 1 of this article serves as an introduction. The new developed model with well-posedness, equilibrium points, reproductive potential and sensitivity analysis in Section 2. In Section 3, we verify the existence, uniqueness and stability analysis of the model. In Section 4 we examine numerical scheme using nonstandard finite difference (NSFD) methods with convergence analysis. In Section 5, we develop the optimal control strategies with control parameter and their control results. In Section 6, we discuss the conclusion and future direction of the proposed model for public health.

2 Development of a syphilis disease model

Epidemiological models of infectious diseases, outlined in the literature [36], are vital for understanding disease spread patterns and formulating the different strategies for control. We created the new model, and in this section, we provide a novel mathematical analysis of syphilis, categorizing the entire studied population into six distinct groups. The susceptible group ( S ) comprises individuals at risk of syphilis exposure, while the infected class ( I ) represents those currently affected. The exposed group ( E ) consists of individuals exhibiting symptoms without definitive proof of infection, and the undergoing treatment group ( A ) includes those actively receiving treatment for infections. Recovered individuals ( R ) post-treatment have developed immunity and are free from infection. Additionally, the cautionary/preventative group ( Q ) involves individuals practicing safe sex or taking preventive measures to reduce the risk of infection. This model offers a comprehensive framework for understanding the dynamics of syphilis transmission within a population

(1) d S d t = Δ − β S E − ( ψ ) S + ρ R , d E d t = β S E − ( ω + ψ + λ ) E , d I d t = ω E − ( γ + ψ + d ) I , d A d t = γ I + λ E − ( μ + ψ ) A , d R d t = μ A + d I − ( ψ + ρ + α ) R , d Q d t = α R − ψ Q ,

with the initial conditions

(2) S ( 0 ) ≥ 0 , E ( 0 ) ≥ 0 , I ( 0 ) ≥ 0 , A ( 0 ) ≥ 0 , R ( 0 ) ≥ 0 , Q ( 0 ) ≥ 0 .

The model incorporates several parameters that influence the population behavior over time. The recruitment rate into the susceptible population is denoted by Δ . The transmission rate of syphilis from the exposed class E to the susceptible class S is denoted by β . The natural death rate ψ represents the rate at which individuals from any group die due to causes unrelated to syphilis. The rate at which recovered individuals R lose immunity and return to the susceptible group is represented by ρ . The rate at which exposed individuals E progress to the infected class I is denoted as ω . The rate at which exposed individuals E start treatment and move to the treatment class A is denoted by λ . The rate at which infected individuals I move into the treatment class A is represented by γ . The rate at which infected individuals I recover and move to the recovered class R is represented by d . The rate at which individuals undergoing treatment A recover and move to the recovered class R is denoted by μ . The rate at which recovered individuals R take preventive measures and move to the cautionary group Q is represented by α . The flow diagram of proposed model is shown in Figure 1.

Figure 1

Flow diagram corresponding to the syphilis disease model.

2.1 Positively invariant region

Lemma 1

The region Ξ ∈ R + 6

Ξ = { ( S , E , I , A , R , Q ) ∈ R + 6 : 0 ≤ N }

the solution to system (1) is attracted by every system, and for the specified system in R + 6 , it remains positively invariant when subjected to non-negative initial conditions.

Proof

We will discuss how the systems (1) have positive solutions, and we will give the results as follows:

(3) d S d t S = 0 = Δ ≥ 0 , d E d t E = 0 = 0 ≥ 0 , d I d t I = 0 = ω E ≥ 0 , d A d t A = 0 = γ I ≥ 0 , d R d t R = 0 = μ A + d I ≥ 0 , d Q d t Q = 0 = α R ≥ 0 .

The vector field is purportedly situated within the domain R + 6 on every hyperplane that encompasses the non-negative orthant with the condition that t ≥ 0 in accordance with system (3).□

2.2 Equilibrium points and reproductive number

The proposed model (1) equilibrium states are as follows:

E 0 = ( S + , E + , I + , A + , R + , Q + ) = Δ ψ , 0, 0, 0, 0, 0 , E * = ( S − , E − , I − , A − , R − , Q − ) , S − = λ + ω + ψ β , E − = ( α + ψ + ρ ) ( λ ψ − Δ β + ω ψ + ψ 2 ) ( d μ + d ψ + γ μ + γ ψ + μ ψ + ψ 2 ) M , I − = ( μ ω + ω ψ ) ( α + ψ + ρ ) ( λ ψ − Δ β + ω ψ + ψ 2 ) M , A − = ( α + ψ + ρ ) ( d λ + γ λ + γ ω + λ ψ ) ( λ ψ − Δ β + ω ψ + ψ 2 ) M , R − = M 1 M , Q − = M 2 M 3 ,

where

M = β ψ 4 + β ω ψ 3 + β ψ 3 ρ + α β ψ 3 + β d ψ 3 + β γ ψ 3 + β λ ψ 3 + β μ ψ 3 + α β d ψ 2 + α β γ ψ 2 + α β λ ψ 2 + α β μ ψ 2 + β d λ ψ 2 + α β ω ψ 2 + β d μ ψ 2 + β d ω ψ 2 + β γ λ ψ 2 + β γ μ ψ 2 + β d ψ 2 ρ + β γ ω ψ 2 + β γ ψ 2 ρ + β λ μ ψ 2 + β μ ω ψ 2 + β λ ψ 2 ρ + β μ ψ 2 ρ + β ω ψ 2 ρ + α β d λ μ + α β d λ ψ + α β d μ ω + α β γ λ μ + α β d μ ψ + α β d ω ψ + α β γ λ ψ + α β γ μ ω + α β γ μ ψ + α β γ ω ψ + α β λ μ ψ + β d λ μ ψ + α β μ ω ψ + β d μ ω ψ + β γ λ μ ψ + β d λ ψ ρ + β d μ ψ ρ + β γ μ ω ψ + β γ λ ψ ρ + β γ μ ψ ρ + β γ ω ψ ρ + β μ ω ψ ρ ,

M 1 = d ω ψ 3 + λ μ ψ 3 + d ω 2 ψ 2 + λ 2 μ ψ 2 + d λ μ ψ 2 + d λ 2 μ ψ + d λ ω ψ 2 + d μ ω ψ 2 + d μ ω 2 ψ + γ λ μ ψ 2 + γ λ 2 μ ψ + γ μ ω ψ 2 + γ μ ω 2 ψ + λ μ ω ψ 2 − Δ β d λ μ − Δ β d μ ω − Δ β γ λ μ − Δ β d ω ψ − Δ β γ μ ω − Δ β λ μ ψ + 2 d λ μ ω ψ + 2 γ λ μ ω ψ ,

M 2 = α d ω ψ 3 + α λ μ ψ 3 + α d ω 2 ψ 2 + α λ 2 μ ψ 2 + α d λ μ ψ 2 + α d λ 2 μ ψ + α d λ ω ψ 2 + α d μ ω ψ 2 + α d μ ω 2 ψ + α γ λ μ ψ 2 + α γ λ 2 μ ψ + α γ μ ω ψ 2 + α γ μ ω 2 ψ + α λ μ ω ψ 2 − Δ α β d λ μ − Δ α β d μ ω − Δ α β γ λ μ − Δ α β d ω ψ − Δ α β γ μ ω − Δ α β λ μ ψ ,

M 3 = β ψ 5 + β ω ψ 4 + β ψ 4 ρ + α β ψ 4 + β d ψ 4 + β γ ψ 4 + β λ ψ 4 + β μ ψ 4 + α β d ψ 3 + α β γ ψ 3 + α β λ ψ 3 + α β μ ψ 3 + β d λ ψ 3 + α β ω ψ 3 + β d μ ψ 3 + β d ω ψ 3 + β γ λ ψ 3 + β γ μ ψ 3 + β d ψ 3 ρ + β γ ω ψ 3 + β γ ψ 3 ρ + β λ μ ψ 3 + β μ ω ψ 3 + β λ ψ 3 ρ + β μ ψ 3 ρ + β ω ψ 3 ρ + α β d λ ψ 2 + α β d μ ψ 2 + α β d ω ψ 2 + α β γ λ ψ 2 + α β γ μ ψ 2 + α β γ ω ψ 2 + α β λ μ ψ 2 + β d λ μ ψ 2 + α β μ ω ψ 2 + β d μ ω ψ 2 + β γ λ μ ψ 2 + β d λ ψ 2 ρ + β d μ ψ 2 ρ + β γ μ ω ψ 2 + β γ λ ψ 2 ρ + β γ μ ψ 2 ρ + β γ ω ψ 2 ρ + β μ ω ψ 2 ρ + α β d λ μ ψ + α β d μ ω ψ + α β γ λ μ ψ + α β γ μ ω ψ .

Consider the infectious classes:

d E d t = β S E − ( ω + ψ + λ ) E , d I d t = ω E − ( γ + ψ + d ) E , d A d t = γ I + λ E − ( μ + ψ ) A .

Now, it is imperative for us to calculate the matrices F and V utilizing the next-generation matrix method in the following manner:

F ( E 0 ) = β S + 0 0 0 0 0 0 0 0 V − 1 = 1 ( λ + ω + ψ ) 0 0 ω ( d + γ + ψ ) ( λ + ω + ψ ) 1 ( d + γ + ψ ) 0 d λ + γ λ + γ ω + λ ψ ( μ + ψ ) ( d + γ + ψ ) ( λ + ω + ψ ) γ ( μ + ψ ) ( d + γ + ψ ) 1 ( μ + ψ ) .

Subsequently, the number of reproductive

(4) R 0 = Δ β ψ ( λ + ω + ψ ) .

As widely understood, when the R 0 < 1 , the infection will ultimately go away. If R 0 > 1 , then the disease will propagate across the population. Impact of different parameters on reproductive number is shown in Figures 2 and 3.

Figure 2

Parameters impact on the reproductive number.

Figure 3

Parameters impact on the reproductive number.

2.3 Sensitivity analysis

Sensitivity analysis is applied to investigate how parameters impact the model that was suggested. It is crucial to identify parameters that show high sensitivity, potentially causing disruption in the model dynamics with even a minor numerical alteration. To evaluate the sensitivity of R 0 and results are shown through simulation in Figures 2(a)–(c) and 3(a)–(c).

∂ R 0 ∂ β = Δ ψ ( λ + ω + ψ ) > 0 , ∂ R 0 ∂ Δ = β ψ ( λ + ω + ψ ) > 0 , ∂ R 0 ∂ ψ = − Δ β ψ ( λ + ω + ψ ) 2 − Δ β ψ 2 ( λ + ω + ψ ) < 0 , ∂ R 0 ∂ λ = − Δ β ψ ( λ + ω + ψ ) 2 < 0 , ∂ R 0 ∂ ω = − Δ β ψ ( λ + ω + ψ ) 2 < 0 .

Positive index parameters, especially β and Δ , have a positive effect on R 0 . It can be inferred that a rise in β and Δ values could increase R 0 or initiate an epidemic. On the contrary, increasing the values of parameters that have a negative sensitivity index, such as ψ , λ , and ω , has a detrimental impact on stopping the disease from spreading.

3 Existence, uniqueness, and stability analysis of the system

Consider system (1) as given below:

d S d t = Δ − β S I − ψ S + ρ R = G 1 ( ζ 1 ) , d E d t = β S I − ( ω + ψ + λ ) E = G 2 ( ζ 1 ) , d I d t = ω E − ( γ + ψ + d ) I = G 3 ( ζ 1 ) , d A d t = γ I + λ E − ( μ + ψ ) A = G 4 ( ζ 1 ) , d R d t = μ A + d I − ( ψ + ρ + α ) R = G 5 ( ζ 1 ) , d Q d t = α R − ψ Q = G 6 ( ζ 1 ) ,

where we have ζ 1 = ( t , S , E , I , A , R , Q ) , to ensure the existence and uniqueness of our version, we need to verify the following:

∣ G ( t , S , E , I , R , Q ) ∣ 2 ≤ H i ( 1 + ∣ S ∣ 2 ) , i = 1 , 2, 3, 4, 5, 6.
∀ S , S 1 , ∥ G ( t , S , I , E , A , R , Q ) − G ( t , S 1 , I , E , A , R , Q ) ∥ ∞ 2 ≤ H ¯ i ∥ S − S 1 ∥ ∞ 2 .

We have

∣ G 1 ( ζ 1 ) ∣ 2 = ∣ Δ − ( β E + ψ ) S + ρ R ∣ 2 ≤ 2 Δ 2 + 2 ( β E + ψ ) 2 ∣ S ∣ 2 + 2 ρ 2 ∣ R ∣ 2 ≤ 2 ( Δ 2 + ρ 2 ∣ R ∣ 2 ) + 2 ( β E + ψ ) 2 ∣ S ∣ 2 ≤ 2 ( Δ 2 + ρ 2 ∣ R ∣ 2 ) 1 + ( β E + ψ ) 2 ∣ S ∣ 2 Δ 2 + ρ 2 ∣ R ∣ 2 .

If ( β E + ψ ) 2 ∣ S ∣ 2 Δ 2 + ρ 2 ∣ R ∣ 2 < 1 , then we have ∣ G 1 ( ζ 1 ) ∣ 2 ≤ H 1 ( 1 + ∣ S ∣ 2 ) , where H 1 = 2 ( Δ 2 + ρ 2 ∣ R ∣ 2 ) .

∣ G 2 ( ζ 1 ) ∣ 2 = ∣ β S E − ( ω + ψ + λ ) E ∣ 2 ≤ 2 β 2 ∣ S ∣ 2 + 2 ( ω + ψ + λ ) 2 ∣ E ∣ 2 ≤ 2 ( β 2 ∣ S ∣ 2 ) 1 + ( ω + ψ + λ ) 2 β 2 ∣ S ∣ 2 ∣ E ∣ 2 .

If 1 + ( ω + ψ + λ ) 2 β 2 ∣ S ∣ 2 ∣ E ∣ 2 < 1 , then we have ∣ G 2 ( ζ 1 ) ∣ 2 ≤ H 2 ( 1 + ∣ S ∣ 2 ) , where H 2 = 2 β 2 ∣ S ∣ 2

∣ G 3 ( ζ 1 ) ∣ 2 = ∣ ω E − ( γ + ψ + d ) I ∣ 2 ≤ 2 ω 2 ∣ E ∣ 2 − 2 ( γ + ψ + d ) 2 ∣ I ∣ 2 ≤ 2 ( ω 2 ∣ E ∣ 2 ) 1 + ( γ + ψ + d ) 2 ∣ I ∣ 2 ω 2 ∣ E ∣ 2 .

If 1 + ( γ + ψ + d ) 2 ∣ I ∣ 2 ω 2 ∣ E ∣ 2 < 1 , then we have ∣ G 3 ( ζ 1 ) ∣ 2 ≤ H 3 ( 1 + ∣ S ∣ 2 ) , where H 3 = 2 ( ω 2 ∣ E ∣ 2 )

∣ G 4 ( ζ 1 ) ∣ 2 = ∣ γ I + λ E − ( μ + ψ ) A ∣ 2 ≤ 2 γ 2 ∣ I ∣ 2 + 2 λ 2 ∣ E ∣ 2 + 2 ( μ + ψ ) 2 ∣ A ∣ 2 ≤ 2 ( γ 2 ∣ I ∣ 2 + λ 2 ∣ E ∣ 2 ) 1 + ( μ + ψ ) 2 ∣ A ∣ 2 γ 2 ∣ I ∣ 2 + λ 2 ∣ E ∣ 2 .

If ( μ + ψ ) 2 ∣ A ∣ 2 γ 2 ∣ I ∣ 2 + λ 2 ∣ E ∣ 2 < 1 , then we have ∣ G 4 ( ζ 1 ) ∣ 2 ≤ H 4 ( 1 + ∣ S ∣ 2 ) , where H 4 = 2 ( γ 2 ∣ I ∣ 2 + λ 2 ∣ E ∣ 2 )

∣ G 5 ( ζ 1 ) ∣ 2 = ∣ μ A + d I − ( ψ + ρ + α ) R ∣ 2 ≤ 2 μ 2 ∣ A ∣ 2 + 2 d 2 ∣ I ∣ 2 + 2 ( ψ + ρ + α ) 2 ∣ R ∣ 2 ≤ 2 ( μ 2 ∣ A ∣ 2 + d 2 ∣ I ∣ 2 ) 1 + ( ψ + ρ + α ) 2 ∣ R ∣ 2 μ 2 ∣ A ∣ 2 + d 2 ∣ I ∣ 2 .

If ( ψ + ρ + α ) 2 ∣ R ∣ 2 μ 2 ∣ A ∣ 2 + d 2 ∣ I ∣ 2 < 1 , then we have ∣ G 5 ( ζ 1 ) ∣ 2 ≤ H 5 ( 1 + ∣ S ∣ 2 ) , where H 5 = 2 ( μ 2 ∣ A ∣ 2 + d 2 ∣ I ∣ 2 )

∣ G 6 ( ζ 1 ) ∣ 2 = ∣ α R − ψ Q ∣ 2 ≤ 2 α 2 ∣ R ∣ 2 + 2 ψ 2 ∣ Q ∣ 2 ≤ 2 ( α 2 ∣ R ∣ 2 ) 1 + ψ 2 ∣ Q ∣ 2 α 2 ∣ R ∣ 2 .

If ψ 2 ∣ Q ∣ 2 α 2 ∣ R ∣ 2 < 1 , then we have ∣ G 6 ( ζ 1 ) ∣ 2 ≤ H 6 ( 1 + ∣ S ∣ 2 ) , where H 6 = 2 ( α 2 ∣ R ∣ 2 ) . The growth conditions are consequently satisfied. Subsequently, we demonstrate that the function meets the Lipschitz condition. We obtain

∣ G 1 ( t , S , E , I , A , R , Q ) − G 1 ( t , S 1 , E , I , A , R , Q ) ∣ 2 = ∣ Δ − β S E − ψ S + ρ R − ( Δ − β S 1 E − ψ S 1 + ρ R ) ∣ 2 ≤ ( β E + ψ ) 2 ∣ ( S − S 1 ) ∣ 2 ,

(5) sup t ∈ D s ∣ G 1 ( t , S , E , I , A , R , Q ) − G 1 ( t , S 1 , E , I , A , R , Q ) ∣ 2 ≤ ( β sup t ∈ D s ∣ E ∣ + ψ ) 2 sup t ∈ D s ∣ ( S − S 1 ) ∣ 2 ≤ ( β ∣ E ∣ + ψ ) 2 ∥ S − S 1 ∥ ∞ 2 ≤ H ¯ 1 ∥ S − S 1 ∥ ∞ 2 ,

∣ G 2 ( t , S , E , I , A , R , Q ) − G 2 ( t , S , I 1 , E , A , R , Q ) ∣ 2 = ∣ β S E − ( ω + ψ + λ ) E − ( β S E 1 − ( ω + ψ + λ ) E 1 ) ∣ 2 ≤ ( β ∣ S ∣ + ω + ψ + λ ) 2 ∣ E − E 1 ∣ 2 ,

(6) sup t ∈ D E ∣ G 2 ( t , S , E , I , A , R , Q ) − G 2 ( t , S , E 1 , I , A , R , Q ) ∣ 2 ≤ ( β sup t ∈ D S ∣ S ∣ + ω + ψ + λ ) sup t ∈ D E ∣ E − E 1 ∣ 2 ≤ ( β ∣ S ∣ + ω + ψ + λ ) 2 ∥ E − E 1 ∥ ∞ 2 ≤ H ¯ 2 ∥ E − E 1 ∥ ∞ 2 ,

∣ G 3 ( t , S , E , I , A , R , Q ) − G 3 ( t , S , E , I 1 , A , R , Q ) ∣ 2 = ∣ ω E − ( γ + ψ + d ) I − ω E − ( γ + ψ + d ) I 1 ∣ 2 ≤ ( γ + ψ + d ) 2 ∣ I − I 1 ∣ 2 ,

(7) sup t ∈ D I ∣ G 3 ( t , S , E , I , A , R , Q ) − G 3 ( t , S , E , I 1 , A , R , Q ) ∣ 2 ≤ ( γ + ψ + d ) 2 sup t ∈ D I ∣ I − I 1 ∣ 2 ≤ ( γ + ψ + d ) 2 ∥ I − I 1 ∥ ∞ 2 ≤ H ¯ 3 ∥ I − I 1 ∥ ∞ 2 ,

∣ G 4 ( t , S , E , I , A , R , Q ) − G 4 ( t , S , E , I , A 1 , R , Q ) ∣ 2 = ∣ γ I + λ E − ( μ + ψ ) A − ( γ I + λ E − ( μ + ψ ) A 1 ) ∣ 2 ≤ ( μ + ψ ) 2 ∣ A − A 1 ∣ 2 ,

(8) sup t ∈ D A ∣ G 4 ( t , S , E , I , A , R , Q ) − G 4 ( t , S , E , I , A 1 , R , Q ) ∣ 2 ≤ ( μ + ψ ) 2 sup t ∈ D A ∣ A − A 1 ∣ 2 ≤ ( μ + ψ ) 2 ∥ A − A 1 ∥ ∞ 2 ≤ H ¯ 4 ∥ A − A 1 ∥ ∞ 2 ,

∣ G 5 ( t , S , E , I , A , R , Q ) − G 5 ( t , S , E , I , A , R 1 , Q ) ∣ 2 = ∣ μ A + d I − ( ψ + ρ + α ) R − ( μ A + d I − ( ψ + ρ + α ) R 1 ) ∣ 2 ≤ ( ψ + ρ + α ) 2 ∣ R − R 1 ∣ 2 ,

(9) sup t ∈ D R ∣ G 5 ( t , S , E , I , A , R , Q ) − G 5 ( t , S , E , I , A , R 1 , Q ) ∣ 2 ≤ ( ψ + ρ + α ) 2 sup t ∈ D R ∣ R − R 1 ∣ 2 ≤ ( ψ + ρ + α ) 2 ∥ R − R 1 ∥ ∞ 2 ≤ H ¯ 5 ∥ R − R 1 ∥ ∞ 2 ,

∣ G 6 ( t , S , I , E , A , R , Q ) − G 6 ( t , S , I , E , A , R , Q 1 ) ∣ 2 = ∣ α R − ψ Q − ( α R − ψ Q 1 ) ∣ 2 ≤ ψ 2 ∣ Q − Q 1 ∣ 2 ,

(10) sup t ∈ D Q ∣ G 6 ( t , S , E , I , A , R , Q ) − G 6 ( t , S , E , I , A , R , Q 1 ) ∣ 2 ≤ ψ 2 sup t ∈ D Q ∣ Q − Q 1 ∣ 2 ≤ ψ 2 ∥ Q − Q 1 ∥ ∞ 2 ≤ H ¯ 6 ∥ Q − Q 1 ∥ ∞ 2 .

System (1) therefore has a unique solution under the following condition.

(11) max ( β I + ψ ) 2 ∣ S ∣ 2 Δ 2 + ρ 2 ∣ R ∣ 2 , ( ω + ψ + λ ) 2 ∣ E ∣ 2 β 2 ∣ S ∣ 2 , ( λ + ψ + d ) 2 ∣ I ∣ 2 ω 2 ∣ E ∣ 2 , ( μ + ψ ) 2 ∣ A ∣ γ 2 ∣ I ∣ 2 + λ 2 ∣ E ∣ 2 , ( ψ + ρ + α ) 2 ∣ R ∣ 2 μ 2 ∣ A ∣ 2 + d 2 ( 1 ) ⋅ ∣ I ∣ 2 , ψ 2 ∣ Q ∣ 2 α 2 ∣ R ∣ 2 , < 1 .

3.1 Uniform persistence of equilibrium points

The convergence characteristics of epidemic dynamics are demonstrated in this subsection using the equilibrium points in the system V ( t ) function. The nonlinear function f ( x ) = − ln x + ( x − 1 ) is taken into consideration in the following proof procedure. Then, we can determine that f ( x ) > 0 for any x > 0 and x ≠ 1 . Moreover, f ( x ) achieves a global minimum value of f ( 1 ) = 0 at x = 1 .

Theorem 3.1

The equilibrium point E 0 1 exerts an attractive influence on all solutions of the model delineated in (1) in Ω under the condition that R 0 ≤ 1 .

Proof

Characterize the Lyapunov function V ( t ) that is associated with the specified solution of model (1).

V ( t ) = 1 ⟨ k ⟩ ∑ i = 1 n i P ( i ) S i 0 f S i S i 0 + E i 0 f E i E i 0 + Θ 1 .

The derivative of time is

V ˙ ( t ) = 1 ⟨ k ⟩ ∑ i = 1 n i P ( i ) 1 − S i 0 S i S ˙ i + 1 − E i 0 E i E ˙ i + Θ ˙ 1 .

V ˙ ( t ) = 1 ⟨ k ⟩ ∑ i = 1 n i P ( i ) 1 − S i 0 S i ( Δ − β S i E i − ( ψ ) S i + ρ R i ) + 1 − E i 0 E i ( β S 1 E 1 − ( ω + ψ + λ ) E 1 ) + Θ ˙ 1 .

So,

V ˙ ( t ) = 1 ⟨ k ⟩ ∑ i = 1 n i P ( i ) μ S i 0 2 − S i S i 0 − S i 0 S i + 1 ⟨ k ⟩ ∑ i = 1 n i P ( i ) r S i 0 × 3 − S i 0 S i − E i E i 0 − S i E i 0 S i 0 E i + α 1 Θ 1 ( R 0 − 1 ) .

By establishing a connection between the arithmetic mean and the geometric mean, we can derive the relation

2 − S i S i 0 − S i 0 S i ≤ 0 3 − S i 0 S i − E i E i 0 − S i E i 0 S i 0 E i ≤ 0 .

Considering R 0 ≤ 1 , we have R 0 − 1 ≤ 0 . As a result, V ˙ ( t ) ≤ 0 . At any point, we have S i = S i 0 , E i = E i 0 for i = 1 , 2 , … , n . The above result shows that the system’s limit is reached if the equation is not an autonomous model.□

Theorem 3.2

The equilibrium point E 1 ⋆ of the model delineated in (1) exerts an attractive influence on all trajectories within the region Ω ⁄ E 0 1 when the condition R 0 > 1 is satisfied.

Proof

The following is the Lyapunov function that is being considered:

V ( t ) = 1 ⟨ k ⟩ ∑ i = 1 n i P ( i ) S ˜ i f S i S ˜ i + E ˜ i f E i E ˜ i + Θ ˜ 1 g Θ 1 Θ ˜ 1 .

The derivative of time is

V ˙ ( t ) = 1 ⟨ k ⟩ ∑ i = 1 n i P ( i ) 1 − S ˜ i S i S ˙ i + 1 − E ˜ i E i E ˙ i + Θ ˙ 1 1 − Θ ˜ 1 Θ 1 = 1 ⟨ k ⟩ ∑ i = 1 n i P ( i ) 1 − S ˜ i S i ( Δ − β S E − ( ψ ) S + ρ R ) + 1 − E ˜ i E i ( β S E − ( ω + ψ + λ ) E ) + 1 ⟨ k ⟩ ∑ i = 1 n i 2 P ( i ) ( β 1 S i + σ 1 E i ) Θ ˜ 1 .

If V ˙ ( t ) = 0 , then this implies that the arithmetic mean is greater than or equal to the geometric mean. The result is that the equilibrium point E 1 1 exerts an attractive influence on all trajectories within the region Ω ⁄ E 0 1 . The proof is now complete.□

3.2 Local stability analysis

We qualitatively analyze the suggested model system (1) in order to learn more about its dynamical characteristics and improve our comprehension of how control strategies affect the dynamics of infectious disease transmission. First, the stability properties of the infectious model are examined.

Theorem 3.3

The fractional-order syphilis model suggests that the disease-free equilibrium point is locally asymptotically stable if R 0 < 1 , but it becomes unstable once R 0 > 1 .

Proof

For system (1) at E 0 , the Jacobian J can be expressed as

J ( E 0 ) = − ψ − Δ ψ β 0 0 ρ 0 0 − ω − ψ − λ + Δ ψ β 0 0 0 0 0 ω − d − γ − λ 0 0 0 0 λ γ − μ − ψ 0 0 0 0 d γ + μ − α − ψ − ρ 0 0 0 0 0 α − ψ

The following are the eigenvalues Λ that result from solving the aforementioned matrix: Λ 1 = − ψ , Λ 2 = − ω − ψ − λ + Δ ψ β , Λ 3 = − d − γ − λ , Λ 4 = − μ − ψ , Λ 5 = − ψ − ρ , Λ 6 = − ψ . It is obvious that the point E 0 is locally asympotically stable if R 0 < 1 and all the eigenvalues of J ( E 0 ) are negative.□

3.3 Global stability analysis

If there is an endemic situation, the Lyapunov function { S , I , E , A , R , Q } satisfies d L d t < 0 at the endemic equilibrium E * .

Theorem 3.4

The globally asymptotically stable endemic equilibrium points E * of the syphilis disease model occur when R 0 > 1 .

Proof

The expression for the Lyapunov function is

(12) L ( S * , E * , I * , A * , R * , Q * ) = S − S * − S * log S * S + E − E * − E * log E * E + I − I * − I * log I * I + A − A * − A * log A * A + R − R * − R * log R * R + Q − Q * − Q * log Q * Q

Now, we express as

(13) d L d t = S − S * S ( Δ − β S E − ψ S + ρ R ) + E − E * E ( β S E − ( ω + ψ + λ ) E ) + I − I * I ( ω E − ( γ + ψ + d ) I ) + A − A * A ( γ I + λ E − ( μ + ψ ) A ) + R − R * R ( μ A + d I − ( ψ + ρ + α ) R ) + Q − Q * Q ( α R − ψ Q )

Submitting S = S − S * , I = I − I * , E = E − E * , A = A − A * , R = R − R * , and Q = Q − Q * leads to. After several calculations, in order to minimize complexity, we write

d L d t ≤ Ω − Ω 1

Ω = Δ + β E * ( S − S * ) 2 S + ρ R + ρ R * S * S + β S ( E − E * ) 2 E + ω E + ω E * I * I + γ I + γ I * A * A + λ E + λ E * A * A + μ A + μ A * R * R + d I + d I * R * R + α R + α R * Q * Q ,

Ω 1 = Δ S * S + β E ( S − S * ) 2 S + ψ ( S − S * ) 2 S + ρ R * + ρ R S * S + β S * ( E − E * ) 2 E + ( ω + ψ + λ ) ( E − E * ) 2 E + ω E * + ω E I * I + ( γ + ψ + d ) ( I − I * ) 2 I + γ I * + γ I A * A + λ E * + λ E A * A + ( μ + ψ ) ( A − A * ) 2 A + μ A * + μ A R * R + d I * + d I R * R + ( μ + ρ + α ) ( R − R * ) 2 R + α R * + α R Q * Q + ( Q − Q * ) 2 Q .

Its satisfied that if Ω < Ω 1 , this yields d L d t < 0 , so, its possible when, S = S * , I = I * , E = E * , A = A * , R = R * , and Q = Q * .

(14) 0 = Ω − Ω 1 ⇒ d L d t = 0 .

It is conceivable to provide evidence for the existence of the greatest compact invariant set in the suggested model

(15) { ( S * , I * , E * , A * , R * , Q h * ) } .

It is possible to infer that if a system is stable, then its corresponding energy function E * must likewise be stable by applying LaSalle’s invariance principle. However, this deduction lacks a comprehensive elucidation concerning the indication of the first derivative of the Lyapunov function thereby necessitating additional inquiry. As a result, we shall delve into an alternative approach that capitalizes on equilibrium points to discern between the two aforementioned methods.□

4 Numerical algorithm and analysis of the proposed scheme

This section is dedicated to devising a numerical scheme that replicates the behavior of the continuous model as represented by the set of equation (1) in the system. The work by Anguelov and Lubuma in 2001 and by using the defined criteria in [37], we have h = Δ t > 0 be a step size then

(16) S k + 1 − S k ϕ = Δ − β S k + 1 E k − ψ S k + 1 + ρ R k , E k + 1 − E k ϕ = β S k + 1 E k − ( ω + ψ + λ ) E k + 1 , I k + 1 − I k ϕ = ω E k + 1 − ( λ + ψ + d ) I k + 1 , A k + 1 − A k ϕ = γ I k + 1 + λ E k + 1 − ( μ + ψ ) A k + 1 , R k + 1 − R k ϕ = μ A k + 1 + d I k + 1 − ( ψ + ρ + α ) R k + 1 , Q k + 1 − Q k ϕ = α R k + 1 − ψ Q k + 1 .

The suggested NSFD scheme for the provided model is as follows:

(17) ϕ = 1 − e − ( λ + ψ + d ) h λ + ψ + d .

The discrete approach (16) is formulated in accordance with the standards established by Mickens categorizing it as NSFD scheme as outlined by Mickens in 1994 [38]. The formulation of this scheme follows the methodology described by Anguelov and Lubuma in 2001 [37].

Rule 1: Equation (17) replaces the conventional denominator h = Δ of the discrete derivatives with the complex denominator function, which satisfies the asymptotic connection ϕ ( h ) = h + O ( h 2 ) . This substitution is carried out to ensure that the denominator function J in the continuous model accurately reflects the dynamics of the system taking into account the underlying parameters μ , ρ , ψ , υ . Precise techniques involving such intricate denominator functions find application in a broad spectrum of dynamical systems as demonstrated in works by Lubuma and Patidar [39] and Gumel [40].

Rule 2: The nonlinear variables are incorporated into the system of equations’ right side by the nonlocal approximation (1). For instance, E ( t k ) I ( t k ) ≅ E k + 1 I k we have instead of E ( t k ) I ( t k ) ≅ E k I k .

4.1 Analysis of the proposed novel scheme

Theorem 4.1

The system of equation (1) represents a continuous model system that includes the NSFD scheme (16)within the biologically feasible domain k.

Proof

(18) S k + 1 = ϕ Δ + ρ ϕ R k + S k 1 + ϕ ( β E k + ψ ) , E k + 1 = β ϕ S k + 1 E k + E k 1 + ϕ ( ω + ψ + λ ) , I k + 1 = ϕ ω E k + 1 + I k 1 + ϕ ( λ + ψ + d ) , A k + 1 = ϕ γ I k + 1 + ϕ λ E k + 1 + A k 1 + h ( μ + ψ ) , R k + 1 = μ ϕ A k + 1 + d ϕ I k + 1 + R k 1 + ϕ ( ψ + ρ + α ) , Q k + 1 = α ϕ R k + 1 + Q k 1 + ϕ ψ .

Thus, S k + 1 ≥ 0 , I k + 1 ≥ 0 , E k + 1 ≥ 0 , A k + 1 ≥ 0 , R k + 1 ≥ 0 , and Q k + 1 ≥ 0 whenever, the discrete variables are non-negative from the preceding iteration, k positive invariance must be proved. When the first three equations of (15) are summed, we obtain

(19) [ 1 + ϕ ψ ] H k + 1 = Δ ϕ + H k − [ 1 + ϕ ( λ + ψ + d ) ] I k ≤ Δ ϕ + H k , [ 1 + ϕ ψ ] H k + 1 ≤ Δ h + H k

H k + 1 ≤ Δ ψ if H k ≤ Δ ψ .

The priori connections for Q k + 1 and R k + 1 stem from the fact that Q k + 1 and E k + 1 are less than or equal to H k + 1 and follow a radial pattern. This concludes the proof.□

4.2 Results of the proposed scheme

A nonlinear incidence model for the syphilis disease has been mathematically analyzed and identified. The most sophisticated methods have been used to calculate theoretical results and evaluate their effectiveness. The use of noninteger parametric parameters in the syphilis disease model has led to some fascinating findings. The system under development leverages the following parameter values: β = 0.8 , Δ = 1,00,10,000 , ω = 20 , α = 0.01 , ψ = 0.00709 , μ = 0.9626 , d = 0.0001393 , γ = 0.922 , λ = 0.13 and ρ = 0.000016590 . The susceptible group ( S ) comprises individuals at risk of syphilis exposure while the infected class ( I ) represents those currently affected. The exposed group ( E ) consists of individuals exhibiting symptoms without definitive proof of infection and the undergoing treatment group ( A ) includes those actively receiving treatment for infections. Recovered individuals ( R ) post-treatment have developed immunity, free from infection and the cautionary/preventative group ( Q ). The presented figure illustrates the convergence solutions pertaining to both disease-free and endemic equilibria through the application of the NSFD method, utilizing a step size of h = 0.1 . Also, the convergence solutions for disease-free and endemic equilibria have been effectively represented using the NSFD scheme at ϕ = ϕ ( h ) + O ( h 2 ) . Figures 4, 5, and 6 demonstrate the effectiveness of the strategy used in reducing syphilis disease. It effectively lowers the infection rate while increasing the populations of both susceptible and recovered individuals in both disease-free and endemic states.

Figure 4

Numerical simulation of model with the proposed scheme for DFP.

Figure 5

Numerical simulation of model with the proposed scheme for EEP.

Figure 6

Numerical simulation of model for EEP with changing population.

5 Optimal control strategy for the proposed model

The objective of this section is to examine the relevance of the syphilis model in relation to the optimal control strategy. Preventive control variables, such as education campaigns through social media and health centers targeting population protection ( z 1 ), and a treatment control measure used by patients to diminish the presence of the viral disease in the body ( z 2 ), are considered potential agents in reducing disease prevalence in the community. As the action variable that can be changed to obtain the intended outcome, the control variable establishes the best values for the choice variables. Finding the most effective course of action involves choosing the control variables optimal value. We employ the principles of optimal control theory [41–43] to formulate a control methodology. The adaptation of the control scheme interning people who might come into contact with syphilis class S ( t ) , syphilis-infected class I ( t ) , exposed class exhibiting symptoms but no proof of infection E ( t ) , and people receiving treatment for infection class A ( t ) , recovered class R ( t ) , and people that use caution when having sex Q ( t ) is discussed in the following lines. Now, the model-given equation (1) has a differential formation with the control variable included, which is as follows:

(20) d S d t = Δ − ( 1 − z 2 ) β S ( t ) E ( t ) − ψ S ( t ) + ρ R ( t ) , d E d t = ( 1 − z 2 ) β S ( t ) E ( t ) − ( ω + λ + ψ + z 1 ) E ( t ) , d I d t = ω E ( t ) − ( γ + d + ψ + z 2 ) I ( t ) , d A d t = γ I ( t ) + λ E ( t ) − ( μ + ψ + z 1 ) A ( t ) , d R d t = μ A ( t ) + d I ( t ) − ( ψ + ρ + α ) R ( t ) + z 2 I ( t ) + z 1 A ( t ) , d Q d t = α R ( t ) − ψ Q ( t ) ,

with given initial conditions. Consequently, we examine an optimal control problem wherein the objective function [5] is delineated and articulated in the following manner:

(21) J ( z 1 ( t ) , z 2 ( t ) ) = ∫ 0 t B 1 * E ( t ) + B 2 * I ( t ) + B 3 * A ( t ) + 1 2 C 1 * z 1 2 ( t ) + 1 2 C 2 * z 2 2 ( t ) d t ,

where B 1 * , B 2 * , B 3 * represent weight constants associated with each class, and C i * ≥ 0 , i = 1 , 2 , represent the cost of control. When it comes to control pairs, the objective function is minimized as in

(22) J { z 1 * , z 2 * } = min { J ( z 1 ( t ) , z 2 ( t ) ) , z 1 ( t ) , z 2 ( t ) ∈ Z } .

The control system (20) influences the optimization process, and the control set Z is defined accordingly,

(23) Z = { ( z 1 , z 2 ) \ z i ( t ) is lebesgue measureable on [ 0 , 1 ] , 0 ≤ z i ( t ) ≤ 1 , i = 1 , 2 } .

Furthermore, Pontryagin’s maximum principle is employed to ascertain optimal control strategies in dynamical systems. It furnishes crucial conditions for achieving an optimal solution by directing the selection of controls to minimize or maximize an objective function J . This direction is obtained through an examination at how the Hamiltonian function H behaves at every position in the state variable.

Theorem 5.1

For the optimal control problem (20), there exists z * ( t ) = ( z 1 * ( t ) , z 2 * ( t ) ∈ Z ) such that

min ( z 1 ( t ) , z 2 ( t ) ∈ Z ) J ( z 1 ( t ) , z 2 ( t ) ) = J ( z 1 ( t ) , z 2 ( t ) ) .

Proof

When analyzing the efficacy control scheme [44,45], various methodologies are explored, resulting in nonnegative values observed for both control and state variables. The control system necessary stability is defined by the convexity and closeness of the control variables. Furthermore, it is important to recognize the convex shape of the integrand, which contains an objective function represented as B 1 * E ( t ) + B 2 * I ( t ) + B 3 * A ( t ) + 1 2 C 1 * z 1 2 ( t ) + 1 2 C 2 * z 2 2 ( t ) , serving as a testament to the proof, where B i * , i = 1 , 2, 3, represent weight constants associated with each class, and C i * ≥ 0 , i = 1 , 2, represent the cost of control. For the purpose of resolving our proposed problem, we once again employ the Pontryagin maximum principle [46–48]. The following is an explicit definition of the Hamiltonian function:

(24) H = L ( S ( t ) , E ( t ) , I ( t ) , A ( t ) , R ( t ) , Q ( t ) , z 1 ( t ) , z 2 ( t ) ) + θ 1 d S d t + θ 2 d E d t + θ 3 d I d t + θ 4 d A d t + θ 5 d R d t + θ 6 d Q d t .

Moreover, the presence of nontrivial vector functions, denoted as θ ( t ) = ( θ 1 ( t ) , θ 2 ( t ) , θ 3 ( t ) , θ 4 ( t ) , θ 5 ( t ) , θ 6 ( t ) ) , becomes evident, particularly when evaluating ( V ( t ) , z ( t ) ) as the optimal procedure for the initiated control problem:

(25) d V ( t ) d t = ∂ H ( t , V ( t ) , z ( t ) , θ ( t ) ) ∂ θ k , k = 1 , 2 , … , 6 , 0 = ∂ H ( t , V ( t ) , z ( t ) , θ ( t ) ) ∂ z i , i = 1 , 2 , θ k ′ = ∂ H ( t , V ( t ) , z ( t ) , θ ( t ) ) ∂ V ( t ) . k = 1 , 2 , … , 6 ,

describe the transversally condition as θ k ( t ) = 0 , k = 1 , … , 6 . Moreover, the results are obtained through the application of the necessary condition to the Hamiltonian function.□

5.1 Existence of the optimal control problem

Pontryagin’s maximal principle is used to construct the necessary conditions for an optimal solution, which are persuaded by taking into account the initial time t = 0 . We maintain a control that is Lebesgue measurable and constrained in order to accomplish these goals [45,49]. Initial conditions apply to this control, and it is guaranteed that the system’s solution is upward-bounded. The Lagrangian and Hamiltonian are thoroughly examined in order to address the optimal control problem. The optimal control problem is formulated in the Lagrangian framework, and the relevant equation is expressed as follows: dependent constructed system (20) for control strategy, for defining control set Z

(26) L { S ( t ) , E ( t ) , I ( t ) , A ( t ) , R ( t ) , Q ( t ) , z 1 ( t ) , z 2 ( t ) } = B 1 * E ( t ) + B 2 * I ( t ) + B 3 * A ( t ) + 1 2 C 1 * z 1 2 ( t ) + 1 2 C 2 * z 2 2 ( t ) .

Hamiltonian function H of the model is given in the following form:

(27) H = L ( S ( t ) , E ( t ) , I ( t ) , A ( t ) , R ( t ) , Q ( t ) , z 1 ( t ) , z 2 ( t ) ) + θ 1 d S d t + θ 2 d E d t + θ 3 d I d t + θ 4 d A d t + θ 5 d R d t + θ 6 d Q d t .

By substituting the Lagrangian equation and values of all compartments into equation (27), H can be expressed as follows:

(28) H ( V ( t ) , θ ( t ) , z ( t ) ) = B 1 * E ( t ) + B 2 * I ( t ) + B 3 * A ( t ) + 1 2 C 1 * z 1 2 ( t ) + 1 2 C 2 * z 2 2 ( t ) + θ 1 [ Δ − ( 1 − z 2 ) β S ( t ) E ( t ) − ψ S ( t ) + ρ R ( t ) ] + θ 2 [ ( 1 − z 2 ) β S ( t ) E ( t ) − ( ω + λ + ψ + z 1 ) E ( t ) ] + θ 3 [ ω E ( t ) − ( γ + d + ψ + z 2 ) I ( t ) ] + θ 4 [ γ I ( t ) + λ E ( t ) − ( μ + ψ + z 1 ) A ( t ) ] + θ 5 [ μ A ( t ) + d I ( t ) − ( ψ + ρ + α ) R ( t ) + z 2 I ( t ) + z 1 A ( t ) ] + θ 6 [ α R ( t ) − ψ Q ( t ) ] ,

where V ( t ) = ( S ( t ) , E ( t ) , I ( t ) , A ( t ) , R ( t ) , Q ( t ) ) represent a state variable vector, and θ k ( t ) for k = 1 , … 6 refers to the adjoint variables corresponding to the state equations of the system (20). Furthermore, we establish the essential optimality conditions for Eqs (20) and (21):

(29) d V ( t ) d t = ∂ H ( t , V ( t ) , z ( t ) , θ ( t ) ) ∂ θ k , k = 1 , 2 , … , 6 , 0 = ∂ H ( t , V ( t ) , z ( t ) , θ ( t ) ) ∂ z i , i = 1 , 2 . θ k ′ = ∂ H ( t , V ( t ) , z ( t ) , θ ( t ) ) ∂ V ( t ) , k = 1 , 2 , … , 6 .

The transversality condition can be expressed as θ k ( t ) = 0 , where k = 1 , … , 6 .

The resulting system, involving the adjoint variables θ 1 , θ 2 , θ 3 , θ 4 , θ 5 , and θ 6 , along with the optimal control variables z 1 and z 2 , can be expressed as follows:

(30) θ 1 ′ = ( 1 − z 2 ) β E ( t ) ( θ 1 − θ 2 ) + θ 1 ψ , θ 2 ′ = − B 1 * + ( 1 − z 2 ) β S ( t ) ( θ 1 − θ 2 ) + ω ( θ 2 − θ 3 ) + λ ( θ 2 − θ 4 ) + ( ψ + z 1 ) θ 2 , θ 3 ′ = − B 2 * + γ ( θ 3 − θ 4 ) + d ( θ 3 − θ 5 ) + u 2 ( θ 3 − θ 5 ) + ψ θ 3 , θ 4 ′ = − B 3 * + μ ( θ 4 − θ 5 ) + ψ θ 4 + z 1 ( θ 4 − θ 5 ) , θ 5 ′ = α ( θ 5 − θ 6 ) + ψ θ 5 + ρ ( θ 5 − θ 1 ) , θ 6 ′ = ψ θ 6 , z 1 ( t ) = θ 2 E ( t ) + ( θ 4 − θ 5 ) A ( t ) C 1 * , z 2 ( t ) = ( θ 2 − θ 1 ) β S ( t ) E ( t ) + ( θ 3 − θ 5 ) I ( t ) C 2 * .

Theorem 5.2

Given an optimal control z 1 * , z 2 * and solution S * ( t ) , E * ( t ) , I * ( t ) , A * ( t ) , R * ( t ) , and Q * ( t ) of the corresponding state system (20), there exist adjoint variable θ k ( t ) , k = 1 , … , 6 such that:

(31) θ 1 ′ = ( 1 − z 2 ) β E ( t ) ( θ 1 − θ 2 ) + θ 1 ψ , θ 2 ′ = − B 1 * + ( 1 − z 2 ) β S ( t ) ( θ 1 − θ 2 ) + ω ( θ 2 − θ 3 ) + λ ( θ 2 − θ 4 ) + ( ψ + z 1 ) θ 2 , θ 3 ′ = − B 2 * + γ ( θ 3 − θ 4 ) + d ( θ 3 − θ 5 ) + u 2 ( θ 3 − θ 5 ) + ψ θ 3 , θ 4 ′ = − B 3 * + μ ( θ 4 − θ 5 ) + ψ θ 4 + z 1 ( θ 4 − θ 5 ) , θ 5 ′ = α ( θ 5 − θ 6 ) + ψ θ 5 + ρ ( θ 5 − θ 1 ) , θ 6 ′ = ψ θ 6 ,

by employing the transversality criterion θ k ( t ) , k = 1 , … , 6 .

Proof

Consider the assumption that the state variables S ( t ) = S * ( t ) , E ( t ) = E * ( t ) , I ( t ) = I * ( t ) , A ( t ) = A * ( t ) , R ( t ) = R * ( t ) , Q ( t ) = Q * ( t ) . Additionally, the Hamiltonian is described in terms of the state variables S ( t ) , E ( t ) , I ( t ) , A ( t ) , R ( t ) , and Q ( t ) . Subject to the requirement of transversality, the modified adjoint system is now represented as follows:

(32) θ 1 ′ = ( 1 − z 2 ) β E ( t ) ( θ 1 − θ 2 ) + θ 1 ψ , θ 2 ′ = − B 1 * + ( 1 − z 2 ) β S ( t ) ( θ 1 − θ 2 ) + ω ( θ 2 − θ 3 ) + λ ( θ 2 − θ 4 ) + ( ψ + z 1 ) θ 2 , θ 3 ′ = − B 2 * + γ ( θ 3 − θ 4 ) + d ( θ 3 − θ 5 ) + u 2 ( θ 3 − θ 5 ) + ψ θ 3 , θ 4 ′ = − B 3 * + μ ( θ 4 − θ 5 ) + ψ θ 4 + z 1 ( θ 4 − θ 5 ) , θ 5 ′ = α ( θ 5 − θ 6 ) + ψ θ 5 + ρ ( θ 5 − θ 1 ) , θ 6 ′ = ψ θ 6 .

Theorem 5.3

For the region Z, the optimal control pair ( z 1 * ( t ) , z 2 * ( t ) ) is given as follows:

(33) z 1 * ( t ) = max min θ 2 E ( t ) + ( θ 4 − θ 5 ) A ( t ) C 1 * , 1 , 0 , z 2 * ( t ) = max min ( θ 2 − θ 1 ) β S ( t ) E ( t ) + ( θ 3 − θ 5 ) I ( t ) C 2 * , 1 , 0 .

Proof

When optimal conditions were applied, outcomes like these were obtained:

(34) ∂ H ∂ z 1 = C 1 * z 1 ( t ) − θ 2 E ( t ) − θ 4 A ( t ) + θ 5 A ( t ) , ∂ H ∂ z 2 = C 2 * z 2 ( t ) + θ 1 β E ( t ) S ( t ) − θ 2 β S ( t ) E ( t ) − θ 3 I ( t ) + θ 5 I ( t ) .

The control variable is

(35) z 1 * ( t ) = θ 2 E ( t ) + ( θ 4 − θ 5 ) A ( t ) C 1 * , z 2 * ( t ) = ( θ 2 − θ 1 ) β S ( t ) E ( t ) + ( θ 3 − θ 5 ) I ( t ) C 2 * .

It is evident that the state of control space can still be expressed as follows:

z 1 * ( t ) = 0 , if θ 2 E ( t ) + ( θ 4 − θ 5 ) A ( t ) C 1 * ≤ 0 , θ 2 E ( t ) + ( θ 4 − θ 5 ) A ( t ) C 1 * , if 0 < θ 2 E ( t ) + ( θ 4 − θ 5 ) A ( t ) C 1 * < 1 , 1 , if θ 2 E ( t ) + ( θ 4 − θ 5 ) A ( t ) C 1 * ≥ 1 .

z 2 * ( t ) = 0 , if ( θ 2 − θ 1 ) β S ( t ) E ( t ) + ( θ 3 − θ 5 ) I ( t ) C 2 * ≤ 0 , ( θ 2 − θ 1 ) β S ( t ) E ( t ) + ( θ 3 − θ 5 ) I ( t ) C 2 * , if 0 < ( θ 2 − θ 1 ) β S ( t ) E ( t ) + ( θ 3 − θ 5 ) I ( t ) C 2 * < 1 , 1 , if ( θ 2 − θ 1 ) β S ( t ) E ( t ) + ( θ 3 − θ 5 ) I ( t ) C 2 * ≥ 1 .

In a simplified form, the control variables are

(36) z 1 * ( t ) = max min θ 2 E ( t ) + ( θ 4 − θ 5 ) A ( t ) C 1 * , 1 , 0 , z 2 * ( t ) = max min ( θ 2 − θ 1 ) β S ( t ) E ( t ) + ( θ 3 − θ 5 ) I ( t ) C 2 * , 1 , 0 .

The optimal system is now given as,

(37) d S * d t = Δ − 1 − max min ( θ 2 − θ 1 ) β S ( t ) E ( t ) + ( θ 3 − θ 5 ) I ( t ) C 2 * , 1 , 0 β S * ( t ) E * ( t ) − ψ S * ( t ) + ρ R * ( t ) , d E * d t = 1 − max min ( θ 2 − θ 1 ) β S ( t ) E ( t ) + ( θ 3 − θ 5 ) I ( t ) C 2 * , 1 , 0 β S * ( t ) E * ( t ) − ω + λ + ψ + max min θ 2 E ( t ) + ( θ 4 − θ 5 ) A ( t ) C 1 * , 1 , 0 I * ( t ) , d I * d t = ω E * ( t ) − λ + ψ + γ + d + max min ( θ 2 − θ 1 ) β S ( t ) E ( t ) + ( θ 3 − θ 5 ) I ( t ) C 2 * , 1 , 0 E * ( t ) , d A * d t = γ A * ( t ) + λ E * ( t ) − μ + ψ + max min θ 2 E ( t ) + ( θ 4 − θ 5 ) A ( t ) C 1 * , 1 , 0 A * ( t ) , d R * d t = μ A * ( t ) + d I * ( t ) − ( ψ + ρ + α ) R * ( t ) + max min ( θ 2 − θ 1 ) β S ( t ) E ( t ) + ( θ 3 − θ 5 ) I ( t ) C 1 * , 1 , 0 I * ( t ) + max min θ 2 E ( t ) + ( θ 4 − θ 5 ) A ( t ) C 1 * , 1 , 0 A * ( t ) , d Q * d t = α R * ( t ) − ψ Q * ( t ) .

In the context of optimization, the determination of state and control variables entails the utilization of an adjoint variable in conjunction with the optimal system and initial conditions. It is clear that the objective function second-order derivative with respect to the control variables is positive. As a result, the final Hamiltonian has the following official documentation:

□ H * = B 1 * E * ( t ) + B 2 * I * ( t ) + B 3 * A * ( t ) + 1 2 C 1 * z 1 2 ( t ) + 1 2 C 2 * z 2 2 ( t ) + ∑ k = 1 6 θ k g k ( S * ( t ) , E * ( t ) , I * ( t ) , × A * ( t ) , R * ( t ) , Q * ( t ) ) .

5.2 Numerical simulation and discussion

We start solving the optimal control problem that our SEIARQ model presents numerically in this section. An iterative method is used to determine this optimality system. Initiating the process involves starting with an initial estimation for the control variables, which are represented as z 1 and z 2 . The state variables S , E , I , A , R , and Q are then solved forward in an iterative manner, while concurrently solving the adjoint variables θ i for i = 1, 2, 3, 4, 5, 6 backward at time steps k = 0 and k = t .

In this comprehensive analysis, we conducted numerical simulations to demonstrate our theoretical findings, focusing on the SEIARQ model described by system (20). We used detailed parameter values, setting Δ = 0.5 , β = 2.4 , α = 0.2 , ω = 4 , ψ = 0.9 , d = 0.00016590 , μ = 0.9626 , λ = 0.00709 , γ = 0.13 , and ρ = 0.922 . The initial values were chosen as S ( 0 ) = 1.4117800 , I ( 0 ) = 3.000000 , E ( 0 ) = 3.53534 , A ( 0 ) = 5,000 , R ( 0 ) = 2.95000 , Q ( 0 ) = 1.000 . Moreover, the weight constant is set to be B 1 = 0.6610000 , B 2 = 0.54450 , B 3 = 0.0090030 , C 1 = 0.00044440 , and C 2 * = 0.33550 . Figures 7, 8, 9, 10, 11, 12, 13, 14 demonstrate the dynamics of our proposed optimization strategy through simulations with the effect of control parameters. This thorough analysis, which takes into account all six compartments over a period of days, clearly shows the population-wide advantages of the optimal control approach. Notably, Figure 7 shows a noticeable increase in the susceptible population under the influence of control measures. The effectiveness of the implemented control is evident in Figure 8, illustrating the dynamics of the exposed classes. Additionally, Figures 9 and 10 depict a more rapid reduction in the infected and treatment classes compared to the scenario without control. Figures 11 and 12 emphasize the positive impact of the control scheme on the recovered class and individuals engaging in safe practices. Additionally, Figures 13 and 14 painstakingly show how resilient the control variables the social awareness campaign ( z 1 ( t ) ) and the treatment that helps lower the infected population ( z 2 ( t ) ) are in both the presence and absence of control interventions.

$Figure 7 Evolving dynamics of individuals in the susceptible group S ( t ) S\left(t) , comparing scenarios with and without control measures.$

Figure 7

Evolving dynamics of individuals in the susceptible group S ( t ) , comparing scenarios with and without control measures.

$Figure 8 Evolving dynamics of individuals in the exposed group E ( t ) E\left(t) , delineating distinctions between scenarios with and without implemented control measures.$

Figure 8

Evolving dynamics of individuals in the exposed group E ( t ) , delineating distinctions between scenarios with and without implemented control measures.

$Figure 9 Evolving dynamics of individuals in the infected group I ( t ) I\left(t) , comparing scenarios with and without control measures.$

Figure 9

Evolving dynamics of individuals in the infected group I ( t ) , comparing scenarios with and without control measures.

$Figure 10 Evolving dynamics of individuals in the treatment group A ( t ) A\left(t) , comparing scenarios with and without control measures.$

Figure 10

Evolving dynamics of individuals in the treatment group A ( t ) , comparing scenarios with and without control measures.

$Figure 11 Dynamic behavior of individuals in the recovered group R ( t ) R\left(t) , considering both scenarios with and without control measures.$

Figure 11

Dynamic behavior of individuals in the recovered group R ( t ) , considering both scenarios with and without control measures.

$Figure 12 Dynamic behavior of individuals practicing safe sex in the group Q ( t ) Q\left(t) , considering both scenarios with and without control measures.$

Figure 12

Dynamic behavior of individuals practicing safe sex in the group Q ( t ) , considering both scenarios with and without control measures.

$Figure 13 Dynamics of the control variable z 1 {z}_{1} .$

Figure 13

Dynamics of the control variable z 1 .

$Figure 14 Dynamics of the control variable z 2 {z}_{2} .$

Figure 14

Dynamics of the control variable z 2 .

This comprehensive analysis provides insights into the dynamics of disease transmission and the effectiveness of control strategies, contributing to our understanding of optimal interventions for handling syphilis in a given community.

6 Conclusion

We have examined a mathematical model comprising a set of equations that elucidates the dynamics of the syphilis disease. The investigation of this model has been extensively established and it has been determined that the system exhibits both local and global stability. Furthermore, we have engaged in a thorough sensitivity analysis of the parameters associated with the threshold parameter R 0 . It is imperative to acknowledge that the utilization of NSFD scheme for mathematical models which are predicated on a system of differential equations presents a considerably more potent approach in terms of facilitating the computation of convergent solutions for disease models. Finally, the numerical simulation was executed and all the analytical results were validated through numerical methods employing an NSFD scheme. This study pertains to the development of a novel deterministic mathematical for syphilis infection using optimal control theory, aiming to enhance disease management within human populations. The effective control of any disease necessitates the implementation of recommended measures proposed by researchers and disease control experts. Our model incorporates two crucial control variables: an educational awareness campaign, and treatment protocols intended to reduce the prevalence of infected individuals. We framed an optimal control problem, establishing an objective function to determine optimal values for specific variables, thereby minimizing overall costs. By applying Pontryagin’s maximum principle, we derived essential conditions for optimal solutions. The study explored two optimal control approaches and demonstrated their effectiveness through MATLAB simulations. Graphical results depicted that the simultaneous consideration of both control variables effectively reduced the spread of infection. Discussions encompassed scenarios with and without control measures. Rigorous analyses emphasized the efficacy of the control strategy in achieving optimal states across all six considered compartments. A controlled population increase was observed with control implementation, and a sharper decrease in exposed, infected, and treatment classes occurred with optimal control. Optimal controls also led to a comparatively accelerated recovery, emphasizing the cost-effectiveness of time-dependent controls over time-independent ones. This approach was utilized to decrease the rates of infection for both disease-free and endemic equilibria, thus enabling the regulation of the dissemination of the syphilis disease within the community.

In the future, an intriguing extension of the developed model would involve exploring the anticipated benefits include the enhanced ability of the fractional order scheme to utilize available information, providing further insights into disease transmission dynamics in society. Future studies will delve into various disease control methods, such as vaccines and medications, vector control strategies, and overall public health measures. There will be a particular focus on understanding how these strategies interplay, especially when dealing with multiple diseases concurrently. Also to decipher how these strategies synergize or potentially counteract each other, ultimately paving the way for improved global disease control plans.

Acknowledgments

This study is supported via funding from Prince Sattam bin Abdulaziz University project number (PSAU/2025/R/1446).

Funding information: This study is supported via funding from Prince Sattam bin Abdulaziz University project number (PSAU/2024/R/1445).
Author contributions: M.F.: conceptualization, methodology, formal analysis. writing – original draft, writing review editing. K.J.: conceptualization, methodology, formal analysis, software, visualization, writing original draft. K.S.N.: methodology, formal analysis, software, formatting, writing review, editing. Y.N.A.: formal analysis, software, visualization, writing – review editing. M.U.S.: methodology, formal analysis, software, formatting, writing review, editing. E.H.: formal analysis, supervision, visualization, writing – review, editing. All authors have accepted responsibility for the entire content of this manuscript and approved its submission.
Conflict of interest: The authors state no conflict of interest.
Data availability statement: All data generated or analysed during this study are included in this published article.

References

[1] Aadland D, Finnoff DC, Huang KX. Syphilis cycles. BE J Econ Anal Policy. 2013;13(1):297–348. 10.1515/bejeap-2012-0060Search in Google Scholar

[2] Kiarie J, Mishra CK, Temmerman M, Newman L. Accelerating the dual elimination of mother-to-child transmission of syphilis and HIV: Why now? Int J Gynaecol Obstet. 2015;130:S1–3. 10.1016/j.ijgo.2015.05.002Search in Google Scholar PubMed PubMed Central

[3] World Health Organization. The World health report: 2004: changing history. World Health Organization; 2004. Search in Google Scholar

[4] Brown DL, Frank JE. Diagnosis and management of syphilis. Am Fam Physician. 2003;68(2):283–90. Search in Google Scholar

[5] World Population. WorldoMeters. Website: www.worldometers. 2016. Search in Google Scholar

[6] LaFond RE, Lukehart SA. Biological basis for syphilis. Clin Microbiol Rev. 2006;19(1):29–49. 10.1128/CMR.19.1.29-49.2006Search in Google Scholar PubMed PubMed Central

[7] Saloojee H, Velaphi S, Goga Y, Afadapa N, Steen R, Lincetto O. The prevention and management of congenital syphilis: an overview and recommendations. Bull World Health Organ. 2004;82(6):424–30. Search in Google Scholar

[8] Parran T. Syphilis: a public health problem. Science. 1938;87(2251):147–52. 10.1126/science.87.2251.147Search in Google Scholar PubMed

[9] Oswal S, Lyons G. Syphilis in pregnancy. Cont Educ Anaesth Crit Care Pain. 2008;8(6):224–7. 10.1093/bjaceaccp/mkn042Search in Google Scholar

[10] Brown ST, Zaidi A, Larsen SA, Reynolds GH. Serological response to syphilis treatment: a new analysis of old data. JAMA. 1985;253(9):1296–9. 10.1001/jama.1985.03350330094030Search in Google Scholar

[11] Larsen SA, Steiner BM, Rudolph AH. Laboratory diagnosis and interpretation of tests for syphilis. Clin Microbiol Rev. 1995;8(1):1–21. 10.1128/CMR.8.1.1Search in Google Scholar PubMed PubMed Central

[12] Downing CC. Relapse in a case of early syphilis treated with penicillin. Br J Vener Dis. 1950;26(1):31. 10.1136/sti.26.1.31Search in Google Scholar PubMed PubMed Central

[13] Schoch AG, Alexander LJ. Reinfection and relapse after treatment of early syphilis with penicillin; analysis of 137 cases of infectious failure in a total series of 1,105 cases. Arch Dermatol Syphilol. 1949;60(5, Pt. 1):690–700. 10.1001/archderm.1949.01530050052004Search in Google Scholar PubMed

[14] Garnett GP, Aral SO, Hoyle DV, Willard Cates J, Anderson RM. The natural history of syphilis: implications for the transmission dynamics and control of infection. J Sex Transm Dis. 1997;24(4):185–200. Search in Google Scholar

[15] Fenton KA, Breban R, Vardavas R, Okano JT, Martin T, Aral S, et al. Infectious syphilis in high-income settings in the 21st century. Lancet Infect Dis. 2008;8(4):244–53. 10.1016/S1473-3099(08)70065-3Search in Google Scholar PubMed

[16] Ahmed N, Elsonbaty A, Adel W, Baleanu D, Rafiq M. Stability analysis and numerical simulations of spatiotemporal HIV CD4. T cell model with drug therapy. Chaos. 2020;30(8):083122. 10.1063/5.0010541Search in Google Scholar PubMed

[17] Ghori MB, Naik PA, Zu J, Eskandari Z, Naik MU. Global dynamics and bifurcation analysis of a fractional-order SEIR epidemic model with saturation incidence rate. Math Methods Appl Sci. 2022;45(7):3665–88. 10.1002/mma.8010Search in Google Scholar

[18] Xiao Y, Zhao T, Tang S. Dynamics of an infectious diseases with media/psychology induced non-smooth incidence. Math Biosci Eng. 2013;10(2):445–61. 10.3934/mbe.2013.10.445Search in Google Scholar PubMed

[19] Wang A, Xiao Y. A Filippov system describing media effects on the spread of infectious diseases. Nonlinear Anal Hybrid Syst. 2014;11:84–97. 10.1016/j.nahs.2013.06.005Search in Google Scholar

[20] Jan R, Khan MA, Gómez-Aguilar JF. Asymptomatic carriers in transmission dynamics of dengue with control interventions. Optim Control Appl Methods. 2020;41(2):430–47. 10.1002/oca.2551Search in Google Scholar

[21] Garnett GP, Aral SO, Hoyle DV, Willard Cates J, Anderson J. The natural history of syphilis: implications for the transmission dynamics and control of infection. J Sex Transm Dis. 1997;24(4):185–200. 10.1097/00007435-199704000-00002Search in Google Scholar PubMed

[22] Kimbir RA, Udoo MJ, Aboiyar T. A mathematical model for the transmission dynamics of HIV/AIDS in a two-sex population Counseling and Antiretroviral Therapy (ART). J Math Comput Sci. 2012;2(6):1671–84. Search in Google Scholar

[23] Pourbohloul B, Rekart ML, Brunham RC. Impact of mass treatment on syphilis transmission: a mathematical modeling approach. J Sex Transm Dis. 2003;30(4):297–305. 10.1097/00007435-200304000-00005Search in Google Scholar PubMed

[24] Milner FA, Zhao R. A new mathematical model of syphilis. Math Model Nat Phenom. 2010;5(6):96–108. 10.1051/mmnp/20105605Search in Google Scholar

[25] Verma V, Agarwal M, Verma A. A mathematical model for the novel coronavirus with effect of lockdown. Int J Model Simul Sci Comput. 2023;14(03):2350005. 10.1142/S1793962323500058Search in Google Scholar

[26] Verma V. Optimal control analysis of a mathematical model on smoking. Model Earth Syst Environ. 2020;6(4):2535–42. 10.1007/s40808-020-00847-1Search in Google Scholar

[27] Muatjetjeja B, Khalique CM. First integrals for a generalized coupled Lane–Emden system. Nonlinear Anal Real World Appl. 2011;12(2):1202–12. 10.1016/j.nonrwa.2010.09.013Search in Google Scholar

[28] Muatjetjeja B, Khalique CM. Lie group classification for a generalised coupled Lane-Emden system in dimension one. East Asian J Appl Math 2014;4(4):301–11. 10.4208/eajam.080214.230814aSearch in Google Scholar

[29] Muatjetjeja B. Group classification and conservation laws of the generalized Klein-Gordon-Fock equation. Int J Mod Phys B. 2016;30(28–29):1640023. 10.1142/S0217979216400233Search in Google Scholar

[30] Muatjetjeja B. Coupled Lane-Emden-Klein-Gordon-Fock system with central symmetry: Symmetries and conservation laws. J Differ Equ. 2017;263(12):8322–8. 10.1016/j.jde.2017.08.044Search in Google Scholar

[31] Omame A, Abbas M. The stability analysis of a co-circulation model for COVID-19, dengue, and zika with nonlinear incidence rates and vaccination strategies. Healthc Anal. 2023;3:100151. 10.1016/j.health.2023.100151Search in Google Scholar PubMed PubMed Central

[32] Omame A, Raezah AA, Diala UH, Onuoha C. The optimal strategies to be adopted in controlling the co-circulation of COVID-19, Dengue and HIV: Insight from a mathematical model. Axioms. 2023;12(8):773. 10.3390/axioms12080773Search in Google Scholar

[33] Omame A, Abbas M. Modeling SARS-CoV-2 and HBV co-dynamics with optimal control. Phys A Stat Mech Appl. 2023;615:128607. 10.1016/j.physa.2023.128607Search in Google Scholar PubMed PubMed Central

[34] Omame A, Abbas M, Onyenegecha CP. Backward bifurcation and optimal control in a co-infection model for SARS-CoV-2 and ZIKV. Results Phys. 2022;37:105481. 10.1016/j.rinp.2022.105481Search in Google Scholar PubMed PubMed Central

[35] Nadeem Anjam Y, Shahid I, Emadifar H, Arif Cheema S, ur Rahman M. Dynamics of the optimality control of transmission of infectious disease: a sensitivity analysis. Sci Rep. 2024;14(1):1041. 10.1038/s41598-024-51540-7Search in Google Scholar PubMed PubMed Central

[36] Liu X, ur Rahmamn M, Ahmad S, Baleanu D, Nadeem Anjam Y. A new fractional infectious disease model under the non-singular Mittag-Leffler derivative. Wav Random Complex Media. 2022;35(1):1–27. 10.1080/17455030.2022.2036386Search in Google Scholar

[37] Anguelov R, Lubuma JM. Contributions to the mathematics of the nonstandard finite difference method and applications. Numer Methods Partial Differ Equ. 2001;17(5):518–43. 10.1002/num.1025Search in Google Scholar

[38] Mickens RE, editor. Nonstandard finite difference models of differential equations. World Scientific; 1993. 10.1142/2081Search in Google Scholar

[39] Lubuma JM, Patidar KC. Non-standard methods for singularly perturbed problems possessing oscillatory/layer solutions. Appl Math Comput. 2007;187(2):1147–60. 10.1016/j.amc.2006.09.011Search in Google Scholar

[40] Gumel AB. Mathematics of continuous and discrete dynamical systems. Vol. 618. Providence, Rhode Island: American Mathematical Society; 2014. 10.1090/conm/618Search in Google Scholar

[41] Khajji B, Kouidere A, Balatif O, Rachik M. Mathematical modeling, analysis and optimal control of an alcohol drinking model with liver complication. Commun Math Biol Neurosci. 2020;2020:Article-ID 32. Search in Google Scholar

[42] Khajji B, Labzai A, Kouidere A, Balatif O, Rachik M. A discrete mathematical modeling of the influence of alcohol treatment centers on the drinking dynamics using optimal control. J Appl Math. 2020;2020(1):9284698. 10.1155/2020/9284698Search in Google Scholar

[43] Lee S, Jung E, Castillo-Chavez C. Optimal control intervention strategies in low-and high-risk problem drinking populations. Socio-Econ Plan Sci. 2010;44(4):258–65. 10.1016/j.seps.2010.07.006Search in Google Scholar

[44] Zaman G, Kang YH, Jung IH. Optimal treatment of an SIR epidemic model with time delay. Biosyst. 2009;98(1):43–50. 10.1016/j.biosystems.2009.05.006Search in Google Scholar PubMed

[45] Zakary O, Rachik M, Elmouki I. On the analysis of a multi-regions discrete SIR epidemic model: an optimal control approach. Int J Dyn Control. 2017;5:917–30. 10.1007/s40435-016-0233-2Search in Google Scholar PubMed PubMed Central

[46] Hwang CL, Fan LT. A discrete version of Pontryagin’s maximum principle. Oper Res. 1967;15(1):139–46. 10.1287/opre.15.1.139Search in Google Scholar

[47] Guibout V, Bloch A. A discrete maximum principle for solving optimal control problems. In 2004 43rd IEEE Conference on Decision and Control (CDC) (IEEE Cat. No. 04CH37601)2004. Vol. 2, pp. 1806–11. 10.1109/CDC.2004.1430309Search in Google Scholar

[48] Pontryagin LS. Mathematical theory of optimal processes. London: Routledge. 2018. 10.1201/9780203749319Search in Google Scholar

[49] Todorov E, Jordan MI. Optimal feedback control as a theory of motor coordination. Nat Neurosci. 2002;5(11):1226–35. 10.1038/nn963Search in Google Scholar PubMed

Received: 2024-05-16

Revised: 2024-08-22

Accepted: 2025-01-27

Published Online: 2025-07-01

This work is licensed under the Creative Commons Attribution 4.0 International License.

Articles in the same Issue

https://doi.org/10.1515/nleng-2025-0100

Keywords for this article

syphilis disease model; stability analysis; sensitivity analysis; nonstandard finite difference; optimal control strategy

Creative Commons

BY 4.0