Modeling Crimean-Congo Hemorrhagic fever with behavioral awareness: Mathematical analysis via Chebyshev spectral collocation solutions

Waleed Adel

doi:10.2478/ijmce-2026-0020

Introduction

CCHF is a severe zoonotic viral disease caused by the CCHF virus, characterized by a high case fatality rate and significant outbreak potential. The virus is primarily transmitted to humans through bites of infected ticks or via direct contact with the blood or tissues of infected animals or humans as illustrated in Figure 1. CCHF is endemic across wide geographical regions including Africa, Asia, Eastern Europe, and the Middle East, and has been designated by the World Health Organization (WHO) as a priority disease requiring urgent research and control efforts [1]. Although several tick species may participate in transmission, ticks of the genus Hyalomma are recognized as the principal vectors responsible for maintaining and spreading the virus [2]. Epidemiological studies indicate that the average mortality rate among infected individuals is approximately 30%, underscoring the severity of the disease [3]. The origin of the disease name reflects its historical discovery. The first documented outbreak occurred in 1942 in the Crimean region of the former Soviet Union. Subsequently, in 1956, the causative virus was isolated from a febrile patient in the Democratic Republic of Congo. The recognition that both events were caused by the same pathogen led to the combined designation CCHF, which remains in use today. Serological investigations have demonstrated that domestic animals such as cattle, sheep, and goats frequently become infected through bites from infected ticks. While these animals typically exhibit mild or transient febrile symptoms, they play a critical role as amplifying hosts that facilitate viral circulation in endemic areas. In addition to livestock, the virus infects a wide range of wild animals. In the European regions of the former Soviet Union, rabbits are considered the primary wildlife reservoirs, whereas in many Asian countries, species such as rabbits, rodents, and hedgehogs serve as major sources of viral maintenance and transmission [4,5].

In recent years, Iraq has emerged as one of the most severely affected countries, experiencing recurrent CCHF outbreaks with a marked increase in reported cases across several provinces. This trend highlights the urgent need for strengthened surveillance systems and effective prevention and control strategies. The weekly incidence patterns and spatial distribution of reported cases, as depicted in Figures 2 and 3, demonstrate a clear temporal escalation and pronounced geographical heterogeneity. These patterns reflect the combined influence of ecological conditions, human behavior, livestock management practices, and socio-economic factors on disease transmission [6]. The transmission dynamics of CCHF are inherently complex due to the involvement of multiple interacting populations, including ticks as vectors, livestock as amplifying hosts, and humans as incidental hosts. In endemic settings such as Iraq, most human infections are associated with occupational or environmental exposure to ticks and contact with infected livestock, particularly among agricultural workers and individuals engaged in animal husbandry. Human-to-human transmission occurs less frequently and is generally confined to close-contact environments, including households and healthcare facilities. These epidemiological features strongly motivate the development of multi-compartment mathematical models that explicitly incorporate vector-host interactions and zoonotic spillover processes. Furthermore, growing evidence suggests that public awareness, behavioral changes, and preventive practices substantially influence outbreak dynamics, yet these factors are often inadequately represented in conventional epidemic modeling frameworks.

Mathematical modeling has become an indispensable tool for analyzing the transmission dynamics of infectious diseases, assessing intervention strategies, and supporting evidence-based public health decision-making. Classical compartmental models formulated as systems of nonlinear differential equations have been extensively employed to describe epidemic processes. Nevertheless, accurately representing the nonlinear feedback mechanisms between epidemiological states and human behavior remains a significant challenge. In recent years, substantial progress has been made in the numerical simulation and mathematical analysis of epidemic models, enabling the investigation of a wide range of infectious diseases and complex transmission scenarios. These advances have facilitated the modeling of several important diseases, including multiple formulations for COVID-19 [7], anthrax transmission in animal populations [8], Hepatitis C dynamics [9], infections influenced by environmental persistence [10], mumps virus spread [11], Zika virus transmission [12], mosaic disease dynamics [13], canine distemper virus outbreaks [14], Q fever epidemiology [15], and the coupled co-dynamics of COVID-19 and diabetes [16]. Such studies highlight the flexibility of mathematical frameworks in capturing diverse biological, environmental, and behavioral mechanisms underlying disease propagation. Within the context of CCHF, incorporating awareness-driven behavioral changes into transmission models offers valuable insights into the effectiveness of non-pharmaceutical interventions, including public education programs, personal protective measures, and risk communication strategies. Several recent studies have proposed different compartmental and nonlocal modeling approaches to better understand the transmission dynamics of CCHF. For instance, Sina et al. [17] investigated a fractional-order CCHF model incorporating power-law kernels to capture memory effects in disease transmission. Karrar et al. [18] analyzed a delayed CCHF model in human populations and examined the influence of time delays on system stability and disease persistence. Suman et al. [19] developed a compartmental framework that explicitly accounts for the role of blood-sucking ticks in the transmission cycle. Furthermore, Hakimen et al. [20] employed nonlocal fractional derivatives to describe the full dynamical behavior of the disease, and Juan et al. [21] utilized real-world epidemiological data to construct a multi-compartment model providing a data-driven perspective on CCHF transmission dynamics in Afghanistan. Despite these important contributions, none of the existing CCHF models incorporate public awareness or behavioral response mechanisms, treating human susceptibility as a fixed biological quantity unaffected by risk perception, preventive behavior, or public health communication. This represents a significant gap, given the well-documented role of behavioral changes in shaping outbreak dynamics in endemic settings such as Iraq. More broadly, awareness-based epidemic modeling has received considerable attention in the context of directly transmitted diseases. Classical formulations, such as those reviewed in [22,23], typically incorporate awareness as a simple multiplicative reduction uniformly applied to the transmission rate of a single-host system. In vector-borne disease models, awareness is often represented as a static reduction in the human-vector contact rate, without accounting for the dynamics of awareness generation, decay, or its selective effect on specific transmission pathways. In contrast, the present work introduces a dynamically evolving awareness variable A(t) governed by its own differential equation, driven by the real-time prevalence of infection in both humans and ticks, and subject to natural decay capturing behavioral fatigue. Critically, this awareness variable modulates only the human force of infection through the saturation factor $\frac{1}{1 + k_{A} A (t)}$ n, leaving the enzootic tick-livestock cycle unaffected thus a structurally appropriate and biologically realistic design for a zoonotic disease in which the animal reservoir cycle operates independently of human behavioral responses. To our knowledge, no prior work has combined this awareness architecture with a full three-population tick-livestock-human zoonotic transmission model, nor solved the resulting nonlinear system using a high-order spectral collocation framework. Solving large-scale nonlinear epidemic models over extended time horizons poses substantial numerical challenges, particularly in the presence of strong nonlinear interactions, feedback mechanisms, and multi-compartment coupling. Classical low-order time-stepping schemes often suffer from severe stability restrictions, slow convergence, or significant loss of accuracy when applied to such systems, especially in the context of long-term simulations. Spectral and collocation methods have gained increasing attention owing to their ability to achieve high-order accuracy and rapid often spectral or near-exponential convergence for sufficiently smooth solutions. The use of polynomial bases for the approximation of solutions to various mathematical models has been studied extensively in recent years. For instance, Selcuk et al. [24] developed a robust septic Hermite collocation technique for the heat conduction equation, demonstrating the effectiveness of high-degree Hermite-based discretizations for parabolic problems. Subsequently, the same authors investigated a cubic Hermite B-spline collocation method [25], further establishing the versatility of Hermite-type approximations across different regularity regimes. In both approaches, the state variables are represented by global polynomial expansions, and the governing equations are enforced at carefully selected collocation nodes such as Chebyshev-Gauss-Lobatto points yielding highly accurate differentiation matrices and compact representations of the solution. To handle nonlinearities efficiently, the QLM provides a systematic and robust strategy by transforming the original nonlinear system into a sequence of linear subproblems through successive linearization about the current iterate. This Newton-type iterative procedure preserves the essential structure of the original model while producing rapidly convergent solution sequences under mild regularity conditions. When combined with spectral collocation discretizations, quasilinearization yields a computational framework that maintains high accuracy while significantly enhancing numerical stability and efficiency. Furthermore, the incorporation of domain decomposition techniques reinforces the robustness of the overall scheme by partitioning the computational interval into smaller, overlapping or non-overlapping subdomains. This partitioning enables accurate resolution of long-term dynamics, mitigates error accumulation associated with extended integration horizons, and reduces sensitivity to stiffness without compromising computational efficiency. The resulting integrated quasilinearized spectral collocation framework is therefore particularly well suited for simulating complex epidemic models involving multiple interacting compartments and long-term behavioral dynamics. Motivated by these considerations, the present study develops an awareness-driven CCHF transmission model tailored to the epidemiological and socio-behavioral context of Iraq and proposes a high accuracy numerical solution strategy based on QLM coupled with Chebyshev spectral collocation. The proposed modeling framework explicitly incorporates awareness induced behavioral responses into the transmission dynamics, allowing for a more realistic representation of how risk perception, preventive practices, and public health communication influence disease spread. From a computational perspective, the Chebyshev polynomial basis is selected for the spectral collocation scheme due to its distinctive theoretical and practical advantages over alternative orthogonal bases. In particular, the Chebyshev-Gauss-Lobatto collocation nodes admit the explicit closed-form expression $x_{i} = \cos (i π / L), i = 0, 1, \dots, L$ i enabling efficient computation via the Fast Fourier Transform and avoiding the nonlinear eigenvalue problems required to generate Legendre-Gauss-Lobatto nodes. Moreover, Chebyshev polynomials satisfy the equioscillation property [26, 27], which provides near-optimal uniform approximation and effectively suppresses the Runge phenomenon over long simulation horizons. Also, the associated differentiation matrix possesses a fully explicit and numerically stable structure, and the resulting spectral approximation achieves the optimal convergence rate of $𝒢 (L^{- m})$ j in the weighted L₂ norm for functions in the Sobolev class $H_{w}^{m}$ w_i, matching the theoretical performance of Legendre-based schemes while offering superior computational convenience. Other polynomials such as Legendre polynomials which are orthogonal with respect to the uniform weight and are equally rigorous from a theoretical standpoint represent a viable alternative. However, they lack closed-form node formulas, making Chebyshev polynomials the preferred choice for the high-order framework developed here. The integration of quasilinearization with this Chebyshev spectral collocation yields a robust and efficient numerical scheme capable of resolving strong nonlinearities and long-term dynamics with high accuracy and stability. The framework facilitates a detailed qualitative and quantitative analysis of the model, enables reliable numerical approximations validated against reference solutions, and supports systematic investigation of key epidemiological and behavioral parameters. The numerical results demonstrate the capability of the proposed approach to accurately capture the coupled epidemiological behavioral dynamics of CCHF and highlight the critical role of awareness-based interventions in mitigating outbreaks, particularly in high-risk and resource-limited regions such as Iraq. The novelty of the paper lies within the following points:

We develop an awareness-driven multi-compartment mathematical model to describe the transmission dynamics of CCHF through coupled tick-livestock-human interactions. Unlike existing CCHF models, which treat human susceptibility as fixed and do not account for behavioral responses, the proposed model explicitly integrates a dynamic awareness variable into the transmission process, providing a more realistic representation of non-pharmaceutical interventions in high-risk endemic settings.
In contrast to awareness-based models for directly transmitted or single-host vector-borne diseases, where awareness uniformly reduces the overall transmission rate, the proposed formulation selectively modulates the human force of infection through the saturation factor $\frac{1}{1 + k_{A} A (t)}$ x_i, while preserving the independent enzootic tick-livestock cycle. This distinction reflects the biological reality that human behavioral changes cannot interrupt zoonotic viral circulation among animal reservoirs.
Awareness evolves as a dynamic state variable governed by its own differential equation, driven by the prevalence of infection in both humans and ticks at rates η₁ and η₂ respectively, and subject to natural decay at rate ω_A capturing behavioral fatigue a mathematically richer and more epidemiologically grounded representation than static awareness reductions used in prior works.
To accurately resolve the resulting nonlinear system over long time horizons, we propose a high-order numerical framework coupling the QLM with Chebyshev spectral collocation and domain decomposition a computational approach not previously applied to awareness-driven zoonotic disease models.
We rigorously establish the well-posedness of the model by proving positivity and boundedness of solutions, deriving the disease-free equilibrium, and obtaining the basic reproduction number ℛ₀ via the next-generation matrix approach, with local stability shown to be governed by the reproduction threshold and the sensitivity to the significant parameters.
Numerical simulations validate the theoretical analysis and demonstrate the effectiveness of the proposed spectral scheme, while providing quantitative insights into how awareness-driven behavioral responses reduce infection prevalence and outbreak magnitude even when ℛ₀ > 1.
The proposed framework is sufficiently general to be extended to other zoonotic diseases involving vector host interactions and behavior-mediated transmission, offering a reusable modeling and computational template for endemic disease analysis.

The rest of the paper is organized as follows. Section 2 presents the formulation of the awareness-driven CCHF transmission model and its underlying assumptions. Section 3 provides the qualitative analysis of the model, including positivity, boundedness, disease-free equilibrium, the basic reproduction number, and local stability results. Section 4 introduces the Chebyshev spectral framework and establishes rigorous convergence results. Section 5 develops the proposed quasilinearization-based Chebyshev collocation method with domain decomposition. Section 6 presents numerical simulations illustrating the effectiveness of the proposed approach. Finally, Section 7 concludes the paper and outlines future research directions.

Model formulation

In this section, we develop a deterministic compartmental model that describes the transmission dynamics of CCHF through the coupled interactions among ticks, livestock, and humans. The model explicitly incorporates the zoonotic nature of CCHF, vector-mediated transmission, spillover from animals to humans, limited human-to-human transmission, and the impact of public awareness and behavioral responses on disease spread.

The total system consists of three interacting populations: ticks, livestock, and humans. Each population is subdivided into epidemiologically relevant compartments as follows: susceptible ticks S_T(t) and infected ticks I_T(t). Livestock are classified into susceptible livestock S_L(t), exposed (latent) livestock E_L(t), and infectious livestock I_L(t). The human population is divided into susceptible humans S_H(t), exposed humans E_H(t), infectious humans I_H(t), recovered humans R_H(t), and awareness variable A(t) which represent the level of public awareness and behavioral response toward CCHF prevention where recovered individuals acquire temporary immunity and may return to the susceptible class due to immunity waning. Recruitment into the tick, livestock, and human populations occurs at constant rates Λ_T, Λ_L, and Λ_H, respectively. All populations experience natural mortality at rates μ_T, μ_L, and μ_H. Infected humans may also experience disease-induced mortality at rate δ_H. The CCHF transmission occurs through multiple pathways reflecting the ecological and epidemiological nature of the disease. Susceptible ticks become infected through contact with infectious livestock. Susceptible livestock acquire infection primarily through bites from infected ticks. Humans may become infected via contact with infected ticks, infected livestock, or infectious humans, particularly through exposure to blood or bodily fluids in healthcare or household settings. Let the populations of each category be as follows

N_{T} (t) = S_{T} (t) + I_{T} (t), N_{L} (t) = S_{L} (t) + E_{L} (t) + I_{L} (t), N_{H} (t) = S_{H} (t) + E_{H} (t) + I_{H} (t) + R_{H} (t) .

N_T(t) = S_T(t) + I_T(t), \quad N_L(t) = S_L(t) + E_L(t) + I_L(t), \quad N_H(t) = S_H(t) + E_H(t) + I_H(t) + R_H(t)

denote the total tick, livestock, and human populations, respectively. Under the model assumptions, the total populations remain uniformly bounded for all t > 0. The forces of infection can be defined as the tick infection force is

λ_{T} (t) = β_{L T} \frac{I_{L} (t)}{N_{L} (t)} .

\lambda_T(t) = \beta_{LT} \frac{I_L(t)}{N_L(t)}

and the livestock infection force is

λ_{L} (t) = β_{T L} \frac{I_{T} (t)}{N_{T} (t)} .

\lambda_L(t) = \beta_{TL} \frac{I_T(t)}{N_T(t)}

and the human infection force is

λ_{H} (t) = (β_{T H} \frac{I_{T} (t)}{N_{T} (t)} + β_{L H} \frac{I_{L} (t)}{N_{L} (t)} + β_{H H} \frac{I_{H} (t)}{N_{H} (t)}) \frac{1}{1 + k_{A} A (t)} .

\lambda_H(t) = \left( \beta_{TH} \frac{I_T(t)}{N_T(t)} + \beta_{LH} \frac{I_L(t)}{N_L(t)} + \beta_{HH} \frac{I_H(t)}{N_H(t)} \right) \frac{1}{1 + k_A A(t)}

The human force of infection is modulated by the awareness factor $\frac{1}{1 + k_{A} A (t)}$ f(x_i), which takes the form of a Michaelis-Menten saturation function widely adopted in behavioral epidemiology to capture the diminishing marginal effect of awareness on risk reduction [23,28]. The biological rationale is that at low awareness levels each unit increase in A(t) produces a substantial reduction in risky contacts, whereas at high awareness levels behavioral change saturates, since individuals can only reduce their exposure to a finite minimum regardless of further increases in perceived risk. The parameter k_A > 0 controls the strength of this behavioral response: larger values correspond to stronger modification per unit of awareness, while k_A → 0 recovers the awareness-free model. This functional form is preferred over a linear reduction ( 1 − αA ), which risks producing negative transmission rates for large A, and over an exponential form e^−αA, which is less tractable and harder to interpret epidemiologically. The modulation factor is applied to the full human force of infection, reflecting the fact that awareness-driven protective behaviors recommended by WHO for CCHF including use of protective clothing, insect repellent, avoidance of livestock contact, and improved hygiene act simultaneously across all exposure pathways rather than selectively targeting a single transmission route. Awareness is modeled as a prevalence-driven process rather than an absolute-count process, ensuring scale invariance with respect to population size. The awareness variable A(t) increases proportionally to the prevalence of infection in humans at rate η₁, reflecting media coverage and case-driven risk communication, and to the prevalence of infected ticks at rate η₂, reflecting the active dissemination of tick surveillance data and tick-borne risk warnings by public health authorities in endemic settings such as Iraq [29,30]. Awareness decays naturally at rate ω_A, capturing the well-documented phenomenon of behavioral fatigue, whereby populations tend to relax protective measures as perceived risk diminishes over time. It can be shown that for nonnegative initial conditions all state variables remain nonnegative for all t > 0, ensuring biological feasibility of the model. The transmission terms are formulated using frequency-dependent (proportional) incidence, whereby the force of infection scales with the prevalence I/N rather than the absolute infected count I. This formulation is appropriate when contact rates are approximately independent of population density a reasonable assumption for human and livestock populations within a fixed endemic region with stable social and agricultural contact patterns. The proposed model is developed under the following assumptions:

The populations of ticks, livestock, and humans are homogeneously mixed within each group, so that contact rates depend only on compartment proportions rather than spatial location or individual heterogeneity.
Transmission between populations occurs through effective contact rates formulated as frequency-dependent (proportional) incidence, appropriate for populations with density-independent contact behavior.
Demographic processes (recruitment and natural mortality) occur on a slower time scale compared to disease transmission, so that total population sizes remain approximately constant over the outbreak period.
Recovered humans acquire temporary immunity that wanes at rate ω_R > 0, consistent with clinical evidence indicating the absence of durable long-term immunity in CCHF survivors.
Public awareness reduces effective human exposure through the saturation factor $\frac{1}{1 + k_{A} A (t)}$ A_i applied to the human force of infection, without directly altering biological transmission parameters or affecting the tick-livestock enzootic cycle, which operates independently of human behavioral responses.

Based on the above assumptions, the dynamics of the CCHF transmission model are governed by the following system:

\begin{array}{l} \frac{d S_{L}}{d t} = Λ_{L} - λ_{L} S_{L} - μ_{L} S_{L}, \\ \frac{d E_{L}}{d t} = λ_{L} S_{L} - (σ_{L} + μ_{L}) E_{L}, \\ \frac{d I_{L}}{d t} = σ_{L} E_{L} - (γ_{L} + μ_{L}) I_{L}, \\ \frac{d S_{H}}{d t} = Λ_{H} + ω_{R} H_{} - λ_{h} S_{H} - μ_{H} S_{H}, \\ \frac{d E_{H}}{d t} = λ_{H} S_{H} - (σ_{H} + μ_{H}) E_{H}, \\ \frac{d I_{H}}{d t} = σ_{H} E_{H} - (γ_{H} + μ_{H} + δ_{H}) I_{H}, \\ \frac{d R_{H}}{d t} = γ_{H} I_{H} - (ω_{R} + μ_{H}) R_{H}, \\ \frac{d A}{d t} = η_{1} \frac{I_{H}}{N_{H}} + η_{2} \frac{I_{T}}{N_{T}} - ω_{A} A . \end{array}

\\\begin{array}{ccc}\frac{\displaystyle dS_L}{\displaystyle dt}&=&\Lambda_L-\lambda_LS_L-\mu_LS_L,\\\frac{\displaystyle dE_L}{\displaystyle dt}&=&\lambda_LS_L-{(\sigma_L+\mu_L)}E_L,\\\frac{\displaystyle dI_L}{\displaystyle dt}&=&\sigma_LE_L-{(\gamma_L+\mu_L)}I_{L,}\\\frac{\displaystyle dS_H}{\displaystyle dt}&=&\Lambda_H+\omega_RR_H-\lambda_HS_H-\mu_HS_H,\end{array}

The initial conditions of system are assumed to satisfy nonnegative solutions as follows

S_{i} (0) \geq 0, E_{i} (0) \geq 0, I_{i} (0) \geq 0, R_{H} (0) \geq 0, A (0) \geq 0,

S_i(0) \geq 0, \quad E_i(0) \geq 0, \quad I_i(0) \geq 0, \quad R_H(0) \geq 0, \quad A(0) \geq 0

for all relevant compartments, ensuring biological feasibility of the solutions. The definitions of the state variables and model parameters are summarized in Tables 1 and 2, respectively.

Table 1

Description of the state variables in the CCHF transmission model.

State variable	Description
S_T(t)	Susceptible tick population.
I_T(t)	Infected tick population capable of transmitting CCHF.
S_L(t)	Susceptible livestock population.
E_L(t)	Exposed livestock population in the latent stage.
I_L(t)	Infectious livestock population.
S_H(t)	Susceptible human population.
E_H(t)	Exposed human population during the incubation period.
I_H(t)	Infectious human population.
R_H(t)	Recovered human population with temporary immunity.
A(t)	Level of public awareness and behavioral response.

Table 2

Description of the parameters used in the CCHF transmission model.

Parameter	Description
Λ_T, Λ_L, Λ_H	Recruitment rates of ticks, livestock, and humans.
μ_T, μ_L, μ_H	Natural mortality rates.
β_TL	Transmission rate from infected ticks to livestock.
β_LT	Transmission rate from infected livestock to ticks.
β_TH	Transmission rate from infected ticks to humans.
β_LH	Transmission rate from infected livestock to humans.
β_HH	Human-to-human transmission rate.
σ_L, σ_H	Progression rates from exposed to infectious classes.
γ_L, γ_H	Recovery rates of livestock and humans.
δ_H	Disease-induced mortality rate in humans.
ω_R	Rate of loss of immunity in recovered humans.
η₁, η₂	Awareness generation rates.
ω_A	Natural decay rate of awareness.
k_A	Strength of awareness-induced behavioral response.

Qualitative analysis of model

3.1

Positivity of solutions

We prove that system (1) is epidemiologically well-posed in the sense that all state variables remain nonnegative for all future time whenever they start nonnegative.

Theorem 1.

Let

$x (t) = {(S_{T} (t), I_{T} (t), S_{L} (t), E_{L} (t), I_{L} (t), S_{H} (t), E_{H} (t), I_{H} (t), R_{H} (t), A (t))}^{⊤}$ \mathbf{x}(t) = (S_T(t), I_T(t), S_L(t), E_L(t), I_L(t), S_H(t), E_H(t), I_H(t), R_H(t), A(t))^\top

be the solution of system (1) with initial condition $x (0) \in ℝ_{+}^{10}$ B_i. Assume all parameters are nonnegative and that the forces of infection λ_T, λ_L, λ_H are nonnegative whenever the state variables are nonnegative. Then $x (t) \in ℝ_{+}^{10}$ C_i for all t ≥ 0.

Proof.

The right-hand side of system (1) is locally Lipschitz in $ℝ_{+}^{10}$ D_i; hence, a unique local solution exists for any nonnegative initial condition. To show positivity, it suffices to verify that the vector field is inward pointing on the boundary of the nonnegative orthant, i.e., whenever a compartment equals zero while the other compartments are nonnegative, its derivative is nonnegative as proved as follows:

$\begin{array}{l} {\frac{d S_{T}}{d t} |}_{S_{T} = 0} & = Λ_{T} \geq 0, \\ {\frac{d I_{T}}{d t} |}_{I_{T} = 0} & = λ_{T} S_{T} \geq 0, \\ {\frac{d S_{L}}{d t} |}_{S_{L} = 0} & = Λ_{L} \geq 0, \\ {\frac{d E_{L}}{d t} |}_{E_{L} = 0} & = λ_{L} S_{L} \geq 0, \\ {\frac{d I_{L}}{d t} |}_{L} = 0 & = σ_{L} E_{L} \geq 0, \\ {\frac{d S_{H}}{d t} |}_{S_{H} = 0} & = Λ_{H} + ω_{R} R_{H} \geq 0, \\ {\frac{d E_{H}}{d t} |}_{E_{H} = 0} & = λ_{H} S_{H} \geq 0, \\ {\frac{d I_{H}}{d t} |}_{I_{H} = 0} & = σ_{H} E_{H} \geq 0, \\ {\frac{d R_{H}}{d t} |}_{R_{H} = 0} & = γ_{H} I_{H} \geq 0, \\ {\frac{d A}{d t} |}_{A = 0} & = η_{1} \frac{I_{H}}{N_{H}} + η_{2} \frac{I_{T}}{N_{T}} \geq 0 . \end{array}$ \Omega = \left\{ (S_T, I_T, S_L, E_L, I_L, S_H, E_H, I_H, R_H, A) \in \mathbb{R}_+^{10} : N_T \leq \frac{\Lambda_T}{\mu_T}, N_L \leq \frac{\Lambda_L}{\mu_L}, N_H \leq \frac{\Lambda_H}{\mu_H}, A \leq \frac{\eta_1 + \eta_2}{\omega_A} \right\}

provided N_H > 0 and N_T > 0. Therefore, on every boundary face of $ℝ_{+}^{10}$ E_i, the corresponding derivative is nonnegative, and the solution cannot cross into the negative orthant. Hence $x (t) \in ℝ_{+}^{10}$ F_i for all $t \geq 0$ G_i.

3.2

Boundedness of solutions

We next show that solutions of system (1) are uniformly bounded in a positively invariant region.

Theorem 2.

Letx(t) be a solution of system (1) with $x (0) \in ℝ_{+}^{10}$ H_i. Define the total populations

$\begin{matrix} N_{T} (t) = S_{T} (t) + I_{T} (t), N_{L} (t) = S_{L} (t) + E_{L} (t) + I_{L} (t), \\ N_{H} (t) = S_{H} (t) + E_{H} (t) + I_{H} (t) + R_{H} (t) . \end{matrix}$ \begin{array}{cc}N_T(t)=S_T(t)+I_T(t),\;\;N_L(t)=S_L(t)+E_L(t)+I_L(t),&\\\;\;N_H(t)=S_H(t)+E_H(t)+I_H(t)+R_H(t).&{}\end{array}

Assume that μ_T, μ_L, μ_H > 0 and ω_A > 0. Then, for all t ≥ 0,

$\begin{matrix} 0 \leq N_{T} (t) \leq \max {N_{T} (0), \frac{Λ_{T}}{μ_{T}}}, & 0 \leq N_{L} (t) \leq \max {N_{L} (0), \frac{Λ_{L}}{μ_{L}}}, \\ 0 \leq N_{H} (t) \leq \max {N_{H} (0), \frac{Λ_{H}}{μ_{H}}}, \end{matrix}$ \begin{array}{cc}0\leq N_T(t)\leq{\max\{N_T(0),\frac{\Lambda_T}{\mu_T}\}},&0\leq N_L(t)\leq{\max\{N_L(0),\frac{\Lambda_L}{\mu_L}\}},\\\;\;\;\;\;0\leq N_H(t)\leq{\max\{N_H(0),\frac{\Lambda_H}{\mu_H}\}},&{}\end{array}

and the awareness variable satisfies the bound

$0 \leq A (t) \leq max {A (0), \frac{η_{1} + η_{2}}{ω_{A}}} .$ \mathcal{V}(\mathbf{x}) = \begin{pmatrix} \mu_T I_T \\ (\sigma_L + \mu_L) E_L \\ -\sigma_L E_L + (\gamma_L + \mu_L) I_L \\ (\sigma_H + \mu_H) E_H \\ -\sigma_H E_H + (\gamma_H + \mu_H + \delta_H) I_H \\ -\eta_1 \frac{I_H}{N_H} - \eta_2 \frac{I_T}{N_T} + \omega_A A \end{pmatrix}

Consequently, the set

$Ω = {x \in ℝ_{+}^{10} : N_{T} \leq \frac{Λ_{T}}{μ_{T}}, N_{L} \leq \frac{Λ_{L}}{μ_{L}}, N_{H} \leq \frac{Λ_{H}}{μ_{H}}, A \leq \frac{η_{1} + η_{2}}{ω_{A}}}$ F = \begin{pmatrix} 0 & 0 & \beta_{LT} \frac{\Lambda_T \mu_L}{\Lambda_L \mu_T} & 0 & 0 & 0 \\ \beta_{TL} \frac{\Lambda_L \mu_T}{\Lambda_T \mu_L} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ \beta_{TH} \frac{\Lambda_H \mu_T}{\Lambda_T \mu_H} & 0 & \beta_{LH} \frac{\Lambda_H \mu_L}{\Lambda_L \mu_H} & 0 & \beta_{HH} & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \end{pmatrix}

is positively invariant and attracting for system (1). Positive invariance means that any trajectory of system (1) initiating inside Ω remains in Ω for all future time, i.e.,

$x (0) \in Ω ⟹ x (t) \in Ω for all t \geq 0 .$ V = \begin{pmatrix} \mu_T & 0 & 0 & 0 & 0 & 0 \\ 0 & \sigma_L + \mu_L & 0 & 0 & 0 & 0 \\ 0 & -\sigma_L & \gamma_L + \mu_L & 0 & 0 & 0 \\ 0 & 0 & 0 & \sigma_H + \mu_H & 0 & 0 \\ 0 & 0 & 0 & -\sigma_H & \gamma_H + \mu_H + \delta_H & 0 \\ -\eta_2 \frac{\mu_T}{\Lambda_T} & 0 & 0 & 0 & -\eta_1 \frac{\mu_H}{\Lambda_H} & \omega_A \end{pmatrix}

while attracting means that every trajectory of system (1) originating outside Ω, but within the nonnegative orthant $ℝ_{+}^{10}$ I_i, eventually enters Ω in finite time and remains therein for all subsequent time. More precisely, for any initial condition $x (0) \in ℝ_{+}^{10}$ J_i, there exists a finite time t^∗ ≥ 0 such that

$x (t) \in Ω for all t \geq t^{*} .$ V^{-1} = \begin{pmatrix} \frac{1}{\mu_T} & 0 & 0 & 0 & 0 & 0 \\ 0 & \frac{1}{\sigma_L + \mu_L} & 0 & 0 & 0 & 0 \\ 0 & \frac{\sigma_L}{(\sigma_L + \mu_L)(\gamma_L + \mu_L)} & \frac{1}{\gamma_L + \mu_L} & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1}{\sigma_H + \mu_H} & 0 & 0 \\ 0 & 0 & 0 & \frac{\sigma_H}{(\sigma_H + \mu_H)(\gamma_H + \mu_H + \delta_H)} & \frac{1}{\gamma_H + \mu_H + \delta_H} & 0 \\ \frac{\eta_2}{\omega_A \Lambda_T} & 0 & 0 & \frac{\eta_1 \sigma_H}{\omega_A \Lambda_H (\sigma_H + \mu_H)(\gamma_H + \mu_H + \delta_H)} & \frac{\eta_1}{\omega_A \Lambda_H (\gamma_H + \mu_H + \delta_H)} & \frac{1}{\omega_A} \end{pmatrix}

Proof.

Summing the first two equations in (1) gives

$\frac{d N_{T}}{d t} = \frac{d S_{T}}{d t} + \frac{d I_{T}}{d t} = Λ_{T} - μ_{T} (S_{T} + I_{T}) = Λ_{T} - μ_{T} N_{T} .$ FV^{-1} = \begin{pmatrix} 0 & \frac{\beta_{LT} \Lambda_T \mu_L \sigma_L}{\Lambda_L \mu_T (\sigma_L + \mu_L)(\gamma_L + \mu_L)} & \frac{\beta_{LT} \Lambda_T \mu_L}{\Lambda_L \mu_T (\gamma_L + \mu_L)} & 0 & 0 & 0 \\ \frac{\beta_{TL} \Lambda_L}{\Lambda_T \mu_L} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ \frac{\beta_{TH} \Lambda_H}{\Lambda_T \mu_H} & \frac{\beta_{LH} \Lambda_H \mu_L \sigma_L}{\Lambda_L \mu_H (\sigma_L + \mu_L)(\gamma_L + \mu_L)} & \frac{\beta_{LH} \Lambda_H \mu_L}{\Lambda_L \mu_H (\gamma_L + \mu_L)} & \frac{\beta_{HH} \sigma_H}{(\sigma_H + \mu_H)(\gamma_H + \mu_H + \delta_H)} & \frac{\beta_{HH}}{\gamma_H + \mu_H + \delta_H} & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \end{pmatrix}

Solving this linear equation yields

$N_{T} (t) = N_{T} (0) e^{- μ_{T} t} + \frac{Λ_{T}}{μ_{T}} (1 - e^{- μ_{T} t}),$ \mathcal{R}_0 = \max \left\{ \frac{\beta_{HH} \sigma_H}{(\sigma_H + \mu_H)(\gamma_H + \mu_H + \delta_H)}, \sqrt{\frac{\beta_{LT} \beta_{TL} \sigma_L}{\mu_T (\sigma_L + \mu_L)(\gamma_L + \mu_L)}} \right\}

and hence 0 ≤ N_T(t) ≤ max{N_T(0), Λ_T/μ_T} for all t ≥ 0. Similarly, summing the third, fourth, and fifth equations in (1) gives

$\frac{d N_{L}}{d t} = Λ_{L} - μ_{L} (S_{L} + E_{L} + I_{L}) - γ_{L} I_{L} \leq Λ_{L} - μ_{L} N_{L},$ \mathcal{R}_{01} = \frac{\beta_{HH} \sigma_H}{(\sigma_H + \mu_H)(\gamma_H + \mu_H + \delta_H)}.

which implies by comparison that

$0 \leq N_{L} (t) \leq max {N_{L} (0), \frac{Λ_{L}}{μ_{L}}}, t \geq 0 .$ \mathcal{R}_{02} = \sqrt{\frac{\beta_{LT} \beta_{TL} \sigma_L}{\mu_T (\sigma_L + \mu_L)(\gamma_L + \mu_L)}}.

Summing from the sixth to the ninth equations in (1) yields cancellation of the internal transfer terms (including ω_RR_H) and gives

$\frac{d N_{H}}{d t} = Λ_{H} - μ_{H} (S_{H} + E_{H} + I_{H} + R_{H}) - δ_{H} I_{H} \leq Λ_{H} - μ_{H} N_{H} .$ \Upsilon_p^{\mathcal{R}_0} = \frac{\partial \mathcal{R}_0}{\partial p} \times \frac{p}{\mathcal{R}_0}

Therefore,

$0 \leq N_{H} (t) \leq max {N_{H} (0), \frac{Λ_{H}}{μ_{H}}}, t \geq 0 .$ \Upsilon_{\beta_{LT}}^{\mathcal{R}_{02}} = \frac{\partial \mathcal{R}_{02}}{\partial \beta_{LT}} \times \frac{\beta_{LT}}{\mathcal{R}_{02}} = \frac{1}{2}> 0

Finally, using 0 ≤ I_H/N_H ≤ 1 and 0 ≤ I_T/N_T ≤ 1 (whenever N_H > 0 and N_T > 0 ), we obtain

$\frac{d A}{d t} = η_{1} \frac{I_{H}}{N_{H}} + η_{2} \frac{I_{T}}{N_{T}} - ω_{A} A \leq (η_{1} + η_{2}) - ω_{A} A .$ \Upsilon_{\beta_{TL}}^{\mathcal{R}_{02}} = \frac{\partial \mathcal{R}_{02}}{\partial \beta_{TL}} \times \frac{\beta_{TL}}{\mathcal{R}_{02}} = \frac{1}{2}> 0

By comparison with ẏ = (η₁ + η₂) − wω_Ay, it follows that

$0 \leq A (t) \leq A (0) e^{- ω_{A} t} + \frac{η_{1} + η_{2}}{ω_{A}} (1 - e^{- ω_{A} t}),$ \Upsilon_{\sigma_L}^{\mathcal{R}_{02}} = \frac{\partial \mathcal{R}_{02}}{\partial \sigma_L} \times \frac{\sigma_L}{\mathcal{R}_{02}} = \frac{\mu_L}{2(\sigma_L + \mu_L)}> 0

and hence 0 ≤ A(t) ≤ max{A(0),(η₁ + η₂)/ω_A} for all t ≥ 0.

Combining these bounds shows that every trajectory starting in $ℝ_{+}^{10}$ K_i enters and remains in Ω, which is therefore positively invariant and attracting.

3.3

Disease-free equilibrium (DFE)

The disease-free equilibrium (DFE) of system (1) corresponds to the state where no infection persists in any population (ticks, livestock, or humans) and where public awareness is absent because there is no perceived risk. Setting the infected and latent compartments to zero, namely

I_{T} = E_{L} = I_{L} = E_{H} = I_{H} = R_{H} = A = 0,

\Upsilon_{\mu_T}^{\mathcal{R}_{02}} = \frac{\partial \mathcal{R}_{02}}{\partial \mu_T} \times \frac{\mu_T}{\mathcal{R}_{02}} = -\frac{1}{2} < 0

and imposing steady-state conditions on the remaining susceptible classes in (1), we obtain

0 = Λ_{T} - μ_{T} S_{T}, 0 = Λ_{L} - μ_{L} S_{L}, 0 = Λ_{H} - μ_{H} S_{H} .

\Upsilon_{\mu_L}^{\mathcal{R}_{02}} = \frac{\partial \mathcal{R}_{02}}{\partial \mu_L} \times \frac{\mu_L}{\mathcal{R}_{02}} = -\frac{\mu_L (\sigma_L + \gamma_L + 2\mu_L)}{2(\sigma_L + \mu_L)(\gamma_L + \mu_L)} < 0.

Hence,

S_{T}^{0} = \frac{Λ_{T}}{μ_{T}}, S_{L}^{0} = \frac{Λ_{L}}{μ_{L}}, S_{H}^{0} = \frac{Λ_{H}}{μ_{H}},

\Upsilon_{\gamma_L}^{\mathcal{R}_{02}} = \frac{\partial \mathcal{R}_{02}}{\partial \gamma_L} \times \frac{\gamma_L}{\mathcal{R}_{02}} = -\frac{\gamma_L}{2(\gamma_L + \mu_L)} < 0

and the disease-free equilibrium is given by

ℰ_{0} = (S_{T}^{0}, I_{T}^{0}, S_{L}^{0}, E_{L}^{0}, I_{L}^{0}, S_{H}^{0}, E_{H}^{0}, I_{H}^{0}, R_{H}^{0}, A^{0}) = (\frac{Λ_{T}}{μ_{T}}, 0, \frac{Λ_{L}}{μ_{L}}, 0, 0, \frac{Λ_{H}}{μ_{H}}, 0, 0, 0, 0) .

\Upsilon_{\sigma_H}^{\mathcal{R}_{01}} = \frac{\partial \mathcal{R}_{01}}{\partial \sigma_H} \times \frac{\sigma_H}{\mathcal{R}_{01}} = \frac{\mu_H}{\sigma_H + \mu_H}> 0

Moreover, at the DFE the total populations satisfy

N_{T}^{0} = \frac{Λ_{T}}{μ_{T}}, N_{L}^{0} = \frac{Λ_{L}}{μ_{L}}, N_{H}^{0} = \frac{Λ_{H}}{μ_{H}},

\Upsilon_{\beta_{HH}}^{\mathcal{R}_{01}} = \frac{\partial \mathcal{R}_{01}}{\partial \beta_{HH}} \times \frac{\beta_{HH}}{\mathcal{R}_{01}} = 1> 0

and the awareness-modulated reduction factor in the human force of infection reduces to unity since A⁰ = 0, i.e.,

\frac{1}{1 + k_{A} A^{0}} = 1 .

\Upsilon_{\mu_H}^{\mathcal{R}_{01}} = \frac{\partial \mathcal{R}_{01}}{\partial \mu_H} \times \frac{\mu_H}{\mathcal{R}_{01}} = -\frac{\mu_H (\sigma_H + \gamma_H + 2\mu_H + \delta_H)}{(\sigma_H + \mu_H)(\gamma_H + \mu_H + \delta_H)} < 0

Therefore, ℰ₀ represents the baseline demographic steady state in the absence of CCHF transmission, and it serves as the reference equilibrium for the computation of the basic reproduction number and the subsequent stability analysis.

3.4

Basic reproduction number ℛ₀

The basic reproduction number ℛ₀ is defined as the expected number of secondary infections produced by a single infected individual introduced into a fully susceptible population at the disease-free equilibrium (DFE). We compute ℛ₀ using the next-generation matrix (NGM) approach. We take as infected and latent compartments the vector

z (t) = {(I_{T} (t), E_{L} (t), I_{L} (t), E_{H} (t), I_{H} (t))}^{⊤},

\Upsilon_{\gamma_H}^{\mathcal{R}_{01}} = \frac{\partial \mathcal{R}_{01}}{\partial \gamma_H} \times \frac{\gamma_H}{\mathcal{R}_{01}} = -\frac{\gamma_H}{\gamma_H + \mu_H + \delta_H} < 0

ordered as: infected ticks (1), exposed livestock (2), infectious livestock (3), exposed humans (4), infectious humans (5). System (1) is written as ż = F(z) − V(z), where F collects rates of new infections only and V collects all remaining transition terms (progression, recovery, mortality):

F = (\begin{matrix} β_{L T} \frac{I_{L}}{N_{L}} S_{T} \\ β_{T L} \frac{I_{T}}{N_{T}} S_{L} \\ 0 \\ (β_{T H} \frac{I_{T}}{N_{T}} + β_{L H} \frac{I_{L}}{N_{L}} + β_{H H} \frac{I_{H}}{N_{H}}) S_{H} \\ 0 \end{matrix}), V = (\begin{matrix} μ_{T} I_{T} \\ (σ_{L} + μ_{L}) E_{L} \\ (γ_{L} + μ_{L}) I_{L} - σ_{L} E_{L} \\ (σ_{H} + μ_{H}) E_{H} \\ (γ_{H} + μ_{H} + δ_{H}) I_{H} - σ_{H} E_{H} . \end{matrix})

\Upsilon_{\delta_H}^{\mathcal{R}_{01}} = \frac{\partial \mathcal{R}_{01}}{\partial \delta_H} \times \frac{\delta_H}{\mathcal{R}_{01}} = -\frac{\delta_H}{\gamma_H + \mu_H + \delta_H} < 0

At ℰ₀, the awareness factor satisfies $\frac{1}{1 + k_{A} A^{0}} = 1$ L_i (since A₀ = 0 ), and $S_{T}^{0} / N_{T}^{0} = S_{L}^{0} / N_{L}^{0} = S_{H}^{0} / N_{H}^{0} = 1$ M_i. The Jacobian of F with respect to z, evaluated at ℰ₀, is

F = D F (ℰ_{0}) = (\begin{matrix} 0 & 0 & β_{L T} & 0 & 0 \\ β_{T L} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ β_{T H} & 0 & β_{L H} & 0 & β_{H H} \\ 0 & 0 & 0 & 0 & 0 \end{matrix}),

\frac{dS_T}{dt} = \Lambda_T - (1 - u_1)\lambda_T S_T - \mu_T S_T,

and the Jacobian of V is

V = D V (E_{0}) = (\begin{matrix} μ_{T} & 0 & 0 & 0 & 0 \\ 0 & σ_{L} + μ_{L} & 0 & 0 & 0 \\ 0 & - σ_{L} & γ_{L} + μ_{L} & 0 & 0 \\ 0 & 0 & 0 & σ_{H} + μ_{H} & 0 \\ 0 & 0 & 0 & - σ_{H} & γ_{H} + μ_{H} + δ_{H} \end{matrix}) .

\frac{dI_T}{dt} = (1 - u_1)\lambda_T S_T - \mu_T I_T,

The matrix V is lower block-triangular and non-singular (all diagonal entries are positive), so its inverse V⁻¹ can be computed explicitly. Denoting for brevity

m_{T} = μ_{T}, b_{L} = σ_{L} + μ_{L}, c_{L} = γ_{L} + μ_{L}, d_{H} = σ_{H} + μ_{H}, e_{H} = γ_{H} + μ_{H} + δ_{H}

\frac{dS_L}{dt} = \Lambda_L - (1 - u_2)\lambda_L S_L - \mu_L S_L,

one obtains

V^{- 1} = (\begin{matrix} \frac{1}{m_{T}} & 0 & 0 & 0 & 0 \\ 0 & \frac{1}{b_{L}} & 0 & 0 & 0 \\ 0 & \frac{σ_{L}}{b_{L} c_{L}} & \frac{1}{c_{L}} & 0 & 0 \\ 0 & 0 & 0 & \frac{1}{d_{H}} & 0 \\ 0 & 0 & 0 & \frac{σ_{H}}{d_{H} e_{H}} & \frac{1}{e_{H}} \end{matrix}) .

\frac{dE_L}{dt} = (1 - u_2)\lambda_L S_L - (\sigma_L + \mu_L) E_L,

The NGM is K = FV⁻¹. Carrying out the matrix product row by row gives

K = (\begin{matrix} 0 & \frac{β_{L T} σ_{L}}{b_{L} c_{L}} & \frac{β_{L T}}{c_{L}} & 0 & 0 \\ \frac{β_{T L}}{m_{T}} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ \frac{β_{T H}}{m_{T}} & \frac{β_{L H} σ_{L}}{b_{L} c_{L}} & \frac{β_{L H}}{c_{L}} & \frac{β_{H H} σ_{H}}{d_{H} e_{H}} & \frac{β_{H H}}{e_{H}} \\ 0 & 0 & 0 & 0 & 0 \end{matrix}) .

\frac{dI_L}{dt} = \sigma_L E_L - (\gamma_L + \mu_L) I_L,

The entry K_ij gives the expected number of new infections of type i produced by a single infected individual of type j in a fully susceptible population. Since rows 3 and 5 of K are identically zero, the eigenvalue λ = 0 has multiplicity at least two. For the remaining eigenvalues, we exploit the block structure of K. The submatrix formed by rows and columns {1,2} is

K_{TL} = (\begin{matrix} 0 & \frac{β_{L T} σ_{L}}{b_{L} c_{L}} \\ \frac{β_{T L}}{m_{T}} & 0 \end{matrix}),

\frac{dS_H}{dt} = \Lambda_H + \omega_R R_H - (1 - u_3)\lambda_H S_H - \mu_H S_H,

whose characteristic equation is $λ^{2} = \frac{β_{L T} β_{T L} σ_{L}}{m_{T} b_{L} c_{L}}$ N_i. The positive eigenvalue is therefore

R_{TL} = \sqrt{\frac{β_{L T} β_{T L} σ_{L}}{μ_{T} (σ_{L} + μ_{L}) (γ_{L} + μ_{L})}}

\frac{dE_H}{dt} = (1 - u_3)\lambda_H S_H - (\sigma_H + \mu_H) E_H,

Row 4 does not feed back into rows 1 or 2 (humans are incidental hosts with no effect on the enzootic cycle), and row 5 is zero. The (4,5) submatrix of K is

K_{HH} = (\begin{matrix} \frac{β_{H H} σ_{H}}{d_{H} e_{H}} & \frac{β_{H H}}{e_{H}} \\ 0 & 0 \end{matrix}),

\frac{dI_H}{dt} = \sigma_H E_H - (\gamma_H + \mu_H + \delta_H) I_H,

whose eigenvalues are 0 and

R_{HH} = \frac{β_{H H} σ_{H}}{(σ_{H} + μ_{H}) (γ_{H} + μ_{H} + δ_{H})} .

\frac{dR_H}{dt} = \gamma_H I_H - (\omega_R + \mu_H) R_H,

Because the two blocks do not interact through a feedback loop (the human row receives input from ticks and livestock but does not return infection to them), the spectral radius of K equals the maximum of the two positive eigenvalues:

R_{0} = ρ (K) = max {\sqrt{\frac{β_{L T} β_{T L} σ_{L}}{μ_{T} (σ_{L} + μ_{L}) (γ_{L} + μ_{L})}}, \frac{β_{H H} σ_{H}}{(σ_{H} + μ_{H}) (γ_{H} + μ_{H} + δ_{H})}} .

\frac{dA}{dt} = \eta_1 \frac{I_H}{N_H} + \eta_2 \frac{I_T}{N_T} - \omega_A A

Remark 1. The quantity ℛ_TL represents the geometric-mean reproduction number of the tick-livestock cycle: it equals the square root of the product of the expected number of livestock infected by one tick, and the expected number of ticks infected by one infectious livestock animal. This cycle governs the persistence of CCHF in nature independently of human behaviour. The quantity ℛ_HH captures the contribution of direct human-tohuman transmission; it can amplify outbreaks in healthcare or household settings but cannot sustain endemic transmission on its own in the absence of the enzootic reservoir. The transmission routes β_TH and β_LH (zoonotic spillover from ticks and livestock to humans) do not appear in ℛ₀. This is a well-known feature of multi-host models in which humans are incidental (dead-end) hosts: spillover infections seed human cases but, because humans do not return virus to the tick or livestock populations, they cannot by themselves drive the spectral radius above the threshold. The spillover pathways do, however, directly determine the rate at which human cases are generated once the enzootic cycle is established, and hence critically influence epidemic size and peak timing. Although public awareness does not affect ℛ₀ (since A⁰ = 0 at the DFE), it plays a crucial role in reducing the effective reproduction number during active outbreaks.

3.5

Local stability of the disease-free equilibrium

Theorem 3.

Let ℰ₀ denote the disease-free equilibrium of system (1). Then ℰ₀ is locally asymptotically stable if ℛ₀ < 1 and unstable if ℛ₀ > 1, where ℛ₀ is the basic reproduction number.

Proof.

Linearising the full system (1) at ℰ₀ produces a Jacobian with block upper-triangular structure. The block corresponding to the infection-free (demographic) subsystem has eigenvalues −μ_T, −μ_L, −μ_H (and − ω_A for the awareness equation), all strictly negative. Hence the stability of ℛ₀ is determined entirely by the infected subsystem.

From the standard NGM theory, the infected subsystem has all eigenvalues with strictly negative real parts if and only if ρ(FV⁻¹) < 1, i.e. ℛ₀ < 1. Therefore:

If ℛ₀ < 1, all Jacobian eigenvalues at ℰ₀ have negative real parts, so ℰ₀ is locally asymptotically stable.
If ℛ₀ > 1, at least one eigenvalue has a positive real part, so ℰ₀ is unstable.

3.6

Sensitivity analysis

In this subsection, we evaluate the sensitivity indices of the basic reproduction number ℛ₀ with respect to key epidemiological parameters. Recall from Section 3 that the basic reproduction number is

R_{0} = max {R_{TL}, R_{HH}},

J(u_1, u_2, u_3) = \int_0^T \left( A_1 I_T + A_2 I_L + A_3 I_H + \frac{1}{2} (B_1 u_1^2 + B_2 u_2^2 + B_3 u_3^2) \right) dt

where

R_{TL} = \sqrt{\frac{β_{L T} β_{T L} σ_{L}}{μ_{T} (σ_{L} + μ_{L}) (γ_{L} + μ_{L})}}, R_{HH} = \frac{β_{H H} σ_{H}}{(σ_{H} + μ_{H}) (γ_{H} + μ_{H} + δ_{H})} .

H = A_1 I_T + A_2 I_L + A_3 I_H + \frac{1}{2} B_1 u_1^2 + \frac{1}{2} B_2 u_2^2 + \frac{1}{2} B_3 u_3^2 + \sum_{i=1}^{10} \lambda_i f_i

The quantity ℛ_TL measures the enzootic tick-livestock transmission cycle, which governs the persistence of CCHF in nature, while ℛ_HH captures human-to-human transmission, which may amplify outbreaks but cannot sustain endemic transmission in the absence of the enzootic reservoir. Since ℛ₀ = max{ℛ_TL, ℛ_HH}, its sensitivity indices are those of the dominant branch. Under the parameter values of Table 3, we verify numerically that ℛ_TL ≫ ℛ_HH, so ℛ₀ = ℛ_TL at these parameter values and the sensitivity indices are computed with respect to ℛ_TL. Although public awareness does not influence ℛ₀ directly (since A⁰ = 0 at the DFE), it plays a crucial role in reducing the effective reproduction number during epidemic outbreaks through the modulation factor $\frac{1}{1 + k_{A} A (t)}$ O_i. To quantify the relative influence of each model parameter on ℛ₀, we compute the normalized local sensitivity index

Υ_{p}^{R_{0}} = \frac{\partial R_{0}}{\partial p} \cdot \frac{p}{R_{0}},

\frac{d\lambda_1}{dt} = -\frac{\partial H}{\partial S_T} = (1 - u_1)\lambda_T \lambda_1 - (1 - u_1)\lambda_T \lambda_2 + \mu_T \lambda_1

following the standard approach and under the parameter values of Table 3, the enzootic branch dominates, i.e., ℛ_TL ≫ ℛ_HH, so the sensitivity indices of the active branch are:

\begin{array}{l} Υ_{β_{T L}}^{R_{0}} = + \frac{1}{2}, & Υ_{β_{L T}}^{R_{0}} = + \frac{1}{2}, & Υ_{μ_{T}}^{R_{0}} = - \frac{1}{2}, \\ Υ_{γ_{L}}^{R_{0}} = - \frac{1}{2} \cdot \frac{γ_{L}}{γ_{L} + μ_{L}} \approx - 0.499, & Υ_{σ_{L}}^{R_{0}} = \frac{μ_{L}}{2 (σ_{L} + μ_{L})} \approx + 0.001, & Υ_{k_{A}}^{R_{0}} = 0 \end{array}

\Upsilon_{\beta_{TL}}^{\mathcal{R}_0}=+\frac{1}{2}, \qquad \Upsilon_{\beta_{LT}}^{\mathcal{R}_0}=+\frac{1}{2}, \qquad \Upsilon_{\mu_T}^{\mathcal{R}_0}=-\frac{1}{2}, \\ \Upsilon_{\gamma_L}^{\mathcal{R}_0} =-\frac{1}{2}\cdot\frac{\gamma_L}{\gamma_L+\mu_L} \approx -0.499, \qquad \Upsilon_{\sigma_L}^{\mathcal{R}_0} =\frac{\mu_L}{2(\sigma_L+\mu_L)} \approx +0.001, \qquad \Upsilon_{k_A}^{\mathcal{R}_0}=0.

These results yield three epidemiologically important conclusions. First, the tick-to-livestock and livestock-totick transmission rates β_TL and β_LT are the most influential parameters, each carrying index +1/2, identifying reduction of tick-livestock contact as the highest-priority control target. Second, the tick natural mortality rate μ_T and livestock recovery rate γ_L each carry index −1/2, indicating that acaricide treatment and veterinary intervention are the most effective biological countermeasures for suppressing ℛ₀. Third, and most importantly for the present study, the awareness strength parameter k_A satisfies $Υ_{k_{A}}^{R_{0}} = 0$ P_i exactly, since A⁰ = 0 at the diseasefree equilibrium. This result confirms analytically that awareness-based interventions do not shift the invasion threshold ℛ₀, but instead act by reducing the effective reproduction number during active outbreaks through the saturation factor $\frac{1}{1 + k_{A} A (t)}$ Q_i. This distinction has direct implications for public health strategy: awareness campaigns alone cannot prevent disease invasion when ℛ₀ > 1, but they are critical for mitigating outbreak magnitude and accelerating disease suppression once an outbreak is underway.

Chebyshev functions: convergence analysis

In this section, the Chebyshev polynomial basis is introduced and employed for the numerical approximation of the proposed model. First, some essential properties of Chebyshev polynomials are recalled. Then, their shifted form on a finite interval is constructed, followed by the formulation of the collocation approximation. Finally, a rigorous convergence and error analysis is established in the L²-norm.

4.1

An overview of Chebyshev polynomials

The Chebyshev polynomials of the first kind ${T_{ℓ} (τ)}_{ℓ = 0}^{\infty}$ R_i are defined on the interval τ ∈ [−1, 1] by

T_{ℓ} (τ) = \cos (ℓ \arccos τ), ℓ \in N_{0}

\frac{d\lambda_3}{dt} = -\frac{\partial H}{\partial S_L} = (1 - u_2)\lambda_L \lambda_3 - (1 - u_2)\lambda_L \lambda_4 + \mu_L \lambda_3

The first few Chebyshev polynomials are given by

\begin{array}{l} T_{0} (τ) = 1, \\ T_{1} (τ) = τ, \\ T_{2} (τ) = 2 τ^{2} - 1, \\ T_{3} (τ) = 4 τ^{3} - 3 τ, \\ T_{4} (τ) = 8 τ^{4} - 8 τ^{2} + 1 . \end{array}

\frac{d\lambda_4}{dt} = -\frac{\partial H}{\partial E_L} = (\sigma_L + \mu_L)\lambda_4 - \sigma_L \lambda_5

These polynomials satisfy the recurrence relation

T_{ℓ + 1} (τ) = 2 τ T_{ℓ} (τ) - T_{ℓ - 1} (τ), ℓ \geq 1,

\frac{d\lambda_5}{dt} = -\frac{\partial H}{\partial I_L} = -A_2 + (1 - u_2)\frac{\partial \lambda_L}{\partial I_L} S_L \lambda_3 - (1 - u_2)\frac{\partial \lambda_L}{\partial I_L} S_L \lambda_4 + (\gamma_L + \mu_L)\lambda_5

and form an orthogonal basis with respect to the weight function

w (τ) = \frac{1}{\sqrt{1 - τ^{2}}}, τ \in (- 1, 1) .

\frac{d\lambda_6}{dt} = -\frac{\partial H}{\partial S_H} = (1 - u_3)\lambda_H \lambda_6 - (1 - u_3)\lambda_H \lambda_7 + \mu_H \lambda_6

Let L be a positive integer. We define the vector of Chebyshev polynomials as

T_{L} (τ) := [T_{0} (τ), T_{1} (τ), \dots, T_{L} (τ)] .

\frac{d\lambda_7}{dt} = -\frac{\partial H}{\partial E_H} = (\sigma_H + \mu_H)\lambda_7 - \sigma_H \lambda_8

It is well known that Chebyshev polynomials admit a monomial representation. Let

M_{L} (τ) := [1, τ, τ^{2}, \dots, τ^{L}] .

\frac{d\lambda_8}{dt} = -\frac{\partial H}{\partial I_H} = -A_3 + (1 - u_3)\frac{\partial \lambda_H}{\partial I_H}S_H \lambda_6 - (1 - u_3)\frac{\partial \lambda_H}{\partial I_H}S_H \lambda_7 + (\gamma_H + \mu_H + \delta_H)\lambda_8 - \gamma_H \lambda_9 - \eta_1 \frac{N_H - I_H}{N_H^2}\lambda_{10}

Then, there exists an upper triangular matrix C_L ∈ R^(L+1)×(L+1) such that

T_{L} (τ) = M_{L} (τ) C_{L} .

\frac{d\lambda_9}{dt} = -\frac{\partial H}{\partial R_H} = -\omega_R \lambda_6 + (\omega_R + \mu_H)\lambda_9

Then the matrix C_L is nonsingular and depends only on the polynomial degree L.

4.2

Shifted Chebyshev polynomials

Since the proposed model is defined on a finite interval t ∈ [0, T], we employ the shifted Chebyshev polynomials. Let the linear transformation

τ = \frac{2 t}{T} - 1

\frac{d\lambda_{10}}{dt} = -\frac{\partial H}{\partial A} = \omega_A \lambda_{10}

map the interval [0, T] onto [−1, 1]. The shifted Chebyshev polynomials are defined by

{\tilde{T}}_{ℓ} (t) := T_{ℓ} (\frac{2 t}{T} - 1), ℓ \in N_{0} .

\text{With transversality conditions: } \lambda_i(T) = 0, \quad i = 1, 2, \dots, 10

Accordingly, the vector of shifted Chebyshev polynomials is written as

{\tilde{T}}_{L} (t) := [{\tilde{T}}_{0} (t), {\tilde{T}}_{1} (t), \dots, {\tilde{T}}_{L} (t)] .

\frac{\partial H}{\partial u_1} = B_1 u_1 + \lambda_T S_T \lambda_1 - \lambda_T S_T \lambda_2 = 0

Using the monomial representation, we obtain

{\tilde{T}}_{L} (t) = {\tilde{M}}_{L} (t) C_{L},

\frac{\partial H}{\partial u_2} = B_2 u_2 + \lambda_L S_L \lambda_3 - \lambda_L S_L \lambda_4 = 0

where

{\tilde{M}}_{L} (t) := [1, t, t^{2}, \dots, t^{L}] .

\frac{\partial H}{\partial u_3} = B_3 u_3 + \lambda_H S_H \lambda_6 - \lambda_H S_H \lambda_7 = 0

4.3

Chebyshev approximation and collocation representation

Let u(t) ∈ L²([0, T]) be a sufficiently smooth function. Then u(t) can be expanded in terms of shifted Chebyshev polynomials as

u (t) = \sum_{ℓ = 0}^{\infty} a_{ℓ} {\tilde{T}}_{ℓ} (t), t \in [0, T] .

u_1^* = \max\left(0, \min\left(1, \frac{\lambda_T S_T(\lambda_2 - \lambda_1)}{B_1}\right)\right)

In practical computations, the infinite series is truncated, and the approximate solution is given by

u_{L} (t) := \sum_{ℓ = 0}^{L} a_{ℓ} {\tilde{T}}_{ℓ} (t) .

u_2^* = \max\left(0, \min\left(1, \frac{\lambda_L S_L(\lambda_4 - \lambda_3)}{B_2}\right)\right)

This approximation can be written in compact matrix form as

u_{L} (t) = {\tilde{T}}_{L} (t) a_{L},

u_3^* = \max\left(0, \min\left(1, \frac{\lambda_H S_H(\lambda_7 - \lambda_6)}{B_3}\right)\right)

where

a_{L} := {[a_{0}, a_{1}, \dots, a_{L}]}^{T} .

\Lambda_T = 325.2154, \quad \Lambda_L = 25.1021, \quad \Lambda_H = 15.2145, \quad \mu_T = 0.02, \quad \mu_L = 0.02, \quad \mu_H = 0.02,

4.4

Operational matrices of differentiation

In this subsection, the operational matrices corresponding to the differentiation of shifted Chebyshev polynomials are introduced. These matrices play a crucial role in constructing the collocation scheme and transforming the governing system into an algebraic form.

Lemma 4. Let T̃_L(t) be the vector of shifted Chebyshev polynomials defined on [0, T] as

{\tilde{T}}_{L} (t) = [{\tilde{T}}_{0} (t), {\tilde{T}}_{1} (t), \dots, {\tilde{T}}_{L} (t)] .

\beta_T = 0.00035, \quad \beta_L = 0.000105, \quad \beta_H = 0.00012, \quad \sigma_L = 0.05, \quad \sigma_H = 0.05,

Then, the first-order derivative of T̃_L(t) can be expressed as

\frac{d}{d t} {\tilde{T}}_{L} (t) = {\tilde{T}}_{L} (t) D^{(1)},

\gamma_L = 0.025, \quad \gamma_H = 0.025, \quad \delta_H = 0.05, \quad \omega_R = 0.08, \quad \omega_A = 0.05,

where D(1) ∈ R(L+1)×(L+1) is the Chebyshev operational matrix of first-order differentiation.

Proof. Using the chain rule and the definition of shifted Chebyshev polynomials

{\tilde{T}}_{ℓ} (t) = T_{ℓ} (\frac{2 t}{T} - 1),

\eta_1 = 0.1, \quad \eta_2 = 0.1, \quad A_1 = 120, \quad A_2 = 10, \quad A_3 = 10, \quad B_1 = 50, \quad B_2 = 50, \quad B_3 = 50.

we obtain

\frac{d}{d t} {\tilde{T}}_{ℓ} (t) = \frac{2}{T} T_{ℓ}^{'} (τ), τ = \frac{2 t}{T} - 1 .

\mathcal{R}_{01} = 2.4578> 1

It is well known that the derivative of Chebyshev polynomials satisfies

T_{ℓ}^{'} (τ) = ℓ U_{ℓ - 1} (τ),

\mathcal{R}_{01} = 0.1985 < 1

where U_ℓ(τ) denotes the Chebyshev polynomial of the second kind. Since U_ℓ−1(τ) can be expanded as a finite linear combination of T_k(τ) for k ≤ ℓ −1, the derivative $\frac{d}{d t} {\tilde{T}}_{ℓ} (t)$ S_i can be written as a linear combination of ${{\tilde{T}}_{k} (t)}_{k = 0}^{L}$ T_i. Collecting the coefficients yields the operational matrix D(1).

Lemma 5. The entries of the first-order Chebyshev operational matrix $D^{(1)} = {[d_{i j}^{(1)}]}_{i, j = 0}^{L}$ U_i are given by

d_{i j}^{(1)} = {\begin{cases} \frac{2 j}{T}, & i + j odd and i < j, \\ 0, & otherwise. \end{cases}

d_{ij}^{(1)}={\{\begin{array}{lc}\frac{2j}T,&i+j\textit{odd and}i<j,\\0,&\textit{otherwise.}\end{array}}

Lemma 6. Let D(1) be the first-order operational matrix of differentiation. Then, the higher-order derivatives of T_L(t) satisfy

\frac{d^{n}}{d t^{n}} {\tilde{T}}_{L} (t) = {\tilde{T}}_{L} (t) D^{(n)}, n \geq 1,

\mathcal{R}_{01} = 0.5478 < 1

where

D^{(n)} = {(D^{(1)})}^{n} .

\mathcal{R}_{02} = 1.3458> 1

Proof. The result follows by repeated application of Lemma 4 and linearity.

Lemma 7. Let $u_{L} (t) = {\tilde{T}}_{L} (t) a_{L}$ V_i be the Chebyshev approximation of u(t). Then, its n-th derivative is given by

\frac{d^{n}}{d t^{n}} u_{L} (t) = {\tilde{T}}_{L} (t) D^{(n)} a_{L} .

\mathcal{R}_{02} = 0.0845 < 1

4.5

Convergence analysis of Chebyshev polynomials in the L² norm

In this work, the Chebyshev polynomial basis is employed on the finite interval [0, T], where T > 0. The convergence properties of Chebyshev polynomials are well established and form the theoretical foundation of spectral collocation methods. In this subsection, we rigorously examine the approximation error associated with Chebyshev polynomial expansions in the L₂ norm.

Let N(t) ∈ L²([0, T]) be a given function. Using the shifted Chebyshev polynomials ${{\tilde{T}}_{ℓ} (t)}_{ℓ = 0}^{\infty}$ W_i defined on [0, T], the function n(t) can be represented as

N (t) = \sum_{ℓ = 0}^{\infty} π_{ℓ} {\tilde{T}}_{ℓ} (t), t \in [0, T],

\mathcal{R}_{02} = 0.4215 < 1

where π_ℓ, ℓ ≥ 0, are the Chebyshev coefficients of N(t).

In practical computations, we restrict attention to the finite-dimensional subspace

W_{L} := span {{\tilde{T}}_{0} (t), {\tilde{T}}_{1} (t), \dots, {\tilde{T}}_{L} (t)} \subset L^{2} ([0, T]) .

\mathcal{R}_{02} = 0.3124 < 1

Accordingly, we approximate N(t) by retaining only the first (L+1) Chebyshev modes:

N (t) \approx N_{L} (t) := \sum_{ℓ = 0}^{L} π_{ℓ} {\tilde{T}}_{ℓ} (t), t \in [0, T] .

\mathcal{R}_{03} = 1.8456> 1

For convenience, the truncated approximation (21) can be written in vector form as

N_{L} (t) = {\tilde{T}}_{L} (t) π_{L},

\mathcal{R}_{03} = 0.1124 < 1

where

π_{L} := {[π_{0}, π_{1}, \dots, π_{L}]}^{⊤} .

\mathcal{R}_{03} = 0.5142 < 1

We now define the approximation error

E_{L} (t) := N (t) - N_{L} (t) = \sum_{ℓ = L + 1}^{\infty} π_{ℓ} {\tilde{T}}_{ℓ} (t) .

\mathcal{R}_{03} = 0.4215 < 1

Our objective is to derive an upper bound for ${‖ E_{L} ‖}_{L^{2} ([0, T])}$ X_i. To this end, we recall a classical result concerning Chebyshev polynomial approximation [26]. The following theorem provides an error estimate in the L² norm. Let t_b > 0 and define the affine mapping $x = χ (t) := \frac{2 t}{t_{b}} - 1$ Y_i from [0, t_b] onto [−1, 1]. Let T_ℓ(x) be the Chebyshev polynomials of the first kind and define the shifted Chebyshev polynomials on [0, t_b] by

T_{ℓ}^{*} (t) := T_{ℓ} χ (t), ℓ = 0, 1, 2, \dots .

\mathcal{R}_{0} = 2.4578> 1

We introduce the Chebyshev weight on [0, t_b],

w (t) := \frac{1}{\sqrt{t (t_{b} - t)}},

\mathcal{R}_{0} = 0.1985 < 1

and the weighted inner product and norm

⟨ f, g ⟩_{w} := \int_{0}^{t_{b}} f (t) g (t) w (t) d t, ‖ f ‖_{L_{w}^{2} (0, t_{b})} := \sqrt{⟨ f, f ⟩_{w}} .

\mathcal{R}_{0} = 0.8145 < 1

It is well known that ${T_{ℓ}^{*}}_{ℓ \geq 0}$ Z_i is orthogonal in $L_{w}^{2} (0, t_{b})$ a_i.

Theorem 8.

Assume that $N \in H_{w}^{m} (0, t_{b})$ b_i for some integer m ≥ 1, i.e., $N^{(m)} \in L_{w}^{2} (0, t_{b})$ c_i. Let P_LN be the orthogonal projection of N onto $ℙ_{L} := span {T_{0}^{*}, T_{1}^{*}, \dots, T_{L}^{*}}$ d_i in $L_{w}^{2} (0, t_{b})$ e_i. Then there exists a constant C_m > 0, independent of L and t_b, such that

24 ${‖ N - P_{L} N ‖}_{L_{w}^{2} (0, t_{b})} \leq C_{m} {(\frac{t_{b}}{2})}^{m} L^{- m} {‖ N^{(m)} ‖}_{L_{w}^{2} (0, t_{b})}, L \geq 1 .$ \mathcal{R}_{0} = 0.5478 < 1

Consequently, ${‖ N - P_{L} N ‖}_{L_{w}^{2} (0, t_{b})} \to 0$ f_i as L → ∞.

Proof.

First, we define $\hat{N} (x) := n (\frac{t_{b}}{2} (x + 1))$ g_i for x ∈ [−1, 1]. Using $t = \frac{t_{b}}{2} (x + 1)$ h_i and $d t = \frac{t_{b}}{2} d x$ i_i, we obtain

$w (t) d t = \frac{1}{\sqrt{t (t_{b} - t)}} d t = \frac{1}{\sqrt{\frac{t_{b}}{2} (x + 1) \cdot \frac{t_{b}}{2} (1 - x)}} \cdot \frac{t_{b}}{2} d x = \frac{d x}{\sqrt{1 - x^{2}}} .$ J(u_1, u_2, u_3) = \int_{0}^{T} \left(A_1 I_T + A_2 I_L + A_3 I_H + \frac{1}{2}B_1 u_1^2 + \frac{1}{2}B_2 u_2^2 + \frac{1}{2}B_3 u_3^2\right) dt

Hence,

$‖ N ‖_{L_{w}^{2} (0, t_{b})} = ‖ \hat{N} ‖_{L_{ω}^{2} (- 1, 1)}, ω (x) := \frac{1}{\sqrt{1 - x^{2}}} .$ \[ \|N\|_{L^2_w(0,t_b)}=\|\widehat N\|_{L^2_{\omega}(-1,1)}, \qquad\omega(x):=\frac{1}{\sqrt{1-x^2}}. \]

Moreover, since $\frac{d}{d x} \hat{N} (x) = \frac{t_{b}}{2} N^{'} (t)$ $\frac{d}{dx}\widehat N(x)=\frac{t_b}{2}\,N'(t)$, repeated differentiation gives

${\hat{N}}^{(m)} (x) = {(\frac{t_{b}}{2})}^{m} N^{(m)} (t), \Rightarrow {‖ {\hat{N}}^{(m)} ‖}_{L_{ω}^{2} (- 1, 1)} = {(\frac{t_{b}}{2})}^{m} {‖ N^{(m)} ‖}_{L_{w}^{2} (0, t_{b})} .$ \[ \widehat N^{(m)}(x)=\Big(\frac{t_b}{2}\Big)^{m} N^{(m)}(t), \qquad \Rightarrow\qquad \|\widehat N^{(m)}\|_{L^2_{\omega}(-1,1)}=\Big(\frac{t_b}{2}\Big)^{m}\|N^{(m)}\|_{L^2_w(0,t_b)}. \]

Let ${\hat{P}}_{L} \hat{N}$ $\widehat P_L \widehat N$ be the $L_{ω}^{2} (- 1, 1)$ $L^2_{\omega}(-1,1)$-orthogonal projection of $\hat{N}$ $\widehat N$ onto ${\hat{P}}_{L} := span {T_{0}, T_{1}, \dots, T_{L}}$ $\widehat{\mathbb{P}}_L:=\mathrm{span}\{T_0,T_1,\cdots,T_L\}$. By construction of $T_{ℓ}^{*}$ $T_\ell^{*}$,

$(P_{L} N) (t) = ({\hat{P}}_{L} \hat{N}) (χ (t)), and {‖ N - P_{L} N ‖}_{L_{w}^{2} (0, t_{b})} = {‖ \hat{N} - {\hat{P}}_{L} \hat{N} ‖}_{L_{ω}^{2} (- 1, 1)} .$ \[(P_L N)(t)= (\widehat P_L \widehat N)\big(\chi(t)\big), \qquad\text{and}\qquad \|N-P_L N\|_{L^2_w(0,t_b)}=\|\widehat N-\widehat P_L\widehat N\|_{L^2_\omega(-1,1)}. \]

Since P̂ᴌ is an orthogonal projector onto the finite-dimensional subspace ℙ̂ᴌ in the Hilbert space $L_{ω}^{2} (- 1, 1)$ \textstyle L_\omega^2(-1,1), it yields the best approximation of N̂ from ℙ̂ᴌ in the norm $‖ \cdot ‖_{L_{ω}^{2}}$ $\|\cdot\|_{L^2_\omega}$, and this best approximation is unique. Indeed, since $L_{ω}^{2} (- 1, 1)$ \textstyle L_\omega^2(-1,1) is a Hilbert space and ℙ̂ᴌ is a closed finite-dimensional subspace, the projection theorem guarantees that for every $\hat{N} \in L_{ω}^{2} (- 1, 1)$ $\widehat{N}\in L^2_\omega(-1,1)$ there exists a unique element p̂_LN̂ ∈ ℙ̂_L satisfying

25 ${‖ \hat{N} - {\hat{P}}_{L} \hat{N} ‖}_{L_{ω}^{2} (- 1, 1)} = min_{p \in {\hat{ℙ}}_{L}} ‖ \hat{N} - p ‖_{L_{ω}^{2} (- 1, 1)} .$ \|\widehat{N}-\widehat{P}_L\widehat{N}\|_{L^2_\omega(-1,1)} = \min_{p\in\widehat{\mathbb{P}}_L} \|\widehat{N}-p\|_{L^2_\omega(-1,1)}.

A standard Jackson inequality for Chebyshev-weighted approximation states that for m ≥ 1 and $\hat{N} \in H_{ω}^{m} (- 1, 1)$ $\widehat N\in H^{m}_\omega(-1,1)$ there exists pL ∈ ℙ̂_L such that

26 ${‖ \hat{N} - p_{L} ‖}_{L_{ω}^{2} (- 1, 1)} \leq C_{m} L^{- m} {‖ {\hat{N}}^{(m)} ‖}_{L_{ω}^{2} (- 1, 1)},$ \|\widehat N-p_L\|_{L^2_\omega(-1,1)} \le C_m\,L^{-m}\,\|\widehat N^{(m)}\|_{L^2_\omega(-1,1)},

with C_m independent of L. Combining (25) and (26) yields

${‖ \hat{N} - {\hat{P}}_{L} \hat{N} ‖}_{L_{ω}^{2} (- 1, 1)} \leq C_{m} L^{- m} {‖ {\hat{N}}^{(m)} ‖}_{L_{ω}^{2} (- 1, 1)} .$ \[ \|\widehat N-\widehat P_L\widehat N\|_{L^2_\omega(-1,1)} \le C_m\,L^{-m}\,\|\widehat N^{(m)}\|_{L^2_\omega(-1,1)}. \]

Then we have

${‖ N - P_{L} N ‖}_{L_{w}^{2} (0, t_{b})} = {‖ \hat{N} - {\hat{P}}_{L} \hat{N} ‖}_{L_{ω}^{2} (- 1, 1)} \leq C_{m} L^{- m} {‖ {\hat{N}}^{(m)} ‖}_{L_{ω}^{2} (- 1, 1)} = C_{m} {(\frac{t_{b}}{2})}^{m} L^{- m} {‖ N^{(m)} ‖}_{L_{w}^{2} (0, t_{b})},$ {\Arrowvert N-P_LN\Arrowvert}_{L_w^2{(0,t_b)}}={\Arrowvert\widehat N-{\widehat P}_L\widehat N\Arrowvert}_{L_\omega^2(-1,1)}\leq C_mL^{-m}{\Arrowvert\widehat N^{(m)}\Arrowvert}_{L_\omega^2(-1,1)}=C_m{(\frac{t_b}2)}^mL^{-m}{\Arrowvert N^{(m)}\Arrowvert}_{L_w^2{(0,t_b)}}

which proves (24). Since m ≥ 1, the right-hand side tends to 0 as L → ∞, which completes the proof.

The QLM-Chebyshev matrix technique based on a domain decomposition strategy

5.1

Domain decomposition strategy

In this subsection, we develop an accurate Chebyshev spectral collocation algorithm for solving the proposed CCHF transmission model on the temporal domain [0, t_b], where t_b > 0 is sufficiently large. It is well known that applying classical global collocation techniques on long computational intervals may lead to loss of accuracy or poor convergence. To overcome this difficulty, we adopt a domain decomposition strategy, whereby the global interval is partitioned into several smaller subdomains and the collocation procedure is applied locally in a sequential manner. Let the interval [0, t_b] be partitioned into U ≥ 1 non-overlapping subintervals such that

0 = t_{0} < t_{1} < t_{2} < \dots < t_{U} = t_{b} .

\[0 = t_0 < t_1 < t_2 < \cdots < t_U = t_b. \]

We denote the u-th subdomain by

Q_{u} := [t_{u - 1}, t_{u}], u = 1, 2, \dots, U .

\[Q_u := [t_{u-1},t_u], \qquad u=1,2,\cdots,U. \]

On each subdomain Q_u, we approximate the solution of the model using Chebyshev polynomials. Specifically, the approximate solution on Q_u is expressed as

y_{L}^{u} (t) := \sum_{ℓ = 0}^{L} π_{ℓ}^{u} {\tilde{T}}_{ℓ} (t) = {\tilde{T}}_{L} (t) Π_{L}^{u}, t \in Q_{u},

\mathbf{y}^{\,u}_L(t):= \sum_{\ell=0}^{L} \boldsymbol{\pi}^{\,u}_{\ell}\, \widetilde{T}_{\ell}(t) = \widetilde{\mathbf{T}}_L(t)\,\boldsymbol{\Pi}^{\,u}_L, \qquad t\in Q_u, \label{eq:local_cheb_expansion}

where 𝘛̃_ℓ(𝘵) are the shifted Chebyshev polynomials on $Q_{u}, Π_{L}^{u} := {[π_{0}^{u}, π_{1}^{u}, \dots, π_{L}^{u}]}^{⊤}$ $\boldsymbol{\Pi}^{\,u}_L := [\boldsymbol{\pi}^{\,u}_{0},\boldsymbol{\pi}^{\,u}_{1},\cdots,\boldsymbol{\pi}^{\,u}_{L}]^{\top}$ is the vector of unknown coefficient vectors, and L denotes the polynomial degree. The global approximate solution on [0, t_b] is then given by

y_{L} (t) := \sum_{u = 1}^{U} 1_{Q_{u}} (t) y_{L}^{u} (t), 1_{Q_{u}} (t) = {\begin{matrix} 1, & t \in Q_{u} & , \\ 0, & t \notin Q_{u} . \end{matrix}

\mathbf{y}_L(t):= \sum_{u=1}^{U} \mathbf{1}_{Q_u}(t)\,\mathbf{y}^{\,u}_L(t), \qquad \mathbf{1}_{Q_u}(t)=

To determine the ( L+1 ) unknown coefficient vectors on each subdomain, we employ ( L+1 ) Chebyshev-GaussLobatto collocation points on Q_u, defined as

C_{L}^{u} := {t_{u, l} = \frac{t_{u} - t_{u - 1}}{2} (1 + \cos \frac{l π}{L}) + t_{u - 1} : l = 0, 1, \dots, L} .

\mathcal{C}^{\,u}_L := \left\{ t_{u,l} = \frac{t_u-t_{u-1}}{2}\left(1+\cos\frac{l\pi}{L}\right)+t_{u-1} \,:\, l=0,1,\cdots,L \right\}. \label{eq:collocation_points}

On the first subdomain Q₁, the original initial conditions are imposed. On each subsequent subdomain Q_u(u ≥ 2), the numerical solution obtained at t_u − 1 is used as the initial condition for the local problem, ensuring continuity across subdomains.

Remark 2. The partition above is non-uniform in general, meaning that the subinterval lengths h_u := t_u − t_u − 1 are not required to be equal across subdomains. A uniform partition is recovered as the special case h_u = t_b/U for all u = 1, 2, ⋯, U. In the numerical experiments reported in this work, a uniform partition is adopted for simplicity and reproducibility, as the CCHF model parameters do not exhibit abrupt temporal changes that would necessitate local refinement. The above partition has the advantages of being simple to implement and fully reproducible which requires no prior knowledge of the solution behavior and the Chebyshev-Gauss-Lobatto nodes on each subdomain are generated by the same affine transformation and also reducing implementation complexity which is well suited for problems whose solutions evolve smoothly and uniformly in time, as is the case for the CCHF model considered here. On the other hand, the disadvantages of potentially inefficient when the solution exhibits localized rapid changes or stiff transients in certain time regions, since computational effort is distributed equally across all subdomains regardless of local solution complexity; may require a larger number of subdomains U to maintain accuracy near sharp features.

Next, we will illustrate the application of the QLM-Chebyshev collocation technique.

5.2

Fundamentals of the QLM

The domain decomposition strategy preserves numerical accuracy over long time intervals; however, due to the strong nonlinearity of the underlying epidemiological model, direct collocation may still suffer from slow or unstable convergence. To address this issue, we employ the QLM, which transforms the nonlinear system into a sequence of linear subproblems. The CCHF awareness-behavior model can be written concisely as

\frac{d y (t)}{d t} = H (t, y (t)),

\frac{d\mathbf{y}(t)}{dt}=\mathbf{H}(t,\mathbf{y}(t)),

where

y (t) = (\begin{matrix} S_{T} (t) \\ I_{T} (t) \\ S_{L} (t) \\ E_{L} (t) \\ I_{L} (t) \\ S_{H} (t) \\ E_{H} (t) \\ I_{H} (t) \\ R_{H} (t) \\ A (t) \end{matrix}), H (t, y) = (\begin{matrix} Λ_{T} - λ_{T} (t) S_{T} (t) - μ_{T} S_{T} (t) \\ λ_{T} (t) S_{T} (t) - μ_{T} I_{T} (t) \\ Λ_{L} - λ_{L} (t) S_{L} (t) - μ_{L} S_{L} (t) \\ λ_{L} (t) S_{L} (t) - (σ_{L} + μ_{L}) E_{L} (t) \\ σ_{L} E_{L} (t) - (γ_{L} + μ_{L}) I_{L} (t) \\ Λ_{H} + ω_{R} R_{H} (t) - λ_{H} (t) S_{H} (t) - μ_{H} S_{H} (t) \\ λ_{H} (t) S_{H} (t) - (σ_{H} + μ_{H}) E_{H} (t) \\ σ_{H} E_{H} (t) - (γ_{H} + μ_{H} + δ_{H}) I_{H} (t) \\ γ_{H} I_{H} (t) - (ω_{R} + μ_{H}) R_{H} (t) \\ η_{1} \frac{I_{H} (t)}{N_{H} (t)} + η_{2} \frac{I_{T} (t)}{N_{T} (t)} - ω_{A} A (t) \end{matrix}) .

\[\mathbf{y}(t)= \begin{pmatrix} S_T(t)\\ I_T(t)\\ S_L(t)\\ E_L(t)\\ I_L(t)\\ S_H(t)\\ E_H(t)\\ I_H(t)\\ R_H(t)\\ A(t) \end{pmatrix}, \qquad \mathbf{H}(t,\mathbf{y})=\begin{pmatrix} \Lambda_T-\lambda_T(t)\,S_T(t)-\mu_T S_T(t)\\[2mm] \lambda_T(t)\,S_T(t)-\mu_T I_T(t)\\[2mm] \Lambda_L-\lambda_L(t)\,S_L(t)-\mu_L S_L(t)\\[2mm] \lambda_L(t)\,S_L(t)-(\sigma_L+\mu_L)E_L(t)\\[2mm] \sigma_L E_L(t)-(\gamma_L+\mu_L)I_L(t)\\[2mm] \Lambda_H+\omega_R R_H(t)-\lambda_H(t)\,S_H(t)-\mu_H S_H(t)\\[2mm] \lambda_H(t)\,S_H(t)-(\sigma_H+\mu_H)E_H(t)\\[2mm] \sigma_H E_H(t)-(\gamma_H+\mu_H+\delta_H)I_H(t)\\[2mm] \gamma_H I_H(t)-(\omega_R+\mu_H)R_H(t)\\[2mm] \eta_1\frac{I_H(t)}{N_H(t)}+\eta_2\frac{I_T(t)}{N_T(t)}-\omega_A A(t) \end{pmatrix}. \]

Let y₀(t) denote an initial approximation to the solution of (30). The QLM iteration is defined as

\frac{d y_{q} (t)}{d t} \approx H (t, y_{q - 1} (t)) + H_{y} (t, y_{q - 1} (t)) (y_{q} (t) - y_{q - 1} (t)), q = 1, 2, \dots,

\[ \frac{d\mathbf{y}_{q}(t)}{dt} \approx \mathbf{H}(t,\mathbf{y}_{q-1}(t)) + \mathbf{H}_{\mathbf{y}}(t,\mathbf{y}_{q-1}(t)) \big(\mathbf{y}_{q}(t)-\mathbf{y}_{q-1}(t)\big), \qquad q=1,2,\cdots, \]

where H_y denotes the Jacobian matrix of H with respect to y. After rearrangement, the QLM linearized system is obtained as

\frac{d y_{q} (t)}{d t} + E_{q - 1} (t) y_{q} (t) = r_{q - 1} (t), q = 1, 2, \dots,

\frac{d\mathbf{y}_{q}(t)}{dt} + \mathbf{E}_{q-1}(t)\,\mathbf{y}_{q}(t) = \mathbf{r}_{q-1}(t), \qquad q=1,2,\cdots,

where

E_{q - 1} (t) = - H_{y} (t, y_{q - 1} (t)), R_{q - 1} (t) = H (t, y_{q - 1} (t)) - H_{y} (t, y_{q - 1} (t)) y_{q - 1} (t) .

\[\mathbf{E}_{q-1}(t)= -\mathbf{H}_{\mathbf{y}}(t,\mathbf{y}_{q-1}(t)), \qquad \mathbf{r}_{q-1}(t)= \mathbf{H}(t,\mathbf{y}_{q-1}(t)) - \mathbf{H}_{\mathbf{y}}(t,\mathbf{y}_{q-1}(t))\,\mathbf{y}_{q-1}(t). \]

where the state vector at the q-th QLM iteration is defined as

y_{q} (t) = (\begin{matrix} {(S_{T})}_{q} (t) \\ {(I_{T})}_{q} (t) \\ {(S_{L})}_{q} (t) \\ {(E_{L})}_{q} (t) \\ {(I_{L})}_{q} (t) \\ {(S_{H})}_{q} (t) \\ {(E_{H})}_{q} (t) \\ {(I_{H})}_{q} (t) \\ {(R_{H})}_{q} (t) \\ A_{q} (t) \end{matrix}), H (t, y_{q - 1}) = (\begin{matrix} Λ_{T} - λ_{T} (t) {(S_{T})}_{q - 1} (t) - μ_{T} {(S_{T})}_{q - 1} (t) \\ λ_{T} (t) {(S_{T})}_{q - 1} (t) - μ_{T} {(I_{T})}_{q - 1} (t) \\ Λ_{L} - λ_{L} (t) {(S_{L})}_{q - 1} (t) - μ_{L} {(S_{L})}_{q - 1} (t) \\ λ_{L} (t) {(S_{L})}_{q - 1} (t) - (σ_{L} + μ_{L}) {(E_{L})}_{q - 1} (t) \\ σ_{L} {(E_{L})}_{q - 1} (t) - (γ_{L} + μ_{L}) {(I_{L})}_{q - 1} (t) \\ Λ_{H} + ω_{R} {(R_{H})}_{q - 1} (t) - λ_{H} (t) {(S_{H})}_{q - 1} (t) - μ_{H} {(S_{H})}_{q - 1} (t) \\ λ_{H} (t) {(S_{H})}_{q - 1} (t) - (σ_{H} + μ_{H}) {(E_{H})}_{q - 1} (t) \\ σ_{H} {(E_{H})}_{q - 1} (t) - (γ_{H} + μ_{H} + δ_{H}) {(I_{H})}_{q - 1} (t) \\ γ_{H} {(I_{H})}_{q - 1} (t) - (ω_{R} + μ_{H}) {(R_{H})}_{q - 1} (t) \\ η_{1} \frac{{(I_{H})}_{q - 1} (t)}{{(N_{H})}_{q - 1} (t)} + η_{2} \frac{{(I_{T})}_{q - 1} (t)}{{(N_{T})}_{q - 1} (t)} - ω_{A} A_{q - 1} (t) \end{matrix})

\Lambda_T-\lambda_T(t)\,(S_T)_{q-1}(t)-\mu_T (S_T)_{q-1}(t)\\[2mm] \lambda_T(t)\,(S_T)_{q-1}(t)-\mu_T (I_T)_{q-1}(t)\\[2mm] \Lambda_L-\lambda_L(t)\,(S_L)_{q-1}(t)-\mu_L (S_L)_{q-1}(t)\\[2mm] \lambda_L(t)\,(S_L)_{q-1}(t)-(\sigma_L+\mu_L)(E_L)_{q-1}(t)\\[2mm] \sigma_L (E_L)_{q-1}(t)-(\gamma_L+\mu_L)(I_L)_{q-1}(t)\\[2mm] \Lambda_H+\omega_R (R_H)_{q-1}(t)-\lambda_H(t)\,(S_H)_{q-1}(t)-\mu_H (S_H)_{q-1}(t)\\[2mm] \lambda_H(t)\,(S_H)_{q-1}(t)-(\sigma_H+\mu_H)(E_H)_{q-1}(t)\\[2mm] \sigma_H (E_H)_{q-1}(t)-(\gamma_H+\mu_H+\delta_H)(I_H)_{q-1}(t)\\[2mm] \gamma_H (I_H)_{q-1}(t)-(\omega_R+\mu_H)(R_H)_{q-1}(t)\\[2mm] \eta_1\frac{(I_H)_{q-1}(t)}{(N_H)_{q-1}(t)} +\eta_2\frac{(I_T)_{q-1}(t)}{(N_T)_{q-1}(t)} -\omega_A A_{q-1}(t)

With the ordering (31), the Jacobian ℰ_q−1(t) is computed as the negative of the Jacobian H_y evaluated at the previous iterate. The family of linear systems (31) is supplemented with the initial condition

y^{(q)} (0) = y_{0} = {(S_{T} (0), I_{T} (0), S_{L} (0), E_{L} (0), I_{L} (0), S_{H} (0), E_{H} (0), I_{H} (0), R_{H} (0), A (0))}^{⊤} .

<tex-math><![CDATA[\mathcal{A}_C = 11.3412

which is consistent with the original model.

5.3

Fundamentals of the QLM-Chebyshev algorithm

Our main objective is to solve the sequence of linearized systems of IVPs arising from the QLM formulation on the time interval [0, t_b]. As described previously, the interval [0, t_b] is decomposed into U non-overlapping subdomains Q_u, and the resulting submodels are solved locally in a sequential manner. Therefore, in what follows, we only illustrate the proposed algorithm on a generic local subdomain Q_u, for u = 1, 2, ⋯, U. In each subdomain 𝒬_u, the approximate solution is collocated at the shifted Chebyshev-Gauss-Lobatto (CGL) nodes. These nodes are the mapped images of the extrema of the Chebyshev polynomial T_L(ξ) on [−1, 1] (i.e., the Gauss-Lobatto points, which include both endpoints), in contrast to the Gauss points which are the interior roots of T_L(ξ). On each subdomain 𝒬_u = [t_u−1, t_u], there are exactly L + 1 such nodes, defined explicitly as

t_{i}^{(u)} = \frac{t_{u} - t_{u - 1}}{2} (1 - \cos \frac{i π}{L}) + t_{u - 1}, i = 0, 1, \dots, L

t_i^{(u)} = \frac{t_u - t_{u-1}}{2} \left(1 - \cos\frac{i\pi}{L}\right) + t_{u-1}, \quad i = 0, 1, \cdots, L, \label{eq:CGL_nodes}

obtained by applying the affine transformation $t = \frac{t_{u} - t_{u - 1}}{2} (1 + ξ) + t_{u - 1}$ $t = \frac{t_u - t_{u-1}}{2}(1+\xi) + t_{u-1}$ to the standard CGL points $ξ_{i} = \cos (i π / L) \in [- 1, 1]$ $\xi_i = \cos(i\pi/L) \in[-1,1]$. Note that $t_{0}^{(u)} = t_{u - 1}$ $t_0^{(u)} = t_{u-1}$ and $t_{L}^{(u)} = t_{u}$ $t_L^{(u)} = t_u$ are the left and right endpoints of 𝒬_u, respectively, so these are closed Gauss-Lobatto nodes that include the subdomain boundaries. To proceed, we assume that the solution of the linearized model can be approximated by truncated ( L+1 )-term series expansions. Accordingly, the QLM-approximations at iteration q ≥ 1 on the subdomain Q_u are expressed as

\begin{matrix} {\begin{cases} {(S_{T})}_{L, u}^{q} (t) = \sum_{ℓ = 0}^{L} π_{ℓ, 1}^{q, u} {\tilde{T}}_{ℓ} (t) = {\tilde{T}}_{L} (t) Π_{L, 1}^{q, u}, \\ {(I_{T})}_{L, u}^{q} (t) = \sum_{ℓ = 0}^{L} π_{ℓ, 2}^{q, u} {\tilde{T}}_{ℓ} (t) = {\tilde{T}}_{L} (t) Π_{L, 2}^{q, u}, \\ {(S_{L})}_{L, u}^{q} (t) = \sum_{ℓ = 0}^{L} π_{ℓ, 3}^{q, u} {\tilde{T}}_{ℓ} (t) = {\tilde{T}}_{L} (t) Π_{L, 3}^{q, u}, \\ {(E_{L})}_{L, u}^{q} (t) = \sum_{ℓ = 0}^{L} π_{ℓ, 4}^{q, u} {\tilde{T}}_{ℓ} (t) = {\tilde{T}}_{L} (t) Π_{L, 4}^{q, u}, \\ {(I_{L})}_{L, u}^{q} (t) = \sum_{ℓ = 0}^{L} π_{ℓ, 5}^{q, u} {\tilde{T}}_{ℓ} (t) = {\tilde{T}}_{L} (t) Π_{L, 5}^{q, u}, \\ {(S_{H})}_{L, u}^{q} (t) = \sum_{ℓ = 0}^{L} π_{ℓ, 6}^{q, u} {\tilde{T}}_{ℓ} (t) = {\tilde{T}}_{L} (t) Π_{L, 6}^{q, u}, \\ {(E_{H})}_{L, u}^{q} (t) = \sum_{ℓ = 0}^{L} π_{ℓ, 7}^{q, u} {\tilde{T}}_{ℓ} (t) = {\tilde{T}}_{L} (t) Π_{L, 7}^{q, u}, \\ {(I_{H})}_{L, u}^{q} (t) = \sum_{ℓ = 0}^{L} π_{ℓ, 8}^{q, u} {\tilde{T}}_{ℓ} (t) = {\tilde{T}}_{L} (t) Π_{L, 8}^{q, u}, \\ {(R_{H})}_{L, u}^{q} (t) = \sum_{ℓ = 0}^{L} π_{ℓ, 9}^{q, u} {\tilde{T}}_{ℓ} (t) = {\tilde{T}}_{L} (t) Π_{L, 9}^{q, u}, \\ A_{L, u}^{q} (t) = \sum_{ℓ = 0}^{L} π_{ℓ, 10}^{q, u} {\tilde{T}}_{ℓ} (t) = {\tilde{T}}_{L} (t) Π_{L, 10}^{q, u} . \end{cases} \end{matrix}

(S_T)^{q}_{L,u}(t) =\displaystyle\sum_{\ell=0}^{L}\pi^{q,u}_{\ell,1}\,\widetilde{T}_{\ell}(t) =\widetilde{\mathbf{T}}_{L}(t)\,\boldsymbol{\Pi}^{\,q,u}_{L,1},\\[2mm] (I_T)^{q}_{L,u}(t) =\displaystyle\sum_{\ell=0}^{L}\pi^{q,u}_{\ell,2}\,\widetilde{T}_{\ell}(t) =\widetilde{\mathbf{T}}_{L}(t)\,\boldsymbol{\Pi}^{\,q,u}_{L,2},\\[2mm] (S_L)^{q}_{L,u}(t) =\displaystyle\sum_{\ell=0}^{L}\pi^{q,u}_{\ell,3}\,\widetilde{T}_{\ell}(t) =\widetilde{\mathbf{T}}_{L}(t)\,\boldsymbol{\Pi}^{\,q,u}_{L,3},\\[2mm] (E_L)^{q}_{L,u}(t) =\displaystyle\sum_{\ell=0}^{L}\pi^{q,u}_{\ell,4}\,\widetilde{T}_{\ell}(t) =\widetilde{\mathbf{T}}_{L}(t)\,\boldsymbol{\Pi}^{\,q,u}_{L,4},\\[2mm] (I_L)^{q}_{L,u}(t) =\displaystyle\sum_{\ell=0}^{L}\pi^{q,u}_{\ell,5}\,\widetilde{T}_{\ell}(t) =\widetilde{\mathbf{T}}_{L}(t)\,\boldsymbol{\Pi}^{\,q,u}_{L,5},\\[2mm] (S_H)^{q}_{L,u}(t) =\displaystyle\sum_{\ell=0}^{L}\pi^{q,u}_{\ell,6}\,\widetilde{T}_{\ell}(t) =\widetilde{\mathbf{T}}_{L}(t)\,\boldsymbol{\Pi}^{\,q,u}_{L,6},\\[2mm] (E_H)^{q}_{L,u}(t) =\displaystyle\sum_{\ell=0}^{L}\pi^{q,u}_{\ell,7}\,\widetilde{T}_{\ell}(t) =\widetilde{\mathbf{T}}_{L}(t)\,\boldsymbol{\Pi}^{\,q,u}_{L,7},\\[2mm] (I_H)^{q}_{L,u}(t) =\displaystyle\sum_{\ell=0}^{L}\pi^{q,u}_{\ell,8}\,\widetilde{T}_{\ell}(t) =\widetilde{\mathbf{T}}_{L}(t)\,\boldsymbol{\Pi}^{\,q,u}_{L,8},\\[2mm] (R_H)^{q}_{L,u}(t) =\displaystyle\sum_{\ell=0}^{L}\pi^{q,u}_{\ell,9}\,\widetilde{T}_{\ell}(t) =\widetilde{\mathbf{T}}_{L}(t)\,\boldsymbol{\Pi}^{\,q,u}_{L,9},\\[2mm] A^{q}_{L,u}(t) =\displaystyle\sum_{\ell=0}^{L}\pi^{q,u}_{\ell,10}\,\widetilde{T}_{\ell}(t) =\widetilde{\mathbf{T}}_{L}(t)\,\boldsymbol{\Pi}^{\,q,u}_{L,10}. \end{cases}

For t ∈ Q_u, let

Π_{L, j}^{q, u} = {(π_{0, j}^{q, u} π_{1, j}^{q, u} \dots π_{L, j}^{q, u})}^{T}, j = 1, 2, \dots, 10,

\[ \boldsymbol{\Pi}^{\,q,u}_{L,j} = \begin{pmatrix} \pi^{q,u}_{0,j} & \pi^{q,u}_{1,j} & \cdots & \pi^{q,u}_{L,j} \end{pmatrix}^{T}, \qquad j=1,2,\cdots,10, \]

denote the vectors of unknown Chebyshev coefficients corresponding to the approximate solutions of the state variables. By using the representation given for the shifted Chebyshev basis vector ${\tilde{T}}_{L} (t)$ $\widetilde{\mathbf{T}}_{L}(t)$, the relations in (34) can be recast as

\begin{matrix} \begin{matrix} {(S_{T})}_{L, u}^{q} (t) = {\tilde{T}}_{L} (t) D_{L} Π_{L, 1}^{q, u}, \\ {(S_{L})}_{L, u}^{q} (t) = {\tilde{T}}_{L} (t) D_{L} Π_{L, 3}^{q, u}, \\ {(I_{L})}_{L, u}^{q} (t) = {\tilde{T}}_{L} (t) D_{L} Π_{L, 5}^{q, u}, \\ {(E_{H})}_{L, u}^{q} (t) = {\tilde{T}}_{L} (t) D_{L} Π_{L, 7}^{q, u}, \\ {(R_{H})}_{L, u}^{q} (t) = {\tilde{T}}_{L} (t) D_{L} Π_{L, 9}^{q, u}, \end{matrix} & \begin{matrix} {(I_{T})}_{L, u}^{q} (t) = {\tilde{T}}_{L} (t) D_{L} Π_{L, 2}^{q, u}, \\ {(E_{L})}_{L, u}^{q} (t) = {\tilde{T}}_{L} (t) D_{L} Π_{L, 4}^{q, u}, \\ {(S_{H})}_{L, u}^{q} (t) = {\tilde{T}}_{L} (t) D_{L} Π_{L, 6}^{q, u}, \\ {(I_{H})}_{L, u}^{q} (t) = {\tilde{T}}_{L} (t) D_{L} Π_{L, 8}^{q, u}, \\ A_{L, u}^{q} (t) = {\tilde{T}}_{L} (t) D_{L} Π_{L, 10}^{q, u}, \end{matrix} & t \in Q_{u}, \end{matrix}

\left\{\begin{array}{ccc}\begin{array}{c}{(S_T)}_{L,u}^q(t)={\widetilde{\textbf{T}}}_L(t)D_L\mathbf\Pi_{L,1}^{q,u},\\{(S_L)}_{L,u}^q(t)={\widetilde{\textbf{T}}}_L(t)D_L\mathbf\Pi_{L,3}^{q,u},\\{(I_L)}_{L,u}^q(t)={\widetilde{\textbf{T}}}_L(t)D_L\mathbf\Pi_{L,5}^{q,u},\\{(E_H)}_{L,u}^q(t)={\widetilde{\textbf{T}}}_L(t)D_L\mathbf\Pi_{L,7}^{q,u},\\{({\mathcal R}_H)}_{L,u}^q(t)={\widetilde{\textbf{T}}}_L(t)D_L\mathbf\Pi_{L,9}^{q,u},\end{array}&\begin{array}{c}{(I_T)}_{L,u}^q(t)={\widetilde{\textbf{T}}}_L(t)D_L\mathbf\Pi_{L,2}^{q,u},\\{(E_L)}_{L,u}^q(t)={\widetilde{\textbf{T}}}_L(t)D_L\mathbf\Pi_{L,4}^{q,u},\\{(S_H)}_{L,u}^q(t)={\widetilde{\textbf{T}}}_L(t)D_L\mathbf\Pi_{L,6}^{q,u},\\{(I_H)}_{L,u}^q(t)={\widetilde{\textbf{T}}}_L(t)D_L\mathbf\Pi_{L,8}^{q,u},\\A_{L,u}^q(t)={\widetilde{\textbf{T}}}_L(t)D_L\mathbf\Pi_{L,10}^{q,u},\end{array}&t\;\in\;Q_u;\end{array}\right.

for q = 1, 2, ⋯ and u = 1, 2, ⋯, U.

Since the considered CCHF model is of integer order, the temporal derivatives appearing in the QLM linearized system are classical first-order derivatives. Differentiating the Chebyshev approximation in (35) yields

\frac{d}{d t} y_{L}^{q, u} (t) = \frac{d}{d t} {\tilde{T}}_{L} (t) Π_{L}^{q, u} = {\tilde{T}}_{L} (t) D_{L}^{(1)} Π_{L}^{q, u}, t \in Q_{u},

\frac{d}{dt}\mathbf{y}^{\,q,u}_{L}(t) = \frac{d}{dt}\widetilde{\mathbf{T}}_{L}(t)\,\boldsymbol{\Pi}^{\,q,u}_{L} = \widetilde{\mathbf{T}}_{L}(t)\,\mathbf{D}^{(1)}_{L}\, \boldsymbol{\Pi}^{\,q,u}_{L}, \qquad t\in Q_u,

where $D_{L}^{(1)}$ $\mathbf{D}^{(1)}_{L}$ denotes the first-order shifted Chebyshev operational differentiation matrix.

Substituting the approximations (35) and their derivatives (36) into the QLM linearized system leads to

{\tilde{T}}_{L} (t) D_{L}^{(1)} Π_{L}^{q, u} = E_{q - 1} (t) {\tilde{T}}_{L} (t) Π_{L}^{q, u} + ρ_{q - 1} (t), t \in Q_{u} .

\widetilde{\mathbf{T}}_{L}(t)\,\mathbf{D}^{(1)}_{L}\, \boldsymbol{\Pi}^{\,q,u}_{L} = \mathcal{E}_{q-1}(t)\, \widetilde{\mathbf{T}}_{L}(t)\,\boldsymbol{\Pi}^{\,q,u}_{L} + \boldsymbol{\rho}_{q-1}(t), \qquad t\in Q_u.

Let ${t_{i}^{(u)}}_{i = 0}^{L}$ $\{t^{(u)}_i\}_{i=0}^{L}$ denote the shifted Chebyshev-Gauss-Lobatto collocation points on the subdomain Q_u. By enforcing (37) at these nodes, we obtain the following linear algebraic system:

{\tilde{T}}_{L} (t_{i}^{(u)}) D_{L}^{(1)} Π_{L}^{q, u} = E_{q - 1} (t_{i}^{(u)}) {\tilde{T}}_{L} (t_{i}^{(u)}) Π_{L}^{q, u} + ρ_{q - 1} (t_{i}^{(u)}), i = 0, 1, \dots, L

\mathcal{C}^{\,u}_L := \left\{ t_{u,l} = \frac{t_u-t_{u-1}}{2}\left(1+\cos\frac{l\pi}{L}\right)+t_{u-1} \,:\, l=0,1,\cdots,L \right\}. \label{eq:collocation_points}

On the first subdomain Q₁, the initial condition of the original model is imposed as

y_{L}^{q, 1} (0) = y_{0} .

\mathbf{y}^{\,q,1}_{L}(0)=\mathbf{y}_0.

For u ≥ 2, continuity across adjacent subdomains is enforced through

y_{L}^{q, u} (t_{u - 1}) = y_{L}^{q, u - 1} (t_{u - 1}),

\mathbf{y}^{\,q,u}_{L}(t_{u-1}) = \mathbf{y}^{\,q,u-1}_{L}(t_{u-1}),

which guarantees a globally continuous approximate solution over the entire interval [0, t_b]. The algorithm of the proposed method is summarized in Algorithm 1.

5.4

Residual error functions for testing accuracy

The last part of this section is devoted to the definition of the residual error functions (REFs) associated with the obtained approximations generated by the proposed QLM-spectral collocation scheme. These residuals provide an effective tool for testing the accuracy of the numerical approximations in the absence of exact analytical solutions. The main idea is to substitute the approximate solutions, which satisfy the linearized system at each iteration, into the original nonlinear CCHF model and evaluate the magnitude of the resulting residuals. Let $(\cdot)_{L, u}^{q} (t)$ \gamma_i denote the approximate solutions obtained at the q-th QLM iteration using polynomial degree L. The residual error functions corresponding to the human, livestock, tick, and awareness compartments are defined as follows:

Algorithm 1

QLM-Shifted Chebyshev Collocation with Domain Decomposition.

Partition [0, t_b] into U subdomains Q_u = [t_u-1, t_u], u = 1, ... , U.

For each Q_u, generate the shifted Chebyshev-Gauss-Lobatto nodes ${_{t_{i}}^{(u)}}_{i = 0}^{L}$ \delta_i.

Construct the shifted Chebyshev basis vector $\tilde{T_{L}} (t)$ \epsilon_i and the differentiation matrix D⁽¹⁾_L.

Initialize the QLM iterate on each subdomain and set ${y_{L}}^{q, 1} (t_{0}) = y_{0} .$ \zeta_i

foru = 1 toUdo

ifu ≥ 2 then

Enforce continuity: ${y_{L}}^{q, u} (t_{u - 1}) = {y_{L}}^{q, u - 1} (t_{u - 1})$ \eta_i.

end if

forq = 1 to q_maxdo

Evaluate $E_{q - 1} (t^{u} i_{})$ \theta_i and $ρ_{q - 1} (t_{i}^{(u)})$ \iota_i at all L + 1 shifted Chebyshev-Gauss-Lobatto collocation nodes ${t_{i}^{(u)}}_{i = 0}^{L}$ \kappa_i defined in (33), where $t_{i}^{(u)} = \frac{t_{u} - t_{u - 1}}{2} (1 - \cos \frac{i π}{L}) + t_{u - 1}$ \lambda_i are the L + 1 extremal nodes of the shifted Chebyshev polynomial on ℚu, including both endpoints t_u-1 and t_u.

Assume $y_{L}^{q, u} (t) = {\tilde{T}}_{L} (t) Π_{L}^{q, u}$ \mu_i for t ∈ Q_u.

Collocate the QLM linear system at ${t_{i}^{(u)}}_{i = 0}^{L}$ \nu_i and impose the left-end constraint at t = t_{u • 1}.

Solve for $Π_{L}^{q, u}$ \xi_i and update $y_{L}^{q, u} (t)$ o_i.

if ${max}_{0 \leq i \leq L} {‖ y_{L}^{q, u} (t_{i}^{(u)}) - y_{L}^{q - 1, u} (t_{i}^{(u)}) ‖}_{\infty} \leq Tol$ \pi_ithen

break

end if

end for

Return the assembled approximation y^q(t) on [0, t_b].

{RE}_{1}^{(q, u)} (t) := \frac{d}{d t} S_{H, L, u}^{(q)} (t) - (Λ_{H} + ω_{R} R_{H, L, u}^{(q)} - λ_{H, L, u}^{(q)} S_{H, L, u}^{(q)} - μ_{H} S_{H, L, u}^{(q)}),

\mathrm{RE}^{(q,u)}_{1}(t)&:= \frac{d}{dt} S^{(q)}_{H,L,u}(t) -\Big( \Lambda_H + \omega_R R^{(q)}_{H,L,u} -\lambda^{(q)}_{H,L,u} S^{(q)}_{H,L,u} -\mu_H S^{(q)}_{H,L,u} \Big), \\[4pt] \mathrm{RE}^{(q,u)}_{2}(t) &:= \frac{d}{dt} E^{(q)}_{H,L,u}(t) -\Big( \lambda^{(q)}_{H,L,u} S^{(q)}_{H,L,u} -(\sigma_H+\mu_H)E^{(q)}_{H,L,u} \Big), \\[4pt] \mathrm{RE}^{(q,u)}_{3}(t)&:= \frac{d}{dt} I^{(q)}_{H,L,u}(t) -\Big( \sigma_H E^{(q)}_{H,L,u} -(\gamma_H+\mu_H+\delta_H)I^{(q)}_{H,L,u} \Big), \\[4pt] \mathrm{RE}^{(q,u)}_{4}(t) &:= \frac{d}{dt} R^{(q)}_{H,L,u}(t) -\Big( \gamma_H I^{(q)}_{H,L,u} -(\omega_R+\mu_H)R^{(q)}_{H,L,u} \Big), \\[4pt] \mathrm{RE}^{(q,u)}_{5}(t)&:= \frac{d}{dt} S^{(q)}_{L,L,u}(t) -\Big( \Lambda_L -\lambda^{(q)}_{L,L,u} S^{(q)}_{L,L,u} -\mu_L S^{(q)}_{L,L,u} \Big), \\[4pt] \mathrm{RE}^{(q,u)}_{6}(t) &:= \frac{d}{dt} E^{(q)}_{L,L,u}(t) -\Big( \lambda^{(q)}_{L,L,u} S^{(q)}_{L,L,u} -(\sigma_L+\mu_L)E^{(q)}_{L,L,u} \Big), \\[4pt] \mathrm{RE}^{(q,u)}_{7}(t) &:= \frac{d}{dt} I^{(q)}_{L,L,u}(t) -\Big( \sigma_L E^{(q)}_{L,L,u} -(\gamma_L+\mu_L)I^{(q)}_{L,L,u} \Big), \\[6pt] \mathrm{RE}^{(q,u)}_{8}(t) &:= \frac{d}{dt} S^{(q)}_{T,L,u}(t) -\Big( \Lambda_T -\lambda^{(q)}_{T,L,u} S^{(q)}_{T,L,u} -\mu_T S^{(q)}_{T,L,u} \Big),

{RE}_{2}^{(q, u)} (t) := \frac{d}{d t} E_{H, L, u}^{(q)} (t) - (λ_{H, L, u}^{(q)} S_{H, L, u}^{(q)} - (σ_{H} + μ_{H}) E_{H, L, u}^{(q)}),

\mathrm{RE}^{(q,u)}_{9}(t) &:= \frac{d}{dt} I^{(q)}_{T,L,u}(t) -\Big( \lambda^{(q)}_{T,L,u} S^{(q)}_{T,L,u} -\mu_T I^{(q)}_{T,L,u} \Big), \\[6pt] \mathrm{RE}^{(q,u)}_{10}(t) &:= \frac{d}{dt} A^{(q)}_{L,u}(t) -\Bigg( \eta_1\frac{I^{(q)}_{H,L,u}}{N^{(q)}_{H,L,u}} +\eta_2\frac{I^{(q)}_{T,L,u}}{N^{(q)}_{T,L,u}} -\omega_A A^{(q)}_{L,u} \Bigg).

{RE}_{3}^{(q, u)} (t) := \frac{d}{d t} I_{H, L, u}^{(q)} (t) - (σ_{H} E_{H, L, u}^{(q)} - (γ_{H} + μ_{H} + δ_{H}) I_{H, L, u}^{(q)}),

\mathcal{I}_M = 35.4215

{RE}_{4}^{(q, u)} (t) := \frac{d}{d t} R_{H, L, u}^{(q)} (t) - (γ_{H} I_{H, L, u}^{(q)} - (ω_{R} + μ_{H}) R_{H, L, u}^{(q)}),

{RE}_{5}^{(q, u)} (t) := \frac{d}{d t} S_{L, L, u}^{(q)} (t) - (Λ_{L} - λ_{L, L, u}^{(q)} S_{L, L, u}^{(q)} - μ_{L} S_{L, L, u}^{(q)}),

{RE}_{6}^{(q, u)} (t) := \frac{d}{d t} E_{L, L, u}^{(q)} (t) - (λ_{L, L, u}^{(q)} S_{L, L, u}^{(q)} - (σ_{L} + μ_{L}) E_{L, L, u}^{(q)}),

{RE}_{7}^{(q, u)} (t) := \frac{d}{d t} I_{L, L, u}^{(q)} (t) - (σ_{L} E_{L, L, u}^{(q)} - (γ_{L} + μ_{L}) I_{L, L, u}^{(q)}),

{RE}_{8}^{(q, u)} (t) := \frac{d}{d t} S_{T, L, u}^{(q)} (t) - (Λ_{T} - λ_{T, L, u}^{(q)} S_{T, L, u}^{(q)} - μ_{T} S_{T, L, u}^{(q)}),

{RE}_{9}^{(q, u)} (t) := \frac{d}{d t} I_{T, L, u}^{(q)} (t) - (λ_{T, L, u}^{(q)} S_{T, L, u}^{(q)} - μ_{T} I_{T, L, u}^{(q)}),

{RE}_{10}^{(q, u)} (t) := \frac{d}{d t} A_{L, u}^{(q)} (t) - (η_{1} \frac{I_{H, L, u}^{(q)}}{N_{H, L, u}^{(q)}} + η_{2} \frac{I_{T, L, u}^{(q)}}{N_{T, L, u}^{(q)}} - ω_{A} A_{L, u}^{(q)}) .

Results, simulations and physical interpretation

In this section, we present detailed numerical simulations to illustrate the dynamical behavior of the proposed CCHF transmission model and to assess the accuracy and convergence of the adopted numerical scheme. In addition, the simulations are epidemiologically motivated by real-world data reported for Iraq, a high-incidence region for CCHF. The consistency of the adopted parameters with reported epidemiological characteristics is discussed and validated in the subsequent subsections.

6.1

Epidemiological context and parameter justification

The numerical simulations presented in this section are epidemiologically motivated by reported characteristics of CCHF transmission in high-incidence regions. In particular, Iraq is considered as a representative setting, as it has experienced recurrent CCHF outbreaks in recent years with a substantial number of laboratory-confirmed cases reported by the World Health Organization (WHO) [29,30]. According to WHO situation reports, CCHF transmission in Iraq is predominantly driven by tick-livestock-human interactions, where most human infections arise from direct contact with infected ticks or livestock, while human-to-human transmission occurs less frequently and is mainly associated with close contact in household or healthcare settings [29]. This epidemiological structure directly supports the modeling assumption that tick-mediated transmission constitutes the dominant pathway, followed by livestock-to-human transmission, with comparatively weaker but non-negligible human-to-human transmission. Clinical evidence further indicates that the incubation period of CCHF in humans typically ranges from several days up to approximately two weeks, depending on the route of exposure, while the infectious period generally spans one to two weeks [30]. Accordingly, the progression and recovery rates adopted in the simulations are selected to reflect these reported clinical time scales. These parameters are not intended to reproduce exact case-level trajectories, but rather to capture realistic average disease dynamics at the population level. Demographic parameters, including natural mortality and recruitment rates for ticks, livestock, and humans, are chosen to reflect realistic life expectancies and population turnover on a daily time scale. Recruitment rates are defined proportionally to the corresponding mortality rates in order to maintain approximately constant total population sizes over the simulation horizon. This assumption is reasonable in the Iraq context, as reported CCHF outbreaks are typically seasonal and occur over relatively short to medium time intervals that do not significantly alter overall population sizes [29]. Transmission coefficients are selected to reproduce the relative importance of different infection pathways rather than their exact magnitudes. In particular, higher values are assigned to tick-livestock transmission to reflect the enzootic cycle, while lower values are used for livestock-human and human-to-human transmission in accordance with epidemiological observations. Awareness-related parameters are introduced to represent the impact of public awareness, behavioral change, and preventive measures emphasized in WHO guidelines, allowing the model to assess how awareness-driven responses influence disease spread.

Based on these considerations, Table 3 summarizes the parameter values and initial conditions employed in the simulations. It is important to emphasize that the present study does not aim to perform country-specific parameter estimation or direct fitting to reported case data. Instead, the adopted parameter set is chosen to remain consistent with epidemiological ranges reported for high-burden regions such as Iraq, ensuring that the simulations provide a realistic qualitative representation of CCHF transmission dynamics while enabling systematic investigation of awareness-driven behavioral effects. In addition, the rate of loss of immunity ω_R = 0.05 day⁻¹ is adopted to reflect the short-lived nature of post-infection immunity in CCHF survivors, consistent with clinical observations indicating that long-term protective immunity is not reliably established following infection [29,30]. This value corresponds to an average immunity duration of approximately 20 days, after which recovered individuals return to the susceptible class and may be reinfected upon renewed exposure.

Table 3

Initial conditions and parameter values for model 1.

Parameter	Value	Parameter	Value
$N_{T}^{0}$ \rho_i	100000	$N_{L}^{0}$ \sigma_i	2000
$N_{H}^{0}$ \tau_i	10000	μ_T	0.0027
μ_L	5.48 × 10⁻⁰⁴	μ_H	3.91 × 10⁻⁰⁵
Λ_T	54.79	Λ_L	1.095
Λ_H	0.39	σ_H	1/5
γ_H	1/8	δ_H	0.02
ω_R	0.05	σ_L	1/4
γ_L	1/5	β_TL	0.6
β_LT	0.4	β_TH	0.15
β_LH	0.10	β_HH	0.05
k_A	5.0	η₁	2.0
η₂	0.5	ω_A	0.10

Remark 3. The weekly CCHF infection data for Iraq presented in Figures 2 and 3 are included to provide epidemiological context and to motivate the modeling framework, rather than to serve as a calibration dataset for direct parameter estimation. Direct fitting of the proposed model to raw surveillance data is precluded by several well-recognized challenges: significant underreporting due to limited diagnostic capacity in affected Iraq governorates [29,30]; structural non-identifiability of the tick-livestock-human parameter space from human case data alone; and the unavailability of concurrent entomological and behavioral datasets required for joint inference of awareness-related parameters such as η₁, η₂, and k_A. Instead, parameter values are drawn from the published CCHF literature [19-21] and chosen to be consistent with the epidemiological characteristics of CCHF in Iraq, with the reported data serving as a qualitative reference to confirm the biological plausibility of the simulated outbreak dynamics. A rigorous data-driven calibration incorporating Bayesian inference or nonlinear least-squares fitting, with appropriate correction for underreporting, is left as an important direction for future work.

6.2

Numerical results and validation

The simulations are designed to demonstrate the time evolution of all state variables, the influence of awareness on disease transmission, and the numerical reliability of the spectral collocation method. The results are summarized through graphical illustrations (Figures 4-12) and quantitative error analyses (Tables 4 and 5). First, Figures 4-9 illustrate the temporal evolution of the tick and livestock populations, including the susceptible, exposed, and infected compartments, together with the absolute error computed by comparison with the reference MATLAB solver ODE45. It can be observed that the solutions exhibit smooth and biologically consistent trajectories, which confirms the well-posedness of the proposed model. In particular, the infected tick and livestock populations initially increase due to transmission interactions and subsequently decline as recovery processes and awareness effects become dominant. The long-term behavior of the solutions indicates convergence toward stable equilibrium states, reflecting realistic disease dynamics under sustained control measures. In addition, Figures 7-8 depict the evolution of the human population compartments together with the awareness variable. The results clearly demonstrate that awareness plays a significant role in mitigating disease transmission. As awareness increases, a noticeable reduction in the infected human population is observed, highlighting the effectiveness of awareness-driven behavioral changes. The recovered human population increases accordingly, while disease-induced mortality remains bounded, confirming the stabilizing influence of awareness mechanisms on the overall system dynamics. These figures collectively demonstrate the strong coupling between epidemiological states and awareness dynamics and emphasize the importance of non-pharmaceutical interventions in controlling the spread of CCHF. Furthermore, a focused analysis of the dynamics of both S_H and I_H is illustrated in Figure 9, which demonstrates the effectiveness of the QLM technique in achieving accurate numerical simulations for these two state variables. Next, to quantitatively assess the numerical accuracy of the proposed spectral collocation method, absolute error norms for all state variables are reported in Figures (10-12) and Table 4 for different polynomial degrees L. The results demonstrate a rapid decay of the error as L increases, with several orders of magnitude reduction observed when transitioning from low to moderate values of L. This behavior confirms the high-order accuracy of the proposed numerical method. Table 5 presents the corresponding error analysis for a second simulation scenario. Similar convergence patterns are observed across all state variables, with the errors approaching near machine-precision levels for sufficiently large values of L. The consistency of the convergence behavior across different scenarios highlights the robustness and stability of the numerical scheme. The reported results confirm that the adopted numerical approach is highly efficient and well suited for solving the proposed nonlinear, multi-compartment epidemiological model. Moreover, the results emphasize the critical role of awareness in reducing infection prevalence and validate the spectral collocation method as a powerful tool for simulating complex epidemic models involving multiple interacting populations.

Table 4

Error norms for different state variables at different values of L.

Variable	L = 8	L = 16	L = 32	L = 64
S_T	3.942 × 10⁻⁵	7.661 × 10⁻⁷	2.916 × 10⁻¹⁰	6.411 × 10⁻¹⁰
I_T	3.942 × 10⁻⁵	7.661 × 10⁻⁷	2.849 × 10⁻¹⁰	6.834 × 10⁻¹⁰
S_L	4.536 × 10⁻⁵	9.096 × 10⁻⁷	4.192 × 10⁻¹⁰	7.547 × 10⁻¹⁰
E_L	1.563 × 10⁻⁵	3.058 × 10⁻⁷	4.842 × 10⁻⁹	4.842 × 10⁻⁹
I_L	1.409 × 10⁻⁵	2.827 × 10⁻⁷	5.817 × 10⁻¹⁰	1.502 × 10⁻⁹
S_H	6.783 × 10⁻⁶	1.421 × 10⁻⁷	3.452 × 10⁻⁹	9.442 × 10⁻⁹
E_H	2.558 × 10⁻⁶	4.950 × 10⁻⁸	2.719 × 10⁻¹⁰	2.719 × 10⁻¹⁰
I_H	2.538 × 10⁻⁶	5.004 × 10⁻⁸	5.485 × 10⁻¹¹	1.067 × 10⁻⁹
R_H	2.977 × 10⁻⁶	9.368 × 10⁻⁸	1.921 × 10⁻¹⁰	4.540 × 10⁻¹⁰
A	1.584 × 10⁻⁴	3.094 × 10⁻⁶	1.021 × 10⁻⁸	1.021 × 10⁻⁸

Table 5

Residual Error norms for different state variables at increasing values of L.

Variable	L = 8	L = 16	L = 32	L = 64
S_T	1.100 × 10⁻⁵	2.245 × 10⁻⁷	5.839 × 10⁻¹²	5.843 × 10⁻¹²
I_T	1.100 × 10⁻⁵	2.245 × 10⁻⁷	5.838 × 10⁻¹²	5.840 × 10⁻¹²
S_L	5.670 × 10⁻⁵	1.490 × 10⁻⁶	5.841 × 10⁻¹²	5.843 × 10⁻¹²
E_L	3.768 × 10⁻⁵	1.040 × 10⁻⁶	5.842 × 10⁻¹²	5.843 × 10⁻¹²
I_L	1.631 × 10⁻⁵	3.916 × 10⁻⁷	5.842 × 10⁻¹²	5.840 × 10⁻¹²
S_H	1.250 × 10⁻⁶	2.846 × 10⁻⁸	5.844 × 10⁻¹²	5.841 × 10⁻¹²
E_H	1.165 × 10⁻⁶	2.007 × 10⁻⁸	5.843 × 10⁻¹²	5.841 × 10⁻¹²
I_H	7.326 × 10⁻⁸	8.928 × 10⁻⁹	5.842 × 10⁻¹²	5.893 × 10⁻¹²
R_H	9.888 × 10⁻⁹	4.579 × 10⁻¹⁰	5.840 × 10⁻¹²	5.846 × 10⁻¹²
A	2.373 × 10⁻⁶	4.574 × 10⁻⁸	5.857 × 10⁻¹²	5.897 × 10⁻¹²

Conclusions

In this work, an awareness-driven deterministic compartmental model was developed to investigate the transmission dynamics of CCHF across coupled tick-livestock-human interactions. The well-posedness of the model was rigorously established by proving positivity and uniform boundedness of solutions within a biologically feasible invariant region. The basic reproduction number ℛ₀ was derived via the next-generation matrix approach, with explicit computation of the 5 × 5 next-generation matrix K = FV⁻¹ and its spectral radius, yielding a threshold that decouples into the enzootic tick-livestock cycle and the human-to-human transmission chain. Local asymptotic stability of the disease-free equilibrium was shown to be governed by this threshold in the standard sense: the equilibrium is stable if and only if ℛ₀ < 1. A high-order numerical scheme was proposed by coupling the QLM with Chebyshev spectral collocation and a domain decomposition strategy, and was shown to achieve spectral accuracy for the resulting nonlinear multi-compartment system over extended time horizons. The scheme was validated against the reference MATLAB solver ode45, with absolute pointwise errors consistently in the range 10⁻¹⁰ − 10⁻¹⁴ across all state variables, confirming both the correctness of the implementation and the high-order convergence of the method. Numerical simulations confirmed the theoretical results and demonstrated that, although public awareness does not alter ℛ₀ directly, it significantly reduces outbreak magnitude, lowers peak infection prevalence, and delays disease progression in both the human and livestock populations. These findings highlight the role of awareness-driven behavioural responses as an effective non-pharmaceutical intervention, particularly in high-incidence, resource-limited settings such as Iraq where pharmaceutical options remain constrained.

Modeling Crimean-Congo Hemorrhagic fever with behavioral awareness: Mathematical analysis via Chebyshev spectral collocation solutions

Full Article

Paradigm

My account