Lijie You^*, Sijie Yu
School of Mathematical Sciences, China West Normal University, Nanchong, Sichuan, China.
DOI: 10.4236/jamp.2026.145096   PDF    HTML   XML   27 Downloads   123 Views   Citations

Abstract

In recent years, controlling the spread of rumors has emerged as an issue of significant public concern. Mathematical modeling provides a powerful tool for conducting quantitative analysis of this problem and can offer critical support for decision-making. This paper addresses the control of rumor propagation in online social networks by introducing a novel fractional-order rumor spreading model. The model categorizes the population into seven compartments and simultaneously incorporates three types of intervention measures-educational mechanisms, memory and forgetting mechanisms, and rumor-refutation mechanisms to better reflect real-world rumor mitigation scenarios. We first establish the existence, non-negativity, uniqueness, and boundedness of the solutions to the proposed model. Similar to epidemiological modeling, a basic reproduction number for rumor spread is defined, and the existence and stability of the model’s equilibrium points are analyzed. On this basis, an optimal control problem is formulated with the aim of balancing intervention effectiveness and implementation cost. By applying Pontryagin’s Maximum Principle, we derive the necessary conditions for optimal control and obtain the corresponding optimal strategy. Numerical simulations demonstrate that the synergistic use of the three intervention mechanisms significantly reduces the peak prevalence and final impact of rumors, outperforming any single intervention, which confirms the validity and practical utility of our model.

Keywords

Caputo Fractional Derivative, Basic Reproduction Number, Local Stability, Global Stability, Optimal Control

Share and Cite:

1. Introduction

The rapid development of the internet and social media has dramatically accelerated the speed and expanded the scale of information dissemination. Platforms such as Facebook, Telegram, Instagram, and TikTok have become common channels for spreading misinformation due to their openness and immediacy [1] [2]. While these platforms greatly facilitate the rapid sharing of information, they also create favorable conditions for the propagation of rumors [3]. As a form of false information lacking a factual basis, rumors tend to spread quickly during public emergencies, posing serious threats to social stability, public safety, and economic development [4]. For instance, a widely circulated claim that “Banlangen and Shuanghuanglian oral liquid can effectively prevent viral infection” triggered panic buying in multiple regions, severely disrupting normal epidemic control efforts [5] [6]. Similarly, rumors of radioactive contamination following the Fukushima nuclear accident led to panic-driven salt hoarding and market instability. Such cases illustrate that the harm caused by rumors extends beyond misleading public perception. It can also exert far-reaching impacts on socioeconomic operations. Therefore, investigating the mechanisms of rumor propagation is of critical practical importance. It is essential to adopt appropriate control strategies to curb the spread of rumors and minimize their negative consequences and associated losses.

The modeling of rumor propagation shares significant similarities with epidemiological models in mathematical biology, as both fundamentally describe the process of “state transition” caused by “contact” with “spreaders” within a population [7]. The seminal work on rumor propagation modeling began in the 1960s with Daley and Kendall [8], who proposed the DK (Daley-Kendall) model. This model classifies the population into three compartments: ignorants, spreaders, and stiflers. The model is characterized by a stifling mechanism that operates on state-transition rules distinct from those used in epidemiological models, thus establishing a foundation for subsequent research on rumor propagation. A key refinement to the DK model was introduced by Maki and Thompson [9], who made a fundamental improvement to its stifling mechanism, leading to the model now known as the MT model. Building on the DK and MT models [10] [11], subsequent work has vastly extended these mathematical frameworks [12]-[14]. Recently, Huo et al. [15] developed a novel multi-medium rumor propagation model. This model incorporates diverse dissemination channels, such as social platforms and news websites, and accounts for the differentiated behaviors of spreaders. Their study demonstrated that the cross-media movement of ignorant individuals plays a crucial role in rumor propagation. Moreover, individual spreading behavior, being influenced by media preferences, means that media coverage significantly shapes this process. This insight indicates that effective intervention and control of rumors can be achieved by implementing measures to regulate the impact of media coverage. In addition, Zhao et al. [16] developed a rumor propagation model that incorporates a forgetting mechanism. Their key contribution was the introduction of a constraint relationship between the forgetting rate and the stifling rate, which critically governs the behavior of stiflers and ultimately determines the rumor’s final impact. This formulation revealed new dynamical features of rumor propagation. Numerical simulations further confirmed that the forgetting mechanism can significantly reduce the influence of rumors.

The models established by the aforementioned researchers are predominantly integer-order, characterized by their simplicity and computational efficiency. However, such models oversimplify rumor propagation as a uniform, smooth transition, thereby failing to capture the intrinsic memory effects. Consequently, these model’s fail to accurately capture real-world rumor dynamics, which often exhibit explosive propagation in the initial phase and prolonged attenuation in the later stage. Although some studies have attempted to incorporate memory effects within integer-order differential frameworks, the intrinsic constraints of integer-order calculus continue to impede a faithful representation of such phenomena. Since the 1990s, fractional calculus has garnered significantly increasing interest from the academic and scientific community, evolving from a primarily mathematical curiosity into a vital tool for modeling complex phenomena across diverse fields. The memory and hereditary properties of fractional-order calculus are well-established and have been extensively applied in numerous studies [17]-[19]. For instance, Li et al. [20] extended the traditional Susceptible-Infected-Removed (SIR) model by introducing a new compartment, clarifiers (C), representing individuals who learn the truth and actively debunk rumors. The inclusion of a clarification mechanism increases the model’s fidelity in simulating real-world rumor refutation processes on social media platforms. Moreover, the application of fractional calculus substantially enhances modeling performance. Numerical simulations confirm that the fractional-order framework yields a closer fit to empirical propagation curves, reduces fitting errors, and achieves more accurate predictions of rumor spread dynamics. Very recently, Niu et al. [21] introduced a “doubters (D)” compartment to more precisely characterize the intermediate state of hesitation and verification that individuals undergo from encountering a rumor to deciding whether to spread it. They successfully integrated fractional-order calculus with time delays to construct a rumor propagation model that better reflects real-world complexity. Furthermore, they conducted an in-depth analysis of delay-induced Hopf bifurcation, revealing the intrinsic mechanisms that may lead to periodic oscillations in rumor propagation dynamics. Among the various definitions of fractional derivatives, the Grünwald-Letnikov, Riemann-Liouville, and Caputo operators are the most widely used. In practical modeling, the Caputo derivative is often preferred primarily owing to two key advantages: its compatibility with standard initial conditions and the fact that the derivative of a constant is zero, which align more naturally with physical interpretations and integer-order calculus conventions. Based on these advantages, the Caputo fractional derivative is employed for the subsequent model analysis.

Based on a comprehensive review of existing literature, it is evident that while substantial efforts have been devoted to exploring rumor suppression mechanisms, the majority of studies have focused on examining the isolated effects of individual strategies such as educational campaigns, debunking initiatives, natural forgetting processes, or incentive punishment systems. A critical limitation of these studies lies in their oversight of the complex reality where multiple suppression mechanisms often operate concurrently or sequentially, potentially giving rise to significant synergistic or antagonistic effects. For example, the efficacy of debunking efforts is often contingent on individuals’ prior educational background, whereas the effect of incentive-punishment measures may be undermined by natural forgetting processes. Thus, the failure to account for such interdependencies inherently limits the explanatory and predictive power of models that simply combine individual strategies for real-world rumor intervention. To address this research gap, this paper proposes a novel rumor propagation model termed the Uneducated-Educated-Exposed-Spreader-Hibernated-Debunker-Stifler ($S_{\mu }{S}_{e}ERHDU$ ) model to better capture the dissemination dynamics in Online Social Networks (OSNs). Based on this modeling framework, we design three coordinated intervention strategies-preemptive education, memory-and-forgetting feedback, and active debunking-to more realistically capture multi-dimensional rumor propagation dynamics in real social environments. Furthermore, an optimal control framework is established to minimize both social impact and control costs during rumor containment.

The paper proceeds as follows. Section 2 presents the foundational properties of the Caputo fractional derivative operator and the Mittag-Leffler function, which serve as essential mathematical preliminaries for the subsequent analysis. Section 3 introduces a novel fractional-order rumor propagation model and provides a rigorous theoretical analysis regarding the existence, non-negativity, boundedness, and uniqueness of its solutions. Section 4 performs a stability analysis of both the rumor-free and rumor-spreading equilibria. In Section 5, an optimal control problem is formulated and analytically examined based on the proposed model. Section 6 carries out numerical simulations to validate the theoretical findings. Finally, Section 7 concludes the paper by summarizing the principal findings and discussing their implications.

2. Preliminaries

This section outlines the essential concepts of Caputo fractional calculus that underpin the mathematical model developed in the subsequent sections of this work.

Definition 1 [22] Let $g\left(t\right)\in L^1\left(\left[t_0,\tau \right]\right)$ ($L^1$ is the set of Lebesgue integrable functions), then left and right RL integrals of fractional-order $\alpha$ of the function $g\left(t\right)$ are, respectively, given as follows

${}_{t_0}^{RL}{I}_t^{\alpha }g\left(t\right)=\frac{1}{\Gamma \left(\alpha \right)}{\displaystyle {\int }_{t_0}^t{\left(t-x\right)}^{\alpha -1}g\left(x\right)\text{d}x},$

and

${}_t^{RL}{I}_{\tau }^{\alpha }g\left(t\right)=\frac{1}{\Gamma \left(\alpha \right)}{\displaystyle {\int }_t^{\tau }{\left(x-t\right)}^{\alpha -1}g\left(x\right)\text{d}x},$

where $\Gamma (.)$ is the well-known gamma function.

Definition 2 [22] The left and right RL fractional-order derivatives of the function $g\left(t\right)$ are, respectively, defined as follows

${}_{t_0}^{RL}{D}_t^{\alpha }g\left(t\right)=\frac{1}{\Gamma \left(n-\alpha \right)}{\left(\frac{\text{d}}{\text{d}t}\right)}^n{\displaystyle {\int }_{t_0}^t{\left(t-x\right)}^{n-\alpha -1}g\left(x\right)\text{d}x},$

and

${}_t^{RL}{D}_{\tau }^{\alpha }g\left(t\right)=\frac{1}{\Gamma \left(n-\alpha \right)}{\left(\frac{\text{d}}{\text{d}t}\right)}^n{\displaystyle {\int }_t^{\tau }{\left(x-t\right)}^{n-\alpha -1}g\left(x\right)\text{d}x},$

where $n$ is any positive integer such that $n-1<\alpha \le n$ .

Definition 3 [22] Let $g\left(t\right)\in C^n\left(t_0,\tau \right)$ (i.e. $g^{\left(n-1\right)}\left(t\right)$ is absolutely continuous), then left and right Caputo derivatives of fractional-order $\alpha$ of the function $g\left(t\right)$ are, respectively, defined as follows

${}_{t_0}^C{D}_t^{\alpha }g\left(t\right)=\frac{1}{\Gamma \left(n-\alpha \right)}{\displaystyle {\int }_{t_0}^t{\left(t-x\right)}^{n-\alpha -1}{g}^{\left(n\right)}\left(x\right)\text{d}x},$

and

${}_t^C{D}_{\tau }^{\alpha }g\left(t\right)=\frac{{\left(-1\right)}^n}{\Gamma \left(n-\alpha \right)}{\displaystyle {\int }_t^{\tau }{\left(x-t\right)}^{n-\alpha -1}{g}^{\left(n\right)}\left(x\right)\text{d}x},$

where $n$ is any positive integer such that $n-1<\alpha \le n$ .

Lemma 1 [23] Let $f\left(x\right)\in C\left[a,b\right]$ and ${}_a^c{D}_x^{\alpha }f\left(x\right)\in C\left(a,b\right]$ for $0<\alpha \le 1$ , then

$f\left(x\right)=f\left(a\right)+\frac{1}{\Gamma \left(\alpha \right)}\left({}_a^c{D}_x^{\alpha }g\right)\left(\xi \right){\left(x-a\right)}^{\alpha }$

with $a\le \xi \le x$ , $\forall x\in \left(a,b\right]$ .

Using Lemma 1, we state the following corollary.

Corollary 1 [23] Suppose that $g\left(x\right)\in C\left[0,b\right]$ and ${}_0^c{D}_x^{\alpha }g\left(x\right)\in C\left(0,b\right]$ for $0<\alpha \le 1$ . If ${}_0^c{D}_x^{\alpha }g\left(x\right)\ge 0$ $\forall t\in \left(0,b\right)$ , then the function $g$ is non-decreasing and if ${}_0^c{D}_x^{\alpha }g\left(x\right)\le 0$ $\forall x\in \left(0,b\right)$ , then the function $g$ is non-increasing for all $x\in \left(0,b\right)$ .

Lemma 2 [24] Let $\alpha \in \left(0,1\right)$ and consider a continuous function $x:\left[t_0,\infty \right)\to ℝ$ satisfying the following condition

${}_0^c{D}_t^{\alpha }x\left(t\right)+\mu x\left(t\right)\le \nu ,t\ge t_0,\mu ,\nu \in ℝ,\mu \ne 0.$

Then, we have the inequality

$x\left(t\right)\le \left(x\left(t_0\right)-\frac{\nu }{\mu }\right)E_{\alpha }\left(-\mu {\left(t-t_0\right)}^{\alpha }\right)+\frac{\nu }{\mu },$

for all $t\ge t_0$ , where $E_{\alpha }$ is the Mittag-Leffler function of one parameter defined by

$E_{\alpha }\left(t\right)={\displaystyle \sum }_{k=0}^{\infty }\frac{t^k}{\Gamma \left(\alpha k+1\right)}.$

Definition 4 [25] Let $F\left(s\right)$ is the Laplace transform of the $f\left(t\right)$ . Then,

$ℒ\left\{{}_0^c{D}_t^{\alpha }f\left(t\right),s\right\}=s^{\alpha }F\left(s\right)-{\displaystyle \sum }_{i=0}^{n-1}{s}^{\alpha -i-1}{f}^{\left(i\right)}\left(0\right),\alpha \in \left(n-1,n\right];n\in \mathbb{N}.$

Definition 5 [26] The Mittag-Leffler function $E_{l,m}\left(x\right)$ is given by

$E_{l,m}\left(x\right)={\displaystyle \sum }_{n=0}^{\infty }\frac{x^n}{\Gamma \left(ln+m\right)},x\in ℝ,l>0,m>0,$

and satisfies the property

$E_{l,m}\left(x\right)=x{E}_{l,l+m}\left(x\right)+\frac{1}{\Gamma \left(m\right)},$

and the Laplace transform of $t^{m-1}{E}_{l,m}\left(\pm \lambda t^l\right)$ is given by

$ℒ\left[t^{m-1}{E}_{l,m}\left(\pm \lambda t^l\right)\right]=\frac{s^{l-m}}{s^l\mp \lambda }.$

Lemma 3 [27] Let $0<\alpha <1$ and $g\in C\left[0,T\right]$ be a positive valued function. Then, for all $t\in \left[0,T\right)$ , one has

${}_0^c{D}_t^{\alpha }\left(g\left(t\right)-g^{*}-g^{*}\ln \frac{g\left(t\right)}{g^{*}}\right)\le \left(1-\frac{g^{*}}{g\left(t\right)}\right){}_0^c{D}_t^{\alpha }g\left(t\right),$

for all $g^{\text{*}}\in {ℝ}_{+}$ .

3. Model Derivation

The $S_{\mu }{S}_{e}ERHDU$ Rumor Propagation Model

We developed a fractional-order $S_{\mu }{S}_{e}ERHDU$ model to characterize the dissemination dynamics of rumors in online social networks (OSNs). The total population, denoted by $N$ , was categorized into seven distinct compartments.

1) $S_{\mu }$ (uneducated): people who do not have the ability to identify rumors;

2) $S_e$ (educated): people who have the ability to identify rumors;

3) $E$ (exposed): people who are in contact with rumor spreaders and debunkers;

4) $R$ (spreaders): people who understand the rumor information and begin to spread the rumor;

5) $H$ (hibernated): people who temporarily forget about rumors;

6) $D$ (debunkers): people who know the truth and actively spread corrective information to counter the rumor;

7) $U$ (stiflers): people who do not spread or debunk rumor information.

At any time $t$ , the densities (or proportions) of these categories in the total population are denoted as $S_{\mu }\left(t\right),S_e\left(t\right),E\left(t\right),R\left(t\right),H\left(t\right),D\left(t\right)$ and $U\left(t\right)$ , respectively. Note that the total population is normalized to unity, i.e., $N\left(t\right)=S_{\mu }\left(t\right)+S_e\left(t\right)+E\left(t\right)+R\left(t\right)+H\left(t\right)+D\left(t\right)+U\left(t\right)\equiv 1$ .

We let the total population $N\left(t\right)=S_{\mu }\left(t\right)+S_e\left(t\right)+E\left(t\right)+R\left(t\right)+H\left(t\right)+D\left(t\right)+U\left(t\right)$ . The principles governing rumor propagation can be summarized as follows.

(A-1) The parameter ${\Lambda }^{\alpha }$ , representing the number of new entrants in social communication, enters the compartments $S_{\mu }$ and $S_e$ in distinct proportions: $\left(1-\rho \right)$ and $\rho$ , respectively.

(A-2) The transfer rate ${\eta }^{\alpha }$ from compartment $S_{\mu }$ to $S_e$ quantifies the efficacy of educational interventions in enhancing the ability to identify rumor-based information.

(A-3) The infection rate ${\beta }^{\alpha }$ represents the rate at which susceptible individuals including both $S_{\mu }$ and $S_e$ are influenced by rumor spreaders and debunkers. The parameter $m$ acts as an adjustment coefficient modulating the infection rate associated with rumor debunkers ($D$ ).

(A-4) The transmission rate from the educated susceptible compartment ($S_e$ ) to the exposed compartment ($E$ ) is reduced by a factor of $\left(1-n\right)$ , where $n$ represents the efficacy of education in lowering susceptibility.

(A-5) The exposed individuals ($E$ ) transition to the rumor spreaders ($R$ ) at a rate ${\xi }^{\alpha }$ , and to the rumor debunkers ($D$ ) at a rate ${\sigma }^{\alpha }$ .

(A-6) Rumor spreaders ($R$ ) transition to the rumor debunker compartment ($D$ ) at a rate of ${\delta }^{\alpha }$ by refuting rumors in a timely manner.

(A-7) The forgetting mechanism allows a portion of rumor spreaders ($R$ ) to transition into the dormant compartment ($H$ ) at a rate ${\theta }^{\alpha }$ , representing the loss of active engagement with the rumor. Conversely, the memory mechanism can reactivate individuals in the dormant state ($H$ ), causing them to re-enter the rumor spreader compartment ($R$ ) at a rate $g^{\alpha }$ , reflecting the retrieval or renewed influence of previously encountered information.

(A-8) As times go on, individuals tend to lose interest in both spreading and debunking rumors. As a result, participants gradually transition into the compartment $U$ , which represents those who neither spread nor counteract rumors. This transition is characterized by the following rates: individuals in the rumor spreader compartment ($R$ ) enter $U$ at a rate ${\omega }^{\alpha }$ , those in the hibernator compartment ($H$ ) transition at a rate ${\gamma }^{\alpha }$ , and individuals in the debunker compartment ($D$ ) move to $U$ at a rate ${\mu }^{\alpha }$ .

Based on the transmission mechanisms described above, the considered fractional-order $S_{\mu }{S}_{e}ERHDU$ rumor propagation model is formally defined by the following system of equations:

$\{\begin{array}{l}{}_0^c{D}_t^{\alpha }{S}_{\mu }\left(t\right)=\left(1-\rho \right){\Lambda }^{\alpha }-{\beta }^{\alpha }\widehat{\kappa }{S}_{\mu }\left(t\right)\left(R\left(t\right)+mD\left(t\right)\right)-\left({\eta }^{\alpha }+d^{\alpha }\right)S_{\mu }\left(t\right),\\ {}_0^c{D}_t^{\alpha }{S}_e\left(t\right)=\rho {\Lambda }^{\alpha }+{\eta }^{\alpha }{S}_{\mu }\left(t\right)-\left(1-n\right){\beta }^{\alpha }\widehat{\kappa }{S}_e\left(t\right)\left(R\left(t\right)+mD\left(t\right)\right)\\ -\left({\lambda }^{\alpha }+d^{\alpha }\right)S_e\left(t\right),\\ {}_0^c{D}_t^{\alpha }E\left(t\right)={\beta }^{\alpha }\widehat{\kappa }{S}_{\mu }\left(t\right)\left(R\left(t\right)+mD\left(t\right)\right)+\left(1-n\right){\beta }^{\alpha }\widehat{\kappa }{S}_e\left(t\right)\left(R\left(t\right)+mD\left(t\right)\right)\\ -\left({\xi }^{\alpha }+{\sigma }^{\alpha }+{\epsilon }^{\alpha }+d^{\alpha }\right)E\left(t\right),\\ {}_0^c{D}_t^{\alpha }R\left(t\right)={\xi }^{\alpha }E\left(t\right)+g^{\alpha }H\left(t\right)-{\theta }^{\alpha }R\left(t\right)-{\delta }^{\alpha }R\left(t\right)-\left({\omega }^{\alpha }+d^{\alpha }\right)R\left(t\right),\\ {}_0^c{D}_t^{\alpha }H\left(t\right)={\theta }^{\alpha }R\left(t\right)-g^{\alpha }H\left(t\right)-\left({\gamma }^{\alpha }+d^{\alpha }\right)H\left(t\right),\\ {}_0^c{D}_t^{\alpha }D\left(t\right)={\sigma }^{\alpha }E\left(t\right)+{\delta }^{\alpha }R\left(t\right)-\left({\mu }^{\alpha }+d^{\alpha }\right)D\left(t\right),\\ {}_0^c{D}_t^{\alpha }U\left(t\right)={\lambda }^{\alpha }{S}_e\left(t\right)+{\epsilon }^{\alpha }E\left(t\right)+{\omega }^{\alpha }R\left(t\right)+{\gamma }^{\alpha }H\left(t\right)+{\mu }^{\alpha }D\left(t\right)-d^{\alpha }U\left(t\right),\end{array}$(1)

with the following non-negative initial conditions:

$\begin{array}{l}{S}_{\mu }\left(0\right)=S_{\mu 0},S_e\left(0\right)=S_{e0},E\left(0\right)=E_0,R\left(0\right)=R_0,\\ H\left(0\right)=H_0,D\left(0\right)=D_0,U\left(0\right)=U_0,\end{array}$ (2)

where ${}_0^c{D}_t^{\alpha }$ is the Caputo fractional-order derivative with $0<\alpha \le 1$ . The schematic diagram of rumor propagation process is presented in Figure 1.

Figure 1. Schematic diagram of the $S_{\mu }{S}_{e}ERHDU$ rumor spreading model.

The parameter description in the rumor propagation model (1) is shown in Table 1. The units of all parameters in Table 1 are hour⁻¹.

Table 1. Parameters meaning of model (1).

Parameters	Parameters meaning
$\Lambda$	The entry rate of individuals into the social communication system
$\rho$	The proportion of new entrants joining the $S_e$ compartment
$\beta$	The probability of successful information transmission per effective contact
$\eta$	The transition rate from $S_{\mu }$ to $S_e$ , induced by exposure to a rumor event
$\widehat{\kappa }$	Contact efficiency of information disseminators
$m$	The relative transmission efficacy of debunkers $\left(D\right)$ compared to rumor spreaders $\left(R\right)$
$d$	Natural loss rate of each compartment
$n$	The transition attenuation rate from $S_e$ to $E$
$\lambda$	The transition rate from $S_e$ to $U$
$\xi$	The transition rate from $E$ to $R$
$\sigma$	The transition rate from $E$ to $D$
$\epsilon$	The transition rate from $E$ to $U$
$g$	The transition rate from $H$ to $R$
$\theta$	The transition rate from $R$ to $H$
$\delta$	The transition rate from $R$ to $D$
$\omega$	The transition rate from $R$ to $U$
$\gamma$	The transition rate from $H$ to $U$
$\mu$	The transition rate from $D$ to $U$

4. Qualitative Analysis

4.1. Existence and Uniqueness

In this subsection, we investigate the existence and uniqueness of the solutions to system (1) within the domain:

$\Omega =\left\{\left(S_{\mu },S_e,E,R,H,D,U\right)\in {ℝ}^7:\max \left\{\left|S_{\mu }|,|S_e\right|,\left|E\right|,\left|R\right|,\left|H\right|,\left|D\right|,\left|U\right|\right\}\le M\right\}.$ (3)

Theorem 1 For each initial condition

$X\left(0\right)=\left(S_{\mu }\left(0\right),S_e\left(0\right),E\left(0\right),R\left(0\right),H\left(0\right),D\left(0\right),U\left(0\right)\right)\in \Omega$

in the system (1), there always exists a unique solution

$X=\left(S_{\mu },S_e,E,R,H,D,U\right)\in \Omega$

for all $t\ge 0$ .

Proof. To establish the existence and uniqueness of solutions, we follow the methodology employed in [28]. Let

$X\left(t\right)=\left(S_{\mu }\left(t\right),S_e\left(t\right),E\left(t\right),R\left(t\right),H\left(t\right),D\left(t\right),U\left(t\right)\right),$

$\widehat{X}\left(t\right)=\left({\widehat{S}}_{\mu }\left(t\right),{\widehat{S}}_e\left(t\right),\widehat{E}\left(t\right),\widehat{R}\left(t\right),\widehat{H}\left(t\right),\widehat{D}\left(t\right),\widehat{U}\left(t\right)\right).$

For brevity, let $X\left(t\right)=X$ and $\widehat{X}\left(t\right)=\widehat{X}$ . Consider $F\left(X\right)=\left(F_1\left(X\right),\cdots ,F_7\left(X\right)\right)$ , with

$\begin{array}{l}{F}_1\left(X\right)=\left(1-\rho \right){\Lambda }^{\alpha }-{\beta }^{\alpha }\widehat{\kappa }{S}_{\mu }\left(R+mD\right)-\left({\eta }^{\alpha }+d^{\alpha }\right)S_{\mu },\\ F_2\left(X\right)=\rho {\Lambda }^{\alpha }+{\eta }^{\alpha }{S}_{\mu }-\left(1-n\right){\beta }^{\alpha }\widehat{\kappa }{S}_e\left(R+mD\right)-\left({\lambda }^{\alpha }+d^{\alpha }\right)S_e,\\ F_3\left(X\right)={\beta }^{\alpha }\widehat{\kappa }{S}_{\mu }\left(R+mD\right)+\left(1-n\right){\beta }^{\alpha }\widehat{\kappa }{S}_e\left(R+mD\right)-\left({\xi }^{\alpha }+{\sigma }^{\alpha }+{\epsilon }^{\alpha }+d^{\alpha }\right)E,\\ F_4\left(X\right)={\xi }^{\alpha }E+g^{\alpha }H-\left({\theta }^{\alpha }+{\delta }^{\alpha }+{\omega }^{\alpha }+d^{\alpha }\right)R,\\ F_5\left(X\right)={\theta }^{\alpha }R-\left(g^{\alpha }+{\gamma }^{\alpha }+d^{\alpha }\right)H,\\ F_6\left(X\right)={\sigma }^{\alpha }E+{\delta }^{\alpha }R-\left({\mu }^{\alpha }+d^{\alpha }\right)D,\\ F_7\left(X\right)={\lambda }^{\alpha }{S}_e+{\epsilon }^{\alpha }E+{\omega }^{\alpha }R+{\gamma }^{\alpha }H+{\mu }^{\alpha }D-d^{\alpha }U.\end{array}$ (4)

Now for any $X,\widehat{X}\in \Omega$ , we have:

$\begin{array}{l}\left|F\left(X\right)-F\left(\widehat{X}\right)\right|={\displaystyle \sum }_{i=1}^7\left|F_i\left(X\right)-F_i\left(\widehat{X}\right)\right|\\ =\left|-{\beta }^{\alpha }\widehat{\kappa }\left(S_{\mu }R-{\widehat{S}}_{\mu }\widehat{R}\right)-m{\beta }^{\alpha }\widehat{\kappa }\left(S_{\mu }D-{\widehat{S}}_{\mu }\widehat{D}\right)-\left({\eta }^{\alpha }+d^{\alpha }\right)\left(S_{\mu }-{\widehat{S}}_{\mu }\right)\right|\\ +|{\eta }^{\alpha }\left(S_{\mu }-{\widehat{S}}_{\mu }\right)-\left(1-n\right){\beta }^{\alpha }\widehat{\kappa }\left(S_{e}R-{\widehat{S}}_e\widehat{R}\right)\\ -m\left(1-n\right){\beta }^{\alpha }\widehat{\kappa }\left(S_{e}D-{\widehat{S}}_e\widehat{D}\right)-\left({\lambda }^{\alpha }+d^{\alpha }\right)\left(S_e-{\widehat{S}}_e\right)|\\ +|{\beta }^{\alpha }\widehat{\kappa }\left(S_{\mu }R-{\widehat{S}}_{\mu }\widehat{R}\right)+m{\beta }^{\alpha }\widehat{\kappa }\left(S_{\mu }D-{\widehat{S}}_{\mu }\widehat{D}\right)+\left(1-n\right){\beta }^{\alpha }\widehat{\kappa }\left(S_{e}R-{\widehat{S}}_e\widehat{R}\right)\\ +m\left(1-n\right){\beta }^{\alpha }\widehat{\kappa }\left(S_{e}D-{\widehat{S}}_e\widehat{D}\right)-\left({\xi }^{\alpha }+{\sigma }^{\alpha }+{\epsilon }^{\alpha }+d^{\alpha }\right)\left(E-\widehat{E}\right)|\\ +\left|{\xi }^{\alpha }\left(E-\widehat{E}\right)+g^{\alpha }\left(H-\widehat{H}\right)-\left({\theta }^{\alpha }+{\delta }^{\alpha }+{\omega }^{\alpha }+d^{\alpha }\right)\left(R-\widehat{R}\right)\right|\end{array}$

$\begin{array}{l}+\left|{\theta }^{\alpha }\left(R-\widehat{R}\right)-\left(g^{\alpha }+{\gamma }^{\alpha }+d^{\alpha }\right)\left(H-\widehat{H}\right)\right|\\ +\left|{\sigma }^{\alpha }\left(E-\widehat{E}\right)+{\delta }^{\alpha }\left(R-\widehat{R}\right)-\left({\mu }^{\alpha }+d^{\alpha }\right)\left(D-\widehat{D}\right)\right|\\ +|{\lambda }^{\alpha }\left(S_e-{\widehat{S}}_e\right)+{\epsilon }^{\alpha }\left(E-\widehat{E}\right)+{\omega }^{\alpha }\left(R-\widehat{R}\right)\\ +{\gamma }^{\alpha }\left(H-\widehat{H}\right)+{\mu }^{\alpha }\left(D-\widehat{D}\right)-d^{\alpha }\left(U-\widehat{U}\right)|\\ \le \left(2{\eta }^{\alpha }+d^{\alpha }\right)\left|S_{\mu }-{\widehat{S}}_{\mu }\right|+2{\beta }^{\alpha }\widehat{\kappa }\left|S_{\mu }R-{\widehat{S}}_{\mu }\widehat{R}\right|+2m{\beta }^{\alpha }\widehat{\kappa }\left|S_{\mu }D-{\widehat{S}}_{\mu }\widehat{D}\right|\\ +\left(2{\lambda }^{\alpha }+d^{\alpha }\right)\left|S_e-{\widehat{S}}_e\right|+2\left(1-n\right){\beta }^{\alpha }\widehat{\kappa }\left|S_{e}R-{\widehat{S}}_e\widehat{R}\right|\\ +2m\left(1-n\right){\beta }^{\alpha }\widehat{\kappa }\left|S_{e}D-{\widehat{S}}_e\widehat{D}\right|+\left(2{\xi }^{\alpha }+2{\sigma }^{\alpha }+2{\epsilon }^{\alpha }+d^{\alpha }\right)\left|E-\widehat{E}\right|\\ +\left(2{\mu }^{\alpha }+d^{\alpha }\right)\left|D-\widehat{D}\right|+\left(2{\theta }^{\alpha }+2{\delta }^{\alpha }+2{\omega }^{\alpha }+d^{\alpha }\right)\left|R-\widehat{R}\right|\\ +\left(2g^{\alpha }+2{\gamma }^{\alpha }+d^{\alpha }\right)\left|H-\widehat{H}\right|+d^{\alpha }\left|U-\widehat{U}\right|\\ \le l_1\left|S_{\mu }-{\widehat{S}}_{\mu }\right|+l_2\left|S_e-{\widehat{S}}_e\right|+l_3\left|E-\widehat{E}\right|+l_4\left|R-\widehat{R}\right|\\ +l_5\left|H-\widehat{H}\right|+l_6\left|D-\widehat{D}\right|+l_7\left|U-\widehat{U}\right|\\ \le L\Vert X-\widehat{X}\Vert ,\end{array}$(5)

where

$\begin{array}{l}L=\text{max}\left\{l_1,l_2,l_3,l_4,l_5,l_6,l_7\right\},\\ l_1=2{\eta }^{\alpha }+d^{\alpha }+2{\beta }^{\alpha }\widehat{\kappa }M+2m{\beta }^{\alpha }\widehat{\kappa }M,\\ l_2=2{\lambda }^{\alpha }+d^{\alpha }+2\left(1-n\right){\beta }^{\alpha }\widehat{\kappa }M+2m\left(1-n\right){\beta }^{\alpha }\widehat{\kappa }M,\\ l_3=2{\xi }^{\alpha }+2{\sigma }^{\alpha }+2{\epsilon }^{\alpha }+d^{\alpha },\\ l_4=2{\theta }^{\alpha }+2{\delta }^{\alpha }+2{\omega }^{\alpha }+d^{\alpha }+2{\beta }^{\alpha }\widehat{\kappa }M+2\left(1-n\right){\beta }^{\alpha }\widehat{\kappa }M,\\ l_5=2g^{\alpha }+2{\gamma }^{\alpha }+d^{\alpha },\\ l_6=2{\mu }^{\alpha }+d^{\alpha }+2m{\beta }^{\alpha }\widehat{\kappa }M+2m\left(1-n\right){\beta }^{\alpha }\widehat{\kappa }M,\\ l_7=d^{\alpha }.\end{array}$(6)

Hence, $F\left(X\right)$ satisfies the Lipschitz condition. Therefore, the fractional-order rumor model (1) always possesses a unique solution.

4.2. Non-Negativity and Uniform Boundedness

Theorem 2 The solution of the fractional-order system (1)-(2) is nonnegative.

Proof. From the model (1), we reach

$\begin{array}{l}{{}_0^c{D}_t^{\alpha }{S}_{\mu }\left(t\right)|}_{S_{\mu }=0}=\left(1-\rho \right){\Lambda }^{\alpha }\ge 0,\\ {{}_0^c{D}_t^{\alpha }{S}_e\left(t\right)|}_{S_e=0}=\rho {\Lambda }^{\alpha }+{\eta }^{\alpha }{S}_{\mu }\left(t\right)\ge 0,\\ {{}_0^c{D}_t^{\alpha }E\left(t\right)|}_{E=0}={\beta }^{\alpha }\widehat{\kappa }{S}_{\mu }\left(t\right)\left(R\left(t\right)+mD\left(t\right)\right)+\left(1-n\right){\beta }^{\alpha }\widehat{\kappa }{S}_e\left(t\right)\left(R\left(t\right)+mD\left(t\right)\right)\ge 0,\\ {{}_0^c{D}_t^{\alpha }R\left(t\right)|}_{R=0}={\xi }^{\alpha }E\left(t\right)+g^{\alpha }H\left(t\right)\ge 0,\\ {{}_0^c{D}_t^{\alpha }H\left(t\right)|}_{H=0}={\theta }^{\alpha }R\left(t\right)\ge 0,\\ {{}_0^c{D}_t^{\alpha }D\left(t\right)|}_{D=0}={\sigma }^{\alpha }E\left(t\right)+{\delta }^{\alpha }R\left(t\right)\ge 0,\\ {{}_0^c{D}_t^{\alpha }U\left(t\right)|}_{U=0}={\lambda }^{\alpha }{S}_e\left(t\right)+{\epsilon }^{\alpha }E\left(t\right)+{\omega }^{\alpha }R\left(t\right)+{\gamma }^{\alpha }H\left(t\right)+{\mu }^{\alpha }D\left(t\right)\ge 0.\end{array}$ (7)

From Corollary 1, we can deduce that the feasible solutions of system (1) with the initial condition (2) are non-negative.

Now we prove the boundedness of the fractional-order model (1).

Theorem 3 The set

$\begin{array}{l}\Theta =\{\left(S_{\mu }\left(t\right),S_e\left(t\right),E\left(t\right),R\left(t\right),H\left(t\right),D\left(t\right),U\left(t\right)\right)\in R_{+}^7:\underset{}{\overset{}{}}\\ S_{\mu }\left(t\right)\le \frac{\left(1-\rho \right){\Lambda }^{\alpha }}{W_1},S_e\left(t\right)\le \frac{\left({\eta }^{\alpha }+\rho d^{\alpha }\right){\Lambda }^{\alpha }}{W_1{W}_2},N\left(t\right)\le \frac{{\Lambda }^{\alpha }}{d^{\alpha }}\}\end{array}$

is positively invariant with respect to model (1), where $W_1={\eta }^{\alpha }+d^{\alpha }$ and $W_2={\lambda }^{\alpha }+d^{\alpha }$ .

Proof. Summing up the seven equations of model (1) yields the governing equation for the total population dynamics as follows:

${}_0^c{D}_t^{\alpha }N\left(t\right)={\Lambda }^{\alpha }-d^{\alpha }N\left(t\right)$ (8)

Using the standard comparison Lemma 2.2, we have

$N\left(t\right)\le \left(N\left(0\right)-\frac{{\Lambda }^{\alpha }}{d^{\alpha }}\right)E_{\alpha }\left(-d^{\alpha }{t}^{\alpha }\right)+\frac{{\Lambda }^{\alpha }}{d^{\alpha }},$ (9)

where $N\left(0\right)=S_{\mu }\left(0\right)+S_e\left(0\right)+E\left(0\right)+R\left(0\right)+H\left(0\right)+D\left(0\right)+U\left(0\right)$ . It follows that $N\left(t\right)\to \frac{{\Lambda }^{\alpha }}{d^{\alpha }}$ as $t\to +\infty$ .

We now turn our attention to the first equation of model (1).

${}_0^c{D}_t^{\alpha }{S}_{\mu }\left(t\right)=\left(1-\rho \right){\Lambda }^{\alpha }-{\beta }^{\alpha }\widehat{\kappa }{S}_{\mu }\left(t\right)\left(R\left(t\right)+mD\left(t\right)\right)-\left({\eta }^{\alpha }+d^{\alpha }\right)S_{\mu }\left(t\right).$ (10)

Since

${\beta }^{\alpha }\widehat{\kappa }{S}_{\mu }\left(t\right)\left(R\left(t\right)+mD\left(t\right)\right)\ge 0,$ (11)

we obtain

${}_0^c{D}_t^{\alpha }{S}_{\mu }\left(t\right)\le \left(1-\rho \right){\Lambda }^{\alpha }-\left({\eta }^{\alpha }+d^{\alpha }\right)S_{\mu }\left(t\right),$ (12)

which consequently yields

$S_{\mu }\left(t\right)\le \left(S_{\mu }\left(0\right)-\frac{\left(1-\rho \right){\Lambda }^{\alpha }}{{\eta }^{\alpha }+d^{\alpha }}\right)E_{\gamma }\left(-\left({\eta }^{\alpha }+d^{\alpha }\right)t^{\alpha }\right)+\frac{\left(1-\rho \right){\Lambda }^{\alpha }}{{\eta }^{\alpha }+d^{\alpha }}.$ (13)

We may conclude that $S_{\mu }\left(t\right)\to \frac{\left(1-\rho \right){\Lambda }^{\alpha }}{{\eta }^{\alpha }+d^{\alpha }}$ as $t\to +\infty$ .

Similarly

$\begin{array}{l}\underset{t\to +\infty }{\lim }\sup S_e\left(t\right)\le \frac{\left({\eta }^{\alpha }+\rho d^{\alpha }\right){\Lambda }^{\alpha }}{W_1{W}_2}.\hfill \end{array}$ (14)

Hence every solution of the model (1) are lying within in Θ.

4.3. Existence of Equilibria

The equilibrium points are essential for the dynamical analysis of a system. This section focuses on the existence analysis of the two equilibrium points of the proposed model. Since the first six equations of the system (1) are independent of the compartment U, for simplicity of analysis, we consider the following subsystem:

$\begin{array}{l}{}_0^c{D}_t^{\alpha }{S}_{\mu }\left(t\right)=\left(1-\rho \right){\Lambda }^{\alpha }-{\beta }^{\alpha }\widehat{\kappa }{S}_{\mu }\left(t\right)\left(R\left(t\right)+mD\left(t\right)\right)-\left({\eta }^{\alpha }+d^{\alpha }\right)S_{\mu }\left(t\right),\\ {}_0^c{D}_t^{\alpha }{S}_e\left(t\right)=\rho {\Lambda }^{\alpha }+{\eta }^{\alpha }{S}_{\mu }\left(t\right)-\left({\lambda }^{\alpha }+d^{\alpha }\right)S_e\left(t\right)-\left(1-n\right){\beta }^{\alpha }\widehat{\kappa }{S}_e\left(t\right)\left(R\left(t\right)+mD\left(t\right)\right),\\ {}_0^c{D}_t^{\alpha }E\left(t\right)={\beta }^{\alpha }\widehat{\kappa }{S}_{\mu }\left(t\right)\left(R\left(t\right)+mD\left(t\right)\right)+\left(1-n\right){\beta }^{\alpha }\widehat{\kappa }{S}_e(t)\left(R\left(t\right)+mD\left(t\right)\right)\\ -\left({\xi }^{\alpha }+{\sigma }^{\alpha }+{\epsilon }^{\alpha }+d^{\alpha }\right)E\left(t\right),\\ {}_0^c{D}_t^{\alpha }R\left(t\right)={\xi }^{\alpha }E\left(t\right)+g^{\alpha }H\left(t\right)-\left({\theta }^{\alpha }+{\delta }^{\alpha }+{\omega }^{\alpha }+d^{\alpha }\right)R\left(t\right),\\ {}_0^c{D}_t^{\alpha }H\left(t\right)={\theta }^{\alpha }R\left(t\right)-\left(g^{\alpha }+{\gamma }^{\alpha }+d^{\alpha }\right)H\left(t\right),\\ {}_0^c{D}_t^{\alpha }D\left(t\right)={\sigma }^{\alpha }E\left(t\right)+{\delta }^{\alpha }R\left(t\right)-\left({\mu }^{\alpha }+d^{\alpha }\right)D\left(t\right),\end{array}$ (15)

with the non-negative initial conditions

$S_{\mu }\left(0\right)=S_{\mu 0},S_e\left(0\right)=S_{e0},E\left(0\right)=E_0,R\left(0\right)=R_0,H\left(0\right)=H_0,D\left(0\right)=D_0.$ (16)

The equilibrium points are determined by setting the right-hand side of the fractional-order equations in system (1) to zero, which yields:

$\{\begin{array}{l}\left(1-\rho \right){\Lambda }^{\alpha }-{\beta }^{\alpha }\widehat{\kappa }{S}_{\mu }\left(R+mD\right)-\left({\eta }^{\alpha }+d^{\alpha }\right)S_{\mu }=0,\\ \rho {\Lambda }^{\alpha }+{\eta }^{\alpha }{S}_{\mu }-\left(1-n\right){\beta }^{\alpha }\widehat{\kappa }{S}_e\left(R+mD\right)-\left({\lambda }^{\alpha }+d^{\alpha }\right)S_e=0,\\ {\beta }^{\alpha }\widehat{\kappa }{S}_{\mu }\left(R+mD\right)+\left(1-n\right){\beta }^{\alpha }\widehat{\kappa }{S}_e\left(R+mD\right)-\left({\xi }^{\alpha }+{\sigma }^{\alpha }+{\epsilon }^{\alpha }+d^{\alpha }\right)E=0,\\ {\xi }^{\alpha }E+g^{\alpha }H-\left({\theta }^{\alpha }+{\delta }^{\alpha }+{\omega }^{\alpha }+d^{\alpha }\right)R=0,\\ {\theta }^{\alpha }R-\left(g^{\alpha }+{\gamma }^{\alpha }+d^{\alpha }\right)H=0,\\ {\sigma }^{\alpha }E+{\delta }^{\alpha }R-\left({\mu }^{\alpha }+d^{\alpha }\right)D=0.\end{array}$ (17)

Solving this system allows one to identify the rumor-free equilibrium (RFE) and the rumor-spreading equilibrium (RSE).

4.3.1. RFE Equilibrium

The RFE is defined as the state in which the densities of all rumor-related compartments (e.g. spreaders, exposeed individuals) are zero, leaving only the susceptible population and possibly other non-infected groups. By solving the system (17) with no infected classes, it is obvious that the system (15) always exists a RFE point

$P^0=\left(S_{\mu }^0,S_e^0,E^0,R^0,H^0,D^0\right)=\left(\frac{\left(1-\rho \right){\Lambda }^{\alpha }}{W_1},\frac{\left({\eta }^{\alpha }+\rho d^{\alpha }\right){\Lambda }^{\alpha }}{W_1{W}_2},0,0,0,0\right).$ (18)

4.3.2. Basic Reproduction Number

In the study of rumor propagation dynamics, the basic reproduction number, $R_0^{\alpha }$ , serves as a crucial threshold parameter to quantify the spreading potential. Following the methodology in [29], we derive $R_0^{\alpha }$ using the next-generation matrix approach. Set $Y={\left(E,R,H,D\right)}^{\text{T}}$ , from system (15), we have ${}_0^c{D}_t^{\alpha }Y=ℱ-\mathcal{V}$ , where

$ℱ=\left[\begin{array}{c}{\beta }^{\alpha }\widehat{\kappa }\left(R+mD\right)\left(S_{\mu }+\left(1-n\right)S_e\right)\\ 0\\ 0\\ 0\end{array}\right],\mathcal{V}=\left[\begin{array}{c}{W}_{3}E\\ -{\xi }^{\alpha }E-g^{\alpha }H+W_{4}R\\ -{\theta }^{\alpha }R+W_{5}H\\ -{\sigma }^{\alpha }E-{\delta }^{\alpha }R+W_{6}D\end{array}\right].$

Here, the auxiliary parameters are defined as:

$W_3={\xi }^{\alpha }+{\sigma }^{\alpha }+{\epsilon }^{\alpha }+d^{\alpha },W_4={\theta }^{\alpha }+{\delta }^{\alpha }+{\omega }^{\alpha }+d^{\alpha },$

$W_5=g^{\alpha }+{\gamma }^{\alpha }+d^{\alpha },W_6={\mu }^{\alpha }+d^{\alpha }.$

The Jacobian matrices $F$ and $V$ at the RFE $P^0$ are

$\begin{array}{l}F=\left[\begin{array}{cccc}0& {\beta }^{\alpha }\widehat{\kappa }\left(S_{\mu }^0+\left(1-n\right)S_e^0\right)& 0& m{\beta }^{\alpha }\widehat{\kappa }\left(S_{\mu }^0+\left(1-n\right)S_e^0\right)\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right],\\ V=\left[\begin{array}{cccc}{W}_3& 0& 0& 0\\ -{\xi }^{\alpha }& W_4& -g^{\alpha }& 0\\ 0& -{\theta }^{\alpha }& W_5& 0\\ -{\sigma }^{\alpha }& -{\delta }^{\alpha }& 0& W_6\end{array}\right].\end{array}$

Now, let $\Delta =W_3{W}_6\left(W_4{W}_5-g^{\alpha }{\theta }^{\alpha }\right)$ . The inverse matrix $V^{-1}$ is

$V^{-1}=\frac{1}{\Delta }\left[\begin{array}{cccc}\frac{\Delta }{W_3}& 0& 0& 0\\ {\xi }^{\alpha }{W}_5{W}_6& W_3{W}_5{W}_6& g^{\alpha }{W}_3{W}_6& 0\\ {\xi }^{\alpha }{\theta }^{\alpha }{W}_6& {\theta }^{\alpha }{W}_3{W}_6& W_3{W}_4{W}_6& 0\\ {\Sigma }_{41}& {\delta }^{\alpha }{W}_3{W}_5& {\delta }^{\alpha }{g}^{\alpha }{W}_3& \frac{\Delta }{W_6}\end{array}\right],$

where ${\Sigma }_{41}={\sigma }^{\alpha }\left(W_4{W}_5-g^{\alpha }{\theta }^{\alpha }\right)+{\xi }^{\alpha }{\delta }^{\alpha }{W}_5$ .

Thus, $F{V}^{-1}=\frac{{\beta }^{\alpha }\widehat{\kappa }\left(S_{\mu }^0+\left(1-n\right)S_e^0\right)}{\Delta }\cdot M$ , where $M$ is a matrix with non-zero first row elements:

$M_{11}={\xi }^{\alpha }{W}_5{W}_6+m{\Sigma }_{41},M_{12}=W_3{W}_5{W}_6+m{\delta }^{\alpha }{W}_3{W}_5,$

$M_{13}=g^{\alpha }{W}_3{W}_6+m{\delta }^{\alpha }{g}^{\alpha }{W}_3,M_{14}=m\frac{\Delta }{W_6}.$

Finally, the basic reproduction number $R_0^{\alpha }$ is obtained as:

$R_0^{\alpha }=\rho \left(F{V}^{-1}\right)=\frac{{\beta }^{\alpha }\widehat{\kappa }\left(S_{\mu }^0+\left(1-n\right)S_e^0\right)W_0}{W_3{W}_6\left(W_4{W}_5-g^{\alpha }{\theta }^{\alpha }\right)},$

where $W_0={\xi }^{\alpha }{W}_5{W}_6+m{\sigma }^{\alpha }\left(W_4{W}_5-g^{\alpha }{\theta }^{\alpha }\right)+m{\xi }^{\alpha }{\delta }^{\alpha }{W}_5$ .

4.3.3. Existence and Uniqueness of RSE

The existence of a RSE in rumor propagation model is mathematically determined by specific conditions related to the basic reproduction number and model parameters. From model (17), the existence of RSE $P^{\text{*}}=\left(S_{\mu }^{\text{*}},S_e^{\text{*}},E^{\text{*}},R^{\text{*}},H^{\text{*}},D^{\text{*}}\right)$ depends on the positive root $R^{\text{*}}$ of the following equation

$G\left(R\right)=b_2{\left({\beta }^{\alpha }\widehat{\kappa }{W}_{0}R\right)}^2+b_1\left({\beta }^{\alpha }\widehat{\kappa }{W}_{0}R\right)+b_0=0,$ (19)

where

$\begin{array}{l}{S}_{\mu }^{*}=\frac{\left(1-\rho \right){\Lambda }^{\alpha }}{{\beta }^{\alpha }\widehat{\kappa }\left(R^{*}+m{D}^{*}\right)+W_1},S_e^{*}=\frac{\rho {\Lambda }^{\alpha }+{\eta }^{\alpha }{S}_{\mu }^{*}}{\left(1-n\right){\beta }^{\alpha }\widehat{\kappa }\left(R^{*}+m{D}^{*}\right)+W_2}\\ E^{*}=\frac{W_4{W}_5-g^{\alpha }{\theta }^{\alpha }}{{\xi }^{\alpha }{W}_5}{R}^{*},H^{*}=\frac{{\theta }^{\alpha }}{W_5}{R}^{*},D^{*}=\frac{{\sigma }^{\alpha }\left(W_4{W}_5-g^{\alpha }{\theta }^{\alpha }\right)+{\xi }^{\alpha }{\delta }^{\alpha }{W}_5}{{\xi }^{\alpha }{W}_5{W}_6}{R}^{*}\end{array}$ (20)

and

$\begin{array}{l}{b}_2=\left(1-n\right)W_3\left(W_4{W}_5-g^{\alpha }{\theta }^{\alpha }\right),\\ b_1={\xi }^{\alpha }\left(\left(1-n\right)W_1+W_2\right)W_3\left(W_4{W}_5-g^{\alpha }{\theta }^{\alpha }\right)W_5{W}_6-\left(1-n\right){\beta }^{\alpha }\widehat{\kappa }{\Lambda }^{\alpha }{\xi }^{\alpha }{W}_5{W}_0,\\ b_0={\xi }^{2\alpha }{W}_1{W}_2{W}_3\left(W_4{W}_5-g^{\alpha }{\theta }^{\alpha }\right)W_5^2{W}_6^2\left(1-R_0^{\alpha }\right).\end{array}$ (21)

We are going to prove it in two cases.

If $R_0^{\alpha }<1$ , then $\left(\left(1-n\right)W_1+W_2\right)W_3{W}_6\left(W_4{W}_5-g^{\alpha }{\theta }^{\alpha }\right)>\left(1-n\right){\beta }^{\alpha }\widehat{\kappa }{\Lambda }^{\alpha }{W}_0$ . It follows that $b_2>0$ , $b_1>0$ , $b_0>0$ . Thus, all the coefficients of $G\left(R\right)$ are positive, which suggests $G\left(R\right)$ has no positive real roots.

If $R_0^{\alpha }>1$ , then $b_2>0$ and $b_0<0$ . Under this assumption, the two roots of the equation $G\left(R\right)$ are given by:

$R_1^{*}=\frac{-b_1+\sqrt{b_1^2-4b_2{b}_0}}{2b_2{\beta }^{\alpha }\widehat{\kappa }{W}_0},R_2^{*}=\frac{-b_1-\sqrt{b_1^2-4b_2{b}_0}}{2b_2{\beta }^{\alpha }\widehat{\kappa }{W}_0}.$(22)

This shows $G\left(R\right)$ has a unique positive real root.

4.4. Stability Analysis of Equilibria

Practically, the stability properties of system (1) are of great importance in characterizing its dynamical behavior. This section focuses on analyzing the stability of the equilibrium points to reveal the underlying trends of rumor propagation.

4.4.1. Local Stability of RFE and RSE in Model (15)

The local asymptotic stability (LAS) of the RFE $P^0$ and the RSE $P^{\text{*}}$ of model (15) is established in this subsection.

Theorem 4 The RFE $P^0$ of model (15) is LAS if $R_0^{\alpha }<1$ and unstable if $R_0^{\alpha }>1$ .

Proof. The Jacobian matrix of model (15) at $P^0$ is:

$J\left(P^0\right)=\left[\begin{array}{cccccc}-W_1& 0& 0& -{\beta }^{\alpha }\widehat{\kappa }{S}_{\mu }^0& 0& -m{\beta }^{\alpha }\widehat{\kappa }{S}_{\mu }^0\\ {\eta }^{\alpha }& -W_2& 0& -\left(1-n\right){\beta }^{\alpha }\widehat{\kappa }{S}_e^0& 0& -m\left(1-n\right){\beta }^{\alpha }\widehat{\kappa }{S}_e^0\\ 0& 0& -W_3& ℬ& 0& mℬ\\ 0& 0& {\xi }^{\alpha }& -W_4& g^{\alpha }& 0\\ 0& 0& 0& {\theta }^{\alpha }& -W_5& 0\\ 0& 0& {\sigma }^{\alpha }& {\delta }^{\alpha }& 0& -W_6\end{array}\right],$ (23)

where $ℬ={\beta }^{\alpha }\widehat{\kappa }\left(S_{\mu }^0+\left(1-n\right)S_e^0\right)$ . The characteristic polynomial is:

${\Delta }^1\left(s^{\alpha }\right)=\left|s^{\alpha }I-J\left(P^0\right)\right|=\left(s^{\alpha }+W_1\right)\left(s^{\alpha }+W_2\right){\Delta }^0\left(s^{\alpha }\right).$ (24)

Obviously, ${\text{Δ}}^1\left(s^{\alpha }\right)$ has two negative real roots: $-W_1$ and $-W_2$ . We now examine the distribution of the roots of ${\text{Δ}}^0\left(s^{\alpha }\right)$ . In fact, ${\text{Δ}}^0\left(s^{\alpha }\right)=0$ satisfies

$\begin{array}{l}\left(s^{\alpha }+W_3\right)\left(s^{\alpha }+W_6\right)\left[\left(s^{\alpha }+W_4\right)\left(s^{\alpha }+W_5\right)-g^{\alpha }{\theta }^{\alpha }\right]\\ =ℬ[{\xi }^{\alpha }\left(s^{\alpha }+W_5\right)\left(s^{\alpha }+W_6\right)+m{\sigma }^{\alpha }\left(\left(s^{\alpha }+W_4\right)\left(s^{\alpha }+W_5\right)-g^{\alpha }{\theta }^{\alpha }\right)\\ +m{\xi }^{\alpha }{\delta }^{\alpha }\left(s^{\alpha }+W_5\right)].\end{array}$ (25)

To prove that all remaining eigenvalues have negative real parts when $R_0^{\alpha }<1$ , we proceed by contradiction. Suppose $s^{\alpha }$ is an eigenvalue with a nonnegative real part, that is, $Re\left(s^{\alpha }\right)\ge 0$ . Then, we find that it satisfies the inequality:

$\begin{array}{c}1=\frac{ℬ[\dots ]}{\left(s^{\alpha }+W_3\right)\left(s^{\alpha }+W_6\right)\left[\left(s^{\alpha }+W_4\right)\left(s^{\alpha }+W_5\right)-g^{\alpha }{\theta }^{\alpha }\right]}\\ \le \frac{ℬ\left({\xi }^{\alpha }{W}_5{W}_6+m{\sigma }^{\alpha }\left(W_4{W}_5-g^{\alpha }{\theta }^{\alpha }\right)+m{\xi }^{\alpha }{\delta }^{\alpha }{W}_5\right)}{W_3{W}_6\left(W_4{W}_5-g^{\alpha }{\theta }^{\alpha }\right)}\\ =R_0^{\alpha },\end{array}$ (26)

which is a contradiction. Therefore, if $R_0^{\alpha }<1$ , all eigenvalues of ${\text{Δ}}^0\left(s^{\alpha }\right)$ possess negative real parts, which implies that the RFE $P^0$ of system (15) is LAS. Conversely, if $R_0^{\alpha }>1$ , then

${\text{Δ}}^0\left(s^{\alpha }\right)=s^{4\alpha }+A_3{s}^{3\alpha }+A_2{s}^{2\alpha }+A_1{s}^{\alpha }+A_0,$ (27)

where

$\begin{array}{l}{A}_3=W_3+W_4+W_5+W_6,\\ A_2=W_3\left(W_4+W_5+W_6\right)+W_4\left(W_5+W_6\right)+W_5{W}_6-g^{\alpha }{\theta }^{\alpha }\\ -{\beta }^{\alpha }\widehat{\kappa }\left(S_{\mu }^0+\left(1-n\right)S_e^0\right)\left({\xi }^{\alpha }+m{\sigma }^{\alpha }\right),\\ A_1=W_3\left(W_4{W}_5+W_4{W}_6+W_5{W}_6\right)+W_4{W}_5{W}_6-g^{\alpha }{\theta }^{\alpha }\left(W_3+W_6\right)\\ -{\beta }^{\alpha }\widehat{\kappa }\left(S_{\mu }^0+\left(1-n\right)S_e^0\right)\left({\xi }^{\alpha }{W}_5+{\xi }^{\alpha }{W}_6+m{\xi }^{\alpha }{\delta }^{\alpha }+m{\sigma }^{\alpha }{W}_4+m{\sigma }^{\alpha }{W}_5\right),\\ A_0=W_3{W}_6\left(W_4{W}_5-g^{\alpha }{\theta }^{\alpha }\right)-{\beta }^{\alpha }\widehat{\kappa }\left(S_{\mu }^0+\left(1-n\right)S_e^0\right)\\ \cdot \left({\xi }^{\alpha }{W}_5{W}_6+m{\sigma }^{\alpha }\left(W_4{W}_5-g^{\alpha }{\theta }^{\alpha }\right)+m{\xi }^{\alpha }{\delta }^{\alpha }{W}_5\right)\\ =W_4{W}_5\left(W_4{W}_5-g^{\alpha }{\theta }^{\alpha }\right)\left(1-R_0^{\alpha }\right).\end{array}$(28)

It is clear that ${\text{Δ}}^0\left(0\right)=A_0<0$ and $\underset{s^{\alpha }\to +\infty }{\lim }{\Delta }^0\left(s^{\alpha }\right)=+\infty$ . Consequently, the equation ${\text{Δ}}^0\left(s^{\alpha }\right)=0$ possesses at least one positive root, implying the instability of $P^0$ when $R_0^{\alpha }>1$ .

4.4.2. Global Stability of the RFE and RSE

The global stability analysis (GAS) of the RFE $P^0$ and the RSE $P^{\text{*}}$ of model (15) is established in this subsection.

Theorem 5 The RFE $P^0$ is GAS if $R_0^{\alpha }<1$ .

Proof. A Lyapunov function is defined as follows

$\begin{array}{c}{\Phi }_1\left(t\right)=\left(S_{\mu }\left(t\right)-S_{\mu }^0-S_{\mu }^0\ln \frac{S_{\mu }\left(t\right)}{S_{\mu }^0}\right)+\left(S_e\left(t\right)-S_e^0-S_e^0\ln \frac{S_e\left(t\right)}{S_e^0}\right)\\ +E\left(t\right)+k_{1}R\left(t\right)+k_{2}H\left(t\right)+k_{3}D\left(t\right),\end{array}$ (29)

where $k_1=\frac{W_3{W}_5{W}_6+m{\delta }^{\alpha }{W}_3{W}_5}{W_0}$ , $k_2=\frac{g^{\alpha }{W}_3{W}_6+m{\delta }^{\alpha }{g}^{\alpha }{W}_3}{W_0}$ , $k_3=\frac{m{W}_3\left(W_4{W}_5-g^{\alpha }{\theta }^{\alpha }\right)}{W_0}$ .

Taking $\alpha$ order Caputo derivative ${}_0^c{D}_t^{\alpha }$ of ${\Phi }_1\left(t\right)$ , we can get

$\begin{array}{c}{}_0^c{D}_t^{\alpha }{\Phi }_1\left(t\right)\le \left(1-\frac{S_{\mu }^0}{S_{\mu }}\right){}_0^c{D}_t^{\alpha }{S}_{\mu }\left(t\right)+\left(1-\frac{S_e^0}{S_e}\right){}_0^c{D}_t^{\alpha }{S}_e\left(t\right)+{}_0^c{D}_t^{\alpha }E\left(t\right)\\ +k_1{}_0^c{D}_t^{\alpha }R\left(t\right)+k_2{}_0^c{D}_t^{\alpha }H\left(t\right)+k_3{}_0^c{D}_t^{\alpha }D\left(t\right)\\ =\left(1-\frac{S_{\mu }^0}{S_{\mu }}\right)\left[\left(1-\rho \right){\Lambda }^{\alpha }-{\beta }^{\alpha }\widehat{\kappa }{S}_{\mu }\left(R+mD\right)-\left({\eta }^{\alpha }+d^{\alpha }\right)S_{\mu }\right]\\ +\left(1-\frac{S_e^0}{S_e}\right)\left[\rho {\Lambda }^{\alpha }+{\eta }^{\alpha }{S}_{\mu }-\left(1-n\right){\beta }^{\alpha }\widehat{\kappa }{S}_e\left(R+mD\right)-\left({\lambda }^{\alpha }+d^{\alpha }\right)S_e\right]\\ +\left[{\beta }^{\alpha }\widehat{\kappa }\left(S_{\mu }+\left(1-n\right)S_e\right)\left(R+mD\right)-\left({\xi }^{\alpha }+{\sigma }^{\alpha }+{\epsilon }^{\alpha }+d^{\alpha }\right)E\right]\\ +k_1\left[{\xi }^{\alpha }E+g^{\alpha }H-\left({\theta }^{\alpha }+{\delta }^{\alpha }+{\omega }^{\alpha }+d^{\alpha }\right)R\right]\\ +k_2\left[{\theta }^{\alpha }R-\left(g^{\alpha }+{\gamma }^{\alpha }+d^{\alpha }\right)H\right]+k_3\left[{\sigma }^{\alpha }E+{\delta }^{\alpha }R-\left({\mu }^{\alpha }+d^{\alpha }\right)D\right].\end{array}$ (30)

Combining $\left(1-\rho \right){\Lambda }^{\alpha }=\left({\eta }^{\alpha }+d^{\alpha }\right)S_{\mu }^0$ and $\rho {\Lambda }^{\alpha }+{\eta }^{\alpha }{S}_{\mu }^0=\left({\lambda }^{\alpha }+d^{\alpha }\right)S_e^0$ , we can see that

$\begin{array}{c}{}_0^c{D}_t^{\alpha }{\Phi }_1\left(t\right)\le d^{\alpha }{S}_{\mu }^0\left(2-\frac{S_{\mu }}{S_{\mu }^0}-\frac{S_{\mu }^0}{S_{\mu }}\right)+\rho {\Lambda }^{\alpha }\left(2-\frac{S_e}{S_e^0}-\frac{S_e^0}{S_e}\right)\\ +{\eta }^{\alpha }{S}_{\mu }^0\left(3-\frac{S_{\mu }^0}{S_{\mu }}-\frac{S_e}{S_e^0}-\frac{S_{\mu }{S}_e^0}{S_{\mu }^0{S}_e}\right)\\ +{\beta }^{\alpha }\widehat{\kappa }\left(S_{\mu }^0+\left(1-n\right)S_e^0\right)\left(R+mD\right)\\ -\frac{W_3{W}_6\left(W_4{W}_5-g^{\alpha }{\theta }^{\alpha }\right)}{W_0}\left(R+mD\right)\\ =d^{\alpha }{S}_{\mu }^0\left(2-\frac{S_{\mu }}{S_{\mu }^0}-\frac{S_{\mu }^0}{S_{\mu }}\right)+\rho {\Lambda }^{\alpha }\left(2-\frac{S_e}{S_e^0}-\frac{S_e^0}{S_e}\right)\\ +{\eta }^{\alpha }{S}_{\mu }^0\left(3-\frac{S_{\mu }^0}{S_{\mu }}-\frac{S_e}{S_e^0}-\frac{S_{\mu }{S}_e^0}{S_{\mu }^0{S}_e}\right)\\ +\frac{W_3{W}_6\left(W_4{W}_5-g^{\alpha }{\theta }^{\alpha }\right)}{W_0}\left(R_0^{\alpha }-1\right)\left(R+mD\right).\end{array}$(31)

Since ${}_0^c{D}_t^{\alpha }{\Phi }_1\left(t\right)\le 0$ for all $t\ge 0$ when $R_0^{\alpha }<1$ , it follows from LaSalle’s invariance principle that the RFE $P^0$ is GAS.

Theorem 6 The RSE $P^{*}$ is GAS if $R_0^{\alpha }>1$ and $\aleph =S_{\mu }^{*}-\min \left\{\frac{{\lambda }^{\alpha }+d^{\alpha }}{{\eta }^{\alpha }},\frac{\left(1-n\right){\sigma }^{\alpha }\left(W_4{W}_5-g^{\alpha }{\delta }^{\alpha }\right)}{{\xi }^{\alpha }{\delta }^{\alpha }{W}_5}\right\}{S}_e^{*}\le 0$ .

Proof. Let

$\begin{array}{l}x=\frac{S_{\mu }}{S_{\mu }^{\text{*}}},y=\frac{S_e}{S_e^{\text{*}}},z=\frac{E}{E^{\text{*}}},\\ e=\frac{R}{R^{\text{*}}},f=\frac{H}{H^{\text{*}}},v=\frac{D}{D^{\text{*}}}.\end{array}$ (32)

We now transform model (15) into the following form:

$\{\begin{array}{l}{}_0^c{D}_t^{\alpha }x=x\left[\frac{\left(1-\rho \right){\Lambda }^{\alpha }}{S_{\mu }^{*}}\left(\frac{1}{x}-1\right)+{\beta }^{\alpha }\widehat{\kappa }{R}^{*}\left(1-e\right)+m{\beta }^{\alpha }\widehat{\kappa }{D}^{*}\left(1-v\right)\right],\\ {}_0^c{D}_t^{\alpha }y=y[\frac{\rho {\Lambda }^{\alpha }}{S_e^{*}}\left(\frac{1}{y}-1\right)+\frac{{\eta }^{\alpha }{S}_{\mu }^{*}}{S_e^{*}}\left(\frac{x}{y}-1\right)+\left(1-n\right){\beta }^{\alpha }\widehat{\kappa }{R}^{*}\left(1-e\right)\\ \begin{array}{c}\\ \end{array}+m\left(1-n\right){\beta }^{\alpha }\widehat{\kappa }{D}^{*}\left(1-v\right)],\\ {}_0^c{D}_t^{\alpha }z=\frac{z}{E^{*}}[{\beta }^{\alpha }\widehat{\kappa }{S}_{\mu }^{*}{R}^{*}\left(\frac{xe}{z}-1\right)+\left(1-n\right){\beta }^{\alpha }\widehat{\kappa }{S}_e^{*}{R}^{*}\left(\frac{ye}{z}-1\right)\\ +m{\beta }^{\alpha }\widehat{\kappa }{S}_{\mu }^{*}{D}^{*}\left(\frac{xv}{z}-1\right)+m\left(1-n\right){\beta }^{\alpha }\widehat{\kappa }{S}_e^{*}{D}^{*}\left(\frac{yv}{z}-1\right)],\\ {}_0^c{D}_t^{\alpha }e=\frac{e}{R^{*}}\left[{\xi }^{\alpha }{E}^{*}\left(\frac{z}{e}-1\right)+g^{\alpha }{H}^{*}\left(\frac{f}{e}-1\right)\right],\\ {}_0^c{D}_t^{\alpha }f=\frac{f{\theta }^{\alpha }{R}^{*}}{H^{*}}\left(\frac{e}{f}-1\right),\\ {}_0^c{D}_t^{\alpha }v=\frac{v}{D^{*}}\left[{\sigma }^{\alpha }{E}^{*}\left(\frac{z}{v}-1\right)+{\delta }^{\alpha }{R}^{*}\left(\frac{e}{v}-1\right)\right].\end{array}$(33)

Define Lyapunov function ${\Phi }_2\left(t\right)$ as follows:

$\begin{array}{c}{\Phi }_2\left(t\right)=a_1{S}_{\mu }^{*}\left(x-1-\ln x\right)+a_2{S}_e^{*}\left(y-1-\ln y\right)+a_3{E}^{*}\left(z-1-\ln z\right)\\ +a_4{R}^{*}\left(e-1-\ln e\right)+a_5{H}^{*}\left(f-1-\ln f\right)+a_6{D}^{*}\left(v-1-\ln v\right),\end{array}$ (34)

where

$\begin{array}{l}{a}_1=a_2=a_3=1,\\ a_4=\frac{{\beta }^{\alpha }\widehat{\kappa }{R}^{*}\left(S_{\mu }^{*}+\left(1-n\right)S_e^{*}\right)}{{\xi }^{\alpha }{E}^{*}}+\frac{m{\beta }^{\alpha }\widehat{\kappa }{D}^{*}\left(S_{\mu }^{*}+\left(1-n\right)S_e^{*}\right){\delta }^{\alpha }{R}^{*}}{{\xi }^{\alpha }{E}^{*}\left({\sigma }^{\alpha }{E}^{*}+{\delta }^{\alpha }{R}^{*}\right)},\\ a_5=\frac{{\beta }^{\alpha }\widehat{\kappa }{R}^{*}\left(S_{\mu }^{*}+\left(1-n\right)S_e^{*}\right)g^{\alpha }{H}^{*}}{{\xi }^{\alpha }{E}^{*}{\theta }^{\alpha }{R}^{*}}+\frac{m{\beta }^{\alpha }\widehat{\kappa }{D}^{*}\left(S_{\mu }^{*}+\left(1-n\right)S_e^{*}\right){\delta }^{\alpha }{R}^{*}{g}^{\alpha }{H}^{*}}{{\xi }^{\alpha }{E}^{*}{\theta }^{\alpha }{R}^{*}\left({\sigma }^{\alpha }{E}^{*}+{\delta }^{\alpha }{R}^{*}\right)},\\ a_6=\frac{m{\beta }^{\alpha }\widehat{\kappa }{D}^{*}\left(S_{\mu }^{*}+\left(1-n\right)S_e^{*}\right)}{{\sigma }^{\alpha }{E}^{*}+{\delta }^{\alpha }{R}^{*}}.\end{array}$ (35)

By applying the Caputo derivative of order $\alpha$ to ${\Phi }_2\left(t\right)$ , we obtain

$\begin{array}{c}{}_0^c{D}_t^{\alpha }{\Phi }_2\left(t\right)\le a_1{S}_{\mu }^{*}\left(1-\frac{1}{x}\right){}_0^c{D}_t^{\alpha }x+a_2{S}_e^{*}\left(1-\frac{1}{y}\right){}_0^c{D}_t^{\alpha }y+a_3{E}^{*}\left(1-\frac{1}{z}\right){}_0^c{D}_t^{\alpha }z\\ +a_4{R}^{*}\left(1-\frac{1}{e}\right){}_0^c{D}_t^{\alpha }e+a_5{H}^{*}\left(1-\frac{1}{f}\right){}_0^c{D}_t^{\alpha }f+a_6{D}^{*}\left(1-\frac{1}{v}\right){}_0^c{D}_t^{\alpha }v\\ =a_1\left(x-1\right)\left[\left(1-\rho \right){\Lambda }^{\alpha }\left(\frac{1}{x}-1\right)+{\beta }^{\alpha }\widehat{\kappa }{S}_{\mu }^{*}{R}^{*}\left(1-e\right)+m{\beta }^{\alpha }\widehat{\kappa }{S}_{\mu }^{*}{D}^{*}\left(1-v\right)\right]\\ +a_2\left(y-1\right)[\rho {\Lambda }^{\alpha }\left(\frac{1}{y}-1\right)+{\eta }^{\alpha }{S}_{\mu }^{*}\left(\frac{x}{y}-1\right)\\ \stackrel{}{}+\left(1-n\right){\beta }^{\alpha }\widehat{\kappa }{S}_e^{*}{R}^{*}\left(1-e\right)+m\left(1-n\right){\beta }^{\alpha }\widehat{\kappa }{S}_e^{*}{D}^{*}\left(1-v\right)]\\ +a_3\left(z-1\right)[{\beta }^{\alpha }\widehat{\kappa }{S}_{\mu }^{*}{R}^{*}\left(\frac{xe}{z}-1\right)+\left(1-n\right){\beta }^{\alpha }\widehat{\kappa }{S}_e^{*}{R}^{*}\left(\frac{ye}{z}-1\right)\\ +m{\beta }^{\alpha }\widehat{\kappa }{S}_{\mu }^{*}{D}^{*}\left(\frac{xv}{z}-1\right)+m\left(1-n\right){\beta }^{\alpha }\widehat{\kappa }{S}_e^{*}{D}^{*}\left(\frac{yv}{z}-1\right)]\\ +a_4\left(e-1\right)\left[{\xi }^{\alpha }{E}^{*}\left(\frac{z}{e}-1\right)+g^{\alpha }{H}^{*}\left(\frac{f}{e}-1\right)\right]\\ +a_5\left(f-1\right)\left[{\theta }^{\alpha }{R}^{*}\left(\frac{e}{f}-1\right)\right]\\ +a_6\left(v-1\right)\left[{\sigma }^{\alpha }{E}^{*}\left(\frac{z}{v}-1\right)+{\delta }^{\alpha }{R}^{*}\left(\frac{e}{v}-1\right)\right].\end{array}$ (36)

After some simple calculations, we have

$\begin{array}{c}{}_0^c{D}_t^{\alpha }{\Phi }_2\left(t\right)\le h_1\left(2-x-\frac{1}{x}\right)+h_2\left(2-y-\frac{1}{y}\right)+h_3\left(3-\frac{1}{x}-y-\frac{x}{y}\right)\\ +h_4\left(3-\frac{1}{x}-\frac{z}{e}-\frac{xe}{z}\right)+h_5\left(3-\frac{1}{x}-\frac{z}{v}-\frac{xv}{z}\right)+h_6\left(3-\frac{1}{y}-\frac{z}{e}-\frac{ye}{z}\right)\\ +h_7\left(3-\frac{1}{y}-\frac{z}{v}-\frac{yv}{z}\right)+h_8\left(4-\frac{1}{y}-\frac{z}{e}-\frac{e}{v}-\frac{yv}{z}\right)+h_9\left(2-\frac{f}{e}-\frac{e}{f}\right).\end{array}$ (37)

where

$\begin{array}{l}{h}_1=d^{\alpha }{S}_{\mu }^{\text{*}},h_2=\left({\lambda }^{\alpha }+d^{\alpha }\right)S_e^{\text{*}}-{\eta }^{\alpha }{S}_{\mu }^{\text{*}},h_3={\eta }^{\alpha }{S}_{\mu }^{\text{*}},\\ h_4={\beta }^{\alpha }\widehat{\kappa }{S}_{\mu }^{\text{*}}{R}^{\text{*}},h_5=m{\beta }^{\alpha }\widehat{\kappa }{S}_{\mu }^{\text{*}}{D}^{\text{*}},h_6=\left(1-n\right){\beta }^{\alpha }\widehat{\kappa }{S}_e^{\text{*}}{R}^{\text{*}},\\ h_7=\frac{\left(1-n\right)m{\beta }^{\alpha }\widehat{\kappa }{S}_e^{\text{*}}{D}^{\text{*}}{\sigma }^{\alpha }{E}^{\text{*}}-m{\beta }^{\alpha }\widehat{\kappa }{S}_{\mu }^{\text{*}}{D}^{\text{*}}{\delta }^{\alpha }{R}^{\text{*}}}{{\sigma }^{\alpha }{E}^{\text{*}}+{\delta }^{\alpha }{R}^{\text{*}}},\\ h_8=\frac{m{\beta }^{\alpha }\widehat{\kappa }{D}^{\text{*}}\left(S_{\mu }^{\text{*}}+\left(1-n\right)S_e^{\text{*}}\right){\delta }^{\alpha }{R}^{\text{*}}}{{\sigma }^{\alpha }{E}^{\text{*}}+{\delta }^{\alpha }{R}^{\text{*}}},\\ h_9=\frac{{\beta }^{\alpha }\widehat{\kappa }{R}^{\text{*}}\left(S_{\mu }^{\text{*}}+\left(1-n\right)S_e^{\text{*}}\right)g^{\alpha }{H}^{\text{*}}}{{\xi }^{\alpha }{E}^{\text{*}}}+\frac{m{\beta }^{\alpha }\widehat{\kappa }{D}^{\text{*}}\left(S_{\mu }^{\text{*}}+\left(1-n\right)S_e^{\text{*}}\right){\delta }^{\alpha }{R}^{\text{*}}{g}^{\alpha }{H}^{\text{*}}}{{\xi }^{\alpha }{E}^{\text{*}}\left({\sigma }^{\alpha }{E}^{\text{*}}+{\delta }^{\alpha }{R}^{\text{*}}\right)}.\end{array}$(38)

Due to (37), to guarantee $h_i\ge 0$ ($i=1,2,\cdots ,9$ ), we only need:

$S_{\mu }^{\text{*}}\le \text{min}\left\{\frac{{\lambda }^{\alpha }+d^{\alpha }}{{\eta }^{\alpha }},\frac{\left(1-n\right){\sigma }^{\alpha }\left(W_4{W}_5-g^{\alpha }{\delta }^{\alpha }\right)}{{\xi }^{\alpha }{\delta }^{\alpha }{W}_5}\right\}{S}_e^{\text{*}}.$

Therefore, ${}_0^c{D}_t^{\alpha }{\Phi }_2\le 0$ , and the equality holds only for $x=y=1$ , $z=e=f=v$ . That is,

$\left\{\left(x,y,z,e,f,v\right)\in \Omega :{}_0^c{D}_t^{\alpha }{\Phi }_2=0\right\}=\left\{\left(x,y,z,e,f,v\right):x=y=1,z=e=f=v\right\}.$ (39)

Thus,

$\left\{\left(S_{\mu },S_e,E,R,H,D\right):S_{\mu }=S_{\mu }^{*},S_e=S_e^{*},\frac{E^{*}}{E}=\frac{R^{*}}{R}=\frac{H^{*}}{H}=\frac{D^{*}}{D}\right\}\cap \Theta =P^{*}.$ (40)

Based on LaSalle’s invariable principle and asymptotic stability theorem, it can be concluded that the RSE $P^{\text{*}}$ is GAS.

5. Fractional Optimal Control Problem

In this section, we examine two key controlling parameters in the rumor propagation model. Firstly, external intervention measures such as platform content regulation, compulsory dissemination of authoritative information, and social-network traffic restrictions directly promote the transition of rumor spreaders ($R$ ) to the immune group ($U$ ) via the control variable ${\omega }^{\alpha }$ . This process, which is distinct from spontaneous debunking behavior (represented by $\delta$ ), is a compulsory intervention administered through the time-varying control function $\omega \left(t\right)$ . It is designed to rapidly establish herd immunity during the peak phase of rumor propagation. Secondly, preventive interventions and re-education, denoted by the control variable $\gamma$ , apply to dormant potential spreaders ($H$ ). The core effect of this mechanism is to immunize this subpopulation, thus blocking their reactivation and ensuring their direct transition to the immune compartment ($U$ ). The intervention intensity, governed by the time-varying control function $\gamma \left(t\right)$ , is dynamically adjusted. It is proactively intensified during critical periods when rumor resurgence is most likely, aiming to effectively curb the reformation of rumor transmission chains.

Building upon the three intervention mechanisms introduced in the introduction, we now formulate them within an optimal control framework. The control variable $\omega \left(t\right)$ realizes the debunking mechanism by enforcing the direct transfer of active spreaders ($R$ ) into stiflers ($U$ ), while $\gamma \left(t\right)$ operationalizes the educational and memory-forgetting mechanisms by immunizing hibernated individuals ($H$ ) against reactivation, thus blocking their return to the spreading compartment. This mapping allows us to dynamically optimize the intensity of these two control levers to balance rumor suppression with intervention costs.

To incorporate the control strategies, we upgrade the constant parameters ${\omega }^{\alpha }$ and ${\gamma }^{\alpha }$ in system (1) to time-varying control inputs $\omega \left(t\right)$ and $\gamma \left(t\right)$ , respectively. The resulting controlled system is governed by the following equations.

$\{\begin{array}{l}{}_0^c{D}_t^{\alpha }{S}_{\mu }\left(t\right)=\left(1-\rho \right){\Lambda }^{\alpha }-{\beta }^{\alpha }\widehat{\kappa }{S}_{\mu }\left(t\right)\left(R\left(t\right)+mD\left(t\right)\right)-\left({\eta }^{\alpha }+d^{\alpha }\right)S_{\mu }\left(t\right),\\ {}_0^c{D}_t^{\alpha }{S}_e\left(t\right)=\rho {\Lambda }^{\alpha }+{\eta }^{\alpha }{S}_{\mu }\left(t\right)-\left({\lambda }^{\alpha }+d^{\alpha }\right)S_e\left(t\right)-\left(1-n\right){\beta }^{\alpha }\widehat{\kappa }{S}_e\left(t\right)\left(R\left(t\right)+mD\left(t\right)\right),\\ {}_0^c{D}_t^{\alpha }E\left(t\right)={\beta }^{\alpha }\widehat{\kappa }{S}_{\mu }\left(t\right)\left(R\left(t\right)+mD\left(t\right)\right)+\left(1-n\right){\beta }^{\alpha }\widehat{\kappa }{S}_e\left(t\right)\left(R\left(t\right)+mD\left(t\right)\right)\\ -\left({\xi }^{\alpha }+{\sigma }^{\alpha }+{\epsilon }^{\alpha }+d^{\alpha }\right)E\left(t\right),\\ {}_0^c{D}_t^{\alpha }R\left(t\right)={\xi }^{\alpha }E\left(t\right)+g^{\alpha }H\left(t\right)-{\theta }^{\alpha }R\left(t\right)-{\delta }^{\alpha }R\left(t\right)-\omega \left(t\right)R\left(t\right)-d^{\alpha }R\left(t\right),\\ {}_0^c{D}_t^{\alpha }H\left(t\right)={\theta }^{\alpha }R\left(t\right)-g^{\alpha }H\left(t\right)-\gamma \left(t\right)H\left(t\right)-d^{\alpha }H\left(t\right),\\ {}_0^c{D}_t^{\alpha }D\left(t\right)={\sigma }^{\alpha }E\left(t\right)+{\delta }^{\alpha }R\left(t\right)-\left({\mu }^{\alpha }+d^{\alpha }\right)D\left(t\right),\\ {}_0^c{D}_t^{\alpha }U\left(t\right)={\lambda }^{\alpha }{S}_e\left(t\right)+{\epsilon }^{\alpha }E\left(t\right)+\omega \left(t\right)R\left(t\right)+\gamma \left(t\right)H\left(t\right)+{\mu }^{\alpha }D\left(t\right)-d^{\alpha }U\left(t\right).\end{array}$(41)

with the non-negative initial conditions:

$\begin{array}{l}{S}_{\mu }\left(0\right)=S_{\mu 0},S_e\left(0\right)=S_{e0},E\left(0\right)=E_0,R\left(0\right)=R_0,\\ H\left(0\right)=H_0,D\left(0\right)=D_0,U\left(0\right)=U_0.\end{array}$ (42)

Building upon system (41), our goal is to suppress the number of rumor spreaders within a short time horizon while simultaneously minimizing the associated intervention costs (i.e., costs of debunking efforts and education-enhancement measures). To balance these two purposes, we formulate the corresponding optimal control problem. This amounts to finding the control functions $\omega \left(t\right)$ and $\gamma \left(t\right)$ that minimize the objective functional $J$ defined by:

$J\left(\omega \left(t\right),\gamma \left(t\right)\right)={\displaystyle {\int }_0^{\tau }\left[M_{0}R\left(t\right)+\frac{M_1}{2}{\omega }^2\left(t\right)+\frac{M_2}{2}{\gamma }^2\left(t\right)\right]\text{d}t}.$ (43)

where $M_0$ denotes the weight coefficient assigned to $R$ , while $M_1$ and $M_2$ are the weights for the control functions $\omega \left(t\right)$ and $\gamma \left(t\right)$ , respectively. The parameter $\tau$ signifies the total time horizon of the control intervention.

In this objective functional, we focus on minimizing the density of active rumor spreaders $R\left(t\right)$ because they are the primary drivers of rumor propagation and the most direct source of public harm. While the exposed individuals $E\left(t\right)$ and hibernated individuals $H\left(t\right)$ also represent potential risks, they contribute to rumor spread indirectly. Controlling $R\left(t\right)$ effectively reduces the influx into both $E\left(t\right)$ and $H\left(t\right)$ , and this simplified formulation allows for a tractable analytical derivation of the optimal control strategy. The quadratic terms on $\omega \left(t\right)$ and $\gamma \left(t\right)$ penalize excessive intervention effort, reflecting the economic and social costs associated with implementing control measures.

The admissible set for the two control variables, $\omega \left(t\right)$ and $\gamma \left(t\right)$ , is defined as follows:

$\Theta =\left\{\left(\omega ,\gamma \right)|0\le \omega \left(t\right)\le {\omega }_{\max },0\le \gamma \left(t\right)\le {\gamma }_{\max },t\in \left[0,\tau \right]\right\}.$ (44)

where ${\omega }_{\text{max}}$ and ${\gamma }_{\text{max}}$ are fixed positive constants.

The Lagrange function is defined as follows:

$L=M_{0}R\left(t\right)+\frac{M_1}{2}{\omega }^2\left(t\right)+\frac{M_2}{2}{\gamma }^2\left(t\right).$ (45)

The Hamiltonian function is described as follows:

$\begin{array}{c}\mathscr{H}=M_{0}R\left(t\right)+\frac{M_1}{2}{\omega }^2\left(t\right)+\frac{M_2}{2}{\gamma }^2\left(t\right)\\ +{\lambda }_{S\mu }\left[\left(1-\rho \right){\Lambda }^{\alpha }-{\beta }^{\alpha }\widehat{\kappa }{S}_{\mu }\left(t\right)\left(R\left(t\right)+mD\left(t\right)\right)-\left({\eta }^{\alpha }+d^{\alpha }\right)S_{\mu }\left(t\right)\right]\\ +{\lambda }_{Se}\left[\rho {\Lambda }^{\alpha }+{\eta }^{\alpha }{S}_{\mu }\left(t\right)-\left(1-n\right){\beta }^{\alpha }\widehat{\kappa }{S}_e\left(t\right)\left(R\left(t\right)+mD\left(t\right)\right)-\left({\lambda }^{\alpha }+d^{\alpha }\right)S_e\left(t\right)\right]\\ +{\lambda }_E\left[{\beta }^{\alpha }\widehat{\kappa }{S}_{\mu }\left(t\right)\left(R\left(t\right)+mD\left(t\right)\right)+\left(1-n\right){\beta }^{\alpha }\widehat{\kappa }{S}_e\left(t\right)\left(R\left(t\right)+mD\left(t\right)\right)-W_{3}E\left(t\right)\right]\\ +{\lambda }_R\left[{\xi }^{\alpha }E\left(t\right)+g^{\alpha }H\left(t\right)-\left({\theta }^{\alpha }+{\delta }^{\alpha }+\omega \left(t\right)+d^{\alpha }\right)R\left(t\right)\right]\\ +{\lambda }_H\left[{\theta }^{\alpha }R\left(t\right)-\left(g^{\alpha }+\gamma \left(t\right)+d^{\alpha }\right)H\left(t\right)\right]\\ +{\lambda }_D\left[{\sigma }^{\alpha }E\left(t\right)+{\delta }^{\alpha }R\left(t\right)-\left({\mu }^{\alpha }+d^{\alpha }\right)D\left(t\right)\right]\\ +{\lambda }_U\left[{\lambda }^{\alpha }{S}_e\left(t\right)+{\epsilon }^{\alpha }E\left(t\right)+\omega \left(t\right)R\left(t\right)+\gamma \left(t\right)H\left(t\right)+{\mu }^{\alpha }D\left(t\right)-d^{\alpha }U\left(t\right)\right].\end{array}$ (46)

Here ${\lambda }_{S\mu },{\lambda }_{Se},{\lambda }_E,{\lambda }_R,{\lambda }_H,{\lambda }_D,{\lambda }_U$ represent adjoint variables, respectively.

Theorem 5.1. Let $\left(S_{\mu \text{*}},S_{e\text{*}},E_{\text{*}},R_{\text{*}},H_{\text{*}},D_{\text{*}},U_{\text{*}}\right)$ be the optimal state trajectories associated with the optimal controls $\left({\omega }_{\text{*}}\left(t\right),{\gamma }_{\text{*}}\left(t\right)\right)$ that minimize the objective functional $J\left(\omega \left(t\right),\gamma \left(t\right)\right)$ subject to system (43). Then, there exist adjoint variables $\left({\lambda }_{S_{\mu }},{\lambda }_{S_e},{\lambda }_E,{\lambda }_R,{\lambda }_H,{\lambda }_D,{\lambda }_U\right)$ satisfying the following adjoint system:

$\{\begin{array}{l}{}_0^c{D}_t^{\alpha }{\lambda }_{S\mu }={\beta }^{\alpha }\widehat{\kappa }\left(R_{*}+m{D}_{*}\right)\left({\lambda }_{S\mu }-{\lambda }_E\right)+{\eta }^{\alpha }\left({\lambda }_{S\mu }-{\lambda }_{Se}\right)+d^{\alpha }{\lambda }_{S\mu },\\ {}_0^c{D}_t^{\alpha }{\lambda }_{Se}=\left(1-n\right){\beta }^{\alpha }\widehat{\kappa }\left(R_{*}+m{D}_{*}\right)\left({\lambda }_{Se}-{\lambda }_E\right)+{\lambda }^{\alpha }\left({\lambda }_{Se}-{\lambda }_U\right)+d^{\alpha }{\lambda }_{Se},\\ {}_0^c{D}_t^{\alpha }{\lambda }_E={\xi }^{\alpha }\left({\lambda }_E-{\lambda }_R\right)+{\sigma }^{\alpha }\left({\lambda }_E-{\lambda }_D\right)+{\epsilon }^{\alpha }\left({\lambda }_E-{\lambda }_U\right)+d^{\alpha }{\lambda }_E,\\ {}_0^c{D}_t^{\alpha }{\lambda }_R=-M_0+{\beta }^{\alpha }\widehat{\kappa }{S}_{\mu *}\left({\lambda }_{S\mu }-{\lambda }_E\right)+\left(1-n\right){\beta }^{\alpha }\widehat{\kappa }{S}_{e*}\left({\lambda }_{Se}-{\lambda }_E\right)\\ +{\theta }^{\alpha }\left({\lambda }_R-{\lambda }_H\right)+{\delta }^{\alpha }\left({\lambda }_R-{\lambda }_D\right)+{\omega }_{*}\left(t\right)\left({\lambda }_R-{\lambda }_U\right)+d^{\alpha }{\lambda }_R,\\ {}_0^c{D}_t^{\alpha }{\lambda }_H=g^{\alpha }\left({\lambda }_H-{\lambda }_R\right)+{\gamma }_{*}\left(t\right)\left({\lambda }_H-{\lambda }_U\right)+d^{\alpha }{\lambda }_H,\\ {}_0^c{D}_t^{\alpha }{\lambda }_D=m{\beta }^{\alpha }\widehat{\kappa }{S}_{\mu *}\left({\lambda }_{S\mu }-{\lambda }_E\right)+\left(1-n\right)m{\beta }^{\alpha }\widehat{\kappa }{S}_{e*}\left({\lambda }_{Se}-{\lambda }_E\right)\\ +{\mu }^{\alpha }\left({\lambda }_D-{\lambda }_U\right)+d^{\alpha }{\lambda }_D,\\ {}_0^c{D}_t^{\alpha }{\lambda }_U=d^{\alpha }{\lambda }_U.\end{array}$(47)

with the transversality conditions ${\lambda }_{S\mu }\left(\tau \right)=0$ , ${\lambda }_{Se}\left(\tau \right)=0$ , ${\lambda }_E\left(\tau \right)=0$ , ${\lambda }_R\left(\tau \right)=0$ , ${\lambda }_H\left(\tau \right)=0$ , ${\lambda }_D\left(\tau \right)=0$ , ${\lambda }_U\left(\tau \right)=0$ .

Futher, the optimal control ${\omega }_{\text{*}}\left(t\right)$ and ${\gamma }_{\text{*}}\left(t\right)$ are given by

$\begin{array}{l}{\omega }_{*}\left(t\right)=\max \left\{\min \left\{\frac{R\left(t\right)\left({\lambda }_R-{\lambda }_U\right)}{M_1},{\omega }_{\max }\right\},0\right\},\\ {\gamma }_{*}\left(t\right)=\max \left\{\min \left\{\frac{H\left(t\right)\left({\lambda }_H-{\lambda }_U\right)}{M_2},{\gamma }_{\max }\right\},0\right\}.\end{array}$ (48)

Proof. The solution to the fractional optimal control problem, defined for dynamic system (41) with the performance index (43), corresponds to the minimal value of the Lagrangian function at the point $\left(t,S_{\mu \text{*}},S_{e\text{*}},E_{\text{*}},R_{\text{*}},D_{\text{*}},U_{\text{*}},{\omega }_{\text{*}},{\gamma }_{\text{*}}\right)$ . These optimal values are derived by applying the necessary optimality conditions, specifically, the fractional Euler Lagrange equations in the Caputo sense, as established by Agrawal [30]. Accordingly, the following set of equations defines the adjoint system for the problem:

$\begin{array}{l}{}_0^c{D}_t^{\alpha }{\lambda }_{S\mu }=-\frac{\partial \mathscr{H}}{\partial S_{\mu }},{}_0^c{D}_t^{\alpha }{\lambda }_{Se}=-\frac{\partial \mathscr{H}}{\partial S_e},{}_0^c{D}_t^{\alpha }{\lambda }_E=-\frac{\partial \mathscr{H}}{\partial E},\\ {}_0^c{D}_t^{\alpha }{\lambda }_R=-\frac{\partial \mathscr{H}}{\partial R},{}_0^c{D}_t^{\alpha }{\lambda }_H=-\frac{\partial \mathscr{H}}{\partial H},{}_0^c{D}_t^{\alpha }{\lambda }_D=-\frac{\partial \mathscr{H}}{\partial D},{}_0^c{D}_t^{\alpha }{\lambda }_U=-\frac{\partial \mathscr{H}}{\partial U}.\end{array}$ (49)

Since the Hamiltonian function $\mathscr{H}$ is quadratic with respect to the control variables $\omega$ and $\gamma$ , its minimum occurs when the first-order optimality conditions are satisfied, i.e., the partial derivatives of $\mathscr{H}$ with respect to $\omega$ and $\gamma$ equal zero. Applying these conditions, $\frac{\partial \mathscr{H}}{\partial \omega }=0$ and $\frac{\partial \mathscr{H}}{\partial \gamma }=0$ , yields the following expressions for the optimal controls:

$\begin{array}{l}\frac{\partial \mathscr{H}}{\partial \omega }=M_1\omega +R\left({\lambda }_U-{\lambda }_R\right)=0,\\ \frac{\partial \mathscr{H}}{\partial \gamma }=M_2\gamma +H\left({\lambda }_U-{\lambda }_H\right)=0.\end{array}$ (50)

By employing the Pontryagin’s Maximum Principle, the constraints defining the admissible control set are succinctly captured, leading to a compact formulation. The corresponding optimal control values are then derived as follows:

$\begin{array}{l}{\omega }_{\text{*}}\left(t\right)=\text{max}\left\{\text{min}\left\{\frac{R\left(t\right)\left({\lambda }_R-{\lambda }_U\right)}{M_1},{\omega }_{\text{max}}\right\},0\right\},\\ {\gamma }_{\text{*}}\left(t\right)=\text{max}\left\{\text{min}\left\{\frac{H\left(t\right)\left({\lambda }_H-{\lambda }_U\right)}{M_2},{\gamma }_{\text{max}}\right\},0\right\}.\end{array}$ (51)

6. Numerical Simulation

In this section, numerical simulations are conducted to validate the theoretical analysis and examine the effectiveness of the proposed $S_{\mu }{S}_{e}ERHDU$ rumor propagation model. The numerical investigation consists of three main components.

First, the dynamic behaviors of the $S_{\mu }{S}_{e}ERHDU$ model are explored using the Adams–Bashforth–Moulton predictor–corrector (ABMPC) method. Second, different prevention and control strategies are evaluated under a unified simulation framework. Finally, optimal control strategies are implemented and analyzed via the forward–backward sweep method (FBSM). To illustrate the model’s performance under different social interaction scenarios, three representative network settings, referred to Facebook, Twitch RU, and Twitch PT, are considered as benchmark cases. These network designations are utilized primarily for comparative analysis and illustrative objectives, rather than representing empirically derived network structures.

To differ The model parameters used in Scenarios 1 and 2 (Section 6.1), as well as in the optimal control simulations (Section 6.2), are synthetic and chosen for illustrative purposes to satisfy the theoretical threshold conditions ($R_0^{\alpha }<1$ and $R_0^{\alpha }>1$ ) and to demonstrate a range of dynamical behaviors. The natural loss rate $d^{\alpha }$ is set equal to the entry rate ${\Lambda }^{\alpha }$ to maintain a constant total population, facilitating comparison across scenarios.

To differentiate among the various network interaction scenarios, all model parameters remain consistent across the three cases, with the exception of the contact intensity parameter $\widehat{\kappa }$ . The three network cases—referred to as “Facebook”, “Twitch RU”, and “Twitch PT”—serve exclusively as illustrative benchmarks for different social interaction intensities and should not be interpreted as empirical representations of these platforms. This parameter $\widehat{\kappa }$ quantifies the general degree of social interaction and information exposure within a network. Specifically, a value of $\widehat{\kappa }=14.98$ is assigned to the Facebook scenario, whereas $\widehat{\kappa }=32.74$ and $\widehat{\kappa }=43.69$ are used for the Twitch RU and Twitch PT scenarios, respectively. This configuration allows us to isolate the impact of interaction intensity on rumor propagation dynamics. This parameter $\widehat{\kappa }$ quantifies the general degree of social interaction and information exposure within a network. Specifically, a value of $\widehat{\kappa }=14.98$ is assigned to the Facebook scenario, whereas $\widehat{\kappa }=32.74$ and $\widehat{\kappa }=43.69$ are used for the Twitch RU and Twitch PT scenarios, respectively. This configuration allows us to isolate the impact of interaction intensity on rumor propagation dynamics.

6.1. Features of the $S_{\mu }{S}_{e}ERHDU$ Model

We first verify the theoretical findings regarding the basic reproduction number and rumor propagation dynamics via numerical simulations. Since the proposed model is deterministic in a homogeneously mixed population, a single run of the numerical solver is sufficient to obtain the unique solution for each fixed parameter set and initial condition. Below are the numerical findings across varying scenarios of the basic reproduction number.

Scenario 1: ($R_0^{\alpha }<1$ )

To simulate a scenario where rumors naturally die out, the following synthetic parameter values are chosen such that the basic reproduction number remains below unity: $\alpha =0.98$ , $\rho =0.334$ , ${\Lambda }^{\alpha }=0.082$ , $d^{\alpha }=0.082$ , ${\beta }^{\alpha }=0.01$ , $n=0.7$ , $m=1.1$ , ${\eta }^{\alpha }=0.03$ , ${\lambda }^{\alpha }=0.1$ , ${\xi }^{\alpha }=0.1$ , ${\sigma }^{\alpha }=0.1$ , ${ϵ}^{\alpha }=0.12$ , ${\theta }^{\alpha }=0.3$ , ${\delta }^{\alpha }=0.3$ , ${\omega }^{\alpha }=0.12$ , $g^{\alpha }=0.4$ , ${\gamma }^{\alpha }=0.2$ , ${\mu }^{\alpha }=0.4$ . The assumption ${\Lambda }^{\alpha }=d^{\alpha }$ is adopted to simplify the tracking of population transfers between nodes. Based on these parameters, the basic reproduction numbers $R_0^{\alpha }$ are computed as 0.109 for the Facebook network, 0.239 for the Twitch RU network, and 0.319 for the Twitch PT network. Correspondingly, the contact intensity parameter $\widehat{\kappa }$ is set to 14.98, 32.74, and 43.69 for the Facebook, Twitch RU, and Twitch PT scenarios, respectively.

Scenario 2: ($R_0^{\alpha }>1$ )

To simulate a scenario where rumors persist and reach an endemic equilibrium, the following synthetic parameter values are chosen such that the basic reproduction number exceeds unity: $\alpha =0.98$ , $\rho =0.334$ , ${\Lambda }^{\alpha }=0.082$ , $d^{\alpha }=0.082$ , ${\beta }^{\alpha }=0.7$ , $n=0.8$ , $m=1.2$ , ${\eta }^{\alpha }=0.05$ , ${\lambda }^{\alpha }=0.2$ , ${\xi }^{\alpha }=0.3$ , ${\sigma }^{\alpha }=0.4$ , ${ϵ}^{\alpha }=0.15$ , ${\theta }^{\alpha }=0.2$ , ${\delta }^{\alpha }=0.11$ , ${\omega }^{\alpha }=0.2$ , $g^{\alpha }=0.3$ , ${\gamma }^{\alpha }=0.1$ , ${\mu }^{\alpha }=0.25$ , we evaluate the basic reproduction number $R_0^{\alpha }$ for the Facebook, Twitch RU, and Twitch PT networks. The corresponding values are computed as 11.405, 24.893, and 33.218, respectively.

6.1.1. Stability Verification of the RFE and RSE

In this section, we employ the Facebook network as a case study to verify the GAS of both the RFE and the RSE in the proposed model.

Figure 2. Stability of the RFE $P^0$ for different initial values.

The dynamics of individual densities for each node ($S_{\mu },S_e,E,R,H,D$ and $U$ ) over time under different initial conditions are shown in Figure 2. When the basic reproduction number satisfies $R_0^{\alpha }<1$ , the system demonstrates pronounced convergence behavior: $S_{\mu }$ approaches $\frac{\left(1-\rho \right){\Lambda }^{\alpha }}{W_1}$ , $S_e$ tends to $\frac{\left({\eta }^{\alpha }+\rho d^{\alpha }\right){\Lambda }^{\alpha }}{W_1{W}_2}$ , while $E,R,H$ and $D$ all converge to 0. Meanwhile, $U$ approaches $\frac{\left({\eta }^{\alpha }+\rho d^{\alpha }\right){\epsilon }^{\alpha }{\Lambda }^{\alpha }}{d^{\alpha }{W}_1{W}_2}$ . These results not only confirm the GAS of the RFE but also indicate that, regardless of the initial conditions, the system eventually reaches a rumor-extinction equilibrium state when $R_0^{\alpha }<1$ .

Figure 3. Stability of the RSE $P^{\text{*}}$ for different initial values.

Figure 3 illustrates the dynamic evolution of node $R$ under different initial conditions and various values of $\aleph$ . The results demonstrate that when $R_0^{\alpha }>1$ :

1) If $\aleph \le 0$ , the RSE of the model (1) is GAS.

2) If $\aleph >0$ , the RSE remains GAS.

This indicates that, under the condition $R_0^{\alpha }>1$ , the system eventually converges to the RSE regardless of the initial conditions, which leads to the following conjecture.

Conjecture 1. The RSE $P^{\text{*}}$ is GAS if $R_0^{\alpha }>1$ .

Remark. Theorem 4.6 provides $\aleph \le 0$ as a sufficient condition for the global asymptotic stability of $P^{\text{*}}$ , which is established through a rigorous Lyapunov analysis. The numerical results in Figure 3, however, show that $P^{\text{*}}$ remains globally asymptotically stable even when $\aleph >0$ , suggesting that this condition is conservative. Therefore, Conjecture 1, which claims that $R_0^{\alpha }>1$ alone guarantees global stability, is supported by numerical evidence but remains to be proven analytically.

6.1.2. The Dynamic Characteristics of $S_{\mu }{S}_{e}ERHDU$ Model

In this section, we simulate the rumor spreading process on the Facebook, Twitch RU, and Twitch PT networks. We assume the following initial conditions:

We assume a normalized total population ($N=1$ ) with the following initial conditions: a single rumor spreader is introduced into a network where 30% of individuals are initially in the educated susceptible state, and the remainder are uneducated. All other compartments start empty except for the uneducated susceptibles, which constitute the rest of the population.

$\begin{array}{l}{S}_{\mu }\left(0\right)=0.699,S_e\left(0\right)=0.3,E\left(0\right)=0,\\ R\left(0\right)=0.001,H\left(0\right)=0,D\left(0\right)=0,U\left(0\right)=0.\end{array}$

These settings represent a population in which a single initial rumor spreader is introduced into a network with 30% of individuals in the susceptible-exposed class and the remainder in the general susceptible class.

Figure 4. Scenario 1: Evolution of rumor diffusion in different networks when $R_0^{\alpha }<1$ .

Figure 5. Scenario 2: Evolution of rumor diffusion in different networks when $R_0^{\alpha }>1$ .

Figure 4 and Figure 5 visually depict the propagation dynamics of rumors across different networks. When the basic reproduction number satisfies $R_0^{\alpha }>1$ , the density of hesitant individuals increases rapidly, peaks, and then declines before eventually stabilizing. This pattern indicates that rumors propagate swiftly within the three network types. Similarly, the densities of rumor spreaders and hibernators rise to a peak before decreasing, ultimately reaching a steady state. This suggests that, in the absence of external intervention, rumor dissemination naturally tends toward an equilibrium. Concurrently, the density of rumor refuters increases monotonically before stabilizing, reflecting a collective self-correction mechanism. In contrast, when $R_0^{\alpha }<1$ , the density of spreaders remains near zero in all three networks, implying that no rumor outbreak occurs. These results confirm that the basic reproduction number is a key indicator of rumor propagation potential. A comparison of the two propagation scenarios reveals that the negative impact of Scenario 2 is more pronounced; therefore, subsequent mechanism simulations will focus on Scenario 2.

6.1.3. Three Prevention and Control Strategies

In this section, we simulate the effect of three prevention and control mechanisms on rumor propagation:

1) Education mechanism: By enhancing users’ digital literacy and their ability to assess information credibility, the likelihood of accepting and spreading rumors is reduced at the cognitive level.

2) Memory and forgetting mechanism: Drawing on natural forgetting and memory reactivation processes among users, this mechanism captures the gradual attenuation of rumor influence over time.

3) Refutation mechanism: Through the timely release of authoritative information via official or insider channels, this strategy uses factual evidence to clarify rumors and inhibit their further spread.

Figure 6. Density of rumor spreaders over time under different $\rho$ and $\eta$ .

Figure 6 demonstrates that the education mechanism significantly curbs the spread of online rumors. As education intensifies, the public’s ability to screen information improves, which effectively reduces the number of rumor spreaders and slows down rumor propagation. This mechanism offers two key advantages: first, it exhibits strong preventive characteristics, effectively slowing rumor spread in its early stages; second, its impact is cumulative-as the public’s overall cognitive level continues to rise, it builds growing resistance to rumors over time.

Figure 7. Density of rumor spreaders over time under different $\theta$ and $g$ .

As shown in Figure 7, with the enhancement of memory and forgetting mechanisms, the density of rumor spreaders decreased significantly. The internal mechanism of this phenomenon is that the forgetting mechanism promotes the transformation of some rumor spreaders into hibernates, which directly reduces the scale of rumor spread groups. Although the memory mechanism can re-activate some hibernators, its target is still derived from the hibernating group produced by the forgetting mechanism. This dynamic balance eventually leads to a decline in the peak size of rumor spreaders. It is worth noting that this rule is applicable in all three networks, which is consistent with the numerical simulation results of Zhao et al. [16], further confirming that the improvement of memory and forgetting mechanism can effectively inhibit the spread of rumors.

Figure 8. Density of rumor spreaders over time under different $\delta$ .

Figure 8 demonstrates that increasing the rumor refutation rate $\delta$ leads to a significant suppression of rumor propagation across all three networks. Specifically, a higher refutation rate results in a notable reduction in both the peak number of rumor spreaders and the final stabilized density of spreaders, thereby curbing the dissemination of rumors more rapidly. These findings provide strong evidence that the timely release of authoritative information through official channels can effectively restrain the spread of rumors. The simulation results further highlight the crucial role of the refutation mechanism in controlling rumor propagation.

As observed in Figures 8(a)-(c), increasing the transformation rate $\delta$ from 0% to 25% produces a substantial reduction in the peak prevalence of spreaders. In contrast, a further increase in $\delta$ from 25% to 50% produces only about half of the mitigation effect achieved by the initial rise. This pattern suggests that, in the early stage of rumor control, elevating $\delta$ effectively curtails the propagation source while simultaneously competing for a larger share of the susceptible population. Once the intervention reaches a moderate level, the base numbers of both spreaders and susceptible individuals have already been significantly reduced. As a result, the same control effort is applied to a smaller target population, leading to diminishing marginal returns–that is, a weaker mitigating effect per unit increase in $\delta$ .

These findings indicate that although the debunking mechanism is highly effective at the outset, its cost-effectiveness diminishes in later stages, as the additional resource input yields disproportionately smaller outcomes. This phenomenon underscores the importance of rational resource allocation in practical applications. It is not always optimal to maximize the control intensity; instead, resources should be directed toward interventions that deliver the highest marginal benefit to maximize the overall control efficiency.

Figure 9. The evolution condition of rumor spreaders with the change of time for different values of $\rho$ , $\eta$ , $\lambda$ , $\theta$ , $g$ and $\delta$ .

A comparison between Figure 9 and Figures 6-8 reveals that the combined application of the three types of prevention and control mechanisms effectively suppresses sudden rumor outbreaks across all three networks. This integrated approach results in slower propagation speeds and a significantly reduced scale of rumor dissemination. The inhibitory effect of this comprehensive strategy is markedly superior to that achieved by any single mechanism applied alone. Specifically, the combined approach not only delays the arrival of the rumor peak and reduces its magnitude, but also accelerates the decline phase of rumor transmission, thereby significantly shortening its overall lifecycle. Such an outcome cannot be achieved by any individual control mechanism operating in isolation.

6.2. Numerical Simulation of Optimal Control of the $S_{\mu }{S}_{e}ERHDU$ Model

To numerically validate the effectiveness and feasibility of the proposed optimal control strategies, we perform a series of simulations for the fractional-order $S_{\mu }{S}_{e}ERHDU$ rumor-spreading model. Based on the theoretical results obtained in the previous section, the numerical experiments are implemented in the MATLAB environment by employing the fractional Euler method combined with the forward–backward predict–evaluate–correct–evaluate (PECE) scheme [31].

Two time-dependent control functions are incorporated into the system: $\omega \left(t\right)$ , which promotes the direct transition of spreaders to immune individuals, and $\gamma \left(t\right)$ , which suppresses the reactivation of dormant individuals into active spreaders. To enable a quantitative evaluation of intervention effects across different control mechanisms, the following four comparative scenarios are considered to comprehensively examine the impact of individual and combined control strategies.

1) No control;

2) Only $\omega \left(t\right)$ control;

3) Only $\gamma \left(t\right)$ control;

4) Combined $\omega \left(t\right)$ and $\gamma \left(t\right)$ control.

For the optimal control simulations, we adopt a representative synthetic parameter set that satisfies $R_0^{\alpha }>1$ to ensure rumor persistence in the absence of control, thereby allowing a meaningful evaluation of intervention effectiveness. The model parameters are fixed as $\rho =0.334$ , $\Lambda =0.082$ , $d=0.07$ , $\beta =0.7$ , $\widehat{\kappa }=15$ , $n=0.8$ , $m=1.2$ , $\eta =0.05$ , $\lambda =0.2$ , $\xi =0.3$ , $\sigma =0.4$ , $ϵ=0.15$ , $\theta =0.12$ , $\delta =0.07$ , $g=0.35$ , and $\mu =0.18$ , which are chosen to be consistent with the Scenario 2 parameter set. The initial conditions are taken as $S_{\mu }\left(0\right)=0.60$ , $S_e\left(0\right)=0.25$ , $E\left(0\right)=0$ , $R\left(0\right)=0.10$ , $H\left(0\right)=0.05$ , $D\left(0\right)=0$ , and $U\left(0\right)=0$ .

Figure 10 illustrates the temporal evolution of the densities of spreaders $R\left(t\right)$ and hibernators $H\left(t\right)$ in the absence of control measures under different fractional orders $\alpha$ . The comparison highlights the influence of the memory effect embedded in the fractional-order model.

(a) (b)

Figure 10. Densities of spreaders ($R$ ) and hibernators ($H$ ) in the absence of control measures. (a) Spreader density $R\left(t\right)$ for different $\alpha$ , (b) Hibernator density $H\left(t\right)$ for different $\alpha$ .

As illustrated in Figure 10(a), the spreader density $R\left(t\right)$ rises sharply in the initial phase across all considered values of $\alpha$ , which signals the onset of rumor outbreak. However, distinct discrepancies emerge in terms of both the peak magnitude and the subsequent decay rate. Specifically, a larger fractional order $\alpha$ yields a higher peak of $R\left(t\right)$ followed by a more rapid decline, whereas smaller values of $\alpha$ correspond to a lower peak and a more protracted spreading process. This phenomenon underscores that a stronger memory effect (associated with smaller $\alpha$ ) decelerates the system response and prolongs the impact of historical states on the current dynamics of rumor propagation.

Figure 10(b) depicts the temporal evolution of the hibernator population $H\left(t\right)$ across varying values of $\alpha .$ In contrast to the spreader population, $H\left(t\right)$ exhibits a more gradual and smooth variation over time. Nevertheless, the fractional order still exerts a pronounced influence on its long-term dynamical behavior. For smaller values of $\alpha ,$ the accumulation of hibernators proceeds at a slower pace, and the convergence toward the steady-state level is significantly delayed, which suggests that memory effects impede the transition of individuals into the dormant state. Conversely, larger $\alpha$ values expedite this transition process and thus facilitate the earlier stabilization of $H\left(t\right).$

In summary, Figure 10 reveals that the fractional order $\alpha$ significantly influences the uncontrolled rumor dynamics. A smaller $\alpha$ prolongs rumor persistence and retards system stabilization, whereas a larger $\alpha ,$ associated with weaker memory effects, leads to more rapid attenuation. These findings underscore the need for implementing control strategies to effectively curb rumor propagation, particularly in systems exhibiting pronounced memory characteristics.

Figure 11. Impact of the single control strategy $\omega \left(t\right)$ on the densities of spreaders ($R$ ) and hibernators ($H$ ).

Conversely, when compared with the uncontrolled scenario, the temporal evolution of the hibernator population $H\left(t\right)$ also undergoes a marked reduction under the control strategy that only employs $\omega \left(t\right)$ . While $\omega \left(t\right)$ does not exert a direct regulatory effect on hibernators, the effective suppression of spreaders serves to curtail the influx of individuals transitioning into the hibernating state.

Consequently, the peak of $H\left(t\right)$ is reduced, and its overall magnitude throughout the simulation period declines correspondingly. This finding indicates that the impact of $\omega \left(t\right)$ on the hibernator population, while indirect, remains significant and primarily arises from the coupling between spreaders and hibernators within the rumor propagation process.

Collectively, the results presented in Figure 11 confirm that the single $\omega \left(t\right)$ control strategy is highly effective in mitigating both the intensity and duration of rumor spreading by directly suppressing active spreaders, while also producing an indirect yet discernible effect on the hibernating population.

Figure 12. Impact of the single control strategy $\gamma \left(t\right)$ on the densities of spreaders ($R$ ) and hibernators ($H$ ).

Figure 12 shows the impacts of the standalone $\gamma \left(t\right)$ control strategy on the temporal evolution of the spreader density $R\left(t\right)$ and hibernator density $H\left(t\right)$ . In comparison with the uncontrolled scenario, the implementation of $\gamma \left(t\right)$ induces a marked reduction in the hibernator density $H\left(t\right)$ across the entire simulation period. Notably, both the peak magnitude and the long-term equilibrium level of $H\left(t\right)$ are diminished substantially, which indicates that the $\gamma \left(t\right)$ control strategy can effectively inhibit the reactivation of dormant individuals.

With respect to the spreader population $R\left(t\right)$ , the impact exerted by $\gamma \left(t\right)$ is relatively moderate. Although $R\left(t\right)$ exhibits a noticeable reduction compared with the uncontrolled scenario, the magnitude of this decline is smaller than that induced by the $\omega \left(t\right)$ -only control strategy. This phenomenon is consistent with theoretical expectations, since $\gamma \left(t\right)$ does not act directly on spreaders; instead, it curtails the replenishment of the spreader class by restricting the transition of individuals from the hibernating state.

On the whole, Figure 12 demonstrates that the single $\gamma \left(t\right)$ control strategy plays a crucial role in reducing the size of the hibernator population and mitigating the potential for rumor resurgence, while its influence on active spreaders remains indirect. These results highlight the complementary nature of $\gamma \left(t\right)$ and $\omega \left(t\right)$ in the rumor suppression process.

Figure 13. Dynamics of spreaders $R$ , hibernators $H$ , and corresponding optimal control intensities under the combined $\omega \left(t\right)$ and $\gamma \left(t\right)$ strategy.

Figure 13 depicts the combined effects of the optimal control strategies $\omega \left(t\right)$ and $\gamma \left(t\right)$ on the dynamics of rumor propagation. As illustrated in Figure 13(a), the concurrent implementation of these two control measures yields the lowest peak value and the most rapid decline in the spreader population R(t), in contrast to the uncontrolled scenario and the individual single-control strategies. This finding demonstrates a distinct synergistic effect between $\omega \left(t\right)$ and $\gamma \left(t\right)$ in curbing the activity of rumor spreaders.

Figure 13(b) presents the evolutionary trend of the hibernator population $H\left(t\right)$ across distinct control scenarios. Under the combined control strategy, the density of hibernators maintains a consistently lower level over the entire simulation period. This trend indicates that $\gamma \left(t\right)$ effectively curbs the reactivation of dormant individuals, while $\omega \left(t\right)$ indirectly curtails the influx into the hibernating state by suppressing the spreader population.

Figure 13(c) shows the corresponding optimal control profiles. Both $\omega \left(t\right)$ and $\gamma \left(t\right)$ exhibit relatively higher intensities at the early stage of rumor propagation, followed by a gradual decrease as time progresses. This time-dependent pattern is well aligned with the evolution trajectories of $R\left(t\right)$ and $H\left(t\right)$ , reflecting the optimal allocation of control efforts to rapidly curb rumor spreading initially while reducing intervention costs at later stages.

Overall, Figure 13 demonstrates that the combined control strategy outperforms single-control approaches in mitigating both active and dormant rumor dynamics, thereby providing a more effective and balanced intervention mechanism.

Figure 14. Comparison of optimal control profiles under different network settings.

Figure 14 displays the optimal control intensities $\omega \left(t\right)$ and $\gamma \left(t\right)$ for the Facebook, Twitch RU, and Twitch PT networks under the combined control strategy. It can be observed that the temporal profiles of both control functions show highly consistent patterns across the three network configurations. Specifically, the control intensities rise during the initial phase to suppress rapid rumor propagation and then gradually decline as the system approaches a stable state.

While the overall trends show consistency, minor variations can be observed in peak magnitudes, timing of peak occurrences, and initial control levels. These differences primarily stem from the distinct transmission intensity parameter values $\widehat{\kappa }$ adopted for each network. However, the observed variations remain relatively limited, suggesting that the proposed optimal control strategy preserves a stable framework and exhibits robustness against variations in network spreading strength.

7. Conclusions

This paper investigates a fractional-order rumor propagation model of the $S_{\mu }{S}_{e}ERHDU$ type. In addition to the traditional dynamics of rumor spread, the model incorporates three control mechanisms: an education mechanism that enhances public media literacy to reduce susceptibility at the source; a memory-forgetting mechanism that promotes the transition of hibernators into stiflers, thereby blocking their potential for renewed dissemination; and a debunking mechanism that uses authoritative information releases to actively guide public opinion and strengthen societal resilience against rumors.

In the theoretical analysis, we establish the existence, uniqueness, non-negativity, and boundedness of the solutions for the fractional-order system. The basic reproduction number $R_0^{\alpha }$ is defined, and two types of equilibria are examined: the RFE and RSE. It is shown that when $R_0^{\alpha }<1$ , the system is GAS at the RFE, implying that rumors will die out naturally and the system possesses effective self-regulating capacity. When $R_0^{\alpha }>1$ , under certain conditions, the global stability of the RSE is demonstrated, indicating that in the absence of intervention, rumors will persist and continue to spread.

In the optimal control section, educational enhancement and debunking interventions are modeled as time-varying control variables $\gamma \left(t\right)$ and $\omega \left(t\right)$ , leading to a dynamic optimization problem that aims to minimize both the number of spreaders and the control costs. Using Pontryagin’s Maximum Principle, we derive the necessary conditions for the optimal control strategy, which reveal the temporal coordination and functional division of labor between different intervention strategies.

The findings indicate that effective rumor control in practice hinges on rapid early-stage response and adequate resource allocation, supported by a coordinated strategy that simultaneously suppresses active spreading and clears potential risks. This underscores the importance of enhancing public media literacy, establishing authoritative debunking channels, and designing intelligent dynamic intervention schemes to build a multi-layered and adaptive immune system against online rumors.

Moreover, the framework developed in this study integrating fractional-order modeling, stability analysis, and optimal control can be readily extended to other socio-dynamic processes with similar diffusion structures and intervention principle. Examples include the spread of computer viruses in cybersecurity, panic contagion in financial markets, and behavioral influence in public health campaigns. This approach thus offers a versatile theoretical tool for analyzing and managing a wide range of complex diffusion phenomena.

To address the unresolved issues identified in this work, future studies should focus on the following aspects:

1) Incorporation of time-delay effects: incorporating the time delays inherent in real-world scenarios, such as lags associated with educational effects, debunking responses, and individual behavioral changes. This involves developing fractional-order rumor propagation models with discrete or distributed delays and examining how these delays influence system stability, induce Hopf bifurcation, and affect the efficacy of control strategies.

2) Adoption of generalized abstract incidence functions: replacing the current bilinear incidence term with more realistic nonlinear incidence forms—such as saturated incidence, crowded-contact incidence, or psychologically influenced incidence would allow for a more accurate description of propagation dynamics. This extension would enable deeper analysis of changes in the basic reproduction number threshold, equilibrium stability, and the structure of optimal control strategies under nonlinear transmission effects.

3) Expansion to complex network topologies: moving beyond the homogeneous mixing assumption, future studies could integrate complex network structures (e.g., scale-free networks, small-world networks, multi-layered networks). This would allow investigation of how network heterogeneity, degree distributions, and community structure influence rumor dynamics and the effectiveness of intervention measures. Such analyses could support the design of targeted control strategies tailored to heterogeneous or multi-layered social networks.

Acknowledgements

This work was supported by Fundamental Research Funds of China West Normal University(24kc003) and the Research Project on Graduate Education Reform of China West Normal University (2022XM24, 2024XM05).