Computational Framework for Fractional Order Neurological Disorder Model under Interpreting Transmission Patterns

1  Introduction

The complex mechanisms behind neurodegenerative disorders such as Parkinson’s disease, Alzheimer’s disease, and multiple sclerosis make them major global health concerns. These conditions entail complex interactions between brain networks, immune cells, and abnormal protein aggregation. Conventional mathematical methods provide useful information, but they cannot capture all aspects of neurodegenerative progression, particularly memory-dependent processes. Recent mathematical modeling attempts have revealed insights into how neurodegenerative diseases progress. Zehra et al. [1] examined the physiological and chaotic effects on neurological disorders using fractional operators. Although the precise cause of neurodegenerative diseases remains incompletely understood, evidence supports the role of protein aggregation and neuroinflammatory responses [2–6]. Mathematical models have been developed to elucidate the concentric organization of demyelinating lesions and highlight the significance of macrophage activation and mobilization [7–11]. However, these models either rely on integer-order derivatives or do not fully incorporate the coupled dynamics of neurons, microglia, and T-cells under a fractional framework.

Fractional-order derivatives overcome the challenge of capturing time-based linkages in neurodegenerative investigations, such as protein clumping events and synapse loss by replacing classical calculus with operators that incorporate past state conditions [12,13]. Fractional-order models have repeatedly shown favorable outcomes in neurobiology research. Recent fractional-order concepts have been applied to the Schrödinger equation [14], biophysics [15], nonlinear equations [16], fractional-order wave equations [17], time-fractional equations [18], and fluid dynamics [19]. An extensive literature survey indicates that computationally intelligent-based solutions have become a key focus in the current phase of technological advancements. These recommendations highlight the critical role of computational solvers and urge the authors to create a reliable, accurate, and consistent approach for addressing the model. Fractional calculus has many applications across various fields, most commonly in engineering and physics. Among modeling options, fraction order systems are more factual and empirical as they capture the difference between genetic and memory features of mathematical frameworks which is different from classical integer order models [20]. Fractional-order models improve mathematical modeling for complicated problems by successfully addressing dynamical processes under uncertainty (See, [21–23]). Fractional-order models for infections and diseases have more stages of freedom than regular derivatives, according to recent studies like [24,25], suggesting that fractional-order derivations may offer a more accurate representation of physiological processes than classical order models. A hybrid strategy that combines machine learning (ML) and fractional-order dynamical modeling was developed in a recent study [26] to predict Parkinson’s disease (PD) using vocal biomarkers. The interpretability and predictive accuracy of early, non-invasive Parkinson’s disease detection are enhanced by this integration. The framework advances computational neurology by offering chances for improved diagnostic tools and customized monitoring.

This work offers a novel approach to comprehending neurodegeneration by capturing the influence of memory on neuroinflammation using a new fractional differential equation with Mittag-Leffler Kernel. The stability of the endemic and disease-free equilibria, as well as the use of machine learning (NARX-BRBNN) to increase numerical simulation efficiency, are also provided by the analysis. The goal of this study is to use the idea of the fractal-fractional derivative to apply complex non-linear differential equations in order to create a new model for brain pathology. Because fractional order and fractional fractal Mittag-Leffler derivatives can more correctly represent memory strength than earlier models (e.g., [27–30]), they are used to investigate the behavior of brain diseases.

The study uses a methodical approach to present the research findings, starting with an introduction to its suggested format. In Section 2, it develops a revolutionary system for brain pathology, describing specific parameters and basic ideas of the fractal fractional operator. Section 3 qualitatively examines the biological viability of the system of equations that serves as the foundation for fractional calculus. Main analysis includes equilibrium points, well-posedness of the framework, basic reproductive number, sensitivity analysis, the positively invariant set (Ω) of the system, and existence and uniqueness of the model. We discuss theorems for the local stability of equilibria and the corresponding proofs in Section 4. We also utilize the correct Lyapunov function to establish global asymptotic stability. Section 5 has the numerical solution of the fractal fractional brain disease model with a Mittag-Leffler kernel and conclusions in the final Section 6.

2  Formulation of Fractional-Order Brain Model

Over time t, the model divides the whole community N(t) into five different groups. To treat severe cases of brain damage, a fractal fractional order system has been developed by finding the most effective individuals. Each of the five compartments that make up the brain disease model represents unique cellular and extracellular components that are essential.

•   F(t): The functioning neurons’ density;

•   I(t): The infected neurons;

•   Sα(t): Extracellular α-synuclein function;

•   M(t): The number of triggered microglia; and

•   T(t): Active T-cells.

Therefore, the entire population N(t) is built by

N(t)=F(t)+I(t)+Sα(t)+M(t)+T(t).(1)

Important parameters, used in model formulation, are summarized in Table 1.

These parameter values (Table 1) have been assumed using scientific literature [31–33] and biological reasoning because several of these neuroinflammation processes have not been empirically explored. Key biological parameter values, such as neuron generative rate, natural death rate, and infection process, have been established using existing neurodegenerative disease models. The infection and release rates of proteins, such as α-synuclein, were designed to replicate the aggregation rate seen in Parkinson’s and Alzheimer’s disease. Furthermore, immune response characteristics such as the activation of microglia cells and T-cells have been chosen to respond at optimal periods in relation to neural infections.

A more accurate way to represent complicated systems with memory effects and non-integer dimensions is to use fractal fractional derivatives. They are becoming more and more significant in engineering and scientific fields. In this work, a unique framework for comprehending the intricate interactions between immunological response mechanisms and neuronal health in neurodegenerative illnesses is presented. The model represents biological processes such as infection rates, immunological activity, and neural degradation using nonlinear fractional differential equations. Delays in reactions and the influence of prior conditions on the progression of disease are also taken into account. Using the Mittag-Leffler definition, the model below, which is based on the generalized hypothesis mentioned above and the flow chart in Fig. 1, presents the impact as follows

{D0,tη,αFFMF(t)=ΠN−FγSα−μNF,D0,tη,αFFMI(t)=F(t)γSα−α1I−I(t)β1M,D0,tη,αFFMSα(t)=δ1I−dαSα−SαλM,D0,tη,αFFMM(t)=Θ+ψSα−dMM,D0,tη,αFFMT(t)=ψTM+ΠT−μTT,(2)

where the initial conditions are as follows:

F(0)≥0,I(0)≥0,Sα(0)≥0,M(0)≥0,T(0)≥0.(3)

Figure 1: Flow diagram showing the various phases of the brain diseases model.

All of them are biologically feasible.

Some Basic Concepts

We provide some fundamental ideas from the Mittag-Leffler fractional calculus to serve as a basis for the findings in this work.

Definition 1 (Mittag-Leffler Derivative [34]): Suppose that f(t) is a function that may or may not be differentiable. Let 0<η,α≤1, where α is a fractal dimension and an is a fractional order of η. A fractal-fractional derivative of f(t) in context of a Mittag-Leffler kernel can be described as

D0,tη,αFFMf(t)=AB(η)1−ηddtα∫0tf(ξ)Eη{−η1−η(t−ξ)η}dξ.(4)

Definition 2 (Mittag-Leffler Integral [34]): The corresponding integral can be described as

J0,tη,αFFMf(t)=1−ηAB(η)t1−αf(t)+ηAB(η)Γ(η)∫0tf(ξ)(t−ξ)η−1ξ1−αdξ.(5)

3  Qualitative Evaluation of the System’s Biological Viability

By examining the nonlinear differential equations, this part qualitatively examines the system. This analysis focuses on identifying key features of the model and examining its most influential parameters.

3.1 Equilibrium Points

We begin with an analysis of the equilibrium points. To obtain these points, we set the left-hand side of system (2) to zero.

3.1.1 Disease-Free Equilibrium (DFE)

At the disease-free equilibrium, there is no infection (I=0) and no extracellular α-synuclein (Sα=0). Solving the steady-state equations:

0=ΠN−μNF⇒F∗=ΠNμN,(6)

0=−α1I−β1IM⇒I∗=0,(7)

0=δ1I−dαSα−λSαM⇒Sα∗=0,(8)

0=Θ+ψSα−dMM⇒M∗=ΘdM,(9)

0=ΠT+ψTM−μTT⇒T∗=ΠTdM+ψTΘdMμT.(10)

Thus, the disease-free equilibrium (DFE) is

E0=(F0,I0,Sα0,M0,T0)=(ΠNμN, 0, 0, ΘdM, ΠTdM+ψTΘdMμT).(11)

3.1.2 Endemic Equilibrium

The endemic equilibrium E∗=(F∗,I∗,Sα∗,M∗,T∗) with I∗>0 and Sα∗>0 exists when R0>1. Due to the nonlinear coupling between compartments, the endemic equilibrium is obtained by solving the steady-state equations numerically. The explicit closed-form expression is algebraically complex and is therefore omitted here; numerical values are provided in the simulation section.

3.2 Well Posedness of the Framework

In this part, we analyze the boundedness of solutions. Summing all equations in system (2):

D0,tη,αFFMN(t)=ΠN−γFSα−μNF+γFSα−α1I−β1IM+δ1I−dαSα−λSαM+Θ+ψSα−dMM+ΠT+ψTM−μTT.(12)

D0,tη,αFFMN(t)=ΠN−μNF−α1I−β1IM+δ1I−dαSα−λSαM+Θ+ψSα−dMM+ΠT+ψTM−μTT.

Since all variables are nonnegative, we obtain

D0,tη,αFFMN(t)≤ΠN+Θ+ΠT−μmin{F+I+Sα+M+T},

where μmin=min{μN,α1,dα,dM,μT}. Hence, N(t) is bounded and the biologically feasible region is

Ω={(F(t),I(t),Sα(t),M(t),T(t))∈R+5:F+I+Sα+M+T≤ΠN+Θ+ΠTμmin}.(13)

The region Ω is positively invariant, i.e., any solution starting in Ω remains in Ω for all t≥0.

3.3 The Fundamental Reproductive Number

The basic reproduction number R0 is calculated using the next-generation matrix approach [35]. We focus on infected compartments I and Sα (because F, M, and T are not directly implicated in new infections). In the DFE, the Jacobian matrices for new infections (F) and transfer terms (V) are:

F=[0γΠNμN00],V=[α1+β1ΘdM0−δ1dα+λΘdM].(14)

Then, we have

V−1=1(α1dM+β1Θ)(dαdM+λΘ)[dM(dαdM+λΘ)0δ1dMα1dM+β1Θ].

The next-generation matrix is K=FV−1, and R0 is its spectral radius calculates as:

R0=δ1dM2γΠNμN(α1dM+β1Θ)(dαdM+λΘ).(15)

Surfaces representing sensitivity with respect to different parameter pairs are shown in Fig. 2. In general, R0 is observed to be increased by infection rate (γ), secretion rate (δ1), and neuron production rate (ΠN) and decreased by the following clearance or deactivation rates or parameters: microglia clearance parameter (λ), secretion parameter (β1), degradation parameter (dα), and deactivation rate (dM). The surfaces of the pairs (γ,λ) and (γ,dM) represent crucial trade-off relations that imply that the proportionally high clearance or deactivation rate is needed to prevent R0>1. However, the pair (γ,δ1) displays the highest joint impact on the increase of disease severity. Therefore, it would be reasonable to combine the treatment of the infection and secretion in order to obtain the best results.

Figure 2: Examination of the reproductive number based on specific criteria.

3.4 Sensitivity of R0’s Parameters

Sensitivity analysis assesses how changes in parameters affect R0, helping prioritize therapeutic targets. The normalized forward sensitivity index of R0 with respect to a parameter p is defined as

ΓpR0=∂R0∂p⋅pR0.(16)

We have

R0=δ1dM2γΠNμN(α1dM+β1Θ)(dαdM+λΘ).(17)

∂R0∂δ1=dM2γΠNμN(α1dM+β1Θ)(dαdM+λΘ)>0,

∂R0∂γ=δ1dM2ΠNμN(α1dM+β1Θ)(dαdM+λΘ)>0,

∂R0∂ΠN=δ1dM2γμN(α1dM+β1Θ)(dαdM+λΘ)>0,

∂R0∂dM=2δ1dMγΠNμN(α1dM+β1Θ)(dαdM+λΘ)−δ1dM2γΠN[α1(dαdM+λΘ)+dα(α1dM+β1Θ)]μN(α1dM+β1Θ)2(dαdM+λΘ)2>0,

∂R0∂μN=−δ1dM2γΠNμN2(α1dM+β1Θ)(dαdM+λΘ)<0,

∂R0∂α1=−δ1dM3γΠNμN(α1dM+β1Θ)2(dαdM+λΘ)<0,

∂R0∂β1=−Θδ1dM2γΠNμN(α1dM+β1Θ)2(dαdM+λΘ)<0,

∂R0∂Θ=−δ1dM2γΠN{(α1dM+β1Θ)λ+β1(dαdM+λΘ)}μN(α1dM+β1Θ)2(dαdM+λΘ)2<0,

∂R0∂dα=−δ1dM3γΠNμN(α1dM+β1Θ)(dαdM+λΘ)2<0,

∂R0∂λ=−Θδ1dM2γΠNμN(α1dM+β1Θ)(dαdM+λΘ)2<0.

Normalized sensitivity indices are computed in Table 2, evaluated at baseline parameter values.

These normalized sensitivity indices are fuhrer displayed in Fig. 3.

Figure 3: Sensitivity analysis of model parameters (normalized indices).

The sensitivity analysis indicates that the three parameters having a major influence on R0 are ΠN, γ, and δ1 with index values of +0.9981 each, implying that a 1% increment in their values causes a corresponding increase in the value of R0 by roughly 1%. Similarly, an insignificant positive influence is indicated on R0 by dM with an index of +0.1578. In contrast to the above parameters, the influence on R0 due to the parameter μN is highly significant but in the reverse direction, having an index value of −0.9981. The other parameters having moderate influence on R0 include α1 (index value −0.4178), β1 (−0.3956), dα (−0.3652), and λ (−0.3457). The first parameter for the activation rate of microglia, Θ, exhibits the least influence (−0.1262). Based on this information, the therapy will focus on δ1, γ, or μN.

3.5 Positively Invariant Ω Set of the System

Theorem 1: A characteristic of the domain Ω is that it is +ve invariant concerning the brain disease system (2).

Proof: The solutions to the brain disease framework (2) remain physiologically viable (that is, non-negative) under suitable starting conditions, remaining inside the positive orthant of R+5 for all t≥0. We get

{D0,tη,αFFMF(t)|F(t)=0=ΠN≥0,D0,tη,αFFMI(t)|I(t)=0=γF(t)Sα(t)≥0,D0,tη,αFFMSα(t)|Sα(t)=0=δ1I(t)≥0,D0,tη,αFFMM(t)|M(t)=0=Θ+ψSα(t)≥0,D0,tη,αFFMT(t)|T(t)=0=ΠT+ψTM(t)≥0.(18)

The solution to the aforementioned systems will be obtained by using the fractional integral. The result will be nonnegative since there are no negative terms in the system. □

3.6 Existence and Distinctiveness of the Model

The existence and uniqueness of solutions to the fractional-order system (2) follow from standard results in fractional calculus. First, we note that all solutions are bounded due to the following argument:

D0,tη,αFFMN(t)≤ΠN+Θ+ΠT−μminN(t),

where μmin=min{μN,α1,dα,dM,μT}. By the fractional comparison principle [36], this inequality implies N(t)≤max{N(0),(ΠN+Θ+ΠT)/μmin} for all t≥0. Hence, each state variable is bounded, and the sup norms ‖F‖∞, ‖I‖∞, ‖Sα‖∞, ‖M‖∞, ‖T‖∞ are finite. With boundedness established, the right-hand side functions Ui are Lipschitz continuous on the compact domain Ω, and the existence and uniqueness of solutions follow from standard results for fractional differential equations [37]. For the existence and uniqueness of solutions to the system (2), we prove the following theorem.

Theorem 2: Assume that χ¯i and χi are positive constants for the purpose of

A1            |Ui(ui,t)−Ui(ui′,t)|≤χi|ui−ui′|,∀i∈1,2,⋯,5.(19)

A2            |Ui(ui,t)|≤χ¯i{1+|ui|},∀(u,t)∈R3×[0,T].(20)

Proof: We have

{D0,tη,αFFMF(t)=ΠN−γFSα−μNF=U1(t,F),D0,tη,αFFMI(t)=γFSα−α1I−β1IM=U2(t,I),D0,tη,αFFMSα(t)=δ1I−dαSα−λSαM=U3(t,Sα),D0,tη,αFFMM(t)=Θ+ψSα−dMM(t)=U4(t,M),D0,tη,αFFMT(t)=ΠT+ψTM−μTT=U5(t,T).(21)

We start with the function U1(t,F). Next, we will demonstrate that

|U1(F1,t)−U1(F2,t)|2≤χ1|F1−F2|2.(22)

Then, we write

|U1(F1,t)−U1(F2,t)|2=|−γ(F1−F2)Sα−μN(F1−F2)|2.

|U1(F1,t)−U1(F2,t)|2=|{−γSα−μN}(F1−F2)|2.

|U1(F1,t)−U1(F2,t)|2≤{2γ2|Sα|2+2μN2}|(F1−F2)2|2.

|U1(F1,t)−U1(F2,t)|2≤{2γ2sup0≤t≤T|Sα|2+2μN2}|(F1−F2)2|2.

|U1(F1,t)−U1(F2,t)|2≤{2γ2|Sα|∞2+2μN2}|(F1−F2)2|2.

|U1(F1,t)−U1(F2,t)|2≤χ1|F1−F2|2,(23)

where χ1={2γ2|Sα|∞2+2μN2}.

|U2(I1,t)−U2(I2,t)|2=|−α1(I1−I2)−β1(I1−I2)M|2.(24)

|U2(I1,t)−U2(I2,t)|2=|{−α1−β1M}(I1−I2)|2.

|U2(I1,t)−U2(I2,t)|2≤{2α12+β12|M|2}|(I1−I2)2|2.

|U2(I1,t)−U2(I2,t)|2≤{2α12+β12sup0≤t≤T|M|2}|(I1−I2)2|2.

|U2(I1,t)−U2(I2,t)|2≤{2α12+β12|M|∞2}|(I1−I2)2|2.

|U2(I1,t)−U2(I2,t)|2≤χ2|I1−I2|2,(25)

where χ2={2α12+β12|M|∞2}.

|U3(Sα1,t)−U3(Sα,t)|2=|−dα(Sα1−Sα2)−λ(Sα1−Sα2)M|2.(26)

|U3(Sα1,t)−U3(Sα,t)|2=|{−dα−λM}(Sα1−Sα2)|2.

|U3(Sα1,t)−U3(Sα,t)|2≤{2dα2+λ2|M|2}|(Sα1−Sα2)|2.

|U3(Sα1,t)−U3(Sα,t)|2≤{2dα2+λ2sup0≤t≤T|M|2}|(Sα1−Sα2)|2.

|U3(Sα1,t)−U3(Sα,t)|2≤{2dα2+λ2|M|∞2}|(Sα1−Sα2)|2.

|U3(Sα1,t)−U3(Sα,t)|2≤χ3|(Sα1−Sα2)|2,(27)

where χ3={2dα2+λ2|M|∞2}.

|U4(M1,t)−U4(M,t)|2=|−dM(M1−M2)|2.(28)

|U4(M1,t)−U4(M,t)|2≤dM2|(M1−M2)|2.

|U4(M1,t)−U4(M,t)|2≤χ4|(M1−M2)|2,(29)

where χ4=dM2.

|U5(T,t)−U5(T,t)|2=|−μT(T1−T2)|2.(30)

|U5(T,t)−U5(T,t)|2≤μT2|(T1−T2)|2.

|U5(T1,t)−U5(T,t)|2≤χ5|(T1−T2)|2,(31)

where χ4=μT2.

Every function’s starting condition is carefully investigated twice, and our model’s second condition is now been confirmed.

|U1(F,t)|2=|ΠN−γFSα−μNF|2.(32)

|U1(F,t)|2=|ΠN+{−γSα−μN}F|2.

|U1(F,t)|2≤2ΠN2+{2γ2|Sα|2+2μN2}|F|2.

|U1(F,t)|2≤2ΠN2+{2γ2sup0≤t≤T|Sα|2+2μN2}|F|2.

|U1(F,t)|2≤2ΠN2+{2γ2|Sα|∞2+2μN2}|F|2.

|U1(F,t)|2≤2ΠN2[1+{2γ2|Sα|∞2+2μN2}|F|22ΠN2].

|U1(F,t)|2≤χ¯1[1+|F|2],(33)

where χ¯1=2ΠN2 under the condition {2γ2|Sα|∞2+2μN2}2ΠN2<1.

|U2(I,t)|2=|γFSα−α1I−β1IM|2.(34)

|U2(I,t)|2=|γFSα+{−α1−β1M}I|2.

|U2(I,t)|2≤2γ2|F|2|Sα|2+{2α12+2β12|M|2}|I|2.

|U2(I,t)|2≤2γ2sup0≤t≤T|F|2sup0≤t≤T|Sα|2+{2α12+2β12sup0≤t≤T|M|2}|I|2.

|U2(I,t)|2≤2γ2|F|∞2|Sα∞|2+{2α12+2β12|M|∞2}|I|2.

|U2(I,t)|2≤2γ2|F|∞2|Sα∞|2[1+{2α12+2β12|M|∞2}|I|22γ2|F(t)|∞2|Sα∞(t)|2].

|U2(I,t)|2≤χ¯2[1+|I|2],(35)

where χ¯2=2γ2|F(t)|∞2|Sα∞|2 under the condition {2α12+2β12|M|∞2}2γ2|F|∞2|Sα∞|2<1.

|U3(Sα,t)|2=|δ1I−dαSα−λSαM|2.(36)

|U3(Sα,t)|2=|δ1I+{−dα−λM}Sα|2.

|U3(Sα,t)|2≤2δ12|I|2+{2dα2+2λ2|M|2}|Sα|2.

|U3(Sα,t)|2≤2δ12sup0≤t≤T|I|2+{2dα2+2λ2sup0≤t≤T|M|2}|Sα|2.

|U3(Sα,t)|2≤2δ12|I|∞2+{2dα2+2λ2|M|∞2}|Sα|2.

|U3(Sα,t)|2≤2δ12|I|∞2[1+{2dα2+2λ2|M|∞2}|Sα|22δ12|I|∞2].

|U3(Sα,t)|2≤χ¯3[1+|Sα|2],(37)

where χ¯3=2δ12|I|∞2 under the condition {2dα2+2λ2|M|∞2}2δ12|I|∞2<1.

|U4(M,t)|2=|Θ+ψSα−dMM|2.(38)

|U4(M,t)|2≤2Θ2+2ψ2|Sα|2+2dM2|M|2.

|U4(M,t)|2≤2Θ2+2ψ2sup0≤t≤T|Sα|2+2dM2|M|2.

|U4(M,t)|2≤2Θ2+2ψ2|Sα|∞2+2dM2|M|2.

|U4(M,t)|2≤{2Θ2+2ψ2|Sα|∞2}[1+dM2|M|22Θ2+2ψ2|Sα|∞2].

|U4(M,t)|2≤χ¯4[1+|M|2],(39)