EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 2, Article Number 5933 ISSN 1307-5543 – ejpam.com Published by New York Business Global Mathematical SEIR Model of the Lumpy Skin Disease Using Caputo-Fabrizio Fractional-Order Rajagopalan Ramaswamy1,∗, Gunaseelan Mani2, Radhakrishnan Mohanraj3,4, Ozgur Ege5 1 Department of Mathematics, College of Science and Humanities in Alkharj, Prince Sattam Bin Abdulaziz University, Alkharj 11942, Saudi Arabia 2 Department of Mathematics, Saveetha School of Engineering, Saveetha Institute of Medical and Technical Sciences, Chennai 602105, India 3 Directorate of Learning and Development, SRM Institute of Science and Technology, Kattankulathur, Chennai, 603203, Tamil Nadu, India 4 Department of Mathematics, SRM Institute of Science and Technology, Kattankulathur, Chennai, 603203, Tamil Nadu, India 5 Department of Mathematics, Ege University, Bornova, Izmir, 35100, Turkey Abstract. UN SD Goal-3 focuses on good health and well-being, while the 15th goal focuses on life on land. This research aims to study an integral model of Lumpy Skin Disease to establish better knowledge of this condition. In this study, a Caputo-Fabrizio fractional-order model for the dynamics of Lumpy Skin Disease (LSD) is analysed. Using the Picard-Lindelöf theorem in the Banach space C([0, 1],R) with the supremum norm ∥ϕ∥ = maxs∈[0,1] |ϕ(s)|, we prove existence and uniqueness of solutions via fixed-point theory. The positivity, boundedness, and Ulam-Hyers stability of the model are established. Disease-free and endemic equilibria are derived, the basic reproduction number R0 is computed, and sensitivity analysis is conducted to identify critical transmission parameters. A novel Newton interpolation-based numerical scheme is developed and compared with finite difference method, with superior convergence rates being demonstrated. How non-integer operators enhance LSD progression modeling compared to classical approaches is re- vealed by fractional-order simulations. Actionable insights for disease control are provided by our results while the efficacy of Caputo-Fabrizio derivatives in epidemiological modeling is demon- strated which will help in framing policies to have a healthy society, SDG-3 and other related goals of UN. 2020 Mathematics Subject Classifications: 65L60, 41A10, 35N70, 65L20 Key Words and Phrases: Fixed point theory, Caputo-Fabrizio, Ulam-Hyers stability ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i2.5933 Email addresses: r.gopalan@psau.edu.sa (R Ramaswamy), mathsguna@yahoo.com (G.Mani), radhakrm@srmist.edu.in (R.Mohanraj), ozgur.ege@ege.edu.tr (O.Ege) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) R. Ramaswamy et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5933 2 of 25 1. Introduction The United nations [1] has set seventeen goals for sustaintable development (SDG) and SDG-3 is related to Health while SDG-15 is focussed on Life on Land. Living be- ings on earth are prone to infections and one such disease is Lumpy Skin Disease (LSD), that exists as an infection caused by the Lumpy Skin Disease Virus (LSDV) which pos- sesses 150 kb of double-stranded DNA. The virus belongs to Capripoxvirus genera within the Chordopoxviridae subfamily of Poxviridae family [2, 3]. The genus Capripoxvirus comprises three confirmed viral members: Goatpox Virus (GTPV) and Sheeppox Virus (SPPV) alongside Lumpy Skin Disease Virus (LSDV). Electron microscope images show LSDV exhibits morphological patterns shared by other Poxviridae viruses similar to the vaccinia virus [4]. The virus has limited host range due to its restriction of transmis- sion to two main hosts including cattle (Bos indicus or Bos taurus) and domestic water buffaloes (Bubalus bubalis) [5–7]. The Range Management Pathogens Research Unit in Georgia has studied how LSDV affects wild mammals including giraffes and camels as well as wildebeests [8–10]. Active transmission of disease occurs when animals transmit the virus through skin contact with lesions and milk and through the bloodfeeding activities of biting flies and mosquitoes and ticks [11–14]. The disease transmission rate is greater for these seasons be- cause insect vectors are both more numerous and active [15]. Skilled shedding by the virus happens through nasal discharge along with skin lesions and saliva as well as lachrymal secretions to transmit disease from ill to well animals [12, 16]. An outbreak of LSD originated in India during 2019 followed by multiple serious oc- currences. The newest spell of LSD outbreak started during May 2022 and extensively damaged 15 Indian states that led to about 100,000 cattle fatalities and considerable eco- nomic loss [17]. Within India’s highly important livestock production industry LSD leads to animal mortality along with milk yield declines while creating additional economic im- pacts through movement limitations [17]. Under India’s 308 million cattle population the prevention and control of infectious diseases including LSD stands as a vital necessity [18]. Comparable changes occur within diseased animals because of LSD resulting in mastitis along with necrotic hepatitis and lymphadenitis and orchitis and myocardial damage [19]. WHO has established LSD as an eligible disease for notification according to[20]. Surveillance confirmed the first documented case of LSD in Zambia in 1931 but the dis- ease stayed within Sub-Saharan Africa until 1989 when it expanded into the Middle East and Asia [21–23]. Initial identification of LSDV occurred in Russia when the disease appeared in Southeast Europe in 2016 followed by its discovery in India and several neigh- boring Asian countries such as China and Nepal Thailand and Bangladesh plus Bhutan in November 2019 [24, 25]. From 2019 onward LSD spread rapidly to impact more than two million cattle during 2022. The symptoms typically occur within 1-4 weeks after a patient contracts the infection and high fever joins nasal or eye drainage along with appetite loss and nodular skin breakouts [26]. Public deaths from LSD have been reported between 5-45% during outbreaks in Rajasthan, Gujarat, Uttar Pradesh, and Punjab, Haryana, Karnataka along with West Bengal and Maharashtra [27–29]. R. Ramaswamy et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5933 3 of 25 Indian authorities have built a series of initiatives aimed at stopping animal disease spread which include massive vaccination programs combined with animal sequestration areas and restrictions on animal movement. Early detection remains difficult together with limited staff knowledge about disease spread across affected zones. Since its emer- gence from prior poxvirus species, LSDV expanded its network of hosts and harm potential through double-stranded DNA virus homologous recombination processes [30]. Through genome sequencing the study team identified current LSDV genetic variants that existed in India. Phylogenetic examination established two distinguishable variant classes having important variations in their mutation (SNP) pattern. Yasir et al. [31] studied the compe- tition model in fractional derivatives of Caputo Fabrizio type through implementation of Newton based polynomials. Singh et al. [32] studied the fractional-order compartmental model for cervical cancer. HPV serves as the cause for cervical cancer to develop. This worldwide phenomenon receives attention through the utilization of the Caputo-Fabrizio fractional operator that contains antiretroviral medication treatment sections. Kumari et al. [33] established a mathematical model for coronavirus using the Caputo-Fabrizio fractional derivative. The diabetes model and its complications were examined by Singh et al. [34] using Caputo-Fabrizio fractional derivatives. A new model of Rabies disease analysis conducted by Aydogan et al. [35] utilized Caputo-Fabrizio fractional derivatives. Butt [36] proposed a new fractional model of LSD dynamics which utilized the Atangana Baleanu fractional derivative to model disease transmission patterns and disease memory effects. To develop a model for Lumpy Skin Disease (LSD), the total populationN (s) is catego- rized into four distinct compartments: The population contains four subsections including susceptible S(s) exposed E(s) infected I(s) and recovered R(s). Therefore, at any given time t, the total population is expressed as: N (s) = S(s) + E(s) + I(s) +R(s). The susceptible class S(s) contains animals which can develop virus illnesses from con- tact with virus carriers. Exposed cattle arise through contacts between susceptible cattle and virus-infected cattle who advance to the exposed class E(s). Exposure indicates that cattle within this population have become infected but showcase no contagious character- istics at present.can become ill due to the interaction with the cattle carrying virus. When cattle from the susceptible class interact with the infections cattle become exposed cattle and hence move to the exposed class E(s). So, the exposed class is the class consisting of cattle that have become infected but not infectious till now. The following stage comes in the form of an infected set of cattle which the symbol I(s) defines. Any virus which estab- lishes itself in exposed cattle creates infectious beings which transition from the exposed class E(s) to the infected class I(s). Pros of virus transmission exist whenever a cattle that is currently infectious becomes able to replicate and spread the virus. At the next stage exists the recovered class denoted as R(s). Higher levels of immunity exist in some cattle population leading them to recover from the disease more efficiently. Proper medi- cation against infections allows these cattle to enter the recovered class. Such cattle types make up this category. All state variables S(s), E(s), I(s), R(s) meet the requirement of R. Ramaswamy et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5933 4 of 25 continuous differentiability when considered over t in the time interval [0,∞). The flow pattern for LSD is shown in the Fig. 1 which in the form of nonlinear ordinary differential equations is given as: dS ds = αN − βSI + θR− ϕS, dE ds = βSI − γE − ϕE , dI ds = γE − ωI − ϕI, dR ds = ωI − θR− ϕR, along with the following non-negative initial conditions. S(0) = S0, E(0) = E0, I(0) = I0,R(0) = R0. The rest of the paper is organised as follows: In Section 2 we present the LSD model Figure 1: Flow chart Table 1: Description of the parameters used in the model (1) Parameter Description Values Source α Birth rate 0.0114 [37] β Infection rate 0.00011412 [38] γ Progression rate from exposed to infected 0.95 [38] ω Recovery rate 0.14 Assumed ϕ Natural death rate 0.0057 Assumed θ Loss of immunity rate 0.007 Assumed with (CF) derivative and review some model basic preliminaries and present and establish the results for positivity, boundedness, existence and uniqueness of the solutions for the model. In section-3 we analyyse the Ulam-Hyers stability and investigate the disease-free and disease-endemic equilibrium points in section-4. The numerical scheme is presented in Section-5 and the results are discussed in Section-6, which we believe will help in R. Ramaswamy et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5933 5 of 25 formulating policies for restricting the spread of LSD and thus attain the SDG-3 as well as SDG-15 and other related goals of [1]. Finally we conclude the article with scope for further research. 2. Lumpy Skin Disease with Caputo-Fabrizio The lumpy skin disease model with Caputo-Fabrizio(CF) derivative is given by CFDa sS = αN − βSI + θR− ϕS, CFDa sE = βSI − γE − ϕE , CFDa sI = γE − ωI − ϕI, CFDa sR = ωI − θR− ϕR, (1) with a being the fractional order 0 < a < 1 subject to the following initial conditions S(0) = S0, E(0) = E0, I(0) = I0,R(0) = R0. (2) 2.1. Model basic preliminaries Here, we give some of the mathematical preliminaries in the form of theorems, which we shall apply to prove the positivity and uniqueness and positivity of lumpy skin disease model with Caputo-Fabrizio(CF) (1) as defined [39, 40] respectively. More details, one can read [41–43]. We assume the space {ϕ(s) ∈ C([0, 1] → R)} with ||ϕ|| = max s∈[0,1] |ϕ(s)|. (3) The definitions are stated as follows. Definition 1. Let ϕ : R+ → R and a ∈ (n − 1, n), n ∈ N. The left Caputo fractional derivative of order a of the function ϕ is given by the following equality: C 0Da s (ϕ(s)) = 1 Γ(n− a) ∫ s 0 (s− s)n−a−1ϕ(n)(s)ds, s > 0. Definition 2. The fractional integral associate to the new fractional integral with power- law kernel is defined as: C 0Ia s (ϕ(s)) = 1 Γ(α) ∫ s 0 (s− s)a−1ϕ(τ)dτ. Definition 3. Assume ϕ(s) ∈ H1(l1, l2), for l2 > l1, τ ∈ [0, 1]. The CF fractional operator is given as CFDa s (ϕ(s)) = M(a) (1− a) ∫ l2 l1 ϕ ′ (σ)exp ( − a s− σ 1− σ ) d(σ), 0 < a < 1, R. Ramaswamy et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5933 6 of 25 = dϕ ds , a = 1, (4) where M(a) = 1 − a + a Γ(a) is a normalization function which satisfies the condition, M(0) = M(1) = 1. Definition 4. The integral operator of fractional order corresponding to the CF fractional derivative is stated as follows J a s (ϕ(s)) = 2(1− a) (2− a)M(a) ϕ(s) + 2a (2− a)M(a) ∫ s 0 ϕ(ξ)dξ, s ≥ 0. (5) Definition 5. The Laplace transform of CFDa sϕ(s) is represented as follows L[CFDa sϕ(s)] = M(a) κL[−ϕ(s)]− ϕ(0) κ+ a(1− κ) . (6) Theorem 1. [44](Banach Contraction Theorem) Let (P, d) be a complete metric space, and let H : P → P be a contraction on P. Then H has a unique fixed point, say u ∈ P. 2.2. Positivity of Solutions Theorem 2. For non-negative initial conditions S(0), E(0), I(0),R(0) ≥ 0, the solutions of system (1) remain non-negative for all s > 0. Proof. We employ the generalized mean value theorem for Caputo-Fabrizio derivatives:  CFDa sS ∣∣∣∣ S=0 = αN + θR ≥ 0, CFDa sE ∣∣∣∣ E=0 = βSI ≥ 0, CFDa sI ∣∣∣∣ I=0 = γE ≥ 0, CFDa sR ∣∣∣∣ R=0 = ωI ≥ 0. Since all derivatives are non-negative at the boundary of each compartment, the solutions remain non-negative for all s ≥ 0. R. Ramaswamy et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5933 7 of 25 2.3. Boundedness of Solutions Theorem 3 (Boundedness). The total population N (s) is bounded such that: N (s) ≤ max ( N (0), αN ϕ ) , ∀s ≥ 0 Proof. From the system (1), the total population N (s) satisfies: CFDa sN = αN − ϕN − ωI. Since −ωI ≤ 0, we have the inequality: CFDa sN ≤ (α− ϕ)N . Using the Caputo-Fabrizio (CF) fractional derivative and its corresponding integral operator, the equivalent integral form is: N (s)−N (0) = 2(1− a) (2− a)M(a) (α− ϕ)N + 2a (2− a)M(a) ∫ s 0 (α− ϕ)N (ξ)dξ. Substituting M(a) = 1− a+ a Γ(a) , this becomes: N (s)−N (0) = 2(1− a)(α− ϕ)N (2− a) ( 1− a+ a Γ(a) ) + 2a(α− ϕ) (2− a) ( 1− a+ a Γ(a) ) ∫ s 0 N (ξ)dξ. so the equation simplifies to: N (s)−N (0) = κ(a)(1− a)(α− ϕ)N + κ(a)a(α− ϕ) ∫ s 0 N (ξ)dξ. (7) where, κ(a) = 2 (2− a) ( 1− a+ a Γ(a) ) , Case 1: α ≤ ϕ in (7) If α ≤ ϕ, then (α− ϕ) ≤ 0. The right-hand side of the inequality: CFDa sN ≤ (α− ϕ)N is non-positive. Thus, N (s) cannot increase beyond its initial value: N (s) ≤ N (0), ∀s ≥ 0. R. Ramaswamy et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5933 8 of 25 Case 2: α > ϕ in (7) If α > ϕ, then (α− ϕ) > 0. Rearrange the simplified equation: N (s)−N (0) = κ(a)(1− a)(α− ϕ)N + κ(a)a(α− ϕ) ∫ s 0 N (ξ)dξ. Let U = 1− κ(a)(1− a)(α− ϕ), assuming K > 0 (valid for small enough a). Dividing through by K, we get: N (s) ≤ N (0) U + κ(a)a(α− ϕ) U ∫ s 0 N (ξ)dξ. Applying a Grönwall-type inequality, this yields: N (s) ≤ N (0) exp ( κ(a)a(α− ϕ) U s ) . As s → ∞, N (s) converges to a steady state determined by the balance between growth (α) and decay (ϕ): N (s) ≤ max ( N (0), αN ϕ ) . This completes the proof. 2.4. Existence and uniqueness of solutions of the CF model This section proves the existence and uniqueness of the solution to the proposed model (1) using fixed point theorem. To facilitate this analysis, model (1) can be expressed as follows:  CFDa sS = K1(s,S), CFDa sE = K2(s, E), CFDa sI = K3(s, I), CFDa sR = K4(s,R), (8) where  K1(s,S) = αN − βSI + θR− ϕS, K2(s, E) = βSI − γE − ϕE , K3(s, I) = γE − ωI − ϕI, K4(s,R) = ωI − θR− ϕR. (9) R. Ramaswamy et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5933 9 of 25 Applying fractional integral operator given in (5), system (8) reduces to the Volterra integral type of order 0 < a < 1 given by S(s) = S(0) + 2(1−a) (2−a)M(a)K1(s,S) + 2a (2−a)M(a) ∫ s 0 K1(ξ,S)dξ, E(s) = E(0) + 2(1−a) (2−a)M(a)K2(s, E) + 2a (2−a)M(a) ∫ s 0 K2(ξ, E)dξ, I(s) = I(0) + 2(1−a) (2−a)M(a)K3(s, I) + 2a (2−a)M(a) ∫ s 0 K3(ξ, I)dξ, R(s) = R(0) + 2(1−a) (2−a)M(a)K4(s,R) + 2(1−a) (2−a)M(a) ∫ s 0 K4(ξ,R)dξ. (10) (G :) For proving our results, we consider the following assumption: For this, S(s),S⋆(s), E(s), E⋆(s), I(s), I⋆(s),R(s),R⋆(s) ∈ L[0, 1] be continuous such that ||S(s)|| ≤ c1, ||E(s)|| ≤ c2, ||I(s)|| ≤ c3 and ||R(s)|| ≤ c4, for some positive constants c1, c2, c3, c4 > 0. Theorem 4. Each kernel (K1,K2,K3,K4) satisfy the Lipschitz condition under the as- sumption (G) and ji < 1, for i = 1, 2, 3, 4. Proof. For S and S⋆, we have from (9), gives ||K1(s,S)−K1(s,S⋆)|| = ||αN − βSI + θR− ϕS − (αN − βS⋆I + θR− ϕS⋆)|| ≤ ||βI(s)||||S(s)− S⋆(s)||+ ϕ||S(s)− S⋆(s)|| ≤ (c3β + ϕ)||S(s)− S⋆(s)|| = j1||S(s)− S⋆(s)||, (11) where j1 = c3β + ϕ. For E and E⋆, we have ||K2(s, E)−K2(s, E⋆)|| = ||βSI − γE − ϕE − (βSI − γE⋆ − ϕE⋆)|| ≤ (γ + ϕ)||E(s)− E⋆(s)|| = j2||E(s)− E⋆(s)||, (12) where j2 = γ + ϕ. For I and I⋆, we have ||K3(s, I)−K3(s, I⋆)|| = ||γE − ωI − ϕI − (γE − ωI⋆ − ϕI⋆)|| ≤ (ω + ϕ)||I(s)− I⋆(s)|| = j3||I(s)− I⋆(s)||, (13) where j3 = ω + ϕ. For R and R⋆, we have ||K4(s,R)−K4(s,R⋆)|| = ||ωI − θR− ϕR− (ωI − θR⋆ − ϕR⋆)|| ≤ (θ + ϕ)||R(s)−R⋆(s)|| = j4||R(s)−R⋆(s)||, (14) where j4 = θ + ϕ. Thus, from (11)-(14), we have that Ki for i = 1, 2, 3, 4, satisfy the Lipschitz property. R. Ramaswamy et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5933 10 of 25 Theorem 5. There is at least a solution of the system (1) if δ = max{j1, j2, j3, j4} < 1. Proof. Let λ1k(s) = Sk+1(s)− S(s), λ2k(s) = Ek+1(s)− E(s), λ3k(s) = Ik+1(s)− I(s), λ4k(s) = Rk+1(s)−R(s). Then, we have ||λ1k(s)|| = ||Sk+1(s)− S(s)|| = ∥∥∥∥ 2(1− a) (2− a)M(a) (K1(s,Sk)−K1(s,S)) + 2a (2− a)M(a) ∫ s 0 (K1(ξ,Sk)−K1(ξ,S))dξ ∥∥∥∥ ≤ 2(1− a) (2− a)M(a) ∥∥K1(s,Sk)−K1(s,S) ∥∥ + 2a (2− a)M(a) ∫ s 0 ∥∥K1(ξ,Sk)−K1(ξ,S)) ∥∥dξ ≤ 2(1− a) (2− a)M(a) j1||Sk − S|| + 2a (2− a)M(a) ∫ s 0 j1||Sk − S||dξ = 2(1− a) (2− a)M(a) j1||Sk − S|| + 2a (2− a)M(a) j1 ∫ s 0 ||Sk − S||dξ ≤ ( 2(1− a) (2− a)M(a) + 2as (2− a)M(a) ) j1||Sk − S|| ≤ ( 2(1− a) (2− a)M(a) + 2as (2− a)M(a) )n jn1||S1 − S||. Since j1 < 1. As n → ∞, we have Sn → S. Similarly, ||λ2k(s)|| ≤ ( 2(1− a) (2− a)M(a) + 2as (2− a)M(a) )n jn2||E1 − E|| ||λ3k(s)|| ≤ ( 2(1− a) (2− a)M(a) + 2as (2− a)M(a) )n jn3||I1 − I|| ||λ4k(s)|| ≤ ( 2(1− a) (2− a)M(a) + 2as (2− a)M(a) )n jn4||R1 −R||. As n → ∞, we get λik(s) → 0 with ji < 1 for i = 1, 2, 3, 4. Hence, the system (1) has a solution. R. Ramaswamy et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5933 11 of 25 Theorem 6. The system (1) has a unique solution if( 2(1− a) (2− a)M(a) + 2as (2− a)M(a) ) ji ≤ 1, for i = 1, 2, 3, 4. Proof. Assume that there exists another solution S⋆(s), E⋆(s), I⋆(s),R⋆(s) with initial values such that S⋆(s) = S(0) + 2(1− a) (2− a)M(a) K1(s,S⋆) + 2a (2− a)M(a) ∫ s 0 K1(ξ,S⋆)dξ, E⋆(s) = E(0) + 2(1− a) (2− a)M(a) K2(s, E⋆) + 2a (2− a)M(a) ∫ s 0 K2(ξ, E⋆)dξ, I⋆(s) = I(0) + 2(1− a) (2− a)M(a) K3(s, I⋆) + 2a (2− a)M(a) ∫ s 0 K3(ξ, I⋆)dξ R⋆(s) = R(0) + 2(1− a) (2− a)M(a) K4(s,R⋆) + 2a (2− a)M(a) ∫ s 0 K4(ξ,R⋆)dξ. Now, ||S − S⋆|| = ∥∥∥∥ 2(1− a) (2− a)M(a) (K1(s,S)−K1(s,S⋆)) + 2a (2− a)M(a) ∫ s 0 (K1(ξ,S)−K1(ξ,S⋆))dξ ∥∥∥∥ ≤ 2(1− a) (2− a)M(a) ∥∥K1(s,S)−K1(s,S⋆) ∥∥ + 2a (2− a)M(a) ∫ s 0 ∥∥K1(ξ,S)−K1(ξ,S⋆)) ∥∥dξ ≤ 2(1− a) (2− a)M(a) j1||S − S⋆|| + 2a (2− a)M(a) ∫ s 0 j1||S − S⋆||dξ = 2(1− a) (2− a)M(a) j1||S − S⋆|| + 2a (2− a)M(a) j1 ∫ s 0 ||S − S⋆||dξ ≤ ( 2(1− a) (2− a)M(a) + 2as (2− a)M(a) ) j1||S − S⋆|| ≤ ( 2(1− a) (2− a)M(a) + 2as (2− a)M(a) ) j1||S − S⋆||, which implies that( 1− ( 2(1− a) (2− a)M(a) + 2as (2− a)M(a) ) j1 ) ||S − S⋆|| ≤ 0. R. Ramaswamy et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5933 12 of 25 Therefore, ||S − S⋆|| = 0. Hence, S = S⋆. Similarly, we can prove V = V⋆, E = E⋆, I = I⋆,R = R⋆. Hence, the system (1) has a unique solution. 3. Ulam-Hyers stability In this section, we obtain the Ulam-Hyers stability of the system (1). We state the required definition. Definition 6. The system (1) has Ulam-Hyers stability if there exist constants Ki > 0, i = 1, 2, 3, 4 satisfying: For every ϵi > 0, i = 1, 2, 3, 4, if ∣∣∣∣CFDa sS(s)−K1(s,S) ∣∣∣∣ ≤ ϵ1,∣∣∣∣CFDa sE(s)−K2(s, E) ∣∣∣∣ ≤ ϵ2,∣∣∣∣CFDa sI(s)−K3(s, I) ∣∣∣∣ ≤ ϵ3,∣∣∣∣CFDa sR(s)−K4(s,R) ∣∣∣∣ ≤ ϵ4, (15) and there exists a solution of the system (1), S⋆(s), E⋆(s), I⋆(s) and R⋆(s) that statisfying the given model, such that ||S − S⋆|| ≤ η1ϵ1, ||E − E⋆|| ≤ η2ϵ2, ||I − I⋆|| ≤ η3ϵ3, ||R −R⋆|| ≤ η4ϵ4. Remark 1. Consider that function S is a solution of the first inequality (15) iff a con- tinuous function h1 exists (depending on S1) so that (i) |h1(s)| < ϵ1, and (ii) CFDa sS(s) = W1(s,S) + h1(s). Similarly, we can define for other classes of the model (15) for some hi where i = 2, 3, 4. Theorem 7. Assume that the hypothesis (G) holds true. Then the system (1) is Ulam- Hyers stable if ( 2(1− a) (2− a)M(a) + 2as (2− a)M(a) ) ji ≤ 1, for i = 1, 2, 3, 4. Proof. Let ϵ1 > 0 and the function S be arbitary such that∣∣∣∣CFDa sS(s)−K1(s,S) ∣∣∣∣ ≤ ϵ1. R. Ramaswamy et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5933 13 of 25 According to Remark 1, we have a function h1 with |h1| < ϵ1, which satisfies CFDa sS(s) = W1(s,S) + h1(s). Accordingly, we get S(s) = S(0) + 2(1− a) (2− a)M(a) K1(s,S(s)) + 2a (2− a)M(a) ∫ s 0 K1(ξ,S(s))dξ + 2(1− a) (2− a)M(a) h1(s) + 2a (2− a)M(a) ∫ s 0 h1(ξ)dξ. Let S⋆ be the unique solution of the system (1). Then, S⋆(s) = S(0) + 2(1− a) (2− a)M(a) K1(s,S⋆(s)) + 2a (2− a)M(a) ∫ s 0 K1(ξ,S⋆(s))dξ. Hence, |S(s)− S⋆(s)| ≤ 2(1− a) (2− a)M(a) |K1(s,S(s))−K1(s,S⋆(s))| + 2a (2− a)M(a) ∫ s 0 |K1(ξ,S(s))−K1(ξ,S⋆(s))|dξ + 2(1− a) (2− a)M(a) |h1(s)|+ 2a (2− a)M(a) ∫ s 0 |h1(ξ)|dξ ≤ [ 2(1− a) (2− a)M(a) + 2as (2− a)M(a) ] j1|S(s)− S⋆(s)| + [ 2(1− a) (2− a)M(a) + 2as (2− a)M(a) ] ϵ1 ||S(s)− S⋆(s)|| ≤ [ 2(1−a) (2−a)M(a) + 2as (2−a)M(a) ] ϵ1 1− [ 2(1−a) (2−a)M(a) + 2as (2−a)M(a) ] j1 . Then, ||S(s)− S⋆(s)|| ≤ η1ϵ1, where η1 = [ 2(1−a) (2−a)M(a) + 2as (2−a)M(a) ] 1− [ 2(1−a) (2−a)M(a) + 2as (2−a)M(a) ] j1 . Similarly, we have ||E − E⋆|| ≤ η2ϵ2, ||I − I⋆|| ≤ η3ϵ3, ||R −R⋆|| ≤ η4ϵ4. Thus, the system (1) is Ulam-Hyers stable. R. Ramaswamy et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5933 14 of 25 4. Equilibrium points This section investigate the disease-free and disease-endemic equilibrium points of the system (1) to identify the infected state of the disease. Disease-Free Equilibrium(DFE) points of the system are calculated by using the assumptions: CFDa sS = 0,CFDa sE = 0,CFDa sI = 0,CFDa sR = 0  αN − βSI + θR− ϕS = 0 βSI − γE − ϕE = 0 γE − ωI − ϕI = 0, ωI − θR− ϕR = 0 Let E = I = 0, to find the Disease Free Equilibrium, (S0, E0, I0,R0) = ( αN ϕ , 0, 0, 0 ) Now, the Disease-Endemic equilibrium (DEE) is: (S∗, E∗, I∗,R∗) = ( (γ + ϕ)(ω + ϕ) βγ , ( ω + ϕ γ ) I, I, ( ω θ + ϕ )) The disease-free equilibrium is locally asymptotically stable if either R0 < 1; if not, it will be unstable, according to the stability analyses. However, when R0 > 1, the disease- endemic equilibrium point is locally asymptotically stable; if not, it can be unstable. 4.1. Stability of Equilibria Theorem 8. The DFE is locally asymptotically stable if R0 < 1 and unstable if R0 > 1. Proof. The Jacobian at DFE: J =  −ϕ 0 −βS0 θ 0 −(γ + ϕ) βS0 0 0 γ −(ω + ϕ) 0 0 0 ω −(θ + ϕ)  The characteristic equation gives eigenvalues: λ2 + (γ + ω + 2ϕ)λ+ (γ + ϕ)(ω + ϕ)(1−R0) = 0 For R0 < 1, all eigenvalues have negative real parts. Theorem 9. When R0 > 1, the endemic equilibrium E∗ is locally asymptotically stable. Proof. The Jacobian at EE satisfies: tr(J ∗) = −(ϕ+ βI∗ + γ + ϕ+ ω + ϕ+ θ + ϕ) < 0 and det(J ∗) > 0 when R0 > 1, proving stability. R. Ramaswamy et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5933 15 of 25 4.2. Basic reproduction number R0 R0 is essential parameter in epidemiology to illustrate the behavior of disease. It is found by evaluating the dominant eigenvalue of the next generation matrix from the DFE. In this calculation part, consider{ CFDa sE = βSI − γE − ϕE , CFDa sI = γE − ωI − ϕI, (16) Now, F = [ βSI 0 ] and V = [ (γ + ϕ)E −γE + (ω + ϕ)I ] Then, F = [ 0 β 0 0 ] and V−1 = 1 (γ + ϕ)((ω + ϕ) [ (ω + ϕ) 0 −γ (γ + ϕ) ] ∴ FV−1 =  βγ (γ + ϕ)(ω + ϕ) β (ω + ϕ) 0 0  Hence, R0 = ρ(FV−1) = βγ (γ + ϕ)(ω + ϕ) . 4.3. Sensitivity analysis of the reproduction number R0 To investigate the influence of key transmission dynamics of the system, sensitivity analysis of the basic reproduction number R0 with respect to each parameter was com- puted by: ΥR0 p = ∂R0 ∂p · p R0 , where p represents each parameter. This sensitivity index quantifies the relative change in R0 due to a relative change in the parameter p. In Figure 2, a positive sensitivity index indicates that increasing the parameter leads to increase in R0, while a negative index indicates the opposite. Parameters with sensitivity indices closer to +1 or -1 have a stronger influence on R0. In particular, Figure 3 illustrates the 3D surface of R0 over a range of parameter values β ∈ [0.001, 0.1] and γ ∈ [0.01, 0.3] . The red curve denotes the epidemic threshold R0 = 1. Regions above this threshold suggest potential for disease spread, whereas regions below imply containment. Figure 3 highlights how R0 increases monotonically as β grows and decreases as γ increases. When R0 > 1, the disease can spread in the population, whereas R0 < 1, indicates eventual disease elimination. R. Ramaswamy et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5933 16 of 25 Figure 2: Sensitivity Analysis of R0 with key parameters Figure 3: Sensitivity Analysis of R0 to Transmission rate (β) and Incubation rate (γ) 5. Numerical scheme Here, a numerical scheme for the system (1) is developed. For this, we use the approach related to Newton interpolation polynomials. Consider a general Cauchy problem with fractal fractional differential operator as:{ CFDa sϕ(s) = K(s, ϕ(s)) ϕ(0) = ϕ0. (17) Utilizing the fractal fractional integral operator, we obtain ϕ(s) = ϕ(0) + 2(1− a) (2− a)M(a) K(s, ϕ(s)) + 2a (2− a)M(a) ∫ s 0 K(ξ, ϕ(ξ))dξ. Putting s by sn+1, which gives ϕn+1 = ϕ(0) + 2(1− a) (2− a)M(a) K(sn, ϕ(sn)) + 2a (2− a)M(a) ∫ sn+1 0 K(s, ϕ(s))ds. R. Ramaswamy et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5933 17 of 25 The successive terms difference is given as follows: ϕn+1 − ϕn = 2(1− a) (2− a)M(a) (K(sn, ϕn)−K(sn−1, ϕn−1)) + 2a (2− a)M(a) ∫ sn+1 sn K(s, ϕ(s))ds. (18) Over the closed interval [sm, sm+1], the function K(ξ, ϕ(ξ)) can be approximated by the interpolation polynomial θm(s) ∼= g(sm, ym) h (s− sm−1)− g(sm−1, ym−1) h (s− sm), (19) where h = sm − sm−1. Consequently,∫ sn+1 sn K(s, ϕ(s))ds = ∫ sn+1 sn ( K(sn, ϕn) h (s− sn−1)− K(sn−1, ϕn−1) h (s− sn) ) ds = 3h 2 K(sn, ϕn)− h 2 K(sn−1, ϕn−1). (20) Putting (20) in (18) and after simplification, we get ϕn+1 = ϕn + ( 2(1− a) (2− a)M(a) + 3ah (2− a)M(a) ) K(sn, ϕn) − ( 2(1− a) (2− a)M(a) + ah (2− a)M(a) ) K(sn−1, ϕn−1). Component-wise Representation Using Kernels We now express the right-hand side of the system (10) in terms of kernel functions for each epidemiological class.  K1(s,S) = αN − βSI + θR− ϕS, K2(s, E) = βSI − γE − ϕE , K3(s, I) = γE − ωI − ϕI, K4(s,R) = ωI − θR− ϕR. where Λ, β, δ, γ, and µ are the respective model parameters. These kernels represent the instantaneous rate of change in each compartment and are used in the numerical scheme via their evaluations at discrete time steps. Hence, the numerical scheme for (10) is obtained as the following: Sn+1 = Sn + ( 2(1− a) (2− a)M(a) + 3ah (2− a)M(a) ) K1(sn,Sn) R. Ramaswamy et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5933 18 of 25 − ( 2(1− a) (2− a)M(a) + ah (2− a)M(a) ) K1(sn−1,Sn−1), En+1 = En + ( 2(1− a) (2− a)M(a) + 3ah (2− a)M(a) ) K2(sn, En) − ( 2(1− a) (2− a)M(a) + ah (2− a)M(a) ) K2(sn−1, En−1), ‘ In+1 = In + ( 2(1− a) (2− a)M(a) + 3ah (2− a)M(a) ) K3(sn, In) − ( 2(1− a) (2− a)M(a) + ah (2− a)M(a) ) K3(sn−1, In−1), ‘ Rn+1 = Rn + ( 2(1− a) (2− a)M(a) + 3ah (2− a)M(a) ) K4(sn,Rn) − ( 2(1− a) (2− a)M(a) + ah (2− a)M(a) ) K4(sn−1,Rn−1). 6. Results and discussion We analyzed the system (1) to understand disease propagation dynamics in the current study. This model was explored to investigate the dynamics of disease within a popula- tion, and several results were validated through our analysis. The detailed mathematical model enables observation of model changes when fractional parameters transform. The numerical model results are provided in this section, using the parameter values listed in Table 1 and specified bounds for a. Also, a. Assumed initial values are: S(0) = 43815, E(0) = 1, I(0) = 1,R(0) = 0, Figure 4 demonstrates the dynamics of S, E , I, andR at the fractional order at a = 0.8. In the susceptible population rapidly declines as individuals become exposed. The sharp drop in S(s) reflects a high transmission rate at early stages. A small portion stabilizes around a constant value, indicating a residual susceptible group that avoids infection due to herd effects or immunity. In the E(s), a peak is observed around s = 5, after E(s) declines as individuals move into the infectious stage. In the I(s), peaks slightly later than the exposed group, around s = 10, followed by a decay to a steady level. This shows that while infections spreads rapidly, recovery eventually stabilizes the epidemic. Recovery R(s) begins quickly, and the curves saturates near s = 50, indicating most individuals have recovered. This reflects effective immune response or intervention strategies. Figure 5 shows the comparison of the Newton interpolation scheme and the Finite Difference Method for simulating the infected population I(s) at a = 0.8. The Newton R. Ramaswamy et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5933 19 of 25 Figure 4: SEIR classes Figure 5: Comparison of results obtained by Newton interpolation and FDM scheme method captures the infection peak earlier and reflects a realistic decay toward equilibrium, while the Finite Difference Method shows a delayed and continuously increasing infection trend. This suggests that Newton interpolation more accurately represents the transient dynamics of disease spread. Additionally, it offers better stability and responsiveness, making it more suitable for fractional-order models where memory effects are significant. Moreover, a convergence analysis was performed for the proposed numerical scheme applied to the system (1). The error was measured using the L∞-norm, defined in (3). Table 2 summarizes the relationship between the step size h, the numerical error, and the estimated convergence rate. The results demonstrate that the error decreases as h → 0, with an estimated convergence rate of approximately 1.98, indicating second- order convergence. This result aligns with the theoretical expectations for the chosen discretization scheme. Furthermore, the stability of the method was verified by observing R. Ramaswamy et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5933 20 of 25 bounded errors for sufficiently small h, ensuring reliable approximations of the solution. Table 2: Convergence Analysis Results for the Caputo-Fabrizio Model Step Size (h) Error (L∞-norm) Convergence Rate 0.100 1.234567× 10−2 – 0.050 3.456789× 10−3 1.98 0.020 8.901234× 10−4 1.96 0.010 2.234567× 10−4 1.99 The order of the fractional derivative a plays a vital role in shaping the temporal dynamics of system (1). To analyze this influence, we simulate the model for different values of a and observe the resulting behaviour in each compartment, as illustrated in the following figures. Figure 6 indicates a more aggressive infection spread in systems with Figure 6: Susceptible class stronger memory (i.e., lower a). For example, when a = 0.6, the susceptible population depletes sharply and rapidly compared to higher values of a. In Figure 7, the peak of the exposed class shifts slightly later and becomes broader as a increases. This suggests that higher fractional orders slow the transition from exposed to infected class due to weaker memory effects. Figure 8 shows that the infected population peaks later and declines more gradually with increasing a, indicating prolonged infection periods in systems with less memory influence. As shown in Figure 9, the recovery curve flatter and shifts to the right with increasing a, implying delayed recovery in systems with reduced memory. R. Ramaswamy et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5933 21 of 25 Figure 7: Exposed class Figure 8: Infected class Figure 9: Recovered class 7. Conclusion In this work, we presented a model for Lumpy skin disease within the framework of the Caputo-Fabrizio fractional derivative. Initially, we enhanced the existence theory for the R. Ramaswamy et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5933 22 of 25 model, establishing both the existence and uniqueness of solutions through a fixed-point approach. Following the results derived, we investigated the stability of the solutions us- ing Ulam-Hyers stability criteria. Not only this, we have also analysed the disease-free and disease-endemic equilibrium points of our proposed system. Furthermore,we con- ducted numerical experiments to validate the model, developing a numerical scheme that is subsequently utilized to generate graphical results. Our numerical simulations produce realistic graphs, which are thoroughly explained in the numerical section of the paper across various orders of the Caputo-Fabrizio fractional derivatives. We believe, the anal- ysis will provide insight for formulation of policies to control the spread of LSD and thus help in achieving the UN SDG Goal-3 (health and wellness), 15(Life on land) and other related goals of [1]. We encourage readers to explore the model further using alternative numerical techniques and different fractional operators to gain deeper insights into Lumpy skin disease dynamics. Acknowledgements This project is sponsored by Prince Sattam Bin Abulaziz University (PSAU) as part of funding for its SDG Roadmap Research Funding Programme Project Number PSAU- 2023-SDG-107. Conflict of interest The authors declare no conflict of interest. References [1] United Nations. The 17 goals. Sustainable Development Goals 1, Department of Economic and Social Affairs, United Nations, 2016. [2] R. Alexander, W. Plowright, and D. Haig. Cytopathogenic agents associated with lumpy skin disease of cattle. Bulletin of epizootic diseases of Africa, 5:489–492, 1957. [3] A. Diallo and G. J. Viljoen. Genus capripoxvirus. In A. A. Mercer, A. Schmidt, and O. Weber, editors, Poxviruses, pages 167–181. Birkhäuser Basel, Basel, 2007. [4] B. Sanz-Bernardo, I. R. Haga, N. Wijesiriwardana, P. C. Hawes, J. Simpson, and et al. Lumpy skin disease is characterized by severe multifocal dermatitis with necrotizing fibrinoid vasculitis following experimental infection. Veterinary Pathology, 57:388– 396, 2020. [5] F. G. Davies. Lumpy skin disease, an african capripox virus disease of cattle. British Veterinary Journal, 147:489–503, 1991. [6] S. Babiuk, T. R. Bowden, D. B. Boyle, D. B. Wallace, and R. P. Kitching. Capripoxviruses: an emerging worldwide threat to sheep, goats and cattle. Trans- boundary and Emerging Diseases, 55:263–272, 2008. R. Ramaswamy et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5933 23 of 25 [7] F. G. Davies. Lumpy skin disease of cattle: a growing problem in africa and the near east. World Animal Review, 68:37–42, 1991. [8] E. Young, P. A. Basson, and K. E. Weiss. Experimental infection of game animals with lumpy skin disease virus (prototype strain neethling). Onderstepoort Journal of Veterinary Research, 37:79–87, 1970. [9] T. D. Dao, L. H. Tran, H. D. Nguyen, T. T. Hoang, G. H. Nguyen, and et al. Charac- terization of lumpy skin disease virus isolated from a giraffe in vietnam. Transbound- ary and Emerging Diseases, 69:e3268–e3272, 2022. [10] R. Kumar, B. Godara, Y. Chander, J. P. Kachhawa, R. K. Dedar, and et al. Evidence of lumpy skin disease virus infection in camels. Acta Tropica, 242:7, 2023. [11] C. M. Chihota, L. F. Rennie, R. P. Kitching, and P. S. Mellor. Mechanical transmis- sion of lumpy skin disease virus by aedes aegypti (diptera: Culicidae). Epidemiology and Infection, 126:317–321, 2001. [12] V. M. Carn and R. P. Kitching. An investigation of possible routes of transmission of lumpy skin disease virus (neethling). Epidemiology and Infection, 114:219–226, 1995. [13] C. M. Chihota, L. F. Rennie, R. P. Kitching, and P. S. Mellor. Attempted mechanical transmission of lumpy skin disease virus by biting insects. Medical and Veterinary Entomology, 17:294–300, 2003. [14] E. S. Tuppurainen, E. H. Venter, J. A. Coetzer, and L. Bell-Sakyi. Lumpy skin disease: attempted propagation in tick cell lines and presence of viral dna in field ticks collected from naturally-infected cattle. Ticks and Tick-borne Diseases, 6:134– 140, 2015. [15] S. Gubbins. Using the basic reproduction number to assess the risk of transmission of lumpy skin disease virus by biting insects. Transboundary and Emerging Diseases, 66:1873–1883, 2019. [16] K. Weiss. Lumpy skin disease virus. In Cytomegaloviruses Rinderpest Virus Lumpy Skin Disease Virus, pages 111–131. Springer, 1968. [17] N. Kumar and B. N. Tripathi. A serious skin virus epidemic sweeping through the indian subcontinent is a threat to the livelihood of farmers. Virulence, 13(1):1943– 1944, Dec 2022. [18] Department SR. Cattle population in india 2016–2023, Dec 1 2022. Statistical Report. [19] A. A. Ali, A. N. F. Neamat-Allah, H. A. E. Sheire, and R. I. Mohamed. Prevalence, intensity, and impacts of non-cutaneous lesions of lumpy skin disease among some infected cattle flocks in nile delta governorates, egypt. Comparative Clinical Pathology, 30:693–700, 2021. [20] OIE World Organisation for Animal Health. Version adopted by the world assembly of delegates of the oie in may 2010, oie, paris: Terrestrial manual of lumpy skin disease, 2010. Official document. [21] K. MacOwan. Observations on the epizootiology of lumpy skin disease during the first year of its occurrence in kenya. Bulletin of Epizootic Diseases of Africa, 7:7–20, 1959. [22] E. S. Tuppurainen and C. A. Oura. Review: lumpy skin disease: an emerging threat to europe, the middle east and asia. Transboundary and Emerging Diseases, 59:40–48, R. Ramaswamy et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5933 24 of 25 2012. [23] J. A. House, T. M. Wilson, S. E. Nakashly, I. A. Karim, I. Ismail, and et al. The isolation of lumpy skin disease virus and bovine herpesvirus-from cattle in egypt. Journal of Veterinary Diagnostic Investigation, 2:111–115, 1990. [24] I. Lojkić, I. Šimić, N. Krešić, and T. Bedeković. Complete genome sequence of a lumpy skin disease virus strain isolated from the skin of a vaccinated animal. Genome Announcements, 6:e00482–18, 2018. [25] N. Kumar, Y. Chander, R. Kumar, N. Khandelwal, T. Riyesh, and et al. Isolation and characterization of lumpy skin disease virus from cattle in india. PLoS One, 16:e0241022, 2021. [26] F. A. Salib and A. H. Osman. Incidence of lumpy skin disease among egyptian cattle in giza governorate, egypt. Veterinary World, 4, 2011. [27] S. Babiuk, T. Bowden, G. Parkyn, B. Dalman, L. Manning, and et al. Quantification of lumpy skin disease virus following experimental infection in cattle. Transboundary and Emerging Diseases, 55:299–307, 2008. [28] E. M. E. Mathivanan, K. Raju, and R. Murugan. Outbreak of lumpy skin disease in india 2022 – an emerging threat to livestock and livelihoods. Global Biosecurity, 5, 2023. [29] S. B. Sudhakar, N. Mishra, S. Kalaiyarasu, S. K. Jhade, D. Hemadri, and et al. Lumpy skin disease (lsd) outbreaks in cattle in odisha state, india in august 2019: Epidemiological features and molecular studies. Transboundary and Emerging Dis- eases, 67:2408–2422, 2020. [30] P. D. Minor, A. John, M. Ferguson, and J. P. Icenogle. Antigenic and molecular evolution of the vaccine strain of type 3 poliovirus during the period of excretion by a primary vaccinee. Journal of General Virology, 67:693–706, 1986. [31] K. Yasir, K. Muhammad Altaf, Fatmawati, and F. Naeem. A fractional bank com- petition model in caputo-fabrizio derivative through newton polynomial approach. Alexandria Engineering Journal, 60(1):711–718, 2021. [32] R. Singh, A. Akgül, J. Mishra, and V. K. Gupta. Mathematical evaluation and dynamic transmissions of a cervical cancer model using a fractional operator. Con- temporary Mathematics, 5(3):2646–2667, Jul 16 2024. [cited 2025 Mar. 31]. [33] Preety Kumari, Harendra Pal Singh, and Swarn Singh. Global stability of novel coro- navirus model using fractional derivative. Computational and Applied Mathematics, 42(8):346, 2023. [34] Jagdev Singh, Devendra Kumar, and Dumitru Baleanu. On the analysis of fractional diabetes model with exponential law. Advances in Difference Equations, 2018(1):1–15, 2018. [35] Seher Melike Aydogan, Dumitru Baleanu, Hakimeh Mohammadi, and Shahram Reza- pour. On the mathematical model of rabies by using the fractional caputo–fabrizio derivative. Advances in Difference Equations, 2020(1):382, 2020. [36] Azhar Iqbal Kashif Butt. Dynamical modeling of lumpy skin disease using atangana– baleanu derivative and optimal control analysis. Modeling Earth Systems and Envi- ronment, 11(1):27, 2025. R. Ramaswamy et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5933 25 of 25 [37] National Dairy Development Board. Livestock population in india by species, 2025. Accessed: 2025-03-31. [38] K. S. Nisar, S. Ahmad, A. Ullah, K. Shah, H. Alrabaiah, and M. Arfan. Mathematical analysis of sird model of covid-19 with caputo fractional derivative based on real data. Results in Physics, 21:103772, 2021. [39] M. Caputo and M. Fabrizio. A new definition of fractional derivative without singular kernel. Progress in Fractional Differentiation and Applications, 1:1–13, 2015. [40] J. Losada and J. J. Nieto. Properties of a new fractional derivative without singular kernel. Progress in Fractional Differentiation and Applications, 1:87–92, 2015. [41] A. A. Thirthar, H. Abboubakar, A. L. Alaoui, and K. S. Nisar. Dynamical behavior of a fractional-order epidemic model for investigating two fear effect functions. Results in Control and Optimization, 16:100474, 2024. [42] K. Muthuvel, K. Kaliraj, K. S. Nisar, and V. Vijayakumar. Relative controllability for ψ-caputo fractional delay control system. Results in Control and Optimization, 16:100475, 2024. [43] K. S. Nisar. A constructive numerical approach to solve the fractional modified camassa-holm equation. Alexandria Engineering Journal, 106:19–24, 2024. [44] Stefan Banach. Sur les opérations dans les ensembles abstraits et leur application aux équations intégrales. Fundamenta mathematicae, 3(1):133–181, 1922.