Stability Analysis and Optimal Control of a Ψ-Hilfer Fractional Mathematical Model for HSV-2 Transmission Dynamics

Herpes simplex virus type 2 (HSV-2) is one of the most common sexually transmitted infections globally and is a major health concern relative to female sex workers (FSWs), who are at high risk for acquiring and transmitting Human Immunodeficiency Virus (HIV) themselves. HSV-2 is a lifelong infection with genital ulceration symptoms and asymptomatic viral shedding time making for persistent transmission dynamics [1]. Empirical evidence shows that FSW populations are significantly seroprevalent for HSV-2. For instance, in Guangxi, China, a seroprevalence of 54.9% occurred among those tested, where older age, less education, non-Han ethnicity, and lower tier venues (versus higher tiers) for sex work were significantly associated with increased risk of infection [2]. In a prospective cohort study in Yunnan Province, the incidence of HSV-2 was 21.9 per 100 person-years accrued with 83% of female participants indicating they had never tested positive for any sexually transmitted infections (STIs) previously, which was associated with higher risk factors such as early sexual debut, currently having Neisseria gonorrhoeae infection, and not having a regular sexual partner existed [3]. Biologically, HSV-2 infection renders one more susceptible to HIV by breaking down the genital mucosal barrier in addition to recruiting CD4+ T cells - activated by HIV - into the system [4]. Meta-analyses have revealed that HSV-2 infection increases one's risk of HIV acquisition by three times in the general population and two times in the higher-risk populations [4]. At a general level across populations, HSV-2 prevalence has been found to correlate significantly with HIV population prevalence in a positive linear manner suggesting that HSV-2 could be a biological “marker” for whether a population is on track for an HIV epidemic [5]. At a specific level with FSWs, HIV prevalence increases significantly when HSV-2 prevalence reaches at least 25–50% internationally; thus, the two infections are epidemiologically linked on a universal scale [6]. Therefore, proper intervention strategies to incorporate certain prevention and surveillance efforts for HSV-2 (inclusive of prevention efforts based on education, condom use promotion and then with successful antivirals first and eventual vaccination options down the line) will greatly reduce the dual epidemic burden on vulnerable populations.

The use of mathematical modeling represents an indispensable formulated approach for understanding phenomena in the real world across various fields. It takes from real-life endeavors and forms mathematical abstractions to predict and describe discrete and continuous values from chemical reactions to epidemiological dynamics to engineering responses. Traditionally, calculus-based problem-solving approaches based on integer-order rates of change for formerly continuous phenomena have been used to better understand such situations. However, many dynamically occurring processes inherently expose themselves to memory effects or hereditary responses - features that integer-order forecasts can never support. Therefore, fractional differential equations (FDEs) emerged to study the long-range temporal dependence and nonlocal interactions made possible by variably-dependent considerations when attempting to forecast dynamically [7].

Fractional calculus is the field that extends differentiation beyond integer-order - applying certain derivative principles with an order that is not confined solely to whole number values. Fractional approaches are commonly used now in diverse applications that experience dynamics in a non-uniform or history-dependent way. Over the last few decades fractional modeling has increased relevance in various fields for research to empirically account for anomalous diffusion properties, viscoelastic problems/deformations, heterogeneous heat transfers through solid media and biological factors whose kinetics perform according to delayed responses [7], [8]. Since the kernels of fractional derivatives are built into their constructs by fractional nature, they can simulate and generate relationships over time that would otherwise take a process of relaxation to get to with calculus applied to only integer-order considerations.

Within this context, one of the most general approaches is the Hilfer fractional derivative operator. In 1999 Rainer Hilfer proposed a connection between the two classical operators - the Riemann-Liouville and Caputo derivatives - by introducing a type parameter $\beta$ $\in$ [0, 1]. This type of operator works effectively when viscoelastic qualities are crucial for stability predictions, or if the anomalies fit better into nonlinear relaxation dynamics [7], [8], [9]. The theoretical convenience of an operator that maintains initial and boundary conditions appeals to a real-world based approach that applies fractional order difference-like solutions that should fit necessarily for more predictively meaningful systems.

Yet the classical Hilfer form cannot consider a system where the rate of change depends on another function through time. Kilbas and Almeida bridge this gap by providing context for the fractional derivative relative to another function $\Psi$(t) which implies that integration and differentiation occur based on variable kernels [9]. Through this generalization of variable dependence for integer-like approaches (as opposed to fixed kernels) Sousa and Oliveira showed that a $\Psi$-Hilfer fractional derivative is one which takes a monotonically increasing function $\Psi$ into its definition as a kernel [10]. This generalization consolidates other significant operators - Riemann-Liouville, Caputo, Erdélyi-Kober, Katugampola and Hilfer - from a universalized approach where fractional order developments would be relevant [9], [10], [11].

Recent studies have shown that the $\Psi$-Hilfer operator is especially useful for systems with variable-order behaviors, nonlocal considerations or significant memory effects. It has been used to prove existence/uniqueness for solutions to boundary value problems (inspired through case studies) [9] to assess controllability/reachability of higher-dimensional fractional dynamical systems with multiple delays [8], [9], [10], [11], [12] and to confer several types of stability Ulam-Hyers stability through unique solutions [11]. Thus, the theoretical contributions and applied efforts through $\Psi$-Hilfer derivatives in fractional calculus developments make it a crucial tool for theoretical prediction versus something that can be practically applied through real-world endeavors.

The most novel contribution of this manuscript is the development of an HSV-2 transmission model with respect to the $\Psi$-Hilfer fractional derivative which combines memory, heterogeneous temporal dependence and nonlocality in a single mathematical expression. Such advances over integer-order and Caputo-type fractional derivatives allow infection dynamics to depend explicitly on a kernel function $\psi(t)$ but instead, purely independent and random occurrences devoid of biological relevance for persistence, reactivation from latency, and various exposure history that characterize HSV-2 transmission.

Biological considerations of HSV-2, including latency, recurrence, lifelong infection make fractional-order derivatives the operators of choice as integer-order models are memoryless and cannot appropriately model this type of phenomena. The hereditary aspect of the $\Psi$-Hilfer operator is relevant for stability margins, transient response, and control.

While integer-order derivatives capture HSV-2 dynamics effectively, they do not incorporate historical components of STIs, like a delay in response to behavior, relapse and reinfection. The $\Psi$-Hilfer fractional-order system offers convergence at more gradual rates, stability conditions and persistence from long-memory considerations that integer-order systems neither possess nor require.

The structure of the paper is as follows. Section 2 describes important groundwork needed for the $\Psi$-Hilfer fractional derivative operator and its relevant properties. Section 3 discusses HSV-2 transmission dynamics as functionally determined by behavioral efforts but as dynamics determined by $\Psi$-Hilfer fractional order as a justification for infection persistence - which suggests that it should be approached from an operator perspective instead based on social tendencies. Section 4 tests solutions for existence/uniqueness/positivity/boundedness applicable to an appropriate model existence. Section 5 provides a stability investigation of both the infection-free and endemic equilibrium states based on fractional stability theory. Section 6 develops the fractional optimal control framework, where appropriate control strategies are introduced and the corresponding necessary conditions for optimality are derived. Section 7 includes numerical simulations that demonstrate the theoretical results and show how fractional orders influence system dynamics. Finally, Section 8 summarizes the key findings and offers directions for future research.

Let $a, b \in \mathbb{R}$ with $a \lt b$, and let $\psi \in C^1([a, b], \mathbb{R})$ be a positive, strictly increasing function satisfying $\psi^{\prime}(t) \neq 0$ for all $t \in[a, b]$. For a given $f \in C([a, b])$ and fractional order $\alpha>0$, the $\Psi$-Riemann-Liouville fractional integral of order $\alpha$ is defined as:

$I_{a+}^{a ; \psi} f(t)=\frac{1}{\Gamma(\alpha)} \int_a^t \psi^{\prime}(s)[\psi(t)-\psi(s)]^{\alpha-1} f(s) d s$,

where, $\Gamma(\cdot)$ denotes the Euler gamma function [13, 14].

Let $n-1 \lt \alpha \lt n$ with $n \in \mathbb{N}$ and $\beta \in[0,1]$. The left-sided $\Psi$-Hilfer fractional derivative of order $\alpha$ and type

$\beta$ is given by:

$ { }^H D_{a+}^{\alpha, \beta ; \psi} f(t)=I_{a+}^{\beta(n-\alpha) ; \psi}\left(\frac{1}{\psi^{\prime}(t)} \frac{d}{d t}\right)^n I_{a+}^{(1-\beta)(n-\alpha) ; \psi} f(t), t \in[a, b] $,

which generalizes several well-known operators depending on the choice of $\psi$ and $\beta$ [13], [14], [15]. The special cases of the $\Psi$-Hilfer fractional derivative are summarized in Table 1.

•Inverse property

For $f \in C^1([a, b]),{ }^H D_{a+}^{\alpha, \beta ; \psi} I_{a+}^{\alpha ; \psi} f(t)=f(t)$.

•Link between forms

${ }^H D_{a+}^{\alpha, 0 ; \psi} f(t)=D_{a+}^{\alpha ; \psi} f(t),{ }^H D_{a+}^{\alpha, 1 ; \psi} f(t)={ }^C D_{a+}^{\alpha ; \psi} f(t)$.

•Dependence on the kernel $\psi$

The kernel governs the memory law of the system: $\psi(t)=t^\rho$ (power-law memory), $\psi(t)=\ln t$ (logarithmic memory), and $\psi(t)=e^t$ (exponential-type scaling) [15].

•Composition formula

If $\gamma=\alpha+\beta(1-\alpha)$, then

$I_{a+}^{\gamma ; \psi}{ }_H D_{a+}^{\alpha, \beta ; \psi} f(t)=f(t)-\frac{[\psi(t)-\psi(a)]^{\gamma-1}}{\Gamma(\gamma)} I_{a+}^{1-\gamma ; \psi} f(a)$,

which ensures compatibility and continuity of the operator [13], [14].

The $\Psi$-Hilfer operator unifies Riemann–Liouville, Caputo, Hadamard, Katugampola, and Hilfer derivatives within a single mathematical structure. It has been effectively used in epidemiological, physical, and financial models such as the $\Psi$-Hilfer fractional Black–Scholes system and various fractional epidemic models owing to its adaptable kernel $\psi$, which captures multiple memory effects [15]. The $\Psi$-Hilfer fractional derivative definitions and results employed in this paper are in line with the current advances in fractional calculus and epidemic modeling [10], [16], [17].

The total population is divided into classes according to both the infection status and involvement in sex work. It is assumed that all individuals are homogeneously mixed regarding exposure to the virus. At time $t$, the total population $N(t)$ is defined as:

$ N(t)=S_n(t)+E_n(t)+I_n(t)+S_p(t)+E_p(t)+I_p(t) $

where,

•$S_n(t)$: Susceptible non-sex workers,

•$E_n(t)$: Exposed non-sex workers,

•$I_n(t)$: Infectious non-sex workers,

•$S_p(t)$: Susceptible sex workers,

•$E_p(t)$: Exposed sex workers,

•$I_p(t)$: Infectious sex workers.

The biological rationale for the model's compartment structure relies upon the natural history of the virus (exposure, latency, reactivation, recurrent episodes of infection) and the difference between sex workers and non-sex workers is based upon known heterogeneity of contact rates, access to health care, mental stress, and probability of re-infection, which all influence HSV-2 transmission Individuals are recruited into the population at a constant rate $\Lambda$ and enter the susceptible classes. All individuals experience natural death at a constant rate $\mu$.

Sex workers acquire infection at a rate $\lambda_h$, while non-sex workers contract the disease at a reduced rate $\phi \lambda_p$, where $\phi \in(0,1)$ indicates a lower exposure risk. For sex workers, the force of infection is expressed as:

$ \lambda_h=\frac{\gamma c_h\left(\eta I_n+I_p\right)}{N} $

where, $\gamma$ denotes the probability of transmission per sexual act, $c_h$ is the average number of contacts for sex workers, and $\eta \in(0,1)$ reduces the infectivity of $I_n$ relative to $I_p$.

Non-sex workers acquire infection through:

$ \lambda_p=\frac{\rho c_p\left(S_p+E_p+I_p\right)}{N} $

where, $\rho$ is the rate of interaction with sex workers and $c_p$ represents the average number of contacts for non-sex workers. Because sex workers engage in a higher rate of partnerships, they sustain greater transmission potential and higher prevalence levels than non-sex workers.

Differences in parameter values for sex workers vs. non-sex workers are incorporated into the model because, from an epidemiological standpoint, sex workers are more likely to have HSV-2, be infectious longer, and have more chances to be exposed. Differences like these create a feasible transmission dynamic and a biologically plausible temporal evolution.

Individuals exit sex work at rate $\phi$, except for the infectious class $I_p$. Infectious sex workers may leave the profession faster, at rate $\delta>\phi$, due to health or social constraints.

Latent infection progresses to the active stage at rate $\sigma$ among sex workers and at a reduced rate $\pi \sigma$ among non-sex workers $(\pi \in(0,1))$. Infectious individuals recover to latency at rate $r$. For sex workers, the rate of recovery is reduced to $\varepsilon r(\varepsilon \in(0,1))$, accounting for limited healthcare access and possible drug use.

The recurrence and psychosocial burden associated with HSV-2 increase the likelihood of reactivation among sex workers compared to non–sex workers. Disease-induced mortality is neglected, assuming spontaneous recovery or mild treatment. The resulting nonlinear system of differential equations is summarized schematically in Figure 1.

For $t>0$, let $\psi \in C^1([ 0, T])$ be a strictly increasing function such that $\psi^{\prime}(t)>0$. The $\Psi$-Hilfer fractional derivative of order $0<\alpha \leq 1$ and type $0 \leq \beta \leq 1$ is denoted by:

$ { }^H D_{0^{+}}^{\alpha, \beta ; \psi} x(t) $,

where, $\xi=\alpha+\beta(1-\alpha)$.

$ \begin{aligned} { }^H D_{0^{+}}^{\alpha, \beta ; \psi} S_n(t) & =\Lambda-\left(\varphi \lambda_h+\lambda_p\right) S_n+\phi S_p-\mu S_n, \\ { }^H D_{0^{+}}^{\alpha, \beta ; \psi} E_n(t) & =\varphi \lambda_h S_n+\phi E_p+r I_n-\lambda_p E_n-(\pi \sigma+\mu) E_n, \\ { }^H D_{0^{+}}^{\alpha, \beta ; \psi} I_n(t) & =\pi \sigma E_n+\delta I_p-\omega \lambda_p I_n-(r+\mu) I_n, \\ { }^H D_{0^{+}}^{\alpha, \beta ; \psi} S_p(t) & =\lambda_p S_n-\lambda_h S_p-(\mu+\phi) S_p, \\ { }^H D_{0^{+}}^{\alpha, \beta ; \psi} E_p(t) & =\lambda_p E_n+\lambda_h S_p+e r I_p-(\phi+\sigma+\mu) E_p, \\ { }^H D_{0^{+}}^{\alpha, \beta ; \psi} I_p(t) & =\omega \lambda_p I_n+\sigma E_p-(e r+\delta+\mu) I_p . \end{aligned} $

The initial conditions are expressed through the $\psi$-Riemann-Liouville fractional integral:

$ \begin{array}{lll} \left(I_{0^{+}}^{1-\xi ; \psi} S_n\right)(0)=S_{n 0}, & \left(I_{0^{+}}^{1-\xi ; \psi} E_n\right)(0)=E_{n 0}, & \left(I_{0^{+}}^{1-\xi ; \psi} I_n\right)(0)=I_{n 0} \\ \left(I_{0^{+}}^{1-\xi ; \psi} S_p\right)(0)=S_{p 0}, & \left(I_{0^{+}}^{1-\xi ; \psi} E_p\right)(0)=E_{p 0}, & \left(I_{0^{+}}^{1-\xi ; \psi} I_p\right)(0)=I_{p 0} \end{array} $

where all initial states satisfy $S_{n 0}, E_{n 0}, I_{n 0}, S_{p 0}, E_{p 0}, I_{p 0} \geq 0$.

•When $\beta=1$, the operator reduces to the Caputo- $\psi$ fractional derivative, and the initial conditions become:

$ S_n(0)=S_{n 0}, E_n(0)=E_{n 0}, I_n(0)=I_{n 0}, S_p(0)=S_{p 0}, E_p(0)=E_{p 0}, I_p(0)=I_{p 0} . $

•If $\psi(t)=t$, the model becomes the Hilfer fractional system. If $\psi(t)=t$ and $\beta=1$, it simplifies to the Caputo fractional model.

•Common $\psi$-function forms:

-$\psi(t)=t^\rho$ (power-law kernel),

-$\psi(t)=\ln (1+t)$ (logarithmic kernel),

-$\psi(t)=\sin (t)+t$ (nonlinear kernel).

•The $\Psi$-Hilfer framework provides nonlocal memory effects, enhancing the realism of cross-interactions between the human and population groups.

•Here, $\psi(t)$ represents the cumulative effect of past exposures, behavior persistence, immune memory delay, and psychosocial response. For HSV-2, $\psi(t)$ has a clear biological interpretation, as it operates biologically in this case since latency exists and recurrent reactivation means that what happened yesterday, last month, last year, and years ago makes a difference to what is going on today.

The assumptions taken for this part of the study are commonly adopted in the realm of fractional epidemic modeling such that we possess a locally Lipschitz vector field, non-negative parameters are chosen and the state space is biologically realistic. Such standard assumptions possess mathematical significance without compromising the sensibility of such equations in reality.

Model

Let $J=[a, b]$ with $a \lt b$. Consider the $n$-dimensional Cauchy problem

where, $0<\alpha \leq 1,0 \leq \beta \leq 1$, and $\gamma=\alpha+\beta(1-\alpha) \in(0,1]$. The operator ${ }^H D_{a+}^{\alpha, \beta ; \psi}$ is the $\Psi$-Hilfer fractional derivative, and $I_{a+}^{\mu ; \psi}$ is the $\psi$-Riemann-Liouville integral of order $\mu>0$. We assume $\psi \in C^1(J, \mathbb{R})$ is strictly increasing with $\psi^{\prime}(t) \neq 0$ on $J$.

In the HSV-2 system, $X$ stacks the epidemiological compartments (e.g., $S, E, I, \ldots$ ), and $F(t, X)$ encodes the incidence, progression, recovery, and control terms fitted to your biology.

Function space and norm

Define the weighted space:

$ X=C_{1-\gamma ; \psi}\left(J, \mathbb{R}^n\right)=\left\{X: J \rightarrow \mathbb{R}^n \text { continuous : }\|X\|_{1-\gamma ; \psi}:=\sup _{t \in J}(\psi(t)-\psi(a))^{1-\gamma}\|X(t)\|<\infty\right\} $

This choice matches the natural regularity of $\Psi$-Hilfer problems and makes the integral operator below well defined.

Hypotheses on the right-hand side F

•(H1) Carathéodory regularity

For each fixed $x, t \mapsto F(t, x)$ is measurable on $J$; for a.e. fixed $t, x \mapsto F(t, x)$ is continuous. Moreover, there exists an integrable majorant $m \in L^1(J)$ (with respect to $d \psi$ ) and $c \geq 0$ such that

$ \|F(t, x)\| \leq m(t)+c\|x\| \quad \text { for all }(t, x) \in J \times \mathbb{R}^n. $

•(H2) (Global or local) $\Psi$-Lipschitz property

There exists $L \geq 0$ such that

$ \|F(t, x)-F(t, y)\| \leq L\|x-y\| \quad \text { for all } x, y \in \mathbb{R}^n, \text { a.e. } t \in J. $

•(H3) Positivity/viability

For any $X \in \mathbb{R}_{\geq 0}^n, F(t, X)$ satisfies quasi-positivity (off-diagonal nonnegative production terms and bounded loss terms), ensuring the solution components remain nonnegative on $J$ once they start nonnegative. This standard condition for epidemiological models ensures the positive cone is forward invariant under (1) .

Equivalent Volterra integral formulation

Using the standard $\Psi$-Hilfer calculus, (1) is equivalent to the Volterra equation

Theorem 1.

Assume (H1)–(H2) and let

$ \kappa=\frac{L}{\Gamma(\alpha)}(\psi(b)-\psi(a))^\alpha . $

(a) If $\kappa<1$, then (1) admits a unique solution $X \in \mathcal{X}$ on $J$.

(b) If (H1) holds and $F$ is continuous in $x$, then there exists at least one solution $X \in \mathcal{X}$ on $J$. Under (H3) the solution is nonnegative component-wise on $J$.

(c) If $\kappa<1$, the solution depends Lipschitz-continuously on the initial data $X_0$ and on perturbations in $F$.

Proof Sketch

(1) Operator setting. Define $(\mathcal{T} X)(t)$ by the right-hand side of (2). Under (H1) and standard properties of $I^{\alpha ; \psi}, \mathcal{T}: \mathcal{X} \rightarrow \mathcal{X}$ is well defined and maps bounded sets to equicontinuous sets (Arzelà-Ascoli).

(2) Contraction estimate. For $X, Y \in \mathcal{X}$,

$ \|(\mathcal{T} X-\mathcal{T} Y)(t)\| \leq \frac{L}{\Gamma(\alpha)} \int_a^t \psi^{\prime}(s)(\psi(t)-\psi(s))^{\alpha-1}\|X(s)-Y(s)\| d s $

Multiply by $(\psi(t)-\psi(a))^{1-\gamma}$, take sup over $t \in J$, and get

$ \|\mathcal{T} X-\mathcal{T} Y\|_{1-\gamma ; \psi} \leq \frac{L}{\Gamma(\alpha)}(\psi(b)-\psi(a))^\alpha\|X-Y\|_{1-\gamma ; \psi}=\kappa\|X-Y\|_{1-\gamma ; \psi} . $

Thus $\mathcal{T}$ is a contraction when $\kappa<1$, and Banach's theorem gives a unique fixed point.

(3) Existence without smallness. When $\kappa \geq 1, \mathcal{J}$ is continuous and compact on closed balls in $\mathcal{X}$. Schauder's fixed-point theorem gives existence.

(4) Positivity/viability. If $X_0 \geq 0$ and $F$ is quasi-positive, the kernel in (2) is nonnegative, thus $X(t) \geq 0$ for all $t \in J$.

(5) Continuous dependence. For two data sets $\left(X_0, F\right)$ and $\left(\tilde{X}_0, \tilde{F}\right)$ satisfying the same hypotheses with $\kappa<1$,

$ \|X-\tilde{X}\|_{1-\gamma ; \psi} \leq \frac{1}{1-\kappa}\left[\frac{\left\|X_0-\tilde{X}_0\right\|}{\Gamma(\gamma)}+\frac{(\psi(b)-\psi(a))^\alpha}{\Gamma(\alpha)}\|F-\tilde{F}\|\right] . $

Remarks specific to the HSV-2 system

-Incidence $\lambda(t, X)=\beta(t) \Phi(X)$ and transition terms are locally Lipschitz in $X$ on the biologically feasible region.

-If $\kappa<1$ fails, restrict to a subinterval $[a, a+h]$ to restore uniqueness.

-The initial condition $\left(I_{a+}^{1-\gamma ; \psi} X\right)(a)=X_0$ matches (2) exactly.

-The feasible region $\left\{X \geq 0, \sum_i x_i \leq M\right\}$ is positively invariant.

Model, Setting, and Initial Data

Let $\psi \in C^1([ 0, T])$ be strictly increasing with $\psi^{\prime}(t)>0$. Consider the $\psi$-Hilfer fractional derivative ${ }^H D_{0+}^{\alpha, \beta ; \psi}$ of order $0<\alpha \leq 1$ and type $0 \leq \beta \leq 1$, and put $\xi=\alpha+\beta(1-\alpha)$. The state is

$ X(t)=\left(S_n(t), E_n(t), I_n(t), S_p(t), E_p(t), I_p(t)\right)^{\top} \in \mathbb{R}_{\geq 0}^6, $

governed by the system $\mathcal{S}$ given in the main text (with nonnegative parameters such as birth $\Lambda$, removal $\mu$, and transfer rates $\lambda_h, \lambda_p, \phi, \sigma, \omega, \delta, e, r, \pi$, etc.). The admissible $\psi$-Riemann-Liouville type initial data are

$ \left(I_{0+}^{1-\xi ; \psi} S_n\right)(0)=S_{n 0}, \ldots, \left(I_{0+}^{1-\xi ; \psi} I_p\right)(0)=I_{p 0}, S_{n 0}, \ldots, I_{p 0} \geq 0 $

Notation.

We denote by $E_\alpha(\cdot)$ the one-parameter Mittag-Leffler function and by $K_{\alpha, \beta}^\psi(t, \tau) \geq 0$ the resolvent kernel associated with ${ }^H D_{0+}^{\alpha, \beta ; \psi}$ appearing in the variation-of-constants representation.

Auxiliary Ingredients in the $\psi$-Hilfer Framework

We rely on the following standard tools adapted to ${ }^H D_{0+}^{\alpha, \beta ; \psi}$:

•Comparison / Maximum Principle

If $y$ solves ${ }^H D_{0+}^{\alpha, \beta ; \psi} y(t) \leq f(t, y)$ with $f$ nondecreasing in $y$ and $y(0) \geq 0$, then $y(t) \geq 0$ for $t \geq 0$. In compartmental systems with nonnegative transfers and linear loss in each coordinate, the positive orthant is forward invariant; see the comparison principles for nonlocal derivatives in.

•Linear $\psi$-Hilfer Equation

The solution to ${ }^H D_{0+}^{\alpha, \beta ; \psi} y(t)+\lambda y(t)=g(t)$ with $y(0) \geq 0$ and $g(t) \geq 0$ satisfies

$ 0 \leq y(t) \leq y(0) E_\alpha\left(-\lambda(\psi(t)-\psi(0))^\alpha\right)+\int_0^t K_{\alpha, \beta}^\psi(t, \tau) g(\tau) \mathrm{d} \tau, $

where, $K_{\alpha, \beta}^\psi \geq 0$.

•Fractional Grönwall-Bellman ($\psi$-version)

If ${ }^H D_{0+}^{\alpha, \beta ; \psi} u(t) \leq a-b u(t)$ with $a, b>0$, then

$ u(t) \leq u(0) E_\alpha\left(-b(\psi(t)-\psi(0))^\alpha\right)+\frac{a}{b}\left[ 1-E_\alpha\left(-b(\psi(t)-\psi(0))^\alpha\right)\right] . $

Consequently, $u(t) \leq \max \{u(0), a / b\}$ and $u(t) \rightarrow a / b$ in the fractional sense.

Theorem 2.

For system $\mathcal{S}$ with nonnegative $\psi-R L$ initial data, the nonnegative orthant $\mathbb{R}_{\geq 0}^6$ is forward invariant. In particular, $S_n, E_n, I_n, S_p, E_p, I_p \geq 0$ for all $t \geq 0$.

Proof.

Check each coordinate on the boundary where that coordinate equals zero; all parameters are nonnegative:

$ \begin{aligned} \left.{ }^H D_{0+}^{\alpha, \beta ; \psi} S_n\right|_{S_n=0} & =\Lambda+\phi S_p \geq 0, \\ \left.{ }^H D_{0+}^{\alpha, \beta ; \psi} E_n\right|_{E_n=0} & =\phi E_p+r I_n+\phi \lambda_h S_n \geq 0, \\ \left.{ }^H D_{0+}^{\alpha, \beta ; \psi} I_n\right|_{I_n=0} & =\pi \sigma E_n+\delta I_p \geq 0, \\ \left.{ }^H D_{0+}^{\alpha, \beta ; \psi} S_p\right|_{S_p=0} & =\lambda_p S_n \geq 0, \\ \left.{ }^H D_{0+}^{\alpha, \beta ; \psi} E_p\right|_{E_p=0} & =\lambda_p E_n+\lambda_h S_p+e r I_p \geq 0, \\ \left.{ }^H D_{0+}^{\alpha, \beta ; \psi} I_p\right|_{I_p=0} & =\omega \lambda_p I_n+\sigma E_p \geq 0 . \end{aligned} $

All other terms are outflows proportional to the corresponding coordinate and thus vanish on that boundary. Therefore, by the fractional comparison/maximum principle for nonlocal operators (applied componentwise), the positive orthant is invariant and all components remain nonnegative for $t \geq 0$.

Theorem 3.

Let $N(t)=S_n+E_n+I_n+S_p+E_p+I_p$.

The total population satisfies the linear fractional balance law ${ }^H D_{0+}^{\alpha, \beta ; \psi} N(t)=\Lambda-\mu N(t)$. Consequently, $0 \leq N(t) \leq N(0) E_\alpha \left(-\mu(\psi(t)-\psi(0))^\alpha\right)+\frac{\Lambda}{\mu}\left[ 1-E_\alpha\left(-\mu(\psi(t)-\psi(0))^\alpha\right)\right]$, so $N(t) \leq \max \{N(0), \Lambda / \mu\}$ for all $t \geq 0$ and $N(t) \rightarrow \Lambda / \mu$ in the fractional sense.

Proof.

Summing the six equations of S cancels all inter-compartment transfers, leaving births and natural losses only:

$ { }^H D_{0+}^{\alpha_i \beta ; \psi} N(t)=\Lambda-\mu\left(S_n+E_n+I_n+S_p+E_p+I_p\right)=\Lambda-\mu N(t) . $

By the variation-of-constants formula for ${ }^H D_{0+}^{\alpha, \beta ; \psi}$ with nonnegative kernel $K_{\alpha, \beta}^\psi$ and $g(t) \equiv \Lambda \geq 0$, we obtain the stated estimate using the Mittag-Leffler response to the linear damping $\mu$.

Corollary 1. Define the feasible region $\Omega=\left\{X \in \mathbb{R}_{\geq 0}^6: 0 \leq N \leq \Lambda / \mu\right\}$. Then $\Omega$ is positively invariant for $\mathcal{S}$. Every solution with nonnegative $\psi-R L$ initial data enters and remains in $\Omega$, and each coordinate is uniformly bounded by $\Lambda / \mu$.

Proof. Combine Theorem 2 (nonnegativity) with Theorem 3 (total population bound).

Remark 1. The arguments require only that $\psi$ is $C^1$, strictly increasing with $\psi^{\prime}>0$. Under these conditions, the resolvent kernels associated with ${ }^H D_{0+}^{\alpha_\beta \beta ; \psi}$ are nonnegative and completely monotone, which ensures comparison and Grönwall-type bounds.

Remark 2. This technique matches the standard blueprint in fractional epidemiological modeling, where positivity and boundedness are established before the reproduction number and equilibrium analysis.

We present a self-contained stability analysis package for a six-compartment $\psi$-Hilfer fractional-order epidemic system. The discussion covers: (i) feasible region and positivity; (ii) the disease-free and endemic equilibria; (iii) construction of the next-generation matrix and the basic reproduction number $\mathcal{R}_0$; (iv) Jacobian matrices and local (Mittag-Leffler) stability via the generalized Matignon angle condition for $\psi$-Hilfer dynamics; and (v) global stability of both the disease-free and endemic equilibria via convex Volterra-type Lyapunov functionals and a LaSalle-type invariance principle adapted to the $\psi$-Hilfer setting. The results are stated in a form ready to be integrated into applications and numerical studies. The authors apply relevant procedures from the literature to fractional dynamical systems based on $\psi$-Hilfer operators, Mittag-Leffler stability, generalized Matignon-type stability conditions and other fractional epidemic models and biological systems.

Let $X(t)=\left(S_n, E_n, I_n, S_p, E_p, I_p\right)^{\top}$ denote the state vector. For order $0<\alpha \leq 1$, type $0 \leq \beta \leq 1$, and $\xi=\alpha+\beta(1-\alpha)$, consider the $\psi$-Hilfer fractional system (for $t>0$ ):

$ \begin{aligned} { }^H D_{0+}^{\alpha, \beta ; \psi} S_n & =\Lambda-\left(\varphi \lambda_h+\lambda_p\right) S_n+\phi S_p-\mu S_n, \\ { }^H D_{0+}^{\alpha, \beta ; \psi} E_n & =\varphi \lambda_h S_n+\phi E_p+r I_n-\lambda_p E_n-(\pi \sigma+\mu) E_n, \\ { }^H D_{0+}^{\alpha, \beta ; \psi} I_n & =\pi \sigma E_n+\delta I_p-\omega \lambda_p I_n-(r+\mu) I_n, \\ { }^H D_{0+}^{\alpha, \beta ; \psi} S_p & =\lambda_p S_n-\lambda_h S_p-(\mu+\phi) S_p, \\ { }^H D_{0+}^{\alpha, \beta ; \psi} E_p & =\lambda_p E_n+\lambda_h S_p+e r I_p-(\phi+\sigma+\mu) E_p, \\ { }^H D_{0+}^{\alpha, \beta ; \psi} I_p & =\omega \lambda_p I_n+\sigma E_p-(e r+\delta+\mu) I_p, \end{aligned} $

subject to $\psi$-Riemann-Liouville type initial data $\left(I_{0+}^{1-\xi ; \psi} X\right)(0)=X_0 \geq 0$ componentwise. The forces of infection $\lambda_h, \lambda_p$ are assumed to be nonnegative, locally Lipschitz functions that vanish at the disease-free state and depend on the infectious variables (e.g., $I_n, I_p$ ). All parameters are nonnegative constants with $\Lambda, \mu>0$.

Remark 1. (Feasible region and well-posedness). Let $\Omega=\left\{X \in \mathbb{R}_{\geq 0}^6: X\right.$ bounded $\}$. Under standard local Lipschitz and linear-growth conditions on the right-hand side and for a strictly increasing $\psi$, the system admits a unique nonnegative solution for each nonnegative initial datum, and $\Omega$ is positively invariant.

Remark 2. Positivity and boundedness (and hence existence/uniqueness) follow from standard $\psi$-Hilfer well-posedness arguments and comparison principles.

Disease-Free Equilibrium (DFE)

At an equilibrium with no infection $\left(E_n=I_n=E_p=I_p=0\right)$, we have $\lambda_h=\lambda_p=0$. Solving the remaining algebraic balance:

$ 0=\Lambda+\phi S_p-\mu S_n, 0=-(\mu+\phi) S_p , $

gives $S_p^0=0$ and $S_n^0=\Lambda / \mu$. Therefore the DFE is

$ \varepsilon_0=(\Lambda / \mu, 0,0,0,0,0) . $

Endemic Equilibria

An endemic equilibrium $\varepsilon^*=\left(S_n^*, E_n^*, I_n^*, S_p^*, E_p^*, I_p^*\right) \gg 0$ solves the steady (algebraic) system with $\lambda_h^*=\lambda_h\left(I_n^*, I_p^*\right)$ and $\lambda_p^*=\lambda_p\left(I_n^*, I_p^*\right)>$ 0. Closed forms depend on the explicit incidences $\lambda_h, \lambda_p$; the stability conditions below are stated in terms of the Jacobian evaluated at the equilibrium.

Let $a_n=\frac{\partial \lambda_h}{\partial I_n}, a_p=\frac{\partial \lambda_h}{\partial I_p}, b_n=\frac{\partial \lambda_p}{\partial I_n}, b_p=\frac{\partial \lambda_p}{\partial I_p}$ evaluated at the relevant state. At the DFE, consider the infectious subvector $Z= \left(E_n, I_n, E_p, I_p\right)^{\top}$. The linearization gives

$ { }^H D_{0+}^{\alpha, \beta ; \psi} Z=(F-V) Z, J_{\text {inf }}\left(\varepsilon_0\right)=F-V, $

where the new infection and transition matrices are

$ F=\left(\begin{array}{cccc} \varphi \frac{\Lambda}{\mu} a_n & \varphi \frac{\Lambda}{\mu} a_p & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{array}\right), \quad V=\left(\begin{array}{cccc} \pi \sigma+\mu & -r & -\phi & 0 \\ -\pi \sigma & r+\mu & 0 & -\delta \\ 0 & 0 & \phi+\sigma+\mu & -e r \\ 0 & 0 & -\sigma & e r+\delta+\mu \end{array}\right) . $

Following the van den Driessche–Watmough construction, the basic reproduction number is

$ \mathcal{R}_0=\rho\left(F V^{-1}\right), $

the spectral radius of $F V^{-1}$ evaluated at the DFE. For explicit $\lambda_h, \lambda_p$ (e.g. bilinear or saturated incidence) one obtains a closed-form expression for $\mathcal{R}_0$ by substituting $a_n, a_p, b_n, b_p$ and simplifying $F V^{-1}$.

For $\psi$-Hilfer systems of order $0<\alpha \leq 1$ with strictly increasing $\psi$, the generalized Matignon angle condition asserts: the equilibrium of the linearized system is (Mittag-Leffler) asymptotically stable iff all eigenvalues $\lambda_i$ of the Jacobian satisfy

$ \left|\arg \left(\lambda_i\right)\right|>\alpha \pi / 2 $

Theorem 4. (Local stability of the DFE)

If $\mathcal{R}_0<1$, then all characteristic roots of $J_{\text {inf }}\left(\varepsilon_0\right)=F-V$ satisfy $\left|\arg \left(\lambda_i\right)\right|>\alpha \pi / 2$. Hence the DFE $\varepsilon_0$ is locally asymptotically stable in the $\psi$-Hilfer sense. If $\mathcal{R}_0>1, J_{\mathrm{inf}}\left(\mathcal{E}_0\right)$ admits an eigenvalue with $|\arg (\lambda)| \leq \alpha \pi / 2$ and the $D F E$ is unstable.

Sketch. Using $\mathcal{R}_0<1$ is equivalent to $\rho\left(F V^{-1}\right)<1$, which in turn implies sectorial location of the spectrum of $F-V$ in the region $|\arg (\lambda)|>\alpha \pi / 2$; see the $\psi$-Hilfer linear theory and sectorial resolvent bounds.

Theorem 5. (Local stability of an endemic equilibrium)

Let $\mathcal{E}^* \gg 0$ be an endemic equilibrium. If all eigenvalues of $J\left(\mathcal{E}^*\right)$ satisfy $\left|\arg \left(\lambda_i\right)\right|>\alpha \pi / 2$, then $\mathcal{E}^*$ is locally asymptotically stable.

We utilize convex Volterra-type Lyapunov functionals and a LaSalle-type invariance principle adapted to $\psi$-Hilfer dynamics.

Lemma 1. (Fractional Dini derivative bound)

Let $V: \Omega \rightarrow \mathbb{R}_{\geq 0}$ be $C^1$. Along solutions of ${ }^H D_{0+}^{\alpha, \beta ; \psi} V(X(t)) \leq \nabla V(X(t)) \cdot f(X(t))$, where $f$ denotes the right-hand side of the system.

Global Stability of the DFE for $\mathcal{R}_0<1$

Define $V_0=E_n+I_n+E_p+I_p$. Using Lemma 1, the model equations, and that $\lambda_h, \lambda_p \geq 0$ vanish at the DFE, one obtains the differential inequality

$ { }^H D_{0+}^{\alpha, \beta ; \psi} V_0 \leq-c V_0 $

for some $c>0$ whenever $\mathcal{R}_0<1$. Hence $V_0(t) \rightarrow 0$ with Mittag-Leffler rate in the $\psi$-Hilfer sense, which forces $\left(E_n, I_n, E_p, I_p\right) \rightarrow (0,0,0,0)$ and $\left(S_n, S_p\right) \rightarrow(\Lambda / \mu, 0)$.

Theorem 6. (Global asymptotic stability of the DFE)

If $\mathcal{R}_0<1$, then the DFE $\varepsilon_0$ is globally asymptotically stable in $\Omega$.

Global Stability of the Endemic Equilibrium for $\mathcal{R}_0>1$

Suppose an endemic equilibrium $\varepsilon^* \gg 0$ exists. Consider the convex Volterra functional

$ V^*=\sum_{X \in\left\{S_n, E_n, I_n, S_p, E_p, I_p\right\}} c_X\left(X-X^*-X^* \ln \frac{X}{X^*}\right), $

with positive constants $c_X$ chosen to symmetrize infection-transfer couplings (e.g., $c_{E_n}=\pi \sigma$, $c_{I_n}=r+\mu, c_{E_p}=\sigma, c_{I_p}=e r+\delta+\mu$, while $c_{S_n}, c_{S_p}$ cancel the bilinear infection terms; explicit formulas depend on the chosen $\lambda_h, \lambda_p$ ). From Lemma 1 and equilibrium relations,

$ { }^H D_{0+}^{\alpha, \beta ; \psi} V^*(X(t)) \leq-W\left(X(t)-X^*\right), $

for some positive definite form $W$. Therefore $V^*(t)$ is nonincreasing and fractional LaSalle implies convergence to the largest invariant subset of $\left\{{ }^H D V^*=0\right\}$, namely $\left\{\varepsilon^*\right\}$.

Theorem 7. (Global asymptotic stability of the endemic equilibrium)

If $\mathcal{R}_0>1$ and an endemic equilibrium $\mathcal{E}^* \gg 0$ exists, then $\mathcal{E}^*$ is globally asymptotically stable in $\Omega$.

Practical Notes

•To obtain a closed-form $\mathcal{R}_0$, substitute explicit incidences in $F$ (e.g., bilinear $\lambda_h= \left.\beta_h I_p / N, \lambda_p=\beta_p I_n / N\right)$ and compute $\rho\left(F V^{-1}\right)$.

•The local stability conditions reduce to checking the sector constraint $\left|\arg \left(\lambda_i\right)\right|> \alpha \pi / 2$ for eigenvalues of the Jacobian at the equilibrium of interest.

•The Lyapunov constructions above are robust to common incidence forms (bilinear, standard, saturated), provided $\lambda_h, \lambda_p$ are locally Lipschitz and vanish at the DFE.

Optimal control theory is used here to evaluate biologically feasible intervention control efforts for HSV-2 reduction. The control variables are preventive education levels, therapeutic antiviral treatment and behavioral adjustment— all tangible and frequently used in public health. Optimal control identifies the most cost-effective strategies amidst limited resources.

To mitigate HSV-2 transmission, we formulate an optimal control problem governed by the $\Psi$-Hilfer fractional derivative. Let

$ { }^H D_{t_0^{+}}^{\alpha, \beta ; \Psi} x(t) $

denote the $\Psi$-Hilfer fractional derivative of order $0<\alpha<1$ and type $0 \leq \beta \leq 1$ with respect to the increasing function $\Psi(t)$. For each state variable $x(t)$, the initial condition is defined through the $\Psi$-Riemann-Liouville integral:

$ I_{t_0^{+}}^{(1-\alpha)(1-\beta) ; \Psi} x\left(t_0\right)=x_0 $

Let $S_n, E_n, I_n, S_p, E_p, I_p$ denote the susceptible, exposed, and infected populations of the two interacting subgroups (“normal” and “promiscuous”). Introducing three time-dependent control measures:

•$u_1(t)$: educational/preventive awareness campaigns;

•$u_2(t)$: antiviral treatment for infected individuals;

•$u_3(t)$: behavioral interventions (e.g., condom distribution, counseling).

The controlled $\Psi$-Hilfer fractional system is formulated as

$ \left\{\begin{array}{l} { }^H D_{t_0^{+}}^{\alpha, \beta ; \Psi} S_n=\Lambda_n-\left(1-u_1\right) \beta_n S_n I_p-\mu_n S_n, \\ { }^H D_{t_0^{+}}^{\alpha, \beta ; \Psi} E_n=\left(1-u_1\right) \beta_n S_n I_p-\left(\mu_n+\sigma_n\right) E_n, \\ { }^H D_{t_0^{+}}^{\alpha, \beta ; \Psi} I_n=\sigma_n E_n-\left(\mu_n+\rho_n+u_2\right) I_n, \\ { }^H D_{t_0^{+}}^{\alpha, \beta ; \Psi} S_p=\Lambda_p-\left(1-u_3\right) \beta_p S_p I_n-\mu_p S_p, \\ { }^H D_{t_0^{+}}^{\alpha, \beta ; \Psi} E_p=\left(1-u_3\right) \beta_p S_p I_n-\left(\mu_p+\sigma_p\right) E_p, \\ { }^H D_{t_0^{+}}^{\alpha, \beta ; \Psi} I_p=\sigma_p E_p-\left(\mu_p+\rho_p+u_2\right) I_p . \end{array}\right. $

All parameters are positive and biologically meaningful, and $u_i(t) \in[ 0,1]$ for $i=1,2,3$.

The aim is to minimize both infection prevalence and control costs over a finite horizon $\left[t_0, T\right]:$

$ J\left(u_1, u_2, u_3\right)=\int_{t_0}^T\left[A_n I_n(t)+A_p I_p(t)+1 / 2\left(B_1 u_1^2+B_2 u_2^2+B_3 u_3^2\right)\right] d \Psi(t) $

where, $A_n, A_p$ weight the epidemiological burden and $B_i>0$ denote the unit implementation costs.

The optimal control vector $U^*=\left(u_1^*, u_2^*, u_3^*\right)$ satisfies

$ J\left(U^*\right)=\min _{U \in \mathcal{U}} J(U), \mathcal{U}=\left\{U(t): 0 \leq u_i(t) \leq 1, t \in\left[t_0, T\right]\right\} . $

The existence of $U^*$ is guaranteed since:

•$\mathcal{U}$ is closed, convex, and bounded in $L^{\infty}$;

•the right-hand sides of the state equations are linear in the controls and Lipschitz in the states;

•the integrand

$ L\left(I_n, I_p, u_1, u_2, u_3\right)=A_n I_n+A_p I_p+1 / 2\left(B_1 u_1^2+B_2 u_2^2+B_3 u_3^2\right) $

is convex and coercive in $\left(u_1, u_2, u_3\right)$.

Therefore, at least one optimal control $U^*$ and corresponding state trajectory $X^*$ exist on $\left[t_0, T\right]$

Define adjoint variables $\lambda_i(t)$ corresponding to each state equation and the Hamiltonian

$ H(X, U, \lambda, t)=A_n I_n+A_p I_p+1 / 2\left(B_1 u_1^2+B_2 u_2^2+B_3 u_3^2\right)+\sum_{i=1}^6 \lambda_i f_i(X, U, t) $

where, $f_i(X, U, t)$ are the right-hand sides of the state system. Applying the fractional Pontryagin maximum principle the adjoint equations under the $\Psi$-Hilfer operator are

$ { }^H D_{T^{-}}^{\alpha, \beta ; \Psi} \lambda_i(t)=-\frac{\partial H}{\partial x_i}, \lambda_i(T)=0, i=1, \ldots, 6 . $

For example,

$ \begin{aligned} { }^H D_{T^{-}}^{\alpha, \beta ; \Psi} \lambda_1 & =\lambda_1\left[\left(1-u_1\right) \beta_n I_p+\mu_n\right]-\lambda_2\left(1-u_1\right) \beta_n I_p, \\ { }^H D_{T^{-}}^{\alpha, \beta ; \Psi} \lambda_3 & =-A_n+\lambda_3\left(\mu_n+\rho_n+u_2\right)-\lambda_6\left(1-u_3\right) \beta_p S_p+\cdots, \end{aligned} $

and similarly for the remaining variables.

The transversality conditions arise from the finite control horizon and ensure that any remaining epidemiological or economic benefit is not discounted post-terminal time. They are required for biological interpretation and implementability of the optimal control solutions.

The necessary conditions

$ \frac{\partial H}{\partial u_i}=0 $

yield

$ \begin{aligned} & u_1^*(t)=\max \left\{0, \min \left(\frac{\beta_n S_n I_p\left(\lambda_1-\lambda_2\left(1-p_n\right)-\lambda_4 p_n\right)}{B_1}, 1\right)\right\}, \\ & u_2^*(t)=\max \left\{0, \min \left(\frac{\left(\lambda_3-\lambda_6\right)\left(I_n+I_p\right)}{B_2}, 1\right)\right\}, \\ & u_3^*(t)=\max \left\{0, \min \left(\frac{\beta_p S_p I_n\left(\lambda_4-\lambda_5\right)}{B_3}, 1\right)\right\} . \end{aligned} $

These controls minimize $H$ pointwise almost everywhere on $\left[t_0, T\right]$.

The optimized $u_1^*, u_2^*, u_3^*$ describe the cost-effective balance between prevention, treatment, and behavioral intervention. The fractional order $\alpha$ and type $\beta$ encode memory and hereditary effects of disease transmission, providing a more realistic control strategy than integer-order models.

Consider the $\Psi$-Hilfer fractional system

$ { }^H D_{\Psi}^{\alpha, \beta} x(t)=f(t, x(t)), \quad 0<\alpha \leq 1,0 \leq \beta \leq 1, $

subject to the $\Psi$-Riemann-Liouville integral initial conditions

$ I_{\Psi\left(t_0\right)}^{(1-\alpha)(1-\beta)} x\left(t_0\right)=x_0 $

where, $x(t)=\left[S_n, E_n, I_n, S_p, E_p, I_p\right]^{\top}$ and $\Psi(t)$ is a strictly increasing function as defined in the uploaded document.

Applying the $\Psi$-Riemann-Liouville fractional integral to both sides yields

$ x(t)=x_0+\frac{1}{\Gamma(\alpha)} \int_0^t(\Psi(t)-\Psi(\tau))^{\alpha-1} \Psi^{\prime}(\tau) f(\tau, x(\tau)) d \tau $

This form allows discretization on a uniform $\Psi$–grid.

Let $t_0=0 \lt t_1 \lt … \lt t_N = T$ and define $\tau_n=\Psi\left(t_n\right)$ with uniform step $\Delta_\psi=\tau_{n+1}-\tau_n$. The grid ensures proper weighting of the fractional kernel.

For $n=0,1, \ldots, N-1$, the predictor-corrector steps are:

(a) Predictor (Adams–Bashforth)

$ \hat{x}_{n+1}=x_0+\frac{\left(\Delta_{\Psi}\right)^\alpha}{\Gamma(\alpha+1)} \sum_{k=0}^n b_{n-k}^{(\alpha)} f\left(t_k, x_k\right), $

where, $b_j^{(\alpha)}=(j+1)^\alpha-j^\alpha$.

(b) Evaluate

Compute $\hat{f}_{n+1}=f\left(t_{n+1}, \hat{x}_{n+1}\right)$.

(c) Corrector (Adams-Moulton)

$ x_{n+1}=x_0+\frac{\left(\Delta_{\Psi}\right)^\alpha}{\Gamma(\alpha+2)}\left[\hat{f}_{n+1}+\sum_{k=0}^n a_{n-k}^{(\alpha)} f\left(t_k, x_k\right)\right], $

with $a_j^{(\alpha)}=(j+2)^{\alpha+1}-2(j+1)^{\alpha+1}+j^{\alpha+1}$.

(d) Optional Re-evaluation

Set $f_{n+1}=f\left(t_{n+1}, x_{n+1}\right)$ for the next iteration.

Each compartment evolves as

$ \begin{aligned} S_n^{n+1} & =S_{n, 0}+C_\alpha \sum_{k=0}^n w_{n-k} F_{S_n}\left(t_k, x_k\right) \\ E_n^{n+1} & =E_{n, 0}+C_\alpha \sum_{k=0}^n w_{n-k} F_{E_n}\left(t_k, x_k\right) \\ I_n^{n+1} & =I_{n, 0}+C_\alpha \sum_{k=0}^n w_{n-k} F_{I_n}\left(t_k, x_k\right) \\ S_p^{n+1} & =S_{p, 0}+C_\alpha \sum_{k=0}^n w_{n-k} F_{S_p}\left(t_k, x_k\right) \\ E_p^{n+1} & =E_{p, 0}+C_\alpha \sum_{k=0}^n w_{n-k} F_{E_p}\left(t_k, x_k\right) \\ I_p^{n+1} & =I_{p, 0}+C_\alpha \sum_{k=0}^n w_{n-k} F_{I_p}\left(t_k, x_k\right) \end{aligned} $

where, $C_\alpha=\left(\Delta_{\Psi}\right)^\alpha / \Gamma(\alpha+1)$.

Initialize parameters: $\alpha, \beta, N, T, \Psi(t)$.

Compute $\tau_n=\Psi\left(t_n\right)$ and $\Delta_{\Psi}=(\Psi(T)-\Psi(0)) / N$.

Set $x_0$ and evaluate $f_0=f\left(t_0, x_0\right)$.

Predictor: $\hat{x}_{n+1}=x_0+K_1 \sum_{k=0}^n b_{n-k} f_k$, where $K_1=\left(\Delta_{\Psi}\right)^\alpha / \Gamma(\alpha+1)$.

Evaluate: $\hat{f}_{n+1}=f\left(t_{n+1}, \hat{x}_{n+1}\right)$.

Corrector: $x_{n+1}=x_0+K_2\left(\hat{f}_{n+1}+\sum_{k=0}^n a_{n-k} f_k\right)$, where $K_2=\left(\Delta_{\Psi}\right)^\alpha / \Gamma(\alpha+2)$.

Update: $f_{n+1}=f\left(t_{n+1}, x_{n+1}\right)$.

The scheme converges with global order $\mathcal{O}\left(\left(\Delta_{\Psi}\right)^{1+\alpha}\right)$. Positivity and boundedness of each compartment are preserved for biologically consistent right-hand sides $F$.

The objective of the numerical simulations is to investigate the impact of fractional order $\alpha$, type $\beta$, kernel selection $\psi(t)$ on the dynamics of HSV-2 transmission. The parameters are selected within biologically significant ranges found in the literature which are then cross-examined in simulations to demonstrate deviation from classical integer-order.

The $\Psi$-Hilfer fractional-order transmission model for HSV-2 is implemented and analyzed using the PYTHON application to explore memory and hereditary aspects of disease transmission dynamics. To accommodate the challenges of fractional-order derivatives, the formulated system is solved mumerically by the Adams-Bashforth-Moulton PECE approach, which ensures numerical stability and improved understanding of disease transmission under fractional-order and initial condition variations.

The $\Psi$-Hilfer fractional-order describing HSV-2 transmission model $0<\alpha \leq 1$, type $0 \leq \beta \leq 1$, with kernel $\psi(t)$,

$ \begin{aligned} { }^H D_{0^{+}}^{\alpha, \beta ; \psi} S_n(t) & =\Lambda-\left(\varphi \lambda_h+\lambda_p\right) S_n+\phi S_p-\mu S_n, \\ { }^H D_{0^{+}}^{\alpha, \beta ; \psi} E_n(t) & =\varphi \lambda_h S_n+\phi E_p+r I_n-\lambda_p E_n-(\pi \sigma+\mu) E_n, \\ { }^H D_{0^{+}}^{\alpha, \beta ; \psi} I_n(t) & =\pi \sigma E_n+\delta I_p-\omega \lambda_p I_n-(r+\mu) I_n, \\ { }^H D_{0^{+}}^{\alpha, \beta ; \psi} S_p(t) & =\lambda_p S_n-\lambda_h S_p-(\mu+\phi) S_p, \\ { }^H D_{0^{+}}^{\alpha, \beta ; \psi} E_p(t) & =\lambda_p E_n+\lambda_h S_p+e r I_p-(\phi+\sigma+\mu) E_p, \\ { }^H D_{0^{+}}^{\alpha, \beta ; \psi} I_p(t) & =\omega \lambda_p I_n+\sigma E_p-(e r+\delta+\mu) I_p . \end{aligned} $

with initial conditions:

$ \begin{aligned} & \left(I_{0^{+}}^{1-\xi ; \psi} S_n\right)(0)=17335000, \quad\left(I_{0^{+}}^{1-\xi ; \psi} E_n\right)(0)=1300000, \quad\left(I_{0^{+}}^{1-\xi ; \psi} I_n\right)(0)=80000, \\ & \left(I_{0^{+}}^{1-\xi ; \psi} S_p\right)(0)=40000, \quad\left(I_{0^{+}}^{1-\xi ; \psi} E_p\right)(0)=20000, \quad\left(I_{0^{+}}^{1-\xi ; \psi} I_p\right)(0)=10000 \end{aligned} $

Parameter values employed in the simulations are fixed unless otherwise stated and are provided in Table 2. These values are chosen to reflect a biologically plausible scenario and correspond to findings from contemporary research regarding HSV-2 epidemiology.

Numerical simulations are performed to assess how fractional order $\alpha$ and type $\beta$, in addition to key epidemiological parameters, influence HSV-2 dynamics, whereby shifting $\alpha$ demonstrates how increased memory effects create a system settling delay and shifting transmission/progression rates influence the peak height and duration of infection.

Different $\psi(t)$ correspond to different memory laws-linear, exponential, logarithmic, oscillatory types. These are different orders of persistence related to memory and memory integration which is an important factor in transient behavior and final persistent duration of infection.

Biologically, these figures imply that the larger the memory effect, the longer the length of HSV-2 persistence, which correlates with the chronic nature of this infection and a decreased behavioral response in an increasingly vulnerable population.

Among the key results from the $\Psi$-Hilfer fractional model is the implication of memory driven persistence which implies that even when the basic reproduction number decreases below one, the infection dynamics continue to be influenced by events that have happened before—and from infinitely long ago. Such impulse, via Mittag–Leffler-type decay, cannot be captured by integer-order.

Figure 2.

General Behavior of the Dynamics

For $\psi(t)=t$, the kernel is linear and the temporally established state is constant over time which means $S_n$, $E_n, I_n, S_p, E_p$ and $I_p$ are practically an averaged motion for equilibrium over time. Therefore, it's as if equilibrium can be achieved less quickly the more $\alpha$ decreases (0.9 to 0.7) meaning the less well-defined derivative order is the better remembered derivative of these states within these times. Ultimately, the two infected classes ($I_n, I_p$) keep moving down into equilibrium further as time goes on while the two susceptible classes grow and stabilize which makes sense with the fractional derivative's negative damping and more susceptible to healing/more responsive to it.

Scientific and Interpretive Commentary

This makes sense with the normal Caputo derivative which is appropriate as it's the derivative from which the $\Psi$-Hilfer derivative was derived. This fractional makes sense this way for latency time from infection to visible healing potential - especially since HSV 2 is a recurrent infection. Furthermore, even slightly shifted $\alpha = 1$ promotes dramatic movement in time which suggests that the fractional derivatives of the order have greater meaning in how long over time is best for stabilization of contained infection and recapture of healing rate.

Figure 3.

General Behavior of the Dynamics

Here $\psi(t)=1.3^t$ exponentially proportions what is currently there more and embraces memories of what previously may have existed for transmissive experiences. Therefore, in comparison to this previous entry, stabilization occurs so much more slowly here and the exposed class wavers as does the infected sex-worker class in-between. Additionally, the susceptible populations stabilize over time to levels exponentially equivalent to a delayed stabilization which means more is made good of when there is more of a memory impact of epidemic process.

Scientific and Interpretive Commentary

The exponential kernel applies when intentional exposures add to how intensive transmission currently is. This equals an intentionality of high-risk populations who continue to interact with one another in such a way. Mathematically, this increased memory strength boasts of transience from a value of 1 + t which means that an effective intervention must be in place longer to achieve a bounded equilibrium.

Figure 4.

General Behavior of the Dynamics

For $\psi(t)=0.3 t \ln (5+t)$, nonlinearities dominate the system and create a memory retention speed for saturation over time. Thus, in addition, the populations of the infected classes decrease in slightly regulated fashion and both exposed classes reach equilibrium faster than before. While susceptible classes stabilize much faster than before suggests a less powerful significance from far away memory's decimal time recoveries suggests equilibrium can occur faster than before.

Scientific and Interpretive Commentary

The logarithmic kernel operates as a partial-forgetting form where populations stabilize over time but from predominantly experienced middle interactions. This type of forgetting makes sense with biological realities and real-world observations for HSV-2 over time. The model suggests that types of logarithmic bounds possess solutions where stabilization isn't reality-decayed epidemiologically falsifying.

Figure 5.

General Behavior of the Dynamics

Here $\psi(t)$ is linear and periodic meaning that the contact rate possesses somewhat oscillating characteristics. Thus, small oscillations exist throughout time with the exposed and infected classes pre-stabilization over time. Furthermore, as $\alpha$ becomes decreased, the incrementality of oscillation occurs more and more within this class meaning how fractional memory assists stabilization against transience over time.

Scientific and Interpretive Commentary

This kernel function makes sense for the reality of the sex-worker population if it varies by desire or seasonality (theoretical research suggests sex rates go down when it's cold). Fractional dynamics manage these periodic changes as they're never excessively increased or decreased to destroy positive equity until equilibrium. The $\Psi$-Hilfer derivative allows this mathematical function to make sense that periodic social desirability or environmental desirability are applicable for modeling purposes for transmission purposes.

The relative comparison across Figure 2, Figure 3, Figure 4 and Figure 5 suggests that $\psi(t)$ changes relative change percentages create wildly different relative change percentages in dynamically proportioned transitional evolutions. For $\psi(t)=t$, the average derived linear form comes from sequentially experienced awareness itself; for $\psi(t)=1.3^t$, the average comes from effectively prolonged close memory and experienced awareness for how long HSV2 effectively tries to maintain longevity of epidemiological process - equilibrium successes try to maintain where previous efforts could maintain them longer than control could reach; conversely, $\psi(t)=0.3 t \ln (5+t)$ bounds efforts early to stabilize earlier while $\psi(t)=6 t+4 t \sin (0.02 t)$ reduces hypothetically realistic oscillations that might be most reasonably present (periodic damping). Ultimately these findings support how stabilization can occur under one form of $\Psi$-Hilfer but it's clear that kernel variances and differentiation order of precision impact depth of memory retention optimization effectiveness and relative speed stabilization success of previously worked efforts in time between differentials at best.

The scientific merit of this project is that the Adams-Bashforth-Moulton predictor-corrector scheme employed in this paper demonstrates improved numerical stability and accuracy for $\Psi$-Hilfer fractional systems compared to other more commonly applied numerical approaches and is relevant to memory-based kernels in which such hereditary systems are critical to epidemiology.

Ultimately, the simulations support theoretical existence, uniqueness and positive bounded solutions for developed stability while assessing reality surrounding epidemiological realities surrounding fractional memory for effective disease observations across temporal bounds. Fractional derivatives represent extensive latent periods/cycles and recurrences that define HSV 2's definition; furthermore, how quickly stabilization motivated by $\psi(t)$ function and $\alpha$ order also suggests how quickly equilibrium can be reached with beneficial stabilization potential accounting for oscillation potential. From a prevention perspective, these findings support dynamic parity for educated prevention/treatment strategies in alignment with behavioral insights best aligned against memory inspired transmission responses - and $\Psi$-Hilfer fractional model absorbs all statistical expectations to connect newfound revelations between analytically relevant epidemiological solutions made over time as transmission potential defines itself with increased confidence over time rather than simpler standalone models without research/work behind them.

The $\Psi$-Hilfer fractional system is a biologically relevant and mathematically tunable fractional order model for HSV-2 spread. Fractional elements of memory, a heterogeneous population and optimal control enhance the interpretability and practicality of the model for any potential interventions.

These results indicate that the $\Psi$-Hilfer fractional model is more realistic and flexible to HSV-2 transmission dynamics. These relative qualitative and numerical observations note that memory contributes significantly to stability, infection control and asymptotic behavior which opens new doors to understanding and application for intervention.