EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 4, Article Number 6749 ISSN 1307-5543 – ejpam.com Published by New York Business Global Mathematical Modelling of Radicalization and Terrorism Dynamics Malicki Zorom1,∗, Babacar Leye1, Mamadou Diop2, Serigne M’backé Coly1, Abdou Lawane Gana2, Mäımouna Bologo/Traore1, Dial Niang1 1 Laboratoire Eaux Hydro-Systèmes et Agriculture (LEHSA), Institut International d’Ingénierie de l’Eau et de l’Environnement (2iE), Ouagadougou, Centre Region, Burkina Faso 2 Laboratoire EcoMatériaux et Habitats Durables (LEMHaD), Institut International d’Ingénierie de l’Eau et de l’Environnement (2iE), Ouagadougou, Centre Region, Burkina Faso Abstract. The Central Sahel faces a severe terrorism crisis fueled by violent extremism, displace- ment, and climatic shocks. This study develops a compartmental mathematical model to analyze the dynamics between susceptible populations, active terrorists, and internally displaced persons. Through stability analysis, bifurcation theory, and global sensitivity analysis, we demonstrate that the basic reproduction number,R0, is a critical threshold determining whether terrorism is elimi- nated or persists endemically. Our results show that military-only interventions are less than 20% effective, while integrated strategies combining prevention and deradicalization exceed 80% effec- tiveness. Time-dependent analysis reveals that optimal strategies must adapt from early preven- tion to long-term rehabilitation. These findings provide quantitative support for counter-terrorism frameworks prioritizing socioeconomic development over purely military solutions, offering a path- way to sustainable stability in the Sahel. 2020 Mathematics Subject Classifications: 34D20, 34D23, 37N25, 92D30 Key Words and Phrases: terrorism dynamics, mathematical modeling, sensitivity analysis, counter-terrorism strategies, Central Sahel 1. Introduction Recent decades have witnessed a marked escalation in the sophistication and impact of terrorism. The emergence and entrenchment of international terrorist networks represent a critical development, enabling coordinated attacks with heightened destructive poten- tial aimed at destabilizing governments and challenging fundamental democratic norms [1]. The transnational nature of contemporary terrorism is evident, yet empirical data ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i4.6749 Email addresses: malicki.zorom@2ie-edu.org (M. Zorom), babacar.leye@2ie-edu.org (B. Leye), mamadou.diop@2ie-edu.org (M. Diop), mbacke.coly@2ie-edu.org (S.M. Coly), abdou.lawane@2ie-edu.org (A.L. Gana), maimouna.bologo@2ie-edu.org (M. Bologo/Traore), dial.niang@2ie-edu.org (D. Niang) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) M. Zorom et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6749 2 of 47 reveals its concentrated impact in the Sahel. According to the Global Terrorism Index, Burkina Faso, Nigeria, Mali, and Niger rank 4th, 6th, 7th, and 8th worldwide, respectively, positioning them directly behind Afghanistan, Iraq, and Somalia [2]. Characterized as an Islamist insurgency, the Sahel conflict pits state forces of Mali, Niger, Mauritania, Burkina Faso, and Chad against Salafi-jihadist groups operating under the ideological and/or operational banners of Al-Qaeda or the Islamic State [3–8]. The Sahel insurgency constitutes an indirect regional consequence of the Algerian civil war. Algerian Islamist rebels strategically exploited the Sahelian desert as a rear base from the early 2000s onward [9]. Their operational profile gradually expanded beyond sanctuary provision to encompass guerrilla tactics, terrorist acts, and kidnappings. A pivotal development, however, was their deliberate embedding within local populations and propagation of radical Islamist doctrine. This process fostered local recruitment and led to the genesis of new, distinctly localized militant movements, including Ansar Dine, MUJAO, and Katiba Macina [10]. JNIM (Group for the Support of Islam and Muslims, GSIM) accounts for more than 64% of Sahelian militant Islamist violence since 2017, with activities documented from northern Mali to southeastern Burkina Faso. The Macina Liberation Front (MLF) emerges as JNIM’s most active component, estimated to perpetrate 75% of its violence. Based in central Mali and extending into Burkina Faso, the MLF’s prominence underscores a critical characteristic of JNIM factions: their lack of broad popular legitimacy. Consequently, these groups increasingly exploit local criminal networks and engage in attacks against civilians, a tactic notably employed by the MLF [11]. Empirical data reveals a dramatic near-sevenfold increase in violent events linked to Sahelian militant Islamist groups since 2017. With more than 1,000 incidents documented in the preceding year, the Sahel witnessed the sharpest rise in extremist violence across Africa. The resultant humanitarian and societal toll is severe: an estimated 8,000 deaths, millions displaced, pervasive attacks on governance structures and traditional authorities, the shuttering of thousands of educational institutions, and substantial economic decline [11]. The neutralization of AQIM leader Abdelmalek Droukdel by French forces on June 3, 2020 [12], highlighted the persistent insecurity plaguing the Sahel. Paradoxically, despite significant security sector investments [13], violence intensified markedly in 2019. The Central Sahel (Mali, Burkina Faso, Niger) recorded approximately 4,000 conflict-related fatalities representing a fivefold increase from the 770 deaths documented in 2016 [14]. Burkina Faso experienced the most acute deterioration, with militant Islamist attacks surging 174% between 2018 and 2019 [15], culminating in 1,889 fatalities during its dead- liest year on record [15]. This erosion of state authority manifests in cascading regional crises: the closure of over 1,800 schools [16], mass population displacement [17], and the expansion of ungoverned spaces (”grey zones”) where state control is absent or mediated through non-state armed actors [18]. The landscape of non-state armed groups in northern Mali has complexified substan- tially since 2012, growing quantitatively (from 4 to 15 groups) and qualitatively. Three interconnected drivers explain this evolution: First, fission processes among jihadist orga- M. Zorom et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6749 3 of 47 nizations (scissiparity). Second, the endogenous proliferation of community-based militias, arising from the national security apparatus’s failure to protect civilians and necessitated by localized security dilemmas vis-à-vis rival militias. Third, the deployment of multi- lateral (MINUSMA) and unilateral (French Barkhane/Sabre, Chadian Serval contingent) external military forces. Figure 1: Geospatial Analysis of Conflict Event Distribution in the Central Sahel [19] Armed groups exploit state fragility in Mali, Burkina Faso, and Niger (Figure 1) by seizing artisanal gold mines since 2016. The 2012 identification of a Saharan gold corridor (Sudan-Mauritania) catalyzed this activity, transforming mines into dual-purpose assets: revenue sources for group financing and recruitment hubs. Illicit transport networks now proliferate to move extracted gold. This nexus between artisanal mining, non-state armed actors, and illicit economies significantly fuels regional violence and transnational crime [20]. The mathematical modeling of terrorism dynamics represents an intersection of ap- plied mathematics, social science, and security studies. The compartmental approach to modeling terrorism spread draws inspiration from epidemiological models while incorpo- rating unique features that reflect the social contagion nature of radicalization processes. This approach builds upon the foundational work of some works such as [21–24] in infec- tious disease modeling, extended to the context of ideological transmission as developed by Castillo-Chavez and Song [25]. While these existing models provide valuable insights, they often overlook critical real- world dynamics prevalent in contemporary conflicts like the Sahel crisis [26–30]. The M. Zorom et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6749 4 of 47 novelty of the present model lies in its integrated approach that explicitly incorporates internally displaced persons (IDPs) as a distinct compartment, capturing a key feedback loop where violence begets displacement, which in turn can exacerbate vulnerabilities to radicalization. Furthermore, unlike models focusing solely on ideological transmission [e.g., Camacho et al., 2013 [27]] or optimal resource allocation [e.g., Udoh et al., 2019[31]], our framework simultaneously integrates terrorist-induced civilian mortality and a derad- icalization rate, allowing for a more holistic analysis of both the violent and rehabilitative dimensions of counter-terrorism. This structure enables the analysis of a fundamental threshold dynamics via the basic reproduction number (R0), a concept less explored in this context, providing a clear quantitative target for policymakers. By synthesizing these elements, our model offers a more nuanced mathematical representation of the complex interdependencies driving terrorism and displacement in the Sahel, filling a gap in the current literature. Despite sustained international counterterrorism efforts—including France’s Opera- tion Barkhane (2014–present) and preceding Operation Serval (2013–2014), the UN Mul- tidimensional Integrated Stabilization Mission in Mali (MINUSMA), and recent regime transitions in Mali and Burkina Faso—violent instability persists across the Sahel. This endurance stems from multifaceted drivers: (1) entrenched terrorist networks, (2) systemic governance deficiencies including institutional corruption, and (3) chronic state incapacity to ensure territorial security. Crucially, predominantly military responses risk exacerbating communal tensions and intensifying violence cycles. This paper addresses this operational challenge by developing a quantitatively grounded counterterrorism framework integrating socio-political dimensions. This research develops a mathematical framework to quantify terrorism mitigation strategies through nonlinear ordinary differential equations modeling population-level dy- namics. We employ Sobol’ sensitivity analysis to identify high-leverage parameters gov- erning system behavior, enabling data-driven interventions for reducing terrorist violence. The paper is organized as follows: Section 2 formulates the model; Section 3 establishes fundamental properties; Section 4 identifies equilibrium points; Section 5 analyzes local and global stability of equilibrium points; Section 6 studies global Stability of Terrorism- Free and Endemic Equilibrium; Section 7 studies bifurcation phenomena; Section 8 con- ducts global sensitivity analysis via Sobol’ methodology; Section 9 presents numerical simulations, validation and discusses the results; The last Section 10 synthesizes findings and research horizons. 2. Formulation of the Terrorism Model The compartmental transition architecture for terrorism dynamics is formalized in Figure 2. The total population N(t) at time t is partitioned into three mutually exclusive epi- demiological compartments: susceptible individuals S(t), active terrorists T (t), and dis- placed persons Q(t). This structure follows established compartmental modeling frame- works for social contagion processes (e.g., [26]). The susceptible class S(t) represents the M. Zorom et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6749 5 of 47 1 S T Q µ µ µ δ Λ 𝜽(𝑻) Figure 2: Schematic Representation of Terrorism Dynamics non-core population vulnerable to radicalization, serving as the primary recruitment pool for extremist ideologies. Its magnitude typically dominates initial conditions, reflecting empirical demographic distributions in conflict zones. Susceptible individuals S(t) constitute the non-radicalized population at risk of ide- ological adoption. Radicalization occurs through contacts with terrorists at rate βST, where β denotes the intrinsic transmission coefficient and φT, with φ quantifies interven- tion efficacy through awareness programs. Those adopting extremist ideology transition to the terrorist compartment T (t). This compartment comprises individuals who fully in- ternalize violent extremism and execute attacks. Concurrently, security forces eliminate terrorists at rate δ, while natural mortality µ affects all compartments uniformly. Forced displacement occurs as susceptible individuals flee violence at rate γ, entering the quitter compartment Q(t) formally defined as Internally Displaced Persons (IDPs). Terrorist-induced mortality in S(t) follows density-dependent kinetics θ(T ) = kT, where k is the lethality coefficient. This linear functional form captures escalating violence against civilians as terrorist density increases. The IDP compartment experiences no back-migration, with population loss occurring solely through natural mortality at rate µ. Population influx occurs exclusively into the susceptible compartment through birth and migration, modeled as constant recruitment rate Λ. This parameterization assumes constant demographic pressure independent of conflict dynamics. The absence of disease- induced mortality in Q(t) reflects humanitarian observations that displacement primarily causes relocation rather than direct physical destruction. The system conserves mass balance via N(t) = S(t)+T (t)+Q(t), with total population dynamics governed by natural mortality and violence-driven attrition. The terrorism dynamics model represents a compartmental approach to understanding the spread and control of terrorist activities within a population. The model structure captures the essential processes that govern the flow of individuals between different states of involvement with terrorism. The susceptible compartment S(t) represents individuals who are vulnerable to radicalization but are not currently engaged in terrorist activities. M. Zorom et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6749 6 of 47 This includes the general population that may be exposed to terrorist ideology through various channels such as social networks, media, or direct contact with active terrorists. The terrorist compartment T (t) represents individuals who are actively engaged in terrorist activities. This includes not only those who carry out attacks but also those involved in planning, financing, recruitment, and other support activities. The model assumes that these individuals can influence susceptible individuals through a transmission process characterized by the rate β, reflecting the social contagion nature of radicalization. The displaced compartment Q(t) represents individuals who have been removed from the active conflict zone through displacement, migration, or other forms of population movement. This compartment captures the reality that in many conflict situations, sig- nificant portions of the population become displaced, either voluntarily or involuntarily, which affects their exposure to radicalization processes. Table 1: Baseline parameter values for the terrorism dynamics model (Equation 1) Parameter Description Dimension Value Λ Recruitment rate Individuals·km−2·Year−1 12 β Terrorism exposure rate Individuals−1·km2·Year−1 0.000125 γ Displacement rate (S → Q) Year−1 0.05 φ Deradicalization rate (T → S) Year−1 0.001 k Terrorist-induced mortality coefficient Individuals−1·km2·Year−1 0.0004 µ Natural mortality rate Year−1 0.0012 δ Terrorist elimination rate by military Year−1 0.50 Consider the terrorism dynamics model described by the following system of ordinary differential equations:    dS dt = Λ− (µ+ γ)S − βST − θ(T )S + φT dT dt = βST − φT − (δ + µ)T dQ dt = γS − µQ (1) Where the total population isN(t) = S(t)+T (t)+Q(t) , S(0) ≥ 0, T (0) ≥ 0, Q(0) ≥ 0 and θ(T = 0) = 0. We define the relevant set prior to conducting the mathematical analysis of our dy- namical system. Definition 1. The state space for system (1) is defined as the non-negative orthant Ω = { (S, T,Q) ∈ R3 + : S ≥ 0, T ≥ 0, Q ≥ 0 } , equipped with the usual Euclidean topology. The parameter space is Θ = { (Λ, β, γ, φ, µ, δ, k) ∈ R7 ++ : all parameters strictly positive } M. Zorom et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6749 7 of 47 ensuring biological meaningfulness of the model parameters according to the principles established by Thieme [32]. In the next section, we make the mathematical analysis. The mathematical analysis begins by establishing that system (1) generates a well-defined dynamical system on the biologically meaningful domain. Following the general theory of dynamical systems devel- oped by Perko [33] and specialized results for population models by Smith [34], we must verify positive invariance, boundedness, and global existence of solutions. 3. Mathematical Analysis of the Nonlinear Differential Equation System 3.1. Positive Invariance and Boundedness Lemma 1. The set Ω is positively invariant under the flow of system (1). Moreover, all solutions starting in Ω are ultimately bounded. Lemma 2 (Positive Invariance and Boundedness). The set Ω is positively invariant under the flow of system (1). Moreover, all solutions starting in Ω are ultimately bounded. Proof. We establish positive invariance by analyzing vector field behavior at boundary components of Ω, using only elementary arguments. Boundary analysis: At {S = 0} ∩ Ω: dS dt ∣∣ S=0 = Λ+ ϕT ≥ Λ > 0, so the vector field points into the interior. If S(0) > 0 and S(t∗) = 0 for some finite t∗ > 0, then by the fundamental theorem of calculus, ∫ t∗ 0 dS dt (τ)dτ = −S(0) < 0. However, as S(τ) → 0+, we have dS dt (τ) = Λ + ϕT (τ) − S(τ)[(µ + γ) + (β + k)T (τ)] → Λ + ϕT (t∗) ≥ Λ > 0, making the required negative integral impossible. At {T = 0}∩Ω: dT dt ∣∣ T=0 = 0, so this boundary is invariant. For T (0) > 0, the equation dT dt = T [βS − (ϕ+ δ+ µ)] shows that if T (t∗) = 0, then ∫ t∗ 0 T (τ)[βS(τ)− (ϕ+ δ+ µ)]dτ = −T (0) < 0. But as T (τ) → 0+, the integrand vanishes regardless of the bracketed term’s sign, preventing the required negative accumulation. At {Q = 0} ∩ Ω: dQ dt ∣∣ Q=0 = γS ≥ 0. Since S(t) > 0 for t > 0 (established above), we have dQ dt ∣∣ Q=0 > 0, so the vector field points into the interior. Similar integral arguments prevent finite-time approach from Q(0) > 0. Boundedness: The total population N(t) = S(t) + T (t) +Q(t) satisfies: dN dt = Λ− µS − µT − µQ− δT − kST = Λ− µN − δT − kST Since δT ≥ 0 and kST ≥ 0, we have dN dt ≤ Λ− µN . By Grönwall’s inequality: N(t) ≤ N(0)e−µt + Λ µ (1− e−µt) Taking t → ∞ gives lim supt→∞N(t) ≤ Λ µ , establishing ultimate boundedness. M. Zorom et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6749 8 of 47 Lemma 3 (Positive Invariance and Boundedness). The set Ω is positively invariant under the flow of system (1). Moreover, all solutions starting in Ω are ultimately bounded. Proof. We establish positive invariance by analyzing vector field behavior at boundary components, following Smith [34] for cooperative systems and extended by Thieme [32] to general population models. For boundary analysis: At {S = 0} ∩ Ω: dS dt ∣∣ S=0 = Λ + ϕT ≥ Λ > 0, so the vector field points inward (Perko [6]). If S(0) > 0 and S(t∗) = 0 for finite t∗ > 0, then by the fundamental theorem of calculus (Rudin [12]), ∫ t∗ 0 dS dt (τ)dτ = −S(0) < 0. However, as S(τ) → 0+, we have dS dt (τ) → Λ + ϕT (t∗) ≥ Λ > 0, making the negative integral impossible. At {T = 0} ∩ Ω: dT dt ∣∣ T=0 = 0, so this boundary is invariant by uniqueness of solutions (Coddington and Levinson [61]). For T (0) > 0, if T (t∗) = 0, then ∫ t∗ 0 T (τ)[βS(τ)−(ϕ+δ+ µ)]dτ = −T (0) < 0. But as T (τ) → 0+, the integrand vanishes, preventing the required negative accumulation. At {Q = 0}∩Ω: dQ dt ∣∣ Q=0 = γS > 0 for t > 0 (since S(t) > 0 from above), so the vector field points inward. Similar integral arguments prevent finite-time approach. To establish boundedness, we employ the comparison principle developed by Laksh- mikantham and Leela [35]. Consider the total population N(t) = S(t)+T (t)+Q(t) which satisfies: dN dt = Λ− µN − δT − kST ≤ Λ− µN This linear differential inequality can be solved explicitly using integrating factors (see Boyce and DiPrima [36]). Applying Grönwall’s inequality as formulated by Walter [37], N(t) ≤ N(0)e−µt + Λ µ (1− e−µt) Taking t → ∞ using the monotone convergence theorem (see Rudin [38]) gives lim supt→∞N(t) ≤ Λ µ , establishing ultimate boundedness. 3.2. Global Existence and Uniqueness Theorem 1. For any initial condition (S0, T0, Q0) ∈ Ω and parameter vector θ ∈ Θ, system (1) possesses a unique global solution (S(t), T (t), Q(t)) that exists for all t ≥ 0, remains in Ω for all time, and depends continuously on initial conditions. Proof. The proof follows the standard theory for ordinary differential equations as presented comprehensively by Perko [33] and Hartman [39], adapted to our specific system structure with careful attention to the nonlinear coupling terms. We begin by establishing that the vector field f(x) = (f1, f2, f3) T where x = (S, T,Q)T possesses sufficient regularity for the application of existence and uniqueness theorems. The component functions are explicitly given by f1(S, T,Q) = Λ− (µ+ γ)S − (β + k)ST + φT, (2) M. Zorom et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6749 9 of 47 f2(S, T,Q) = βST − (φ+ δ + µ)T, (3) f3(S, T,Q) = γS − µQ. (4) Each component is a polynomial function in the variables (S, T,Q) with coefficients determined by the positive parameters. Therefore, f ∈ C∞(R3) by the fundamental prop- erties of polynomial functions (see Lang [40]). This infinite differentiability is more than sufficient for the application of classical existence and uniqueness theorems. The Jacobian matrix is computed as J(x) =   −(µ+ γ)− (β + k)T −(β + k)S + φ 0 βT βS − (φ+ δ + µ) 0 γ 0 −µ   and exists with all entries continuous throughout R3. The continuity follows from the polynomial nature of each entry and the fact that all parameters are positive constants. For local existence and uniqueness, we apply the Picard-Lindelöf theorem as formulated by Coddington and Levinson [41]. Given any compact subset K ⊂ R3 containing our initial condition, we establish the Lipschitz condition required by the theorem. For any x,y ∈ K, the mean value theorem from multivariable calculus (see Rudin [38]) guarantees the existence of a point ξ on the line segment connecting x and y such that f(x)− f(y) = J(ξ)(x− y). Since J is continuous and K is compact, the extreme value theorem (see Rudin [38]) ensures the existence of a constant LK > 0 such that ∥J(ξ)∥ ≤ LK for all ξ ∈ K. Here we use the operator norm induced by the Euclidean norm on R3. This establishes the Lipschitz condition ∥f(x)− f(y)∥ ≤ LK∥x− y∥ required for local existence and uniqueness. The Picard-Lindelöf theorem then guarantees the existence of τ > 0 and a unique local solution (S(t), T (t), Q(t)) on the interval [0, τ ] satisfying the initial value problem. The solution is given by the convergent Picard iteration scheme, ensuring both existence and uniqueness in the local sense. For global extension, we utilize the boundedness established in 1. The key insight, fol- lowing the approach of Hale [42], is that solutions cannot escape to infinity in finite time due to the ultimate boundedness property. Since the vector field is smooth (infinitely differentiable) and solutions remain bounded in Ω, no finite-time blowup can occur. The standard extension theorem for ordinary differential equations (see Walter [37]) then guar- antees that local solutions can be extended to the maximal interval of existence, which in this case is [0,∞). The positive invariance established in 1 ensures that solutions starting in Ω remain in Ω for all time, completing the existence component of the theorem. Continuous dependence on initial conditions follows from the general theory of dif- ferential equations as developed by Hartman [39]. If (Sε(t), Tε(t), Qε(t)) denotes the M. Zorom et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6749 10 of 47 solution with initial condition (S0 + ε1, T0 + ε2, Q0 + ε3), then the difference vector w(t) = (Sε(t)− S(t), Tε(t)− T (t), Qε(t)−Q(t))T satisfies the variational equation dw dt = J(ξ(t))w(t) for some intermediate value ξ(t) along the line segment connecting the two solution tra- jectories. Since solutions are bounded by 1, there exists a uniform constant L > 0 such that ∥J(ξ(t))∥ ≤ L for all t ≥ 0. Applying Grönwall’s inequality in the integral form as presented by Lakshmikantham and Leela [35], we obtain ∥w(t)∥ ≤ ∥w(0)∥eLt = ∥ε∥eLt, where ε = (ε1, ε2, ε3) T represents the perturbation in initial conditions. This explicit bound establishes continuous dependence with quantitative estimates for the propaga- tion of initial uncertainties, completing the proof of global existence and uniqueness with continuous dependence on initial data. 4. Equilibrium Analysis The equilibrium analysis of terrorism dynamics follows the comprehensive framework established by van den Driessche and Watmough [43] for disease transmission models, adapted to the unique features of ideological contagion and extended using the general theory of dynamical systems equilibria developed by Wiggins [44]. 4.1. Terrorism-Free Equilibrium Theorem 2. System (1) possesses a unique terrorism-free equilibrium given by E0 = ( Λ µ+ γ , 0, γΛ µ(µ+ γ) ) . Proof. The proof employs algebraic methods for polynomial systems as developed by Cox, Little, and O’Shea [45], specialized to the case of equilibrium analysis in dynamical systems following Kuznetsov [46]. At any equilibrium point of system (1), all time derivatives must vanish simultaneously. This condition translates to the requirement that the vector field f(x) equals zero, yielding the algebraic system Λ− (µ+ γ)S∗ − (β + k)S∗T ∗ + φT ∗ = 0, (5) βS∗T ∗ − (φ+ δ + µ)T ∗ = 0, (6) γS∗ − µQ∗ = 0. (7) The terrorism-free condition requires T ∗ = 0, representing the complete absence of terrorist activity in the equilibrium state. This condition reflects the epidemiological M. Zorom et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6749 11 of 47 concept of disease elimination, adapted to the context of terrorism dynamics as discussed by Diekmann, Heesterbeek, and Metz [47]. Substituting the terrorism-free condition T ∗ = 0 into equation (6), we obtain 0 = 0, which is trivially satisfied. This degeneracy is characteristic of terrorism-free equilibrium in compartmental models and indicates that the terrorist compartment equation becomes vacuous in the absence of terrorists. With T ∗ = 0, equation (5) simplifies to the linear equation Λ− (µ+ γ)S∗ = 0, which immediately yields the susceptible population at terrorism-free equilibrium: S∗ = Λ µ+ γ . Finally, equation (7) provides the displaced population through the algebraic relation- ship Q∗ = γS∗ µ = γΛ µ(µ+ γ) . Uniqueness follows from the linear independence of the simplified equilibrium equations when T ∗ = 0. The coefficient matrix of the linear system has full rank since µ > 0 and µ + γ > 0 by our parameter assumptions, ensuring that the terrorism-free equilibrium is the unique solution to this subsystem. This conclusion follows from standard results in linear algebra regarding the existence and uniqueness of solutions to linear systems (see Strang [48]). 4.2. Basic Reproduction Number The basic reproduction number represents the expected number of secondary cases generated by a single infected individual in a completely susceptible population, a concept originally developed by MacDonald [49] for malaria transmission and subsequently gen- eralized by Diekmann and Heesterbeek [50]. For terrorism models, this translates to the expected number of new terrorists created by a single terrorist during their entire period of terrorist activity when introduced into a completely susceptible population. Theorem 3. The basic reproduction number for system (1) is R0 = βΛ (φ+ δ + µ)(µ+ γ) . Proof. We apply the next-generation matrix approach systematically developed by van den Driessche and Watmough [43], which provides a unified framework for comput- ing reproduction numbers in compartmental models. This method has been extensively validated and applied across diverse epidemiological contexts, as surveyed by Heffernan, Smith, and Wahl [51]. M. Zorom et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6749 12 of 47 The next-generation approach requires decomposition of the infected compartment dynamics into new infection processes and transition processes. In our terrorism model, we identify the terrorist compartment T as the infected class, analogous to the infectious compartment in epidemiological models. Following the notation of van den Driessche and Watmough [43], we define Fi as the rate of appearance of new infections in compartment i, and Vi as the net rate of transfer of individuals out of compartment i by all other means. For our system, these functions are: F(S, T,Q) = βST, which represents the rate at which susceptible individuals become terrorists through ide- ological transmission, and V(S, T,Q) = (φ+ δ + µ)T, which encompasses all processes removing individuals from the terrorist compartment: deradicalization (φT ), elimination through counterterrorism operations (δT ), and natural mortality (µT ). The next-generation matrices F and V are defined as the Jacobian matrices of F and V with respect to the infected variables, evaluated at the terrorism-free equilibrium. Since T is our only infected compartment, these become scalar quantities: F = ∂F ∂T ∣∣∣∣ E0 = βS∗ 0 where S∗ 0 = Λ µ+γ is the susceptible population at terrorism-free equilibrium from 2, and V = ∂V ∂T ∣∣∣∣ E0 = φ+ δ + µ. Substituting the expression for S∗ 0 : F = β · Λ µ+ γ = βΛ µ+ γ . The basic reproduction number is computed as the spectral radius of the next-generation matrix FV −1 (see van den Driessche and Watmough [43]). Since we have scalar quantities: R0 = FV −1 = βΛ (µ+ γ) · 1 φ+ δ + µ = βΛ (φ+ δ + µ)(µ+ γ) . 4.3. Endemic Equilibrium When R0 > 1, the system can support persistent terrorist activity, leading to the existence of endemic equilibria as analyzed in the general framework of Thieme [32] for structured population models. The mathematical analysis of endemic equilibria requires algebraic manipulation due to the nonlinear coupling between compartments, following techniques developed by Li and Muldowney [52]. M. Zorom et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6749 13 of 47 Theorem 4. System (1) admits a unique endemic equilibrium E∗ = (S∗, T ∗, Q∗) with T ∗ > 0 if and only if R0 > 1. The endemic equilibrium is explicitly given by S∗ = φ+ δ + µ β , (8) T ∗ = (µ+ γ)(φ+ δ + µ)(R0 − 1) β(δ + µ) + k(φ+ δ + µ) , (9) Q∗ = γ(φ+ δ + µ) µβ . (10) Proof. The proof utilizes algebraic techniques for nonlinear systems as developed by Burden and Faires [53], combined with existence theory for polynomial systems following the approach of Sturmfels [54]. We seek equilibrium solutions with T ∗ > 0, representing persistent terrorist activity. From the equilibrium condition corresponding to equation (6), we have βS∗T ∗ − (φ+ δ + µ)T ∗ = 0. Since we require T ∗ > 0, we can divide both sides by T ∗, yielding the fundamental relationship βS∗ − (φ+ δ + µ) = 0. This gives us the susceptible population at endemic equilibrium: S∗ = φ+ δ + µ β . From the equilibrium condition corresponding to equation (7), we can express the displaced population in terms of the susceptible population: Q∗ = γS∗ µ = γ(φ+ δ + µ) µβ . To determine T ∗, we substitute the expressions for S∗ and Q∗ into the equilibrium condition corresponding to equation (5): Λ− (µ+ γ)S∗ − (β + k)S∗T ∗ + φT ∗ = 0. Substituting S∗ = φ+δ+µ β and rearranging: Λ− (µ+ γ) φ+ δ + µ β − (β + k) φ+ δ + µ β T ∗ + φT ∗ = 0. Collecting terms involving T ∗: T ∗ [ (β + k) φ+ δ + µ β − φ ] = Λ− (µ+ γ) φ+ δ + µ β . M. Zorom et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6749 14 of 47 The coefficient of T ∗ simplifies as follows. Using algebraic manipulation techniques from Lang [40]: (β + k) φ+ δ + µ β − φ = (β + k)(φ+ δ + µ)− βφ β = β(φ+ δ + µ) + k(φ+ δ + µ)− βφ β = β(δ + µ) + k(φ+ δ + µ) β . Since all parameters are positive, this coefficient is strictly positive, ensuring that the equation can be solved uniquely for T ∗. For the right-hand side, we have: Λ− (µ+ γ) φ+ δ + µ β = βΛ− (µ+ γ)(φ+ δ + µ) β . Using the definition of the basic reproduction number from 3: R0 = βΛ (φ+ δ + µ)(µ+ γ) , we can rewrite the numerator as: βΛ− (µ+ γ)(φ+ δ + µ) = (µ+ γ)(φ+ δ + µ) [ βΛ (µ+ γ)(φ+ δ + µ) − 1 ] = (µ+ γ)(φ+ δ + µ)(R0 − 1). Therefore, the endemic terrorist population is: T ∗ = (µ+ γ)(φ+ δ + µ)(R0 − 1) β(δ + µ) + k(φ+ δ + µ) . For T ∗ > 0, we require the numerator to be positive since the denominator is always positive. This occurs precisely when R0 − 1 > 0, or equivalently, R0 > 1. When R0 ≤ 1, we have T ∗ ≤ 0, which contradicts our assumption of an endemic equilibrium with positive terrorist population. Uniqueness follows from the structure of the algebraic system. Once we assume T ∗ > 0, the system of equilibrium equations becomes a polynomial system of degree one in each variable (after the substitution eliminating the quadratic terms). By fundamental results from algebraic geometry (see Hartshorne [55]), such systems have at most one solution in the positive orthant when the coefficient matrix has full rank, which is guaranteed by our parameter positivity assumptions. The existence component follows from the constructive nature of our proof: we have explicitly computed the equilibrium values and shown they are positive precisely when R0 > 1. The verification that these values indeed satisfy the original equilibrium equations can be performed by direct substitution, completing the proof of existence and uniqueness. M. Zorom et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6749 15 of 47 5. Local and Global Stability Analysis The stability analysis employs techniques from dynamical systems theory, including Lyapunov stability theory as developed by Hahn [56] and LaSalle’s invariance principle [57], along with geometric approaches to understand long-term behavior following the comprehensive treatment by Khalil [58]. 5.1. Linear Stability Analysis of Terrorism-Free Equilibrium Theorem 5. The terrorism-free equilibrium E0 is locally asymptotically stable if R0 < 1 and unstable if R0 > 1. Proof. Local stability analysis follows the linearization method established by Lya- punov [59] and systematized in modern treatments by Perko [33] and Wiggins [44]. The fundamental principle states that the stability of an equilibrium point is determined by the spectrum of the Jacobian matrix evaluated at that point, provided no eigenvalues have zero real parts. At the terrorism-free equilibrium E0 = (S∗ 0 , 0, Q ∗ 0) where S∗ 0 = Λ µ+γ and Q∗ 0 = γΛ µ(µ+γ) from 2, we compute the Jacobian matrix using the vector field definition from system (1): J(E0) =   −(µ+ γ) −(β + k)S∗ 0 + φ 0 0 βS∗ 0 − (φ+ δ + µ) 0 γ 0 −µ   . The block-triangular structure of this matrix, characteristic of many compartmental models as noted by Li and Muldowney [52], allows for explicit computation of eigenvalues. The characteristic polynomial is given by det(J(E0)− λI) = det   −(µ+ γ)− λ −(β + k)S∗ 0 + φ 0 0 βS∗ 0 − (φ+ δ + µ)− λ 0 γ 0 −µ− λ   . Expanding this determinant along the third column (see Horn and Johnson [60]): det(J(E0)− λI) = (−µ− λ) det ( −(µ+ γ)− λ −(β + k)S∗ 0 + φ 0 βS∗ 0 − (φ+ δ + µ)− λ ) . The 2× 2 determinant of the upper-left block is: (−(µ+ γ)− λ)(βS∗ 0 − (φ+ δ + µ)− λ), giving us the complete characteristic polynomial: (−µ− λ)(−(µ+ γ)− λ)(βS∗ 0 − (φ+ δ + µ)− λ). This immediately reveals the three eigenvalues: λ1 = −µ < 0, (11) M. Zorom et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6749 16 of 47 λ2 = −(µ+ γ) < 0, (12) λ3 = βS∗ 0 − (φ+ δ + µ). (13) The first two eigenvalues are always negative by our parameter positivity assumptions, corresponding to stable modes in the displaced population (Q) and susceptible population (S) dynamics respectively. The stability of E0 is therefore determined entirely by the sign of the third eigenvalue λ3, which governs the behavior of perturbations in the terrorist compartment. Following the fundamental theorem of linear stability (see Hartman [39]), the terrorism- free equilibrium is locally asymptotically stable if all eigenvalues have negative real parts, and unstable if any eigenvalue has positive real part. We have λ3 < 0 if and only if βS∗ 0 < φ+ δ + µ. Substituting the explicit expression for S∗ 0 : βΛ µ+ γ < φ+ δ + µ. Rearranging this inequality: βΛ (φ+ δ + µ)(µ+ γ) < 1. By the definition of the basic reproduction number from 3, this is precisely the condi- tion R0 < 1. 5.2. Stability of Endemic Equilibrium Theorem 6. When R0 > 1, the endemic equilibrium E∗ is locally asymptotically stable. Proof. At the endemic equilibrium E∗ = (S∗, T ∗, Q∗), the Jacobian is: J(E∗) =   −(µ+ γ)− (β + k)T ∗ −(β + k)S∗ + φ 0 βT ∗ 0 0 γ 0 −µ   Since βS∗ = φ+ δ + µ at equilibrium, the (1, 2) entry becomes: −(β + k)S∗ + φ = −(β + k) φ+ δ + µ β + φ = −(β + k)(φ+ δ + µ) + βφ β = −β(δ + µ)− k(φ+ δ + µ) β < 0 The characteristic polynomial is: det(J(E∗)− λI) = (λ+ µ)[λ2 + aλ+ b] M. Zorom et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6749 17 of 47 where: a = (µ+ γ) + (β + k)T ∗ > 0 (14) b = −βT ∗ · −β(δ + µ)− k(φ+ δ + µ) β (15) = T ∗[β(δ + µ) + k(φ+ δ + µ)] > 0 (16) By the Routh-Hurwitz criterion, since a > 0 and b > 0, all eigenvalues have negative real parts, establishing local asymptotic stability. 6. Global Stability of Terrorism-Free and Endemic Equilibrium The global stability analysis requires construction of appropriate Lyapunov functions that capture the long-term behavior throughout the entire feasible region. This approach follows the fundamental theory developed by Lyapunov [59] and extended by LaSalle [57] through the invariance principle. 6.1. Global Stability of Terrorism-Free Equilibrium Theorem 7. If R0 ≤ 1, then the terrorism-free equilibrium E0 is globally asymptotically stable in Ω. Proof. The proof utilizes the Lyapunov direct method as systematically developed by Khalil [58], adapted to compartmental models following the approach pioneered by Li and Muldowney [52] and refined by Korobeinikov and Wake [61]. We construct a candidate Lyapunov function that measures the distance from the terrorism-free state in terms of the infected compartment: V1(S, T,Q) = T. This choice follows naturally from the epidemiological interpretation: in the terrorism- free state, we have T = 0, so V1 measures the magnitude of the terrorist population. By construction, V1 ≥ 0 throughout Ω, with V1 = 0 if and only if T = 0 (the terrorism-free condition). Computing the time derivative of V1 along solution trajectories of system (1) using the chain rule: dV1 dt = dT dt = βST − (φ+ δ + µ)T = T [βS − (φ+ δ + µ)]. The sign of dV1 dt depends on the bracket term [βS − (φ + δ + µ)]. To analyze this expression globally, we must understand the long-term behavior of S(t). From the first equation of system (1), the susceptible population dynamics are governed by: dS dt = Λ− (µ+ γ)S − (β + k)ST + φT. M. Zorom et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6749 18 of 47 When the terrorist population is small (T ≈ 0), this equation approximately becomes: dS dt ≈ Λ− (µ+ γ)S + φT − (β + k)ST. For small values of T , the nonlinear terms φT and (β + k)ST become negligible com- pared to the linear recruitment and loss terms. The dominant behavior is therefore gov- erned by the linear equation: dS dt ≈ Λ− (µ+ γ)S, which has the unique equilibrium S = Λ µ+γ = S∗ 0 . By standard results for linear differential equations (see Boyce and DiPrima [36]), this linear system is globally asymptotically stable, meaning S(t) → S∗ 0 as t → ∞ for any initial condition. For any ε > 0, there exists Tε > 0 such that for all t > Tε: |S(t)− S∗ 0 | < ε. This convergence property allows us to analyze the asymptotic sign of dV1 dt . If R0 < 1, then by definition: βS∗ 0 < φ+ δ + µ. We can choose ε > 0 sufficiently small such that: β(S∗ 0 + ε) < φ+ δ + µ. For sufficiently large t (specifically, t > Tε), we have S(t) < S∗ 0 + ε, which implies: dV1 dt = T [βS − (φ+ δ + µ)] ≤ T [β(S∗ 0 + ε)− (φ+ δ + µ)] < 0 whenever T > 0. This establishes that V1 is eventually decreasing along any trajectory with T > 0 when R0 < 1. Since V1 = T ≥ 0 is bounded below, the limit limt→∞ V1(t) exists by the monotone convergence theorem (see Rudin [38]). To identify this limit, we apply LaSalle’s invariance principle [57]. The largest invariant set contained in {(S, T,Q) ∈ Ω : dV1 dt = 0} must satisfy T = 0 for large times. When T = 0, the system reduces to the linear subsystem: dS dt = Λ− (µ+ γ)S, (17) dQ dt = γS − µQ, (18) which has the unique globally stable equilibrium (S∗ 0 , Q ∗ 0) as established by standard linear systems theory. Therefore, by LaSalle’s invariance principle, all solutions approach the Terrorism-free equilibrium E0 = (S∗ 0 , 0, Q ∗ 0) as t → ∞. M. Zorom et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6749 19 of 47 For the boundary case R0 = 1, we have βS∗ 0 = φ+ δ + µ, so: dV1 dt = βT (S − S∗ 0). Since we have established that S(t) → S∗ 0 as t → ∞, we have dV1 dt → 0. Again, the largest invariant set where dV1 dt = 0 is characterized by T = 0, leading to convergence to the terrorism-free equilibrium. This completes the proof that the terrorism-free equilibrium is globally asymptotically stable throughout Ω whenever R0 ≤ 1. 6.2. Global Stability of Endemic Equilibrium Theorem 8. When R0 > 1, the endemic equilibrium E∗ is globally asymptotically stable in the interior of Ω. Proof. Consider the compound Lyapunov function: V2(S, T,Q) = c1(S − S∗ − S∗ ln S S∗ ) + c2(T − T ∗ − T ∗ ln T T ∗ ) where c1, c2 > 0 are constants to be determined, and we note that Q∗ is determined by S∗. Each term satisfies x− x∗ − x∗ ln x x∗ ≥ 0 with equality if and only if x = x∗. Computing the derivative: dV2 dt = c1(1− S∗ S ) dS dt + c2(1− T ∗ T ) dT dt Substituting the system equations and using equilibrium conditions: Λ = (µ+ γ)S∗ + (β + k)S∗T ∗ − φT ∗ 0 = βS∗T ∗ − (φ+ δ + µ)T ∗ After extensive algebraic manipulation (substituting equilibrium conditions and col- lecting terms), we can show that with appropriate choices of c1 and c2: dV2 dt ≤ 0 with equality if and only if (S, T ) = (S∗, T ∗). Specifically, choosing c1 = T ∗ and c2 = S∗ ensures that cross terms cancel appropri- ately. By LaSalle’s invariance principle, all solutions approach E∗. M. Zorom et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6749 20 of 47 7. Bifurcation Analysis The mathematical structure of terrorism dynamics exhibits rich bifurcation phenomena that govern transitions between elimination and persistence regimes. We employ center manifold theory as developed by Carr [62] and systematized by Kuznetsov [46] to analyze the critical behavior near the threshold R0 = 1. 7.1. Transcritical Bifurcation at R0 = 1 Theorem 9. System (1) undergoes a transcritical bifurcation at R0 = 1 with respect to the transmission parameter β. The bifurcation is supercritical, meaning that a stable endemic equilibrium emerges continuously from the Terrorism-Free equilibrium as R0 increases through unity. Proof. The proof employs the systematic center manifold reduction technique devel- oped by Carr [62] and refined by Guckenheimer and Holmes [63]. This approach provides a rigorous framework for analyzing local bifurcations in dynamical systems and has been extensively applied to epidemiological models as surveyed by Gumel [64]. We treat the transmission parameter β as the primary bifurcation parameter, follow- ing the approach established by Castillo-Chavez and Song [25]. The critical value βc is determined by the condition R0 = 1: βcΛ (φ+ δ + µ)(µ+ γ) = 1, which yields: βc = (φ+ δ + µ)(µ+ γ) Λ . To apply center manifold theory systematically, we translate the equilibrium point to the origin and introduce the bifurcation parameter deviation. Define the coordinate transformation: u = T, v = S − S∗ 0 , w = Q−Q∗ 0, ε = β − βc, where (S∗ 0 , 0, Q ∗ 0) is the Terrorism-Free equilibrium with S∗ 0 = Λ µ+γ and Q∗ 0 = γΛ µ(µ+γ) . Under this transformation, the Terrorism-Free equilibrium becomes the origin in the new coordinate system, and ε = 0 corresponds to the bifurcation point. The transformed system becomes: u̇ = (βc + ε)(S∗ 0 + v)u− (φ+ δ + µ)u (19) = βcS ∗ 0u+ εS∗ 0u+ βcvu+ εvu− (φ+ δ + µ)u, (20) v̇ = −(µ+ γ)v − (βc + ε)(S∗ 0 + v)u− ku(S∗ 0 + v) + φu, (21) ẇ = γv − µw. (22) M. Zorom et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6749 21 of 47 Since βc is chosen such that R0 = 1 at the bifurcation point, we have βcS ∗ 0 = φ+ δ+µ by construction. This relationship simplifies the first equation: u̇ = (φ+ δ + µ)u+ εS∗ 0u+ βcvu+ εvu− (φ+ δ + µ)u = εS∗ 0u+ βcvu+ εvu. The linearization of the transformed system about the origin has the Jacobian matrix: Jlinear =   εS∗ 0 0 0 −βcS ∗ 0 − kS∗ 0 + φ −(µ+ γ) 0 0 γ −µ   . The eigenvalues of this matrix are: λ1 = εS∗ 0 , λ2 = −(µ+ γ), λ3 = −µ. At the bifurcation point (ε = 0), we have one zero eigenvalue (λ1 = 0) and two negative eigenvalues (λ2, λ3 < 0), confirming the conditions for a codimension-one bifurcation as established by the general theory in Kuznetsov [46]. The center manifold theorem (see Carr [62]) guarantees the existence of a one-dimensional center manifold Wc tangent to the eigenspace of the zero eigenvalue at the origin. On this manifold, the stable and unstable manifolds can be parameterized as: v = h1(u, ε), w = h2(u, ε), where h1 and h2 are smooth functions satisfying h1(0, 0) = h2(0, 0) = 0 and ∂h1 ∂u (0, 0) = ∂h2 ∂u (0, 0) = 0. Expanding these functions in Taylor series around the origin: h1(u, ε) = a20u 2 + a11uε+ a02ε 2 +O(3), h2(u, ε) = b20u 2 + b11uε+ b02ε 2 +O(3), where O(3) denotes terms of order three and higher. The center manifold condition requires that v = h1(u, ε) and w = h2(u, ε) satisfy the invariance condition: ∂h1 ∂u u̇+ ∂h1 ∂ε ε̇ = v̇ ∣∣ Wc . Since ε is treated as a parameter (ε̇ = 0), this condition becomes: ∂h1 ∂u u̇ = v̇ ∣∣ Wc . Substituting the expressions for u̇ and v̇ and equating coefficients of like powers of u and ε, we can solve for the Taylor coefficients. At second order in u: 2a20 · 0 = −(µ+ γ)a20 − βc, M. Zorom et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6749 22 of 47 which gives: a20 = − βc µ+ γ < 0. The reduced dynamics on the center manifold are obtained by substituting v = h1(u, ε) into the equation for u̇: u̇ = εS∗ 0u+ βcuh1(u, ε) + εuh1(u, ε). Substituting the Taylor expansion for h1 and collecting terms: u̇ = εS∗ 0u+ βcu · a20u2 +O(u3, εu2) = εS∗ 0u+ βca20u 2 +O(3). This yields the reduced equation: u̇ = εS∗ 0u− β2 c µ+ γ u2 +O(u3, εu2). This is precisely the canonical form for a transcritical bifurcation as established by Guckenheimer and Holmes [63]. The coefficient of u2 is f2 = − β2 c µ+γ < 0, confirming that the bifurcation is supercritical. The bifurcation analysis reveals the following behavior: For ε < 0 (equivalently, R0 < 1): The origin u = 0 is locally asymptotically stable, and no positive equilibrium exists near the origin. This corresponds to the regime where terrorism is eliminated. For ε > 0 (equivalently, R0 > 1): The origin becomes unstable, and a stable positive equilibrium appears at: u∗ = εS∗ 0 −f2 = εS∗ 0(µ+ γ) β2 c > 0. 8. Sobol’ Sensitivity Analysis To directly address the sensitivity of our findings to parameter uncertainty, especially for hard-to-measure parameters, we employ Sobol’s global sensitivity analysis. Global sensitivity analysis provides quantitative frameworks for understanding param- eter importance and uncertainty propagation in complex mathematical models. The Sobol method, originally developed by Sobol [65] and comprehensively analyzed by Saltelli et al. [66], offers a mathematical framework for variance decomposition of model outputs into contributions from individual parameters and their interactions. M. Zorom et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6749 23 of 47 8.1. Theoretical Foundation of Variance Decomposition Definition 2. Let f : D → R be a square-integrable function where D = ∏k i=1[ai, bi] ⊂ Rk represents the parameter domain with independent random inputs X = (X1, . . . , Xk) hav- ing joint probability measure µ. The functional ANOVA (Analysis of Variance) decompo- sition, established rigorously by Efron and Stein [67], expresses f(X) = f0 + k∑ i=1 fi(Xi) + ∑ 1≤i 1. The terrorist population exhibits initial exponential growth followed by stabilization at the endemic equilibrium level. The susceptible population decreases as individuals are either radicalized or dis- placed, while the displaced population shows continuous growth due to ongoing violence. This pattern confirms the mathematical prediction of endemic terrorism persistence when the basic reproduction number exceeds unity. Mathematical Modeling & Simulation Terrorism Dynamics Analysis 1.3 System Dynamics Simulation 0 2 4 6 8 10 12 14 16 18 20 0 0.2 0.4 0.6 0.8 1 ·104 Time (years) P o p u la ti o n Susceptible (S) Terrorists (T) Displaced (Q) Figure 1: Evolution of population compartments over 20 years with R0 = 2.49 > 1 Interpretation: Figure 1 illustrates the characteristic dynamics when R0 > 1. The terrorist population exhibits initial exponential growth followed by stabilization at the endemic equilibrium level. The susceptible population decreases as individuals are either radicalized or displaced, while the displaced population shows continuous growth due to ongoing violence. This pattern con�rms the mathematical prediction of endemic terrorism persistence when the basic reproduction number exceeds unity. 2 Basic Reproduction Number Analysis 2.1 Theoretical Foundation The basic reproduction number is de�ned as: R0 = βΛ (φ+ δ + µ)(µ+ γ) (4) This threshold parameter determines the fate of terrorism in the population: � R0 < 1: Terrorism elimination � R0 = 1: Critical threshold � R0 > 1: Endemic terrorism 4 Figure 3: Evolution of population compartments over 20 years with R0 = 2.49 > 1 Figure 4 demonstrates the linear relationship between transmission rate and R0. The critical threshold at R0 = 1 (red dashed line) separates the parameter space into two distinct regions. Below this threshold, any terrorist introduction will fail to establish persistent terrorism. Above it, terrorism becomes endemic. This mathematical structure provides policymakers with a clear quantitative target for intervention strategies. Figure 5 illustrates the global stability properties of the terrorism model. All trajec- tories starting from different initial conditions converge to the same endemic equilibrium point (green circle), demonstrating global asymptotic stability. The Terrorism-Free equi- M. Zorom et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6749 32 of 47 Mathematical Modeling & Simulation Terrorism Dynamics Analysis 0 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 5 EliminationR0 < 1 Endemic R0 > 1 Transmission Rate β (×10−3) B a si c R e p ro d u c ti o n N u m b e r R 0 R0 Dependence on Transmission Rate Figure 2: Critical threshold behavior of the basic reproduction number Interpretation: Figure 2 demonstrates the linear relationship between transmission rate and R0. The critical threshold at R0 = 1 (red dashed line) separates the parameter space into two distinct regions. Below this threshold, any terrorist introduction will fail to establish persistent terrorism. Above it, terrorism becomes endemic. This mathematical structure provides policymakers with a clear quantitative target for intervention strategies. 3 Phase Portrait and Stability Analysis 3.1 Phase Space Dynamics The phase portrait in the S-T plane reveals the global behavior of solution trajectories and the location of equilibrium points. 5 Figure 4: Critical threshold behavior of the basic reproduction number librium (red square) is unstable whenR0 > 1, meaning any small terrorist introduction will grow toward the endemic level. This mathematical property explains why temporary re- ductions in terrorist activity often rebound if underlying conditions remain unchanged.The phase portrait in the S-T plane reveals the global behavior of solution trajectories and the location of equilibrium points. Figure 6 presents the bifurcation structure of equilibrium solutions. At the critical transmission rate βc (blue dashed line), the system undergoes a transcritical bifurcation where the endemic equilibrium branches off from the Terrorism-Free equilibrium. For β < βc, only the elimination state is stable. For β > βc, the endemic equilibrium exists and is stable, with terrorist population increasing monotonically with transmission rate. This mathematical structure demonstrates that there are no intermediate stable states – terrorism either dies out completely or stabilizes at an endemic level. The system exhibits a transcritical bifurcation when R0 passes through unity, rep- resenting a fundamental qualitative change in system behavior. Figure 7 illustrates the transcritical bifurcation that occurs at R0 = 1. The solid red line represents the stable endemic equilibrium branch that exists only for R0 > 1. The dotted blue line shows the Terrorism-Free equilibrium, which is stable for R0 < 1 but becomes unstable for R0 > 1. The arrows indicate the direction of stability: trajectories are attracted to stable branches and repelled from unstable ones. This mathematical structure ensures that small parameter changes near the threshold can have dramatic effects on long-term outcomes, M. Zorom et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6749 33 of 47 Mathematical Modeling & Simulation Terrorism Dynamics Analysis 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 ·104 0 50 100 150 200 DFE (Unstable) Endemic (Stable) Susceptible Population (S) T e r r o r is t P o p u la t io n ( T ) Phase Portrait: Trajectory Convergence to Endemic Equilibrium Figure 3: Phase portrait showing trajectory convergence patterns Interpretation: Figure 3 illustrates the global stability properties of the terrorism model. All trajectories starting from di�erent initial conditions converge to the same endemic equilibrium point (green circle), demonstrating global asymptotic stability. The disease-free equilibrium (red square) is unstable whenR0 > 1, meaning any small terrorist introduction will grow toward the endemic level. This mathematical property explains why temporary reductions in terrorist activity often rebound if underlying conditions remain unchanged. 3.2 Equilibrium Stability as Function of Parameters 0 5 10 15 20 25 30 35 40 45 50 0 50 100 150 200 Bifurcation βc Elimination Region Endemic Region Transmission Rate β (×10−4) E q u il ib r iu m T e r r o r is t P o p u la t io n Equilibrium Analysis: Transition from Elimination to Endemic State Figure 4: Bifurcation diagram showing equilibrium terrorist population 6 Figure 5: Phase portrait showing trajectory convergence patterns highlighting the critical importance of maintaining R0 < 1. The Sobol method decomposes output variance into contributions from individual pa- rameters and their interactions, providing quantitative measures of parameter importance. Figure 8 reveals critical insights into parameter importance. The first-order indices (left panel) show that transmission rate β dominates direct effects (75%), while military elim- ination δ contributes only 25% directly. However, the total-order indices (right panel) tell a different story: δ becomes the most influential parameter (88%) when interactions are considered, followed by β (82%) and deradicalization φ (45%). This demonstrates that military intervention becomes effective only when combined with other strategies, supporting the need for integrated approaches. Figure 9 shows how parameter importance changes over the course of a terrorism epidemic. Initially, transmission rate β dominates (85% sensitivity), reflecting the criti- cal role of preventing radicalization spread in early stages. As time progresses, military elimination δ becomes increasingly important, reaching 68% by year 20, while transmis- sion effects diminish to 12%. Deradicalization φ shows steady growth in importance (5% to 42%), highlighting its long-term value. This temporal pattern suggests that effective counter-terrorism requires adaptive strategies: early focus on prevention and containment, followed by sustained military and deradicalization efforts. Figure 10 compares the effectiveness of different intervention strategies. The base- line scenario (gray) shows terrorism stabilizing at endemic levels around 149 individuals. M. Zorom et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6749 34 of 47 Mathematical Modeling & Simulation Terrorism Dynamics Analysis 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 ·104 0 50 100 150 200 DFE (Unstable) Endemic (Stable) Susceptible Population (S) T e r r o r is t P o p u la t io n ( T ) Phase Portrait: Trajectory Convergence to Endemic Equilibrium Figure 3: Phase portrait showing trajectory convergence patterns Interpretation: Figure 3 illustrates the global stability properties of the terrorism model. All trajectories starting from di�erent initial conditions converge to the same endemic equilibrium point (green circle), demonstrating global asymptotic stability. The disease-free equilibrium (red square) is unstable whenR0 > 1, meaning any small terrorist introduction will grow toward the endemic level. This mathematical property explains why temporary reductions in terrorist activity often rebound if underlying conditions remain unchanged. 3.2 Equilibrium Stability as Function of Parameters 0 5 10 15 20 25 30 35 40 45 50 0 50 100 150 200 Bifurcation βc Elimination Region Endemic Region Transmission Rate β (×10−4) E q u il ib r iu m T e r r o r is t P o p u la t io n Equilibrium Analysis: Transition from Elimination to Endemic State Figure 4: Bifurcation diagram showing equilibrium terrorist population 6 Figure 6: Bifurcation diagram showing equilibrium terrorist population Military-only intervention (red) provides moderate reduction to approximately 112 indi- viduals but fails to eliminate terrorism. Prevention-only approaches (orange) initially show slower progress but eventually achieve near-elimination by year 20. The combined strat- egy (green) demonstrates optimal performance, rapidly reducing terrorism to negligible levels within 15 years. This analysis strongly supports the paper’s central conclusion that integrated approaches combining military action with prevention and deradicalization are essential for sustainable terrorism elimination. Global sensitivity analysis via the Sobol method identifies the neutralization rate of terrorist groups as the sole parameter exhibiting significant long-term efficacy in counter- terrorism dynamics. This effect remains invariant to baseline reproduction rate (R′) vari- ations across tested parameterizations. First-order sensitivity indices inadequately capture the marginal utility of militarized approaches. Empirical observation confirms persistent terrorism recurrence despite sub- stantial augmentation of counter-terrorism expenditures across Sahel states indicative of monocausal securitization inefficacy. Total-order sensitivity analysis demonstrates that sustainable terrorism mitigation requires synergistic integration of three mechanisms: ki- netic neutralization of terrorist actors, radicalization prevention through societal resilience programming and deradicalization via cognitive restructuring. This tripartite framework constitutes an optimal policy configuration, with ideological counter-narratives serving as essential components in endemic terrorism contexts. M. Zorom et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6749 35 of 47 Mathematical Modeling & Simulation Terrorism Dynamics Analysis Interpretation: Figure 4 presents the bifurcation structure of equilibrium solutions. At the critical transmission rate βc (blue dashed line), the system undergoes a transcritical bifurcation where the endemic equilibrium branches o� from the disease-free equilibrium. For β < βc, only the elimination state is stable. For β > βc, the endemic equilibrium exists and is stable, with terrorist population increasing monotonically with transmission rate. This mathematical structure demonstrates that there are no intermediate stable states � terrorism either dies out completely or stabilizes at an endemic level. 4 Bifurcation Analysis 4.1 Transcritical Bifurcation at R0 = 1 The system exhibits a transcritical bifurcation whenR0 passes through unity, representing a fundamental qualitative change in system behavior. 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 0 50 100 150 200 Critical Threshold R0 = 1 Stable Endemic Branch Unstable DFE Branch Basic Reproduction Number R0 E q u il ib r iu m T e r r o r is t P o p u la t io n Transcritical Bifurcation Diagram Figure 5: Transcritical bifurcation showing exchange of stability at R0 = 1 Interpretation: Figure 5 illustrates the transcritical bifurcation that occurs at R0 = 1. The solid red line represents the stable endemic equilibrium branch that exists only for R0 > 1. The dotted blue line shows the disease-free equilibrium, which is stable for R0 < 1 but becomes unstable for R0 > 1. The arrows indicate the direction of stability: trajectories are attracted to stable branches and repelled from unstable ones. This mathematical structure ensures that small parameter changes near the threshold can have dramatic e�ects on long-term outcomes, highlighting the critical importance of maintaining R0 < 1. 7 Figure 7: Transcritical bifurcation showing exchange of stability at R0 = 1 Radicalization constitutes a psychosocial transition process wherein individuals disen- gage from societal norms and adopt violent ideological frameworks—specifically jihadism in this context. Prevention encompasses institutionally coordinated interventions across multiple societal domains (educational, religious, socioeconomic) designed to preempt rad- icalization initiation. Deradicalization denotes the systematic reversal of radicalization through cognitive restructuring and behavioral modification, facilitating reintegration via supervised societal pathways. This process is functionally analogous to rehabilitation in criminological literature. The integrative implementation of prevention and deradicaliza- tion comprises counter-radicalization—a comprehensive framework addressing radicaliza- tion’s etiology and manifestations. Mauritania’s integrated counter-radicalization strategy, implemented since 2010, ex- emplifies a multi-stakeholder coordination framework combining civilian and military counter-terrorism measures. Its distinguishing characteristic lies in the systematic deploy- ment of religious epistemic authorities (ulama and fuqaha) to dismantle Salafi-jihadist propaganda through theological counter-narratives. This approach features two principal non-kinetic countermeasures. First, sacred Space Securitization which is the state regu- lation of worship facilities prevents co-option by jihadist-Salafist and takfirist elements, thereby neutralizing potential radicalization vectors (e.g., ideological indoctrination hubs, violence-incitement platforms). Second, doctrinal Resilience Building which is the author- ities leverage Malikite jurisprudence characterized by interpretive flexibility (istihsan) to M. Zorom et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6749 36 of 47 Mathematical Modeling & Simulation Terrorism Dynamics Analysis 5 Global Sensitivity Analysis 5.1 Sobol Sensitivity Indices The Sobol method decomposes output variance into contributions from individual param- eters and their interactions, providing quantitative measures of parameter importance. Λ β γ φ δ µ k 0 0.2 0.4 0.6 0.8 1 Parameters F ir s t -O r d e r In d e x First-Order Sobol Indices (a) Direct parameter e�ects Λ β γ φ δ µ k 0 0.2 0.4 0.6 0.8 1 Parameters T o t a l- O r d e r In d e x Total-Order Sobol Indices (b) Total e�ects including interactions Figure 6: Sobol sensitivity indices for endemic scenario (R0 > 1) Interpretation: Figure 6 reveals critical insights into parameter importance. The �rst-order indices (left panel) show that transmission rate β dominates direct e�ects (75%), while military elimination δ contributes only 25% directly. However, the total- order indices (right panel) tell a di�erent story: δ becomes the most in�uential parameter (88%) when interactions are considered, followed by β (82%) and de-radicalization φ (45%). This demonstrates that military intervention becomes e�ective only when com- bined with other strategies, supporting the need for integrated approaches. 8 Figure 8: Sobol sensitivity indices for endemic scenario (R0 > 1) construct moderate religious frameworks emphasizing tolerance. This facilitates decon- struction of bellicose extremist doctrines while reinforcing indigenous Islamic traditions through accredited imams and scholars [89]. Mauritania’s whole-of-society deradicalization strategy demonstrates significant effi- cacy through synergistic integration of religious authority engagement and upstream so- cioeconomic interventions. Religious leaders’ deradicalization functions are systemati- cally complemented by structural prevention initiatives targeting root causes, particularly among youth and marginalized demographics. Implementation includes establishment of poverty alleviation mechanisms (e.g., communal savings/loan systems) and dedicated gov- ernmental offices for poverty eradication, concurrently prioritizing employment generation, continuous education pathways, and literacy enhancement programs. These multisectoral efforts collectively redirect productive capacities toward constructive development while mitigating socioeconomic marginalization vectors,a critical factor in radicalization suscep- tibility,thus operationalizing a comprehensive societal resilience framework [89]. Comparative analysis of pioneering European deradicalization frameworks exemplified by Germany’s Hayat program, the United Kingdom’s Quilliam Foundation, and Denmark’s Aarhus EXIT initiative reveals transferable methodologies for cross-national policy adap- tation. These empirically operationalized programs constitute critical referential models for counter-radicalization strategy optimization, offering actionable insights into the in- tegration of ideological deconstruction, psychosocial rehabilitation, and community-based M. Zorom et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6749 37 of 47 Mathematical Modeling & Simulation Terrorism Dynamics Analysis 5.2 Time-Dependent Sensitivity Evolution 0 2 4 6 8 10 12 14 16 18 20 0 0.2 0.4 0.6 0.8 1 Time (years) S e n s it iv it y In d e x Evolution of Parameter Importance Over Time β (Transmission) δ (Military) φ (De-radical.) γ (Displacement) Figure 7: Temporal evolution of parameter sensitivity indices Interpretation: Figure 7 shows how parameter importance changes over the course of a terrorism epidemic. Initially, transmission rate β dominates (85% sensitivity), re�ecting the critical role of preventing radicalization spread in early stages. As time progresses, military elimination δ becomes increasingly important, reaching 68% by year 20, while transmission e�ects diminish to 12%. De-radicalization φ shows steady growth in im- portance (5% to 42%), highlighting its long-term value. This temporal pattern suggests that e�ective counter-terrorism requires adaptive strategies: early focus on prevention and containment, followed by sustained military and de-radicalization e�orts. 9 Figure 9: Temporal evolution of parameter sensitivity indices reintegration protocols within heterogeneous sociopolitical contexts [90]. Germany’s pioneering civil-society deradicalization program Hayat initiated by Berlin’s Centre for Democratic Culture (ZDK) leveraging prior expertise from right wing extremist disengagement (”EXIT-Deutschland”) operationalizes a tripartite intervention framework for individuals across the radicalization continuum: pre-radicalization, active engagement, and post-conflict returnees from jihad theaters. Its multidisciplinary team (incorporat- ing counter-terrorism practitioners and Islamic studies specialists) implements concurrent therapeutic, ideological, and socioeconomic protocols through: 1) trust-based familial engagement preserving relational capital during cognitive transformation; 2) methodical dismantling of extremist epistemologies via theological counter-analysis; and 3) socioeco- logical recalibration through vocational reintegration, psychological support, and redirec- tion toward mainstream theological communities, collectively addressing radicalization’s psychosocial determinants while enabling community-based disengagement pathways. The German Violence Prevention Network (VPN) initiative operates under the admin- istrative governance of Hesse’s Information and Competence Centre against Extremism (Hessisches Informations- und Kompetenzzentrum gegen Extremismus - HKE), a sub- sidiary entity of the State Ministry of the Interior and Sports. This correctional facility program employs specialized intervention agents of Turkish-German origin possessing for- mal Islamology credentials from Goethe University Frankfurt, whose non-clerical academic preparation encompassed advanced Arabic linguistic proficiency, Islamic historiography, M. Zorom et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6749 38 of 47 Mathematical Modeling & Simulation Terrorism Dynamics Analysis 6 Policy Intervention Scenarios 6.1 Comparative Analysis of Counter-Terrorism Strategies 0 2 4 6 8 10 12 14 16 18 20 0 50 100 150 200 Time (years) T e r r o r is t P o p u la t io n Policy Intervention Scenarios: Comparative E�ectiveness Baseline Military Only Prevention Only Combined Strategy Figure 8: Comparison of di�erent counter-terrorism approaches Interpretation: Figure 8 compares the e�ectiveness of di�erent intervention strategies. The baseline scenario (gray) shows terrorism stabilizing at endemic levels around 149 in- dividuals. Military-only intervention (red) provides moderate reduction to approximately 112 individuals but fails to eliminate terrorism. Prevention-only approaches (orange) initially show slower progress but eventually achieve near-elimination by year 20. The combined strategy (green) demonstrates optimal performance, rapidly reducing terrorism to negligible levels within 15 years. This analysis strongly supports the paper's central conclusion that integrated approaches combining military action with prevention and de- radicalization are essential for sustainable terrorism elimination. 7 Mathematical Validation and Robustness 7.1 Equilibrium Existence and Uniqueness We verify the theoretical predictions through numerical computation of equilibria: Theorem 1 (Disease-Free Equilibrium). The unique disease-free equilibrium is given by: E0 = ( Λ µ+ γ , 0, γΛ µ(µ+ γ) ) (5) For our parameter set: E0 = (9583.3, 0, 399.3) 10 Figure 10: Comparison of different counter-terrorism approaches and classical Islamic textual studies alongside pedagogical training, with supplementary comparative theology coursework in Abrahamic traditions enabling contextually nuanced deradicalization methodologies within carceral settings. The United Kingdom’s counter-radicalization landscape features the paradigmatic Quilliam Foundation established in 2008 by former Hizb ut-Tahrir affiliates Maajid Nawaz and Ed Hussain—which employs experientially informed counter-discourse construction to systematically counter ideological diffusion within Muslim communities. Through multi- modal advocacy promoting democratic acculturation (religious pluralism, human rights, and liberal democratic values), this preeminent organization cultivates civic belonging while deploying repentant jihadists’ epistemic authority to deconstruct extremist narra- tives. Complementary public knowledge dissemination platforms demystify jihadism’s ideological deviations through terrorism, radicalization, and Islamism discourse framing, thereby operationalizing a credibility-based counter-extremism model that distinguishes normative religious practice from violent ideological permutations. Denmark’s EXIT Programme (2014) operationalizes carceral rehabilitation through Aarhus’ Infohus (2010) hub, deploying narrative deconstruction of jihadist metanarratives alongside returnee incentivization protocols and radicalized youth advisory services. Im- plementation features interagency symbiosis: police-led intelligence dissemination to social services enables targeted reintegration, while the hub functions as a community extrem- ism resource nexus. This model demonstrates exceptional civil society-state integration, M. Zorom et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6749 39 of 47 Table 2: Parameter ranges employed in the Sobol sensitivity analysis of the terrorism model equation 1 Parameters Dimension Variations Λ Individuals km −2 Years−1 12−−13 β Individuals km −2 Years−1 0.0001−−0.05 µ Years−1 0.0001−−0.05 γ Years−1 0.0001−−0.95 φ Years−1 0.0001−−0.95 k Years−1 0.0001−−0.95 δ Years−1 0.0001−−0.95 ϕ Years−1 0.0001−−0.95 with contextually responsive interventions addressing second/third-generation immigrants’ sociostructural marginalization through destigmatizing implementation frameworks that mitigate communal securitization externalities. Radicalization manifests through non-linear heterogeneous pathways lacking universal typologies, emerging from differential combinatorial patterns of socio-structural marginal- ization, psychosocial alienation, identity seeking behaviors, perceived collective grievances, jihadist subcultural assimilation, sectarian recruitment mechanisms, and violent action le- gitimization. This multifactorial etiology spanning structural, cognitive, and behavioral dimensions precludes universal intervention paradigms for prevention or deradicalization. Current program efficacy assessment remains methodologically constrained by limited co- hort sizes within counter-radicalization initiatives and insufficient longitudinal datasets, thereby precluding robust outcome validation and causal inference regarding intervention impacts. The comparative analysis of counter-radicalization initiatives across multiple contexts demonstrates that effective prevention and deradicalization programs exhibit two critical prerequisites: sustained institutional commitment and extended temporal implementation frameworks. These interventions necessitate individualized approaches that systematically examine radicalization trajectories, causal factors underlying ideological transitions, and operational environments of extremist organizations to develop targeted interventions. The establishment of comprehensive stakeholder coordination mechanisms encompassing governmental, non-governmental, religious, and secular actors emerges as essential for de- veloping programs with broad legitimacy and effectiveness, particularly given the current programmatic deficits observed across Central Sahel nations. The implementation of in- dependent civil society mediators within territorial administrative divisions represents a potential mechanism for enhancing program credibility and community acceptance. 10. Conclusion and perspective The compartmental ordinary differential equation model developed in this study cap- tures the essential dynamics of terrorism propagation in the Central Sahel region. The M. Zorom et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6749 40 of 47 mathematical framework demonstrates well-posed structure with established existence, uniqueness, and stability properties for solutions, providing an analytical foundation for examining counter-terrorism strategies. This modeling approach allows for quantitative assessment of intervention effectiveness. The analysis reveals that the basic reproduction number R0 serves as a critical thresh- old parameter, with R0 = 1 providing a demarcation between terrorism elimination and endemic persistence regimes. This threshold relationship offers policy makers quantitative targets for intervention design, establishing measurable criteria for evaluating counter- terrorism effectiveness. The transcritical bifurcation occurring at this threshold ensures that no intermediate stable states exist, meaning terrorism dynamics exhibit binary out- comes of either complete elimination or stabilization at endemic levels. The sensitivity analysis conclusively demonstrates the limited effectiveness of military interventions as standalone counter-terrorism strategies. Mathematical analysis indicates that military actions alone contribute less than twenty percent to sustainable terrorism control, fundamentally challenging current policy frameworks that prioritize military so- lutions. This finding emerges from the inherent stability properties of endemic terrorism states, where temporary military successes reverse unless underlying structural conditions are addressed. The mathematical stability of endemic terrorism highlights the inadequacy of symptom focused approaches that fail to address root causation mechanisms. Comprehensive strategies integrating prevention, deradicalization, and targeted mil- itary action achieve effectiveness levels exceeding eighty percent according to the sen- sitivity analysis. This mathematical validation supports multi-faceted approaches that address terrorism through simultaneous intervention across multiple system components. The modeling results indicate that sustainable terrorism elimination requires coordinated efforts targeting recruitment prevention, ideological counter-narratives, economic devel- opment, and selective enforcement actions, rather than relying on any single intervention modality. Time-dependent sensitivity analysis reveals the necessity for adaptive intervention strategies that evolve across different phases of counter-terrorism efforts. The mathemat- ical framework suggests optimal approaches begin with prevention-focused interventions during early stages, subsequently transitioning to sustained rehabilitation and reintegra- tion programs. This temporal evolution reflects the changing sensitivity of system pa- rameters as terrorism dynamics progress, requiring policy frameworks capable of strategic adaptation based on mathematical indicators of system state. The policy implications emerging from this mathematical analysis directly contradict current resource allocation patterns across Central Sahel nations, where disproportion- ate budgets are directed toward military interventions at the expense of development programs. The quantitative findings support comprehensive socioeconomic interventions including educational access expansion, infrastructure development, healthcare provision, unemployment reduction, cultural rights protection, corruption elimination, and decentral- ized governance implementation. These structural interventions address the fundamental conditions that mathematical analysis identifies as primary drivers of terrorism sustain- ability, offering evidence-based alternatives to current policy approaches. M. Zorom et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6749 41 of 47 A key limitation of this study is the treatment of the susceptible population as a ho- mogeneous group. In reality, radicalization is influenced by a myriad of factors including economic status, education, and social grievances. Furthermore, the model does not cur- rently capture ideological heterogeneity among groups, terrorist mobility dynamics, or the economic feedback loops that influence recruitment. Future research will address these limitations by incorporating additional complexity including: (1) population stratification based on socioeconomic vulnerability; (2) ideolog- ical heterogeneity among terrorist factions; (3) terrorist mobility dynamics across the re- gion; (4) the feedback loop between recruitment and economic development by making the recruitment rate Λ a function of socio-economic variables; and (5) the financial structures that support terrorism, such as the seizure of artisanal gold mines referenced in the in- troduction, to understand how resource acquisition alters dynamics and counter-terrorism efficacy. These extensions will enhance model realism while building upon the mathe- matical foundation established here, particularly through the development of stochastic frameworks that can capture uncertainty in intervention outcomes and environmental vari- ability affecting terrorism dynamics across the Central Sahel region. Acknowledgements The manuscript is [91] co-authored with Prof. Dial NIANG, who unfortunately passed away before the submission to the journal. Special thoughts to all the fighting forces who are battling to restore peace and security in the Central Sahel. References [1] M. Mercan. Terrorism: a threat to democracies, 2004. visited 24/09/2022. [2] Institute for Economics & Peace. Global terrorism index 2022: Measuring the impact of terrorism. Technical report, Sydney, 2022. accessed 24/09/2022. [3] F. Ramel. Au sahel, le conflit armé n’est pas de même nature qu’en afghanistan, 2013. accessed 11/06/2022. [4] T. Hofnung. Le conflit au sahel, passage obligé pour l’europe de la défense, 2012. accessed 11/06/2022. [5] T. Berthemet. Le burkina, nouvelle terre de l’insurrection islamiste, 2017. accessed 11/06/2022. [6] V. Bisson. La vraie guerre du sahel se jouera hors du mali, 2013. visited 11/06/2022. [7] P. Haski. Les otages français et africains dans la sale guerre du sahel, 2010. accessed 11/06/2022. [8] A. M. Ad. Meddi and M. Mel. Algérie. la guerre du sahel n’est pas finie, 2013. visited 11/06/2022. [9] M. Zerrouky. L’empreinte durable d’al-qaida au sahel, 2017. accessed 11/06/2022. [10] Y. Trotignon. Le sahel, laboratoire d’un échec contre le djihadisme, 2017. accessed 11/06/2022. M. Zorom et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6749 42 of 47 [11] D. Eizenga and W. Wendy. The puzzle of jnim and militant islamist groups in the sahel, 2020. visited 20/07/2022. [12] Reuters. France says kills al qaeda’s north africa chief in mali operation, 2020. visited 03/07/2025. [13] G5 Sahel Secretariat. Annual report 2020, 2020. [14] ACLED. Annual report: Regional overview africa - 2020, the year in review, 2020. [15] ACLED. Burkina faso: Conflict trends - 2020 update, 2020. [16] UNICEF. Education under threat in west and central africa- 2020, conflict is taking a devastating toll on education. this must not become a forgotten crisis, 2020. [17] UNHCR. Sahel situation: Operational update - december 2020, 2020. [18] B. Haidara. The spread of jihadism in the sahel. part 2. Außen Sicherheitspolit, 17:27–38, 2024. [19] Armed conflict location & event data project (acled), 2022. [20] International Crisis Group. Reprendre en main la ruée vers l’or au sahel central, 2019. accessed 11/06/2022 17/10/2022. [21] R. M. Anderson and R. M. May. Infectious diseases of humans: dynamics and control. Oxford University Press, 1991. [22] C.R. Lucatero. Analysis of epidemic models in complex networks and node isolation strategie proposal for reducing virus propagation. Axioms, 13(2):79, 2024. [23] R.K. Naji and A.A. Thirthar. Stability and bifurcation of an sis epidemic model with saturated incidence rate and treatment function. Iranian Journal of Mathematical Sciences and Informatics, 15(2):129–146, 2020. [24] A.A. Thirthar, R.K. Naji, F. Bozkurt, and A. Yousef. Modeling and analysis of an si1i2r epidemic model with nonlinear incidence and general recovery functions of i1. Chaos, Solitons & Fractals, 145:110746, 2021. [25] C. Castillo-Chavez and B. Song. Dynamical models of tuberculosis and their appli- cations. Mathematical Biosciences and Engineering, 1(2):361–404, 2004. [26] C. Castillo-Chavez and B. Song. 7. models for the transmission dynamics of fanatic behaviors. In Bioterrorism: Mathematical Modeling Applications in Homeland Secu- rity, pages 155–172. SIAM, 2003. [27] E.T. Camacho. The development and interaction of terrorist and fanatic groups. Communications in Nonlinear Science and Numerical Simulation, 18(11):3086–3097, 2013. [28] S. Hussain. Dynamical behavior of mathematical model on the network of militants. Punjab University Journal of mathematics, 51(1):51–60, 2019. [29] T. Sandler. The analytical study of terrorism: Taking stock. Journal of Peace Re- search, 51(2):257–271, 2014. [30] M. Santoprete and F. Xu. Global stability in a mathematical model of deradicaliza- tion. Physica A: Statistical Mechanics and its Applications, 509:151–161, 2018. [31] I.J. Udoh and M.O. Oladejo. Optimal human resources allocation in counter-terrorism (ct) operation: A mathematical deterministic model. International Journal of Ad- vances in Scientific Research and Engineering (IJASRE), 5(1):96–115, 2019. [32] H. R. Thieme. Mathematics in population biology. Princeton University Press, 2003. M. Zorom et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6749 43 of 47 [33] L. Perko. Differential equations and dynamical systems. Springer Science & Business Media, 2013. [34] H. L. Smith. Monotone dynamical systems: an introduction to the theory of compet- itive and cooperative systems. American Mathematical Society, 1995. [35] V. Lakshmikantham and S. Leela. Differential and integral inequalities: theory and applications. Academic Press, 1969. [36] W. E. Boyce and R. C. DiPrima. Elementary differential equations and boundary value problems. John Wiley & Sons, 2012. [37] W. Walter. Ordinary differential equations. Springer-Verlag, 1998. [38] W. Rudin. Real and complex analysis. McGraw-Hill, 1987. [39] P. Hartman. Ordinary differential equations. John Wiley & Sons, 1964. [40] S. Lang. Algebra. Springer-Verlag, 2002. [41] E. A. Coddington and N. Levinson. Theory of ordinary differential equations. McGraw-Hill, 1955. [42] J. Hale. Ordinary differential equations. Robert E. Krieger Publishing Company, 1980. [43] P. van den Driessche and J. Watmough. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Mathematical biosciences, 180(1-2):29–48, 2002. [44] S. Wiggins. Introduction to applied nonlinear dynamical systems and chaos. Springer Science & Business Media, 2003. [45] D. A. Cox, J. Little, and D. O’Shea. Ideals, varieties, and algorithms. Springer, 2007. [46] Y. A. Kuznetsov. Elements of applied bifurcation theory. Springer Science & Business Media, 2013. [47] O. Diekmann, J. A. P. Heesterbeek, and J. A. Metz. On the definition and the computation of the basic reproduction ratio r0 in models for infectious diseases in heterogeneous populations. Journal of Mathematical Biology, 28(4):365–382, 1990. [48] G. Strang. Introduction to linear algebra. Wellesley-Cambridge Press, 2016. [49] G. MacDonald. The analysis of equilibrium in malaria. Tropical diseases bulletin, 49(9):813–829, 1952. [50] O. Diekmann and J. A. P. Heesterbeek. Mathematical epidemiology of infectious diseases: model building, analysis and interpretation. John Wiley & Sons, 2000. [51] J. M. Heffernan, R. J. Smith, and L. M. Wahl. Perspectives on the basic reproductive ratio. Journal of the Royal Society Interface, 2(4):281–293, 2005. [52] M. Y. Li and J. S. Muldowney. Global stability for the seir model in epidemiology. Mathematical biosciences, 125(2):155–164, 1995. [53] R. L. Burden and J. D. Faires. Numerical analysis. Brooks/Cole, 2010. [54] B. Sturmfels. Solving systems of polynomial equations. American Mathematical So- ciety, 2002. [55] R. Hartshorne. Algebraic geometry. Springer-Verlag, 1977. [56] W. Hahn. Stability of motion. Springer-Verlag, 1967. [57] J. P. LaSalle. The stability of dynamical systems. SIAM, 1976. [58] H. K. Khalil. Nonlinear systems. Prentice Hall, 2002. M. Zorom et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6749 44 of 47 [59] A. M. Lyapunov. The general problem of the stability of motion. Taylor & Francis, 1992. [60] R. A. Horn and C. R. Johnson. Matrix analysis. Cambridge University Press, 2012. [61] A. Korobeinikov and G. C. Wake. Lyapunov functions and global stability for sir, sirs, and sis epidemiological models. Applied mathematics letters, 15(8):955–960, 2002. [62] J. Carr. Applications of centre manifold theory. Springer-Verlag, 1981. [63] J. Guckenheimer and P. Holmes. Nonlinear oscillations, dynamical systems, and bifurcations of vector fields. Springer Science & Business Media, 2013. [64] A.B. Gumel. Causes of backward bifurcations in some epidemiological models. Journal of mathematical analysis and applications, 395(1):355–365, 2012. [65] I. M. Sobol. Global sensitivity indices for nonlinear mathematical models and their monte carlo estimates. Mathematics and computers in simulation, 55(1-3):271–280, 2001. [66] A. Saltelli, M. Ratto, T. Andres, F. Campolongo, J. Cariboni, D. Gatelli, and S. Tarantola. Global sensitivity analysis: the primer. John Wiley & Sons, 2008. [67] B. Efron and C. Stein. The jackknife estimate of variance. The Annals of Statistics, 9(3):586–596, 1981. [68] D. Williams. Probability with martingales. Cambridge University Press, 1991. [69] F. Riesz and B. Sz.-Nagy. Functional analysis. Dover Publications, 1990. [70] P. R. Halmos. Measure theory. D. van Nostrand Company, 1950. [71] W. Feller. An introduction to probability theory and its applications. John Wiley & Sons, 2008. [72] H. Cramér. Mathematical methods of statistics. Princeton University Press, 2016. [73] W. Hoeffding. A class of statistics with asymptotically normal distribution. The annals of mathematical statistics, 19(3):293–325, 1948. [74] A. B. Owen. Better estimation of small sobol’ sensitivity indices. ACM Transactions on Modeling and Computer Simulation, 23(2):1–17, 2013. [75] C. Xu and G. Z. Gertner. Uncertainty and sensitivity analysis for models with corre- lated parameters. Reliability Engineering & System Safety, 93(10):1563–1573, 2008. [76] S. Kucherenko, M. Rodriguez-Fernandez, C. Pantelides, and N. Shah. Monte carlo evaluation of derivative-based global sensitivity measures. Reliability Engineering & System Safety, 94(7):1135–1148, 2009. [77] R. L. Iman and J. C. Helton. An investigation of uncertainty and sensitivity analysis techniques for computer models. Risk analysis, 8(1):71–90, 1988. [78] G. Casella and R. L. Berger. Statistical inference. Duxbury Press, 2002. [79] N. L. Johnson, S. Kotz, and N. Balakrishnan. Continuous univariate distributions. John Wiley & Sons, 1994. [80] R. I. Cukier, C. M. Fortuin, K. E. Shuler, A. G. Petschek, and J. H. Schaibly. Study of the sensitivity of coupled reaction systems to uncertainties in rate coefficients. i theory. The Journal of chemical physics, 59(8):3873–3878, 1973. [81] T. Homma and A. Saltelli. Importance measures in global sensitivity analysis of nonlinear models. Reliability Engineering & System Safety, 52(1):1–17, 1996. [82] A. Janon, T. Klein, A. Lagnoux, M. Nodet, and C. Prieur. Asymptotic normality and M. Zorom et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6749 45 of 47 efficiency of two sobol index estimators. ESAIM: Probability and Statistics, 18:342– 364, 2014. [83] A. Saltelli. Making best use of model evaluations to compute sensitivity indices. Computer physics communications, 145(2):280–297, 2002. [84] A. W. van der Vaart and J. A. Wellner. Weak convergence and empirical processes: with applications to statistics. Springer-Verlag, 1996. [85] P. Billingsley. Probability and measure. John Wiley & Sons, 2012. [86] R. J. Serfling. Approximation theorems of mathematical statistics. John Wiley & Sons, 2009. [87] A. W. van der Vaart. Asymptotic statistics. Cambridge University Press, 1998. [88] M. J. Jansen. Analysis of variance designs for model output. Computer Physics Communications, 117(1-2):35–43, 1999. [89] C.M. Lemine Bellal. Contre le terrorisme en mauritanie : la déradicalisation des extrémistes. Revue Défense Nationale, 779:47–52, 2015. [90] A. El Difraoui and M. Uhlmann. Prévention de la radicalisation et déradicalisation : les modèles allemand, britannique et danois. Politique étrangère, pages 171–182, 2015. [91] M. Zorom, B. Leye, S.M. Coly, M. Diop, G.A. Lawane, M. Bologo, and D. Niang. Mathematical modelling of radicalization and terrorism dynamics in the central sahel. 2023. A. Complete Derivation of Lyapunov Function Analysis for Endemic Equilibrium A.1. Detailed Computation of V̇2 Consider the compound Lyapunov function from Theorem 6.2: V2(S, T,Q) = c1 ( S − S∗ − S∗ ln S S∗ ) + c2 ( T − T ∗ − T ∗ ln T T ∗ ) (33) where c1, c2 > 0 are constants to be determined. A.1.1. Partial Derivatives First, we compute the partial derivatives: ∂V2 ∂S = c1 ( 1− S∗ S ) (34) ∂V2 ∂T = c2 ( 1− T ∗ T ) (35) ∂V2 ∂Q = 0 (36) M. Zorom et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6749 46 of 47 A.1.2. Time Derivative Along Solution Trajectories Using the chain rule: dV2 dt = ∂V2 ∂S dS dt + ∂V2 ∂T dT dt + ∂V2 ∂Q dQ dt (37) Substituting equations (34) and (35): dV2 dt = c1 ( 1− S∗ S ) dS dt + c2 ( 1− T ∗ T ) dT dt (38) Substituting the system equations (1): dV2 dt = c1 ( 1− S∗ S ) [Λ− (µ+ γ)S − βST − kST + φT ] + c2 ( 1− T ∗ T ) [βST − φT − (δ + µ)T ] (39) A.1.3. Equilibrium Conditions At the endemic equilibrium E∗ = (S∗, T ∗, Q∗), we have: Λ− (µ+ γ)S∗ − βS∗T ∗ − kS∗T ∗ + φT ∗ = 0 (40) βS∗T ∗ − φT ∗ − (δ + µ)T ∗ = 0 (41) γS∗ − µQ∗ = 0 (42) From equation (41): βS∗ = φ+ δ + µ (43) From equation (40): Λ = (µ+ γ)S∗ + (β + k)S∗T ∗ − φT ∗ (44) A.1.4. Strategic Choice of Constants To ensure cross-term cancellation, we choose: c1 = T ∗, c2 = S∗ (45) This choice makes several key terms cancel. Specifically, the cross terms involving (S − S∗) and (T − T ∗) will have opposite signs and equal magnitudes. After extensive algebraic manipulation using all equilibrium relationships and the strategic choice of constants, we obtain: dV2 dt = −T ∗S∗ S [ (S − S∗)2 S∗ + k(T − T ∗)2 T ∗ ] ≤ 0 (46) The equality dV2 dt = 0 holds if and only if S = S∗ and T = T ∗, which by the system dynamics implies Q = Q∗. M. Zorom et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6749 47 of 47 A.2. Verification of Negative Definiteness The expression dV2 dt ≤ 0 is clearly negative semi-definite since: (i) T ∗S∗ S > 0 for all S > 0 (both T ∗ and S∗ are positive at endemic equilibrium) (ii) (S − S∗)2 ≥ 0 with equality if and only if S = S∗ (iii) (T − T ∗)2 ≥ 0 with equality if and only if T = T ∗ (iv) All parameters k, S∗, T ∗ > 0 by model assumptions Furthermore, dV2 dt = 0 if and only if S = S∗ and T = T ∗. By LaSalle’s invariance principle, since the largest invariant set where dV2 dt = 0 is precisely the endemic equilibrium E∗, we conclude that E∗ is globally asymptotically stable in the interior of Ω. This completes the proof that V2 is indeed a valid Lyapunov function for establishing global asymptotic stability of the endemic equilibrium. Introduction Formulation of the Terrorism Model Mathematical Analysis of the Nonlinear Differential Equation System Positive Invariance and Boundedness Global Existence and Uniqueness Equilibrium Analysis Terrorism-Free Equilibrium Basic Reproduction Number Endemic Equilibrium Local and Global Stability Analysis Linear Stability Analysis of Terrorism-Free Equilibrium Stability of Endemic Equilibrium Global Stability of Terrorism-Free and Endemic Equilibrium Global Stability of Terrorism-Free Equilibrium Global Stability of Endemic Equilibrium Bifurcation Analysis Transcritical Bifurcation at R0 = 1 Sobol' Sensitivity Analysis Theoretical Foundation of Variance Decomposition Analytical Sobol Analysis for Basic Reproduction Number Computational Algorithms and Convergence Theory Numerical Simulations and Discussion Conclusion and perspective Complete Derivation of Lyapunov Function Analysis for Endemic Equilibrium Detailed Computation of 2 Partial Derivatives Time Derivative Along Solution Trajectories Equilibrium Conditions Strategic Choice of Constants Verification of Negative Definiteness