A Physics Informed Neural Network (PINN) framework for fractional order modeling of Alzheimer's disease

1 Introduction

Fractional order modeling extends traditional integer-order differential equations by incorporating memory and nonlocal dynamics via Caputo or Riemann-Liouville derivatives. Such models integrate the entire history of a process, not just its current state, enabling accurate depiction of delayed responses, long-term persistence and power-law anomalous behaviors. In epidemiology, this approach has enhanced compartmental models like SIR and SEIR, notably in measles transmission, where Caputo-based formulations improve alignment with known incubation periods and long-tail dynamics (Angstmann et al., 2016, 2015; Akuka et al., 2024). Similarly, fractional SEIRD models for COVID-19 incorporate vaccination, quarantine (Bilgil et al., 2022) and reinfection with memory effects, yielding stronger predictive power and dynamic control insights studeid in Hussain et al. (2025).

In cancer modeling, tumor-immune system dynamics under therapeutic interventions have been enriched using fractal fractional derivatives, capturing immune memory, boundedness and stability in ways inaccessible to integer-order frameworks (Nisar et al., 2024a). Parallel advances in fractional pharmacokinetics highlight how drugs exhibiting irregular accumulation and non-exponential decay require fractional kinetics to model their distribution and clearance effectively, informing dosing strategies and toxicity avoidance (Sopasakis et al., 2019). Viscoelastic behavior in biological tissues–especially in immune cells like macrophages–also benefits from fractional modeling: the fractional Kelvin-Voigt model more accurately mirrors viscoelastic response and drug-induced cytoskeletal changes than integer counterparts, offering enhanced diagnostic and therapeutic characterization (Vo and Ekpenyong, 2022). Together, these applications motivate fractional approaches for chronic, history-dependent diseases such as Alzheimer's, where multi-scale interactions among amyloid aggregation, tau pathology, microglial activation and neuronal degradation unfold over long time horizons. A recent fractal fractional Caputo model for AD dynamics rigorously establishes existence, uniqueness and Ulam-Hyers stability and uses fractional Adams-Bashforth schemes that outperform integer-order simulations in reproducing memory-driven trajectories (Yadav et al., 2024). Some more applications and time fractional effect in sense of memory discuses for system (Zhou and Zhang, 2020) and aortic aneurysm (AAA) phenomena studied (Sumelka et al., 2020). Despite their theoretical appeal, numerical and data-driven modeling of fractional systems remains challenging due to the nonlocality of fractional operators, which limits classical solvers and complicates efficient parameter inference. Traditional discretization methods (e.g., finite differences and spectral methods) become expensive and memory-intensive when capturing nonlocal history terms over long time intervals. To address these limitations, Physics Informed Neural Networks (PINNs) embed known differential equations into neural network training by enforcing governing equations as residuals in the loss function, enabling accurate solutions even with scarce or noisy data while avoiding expensive mesh generation and grid-based methods (Raissi et al., 2019). Recent surveys highlight substantial methodological progress in PINNs such as hybrid optimization schemes, adaptive sampling techniques and multi-PDE frameworks, establishing them as versatile and effective methods for tackling both forward and inverse PDE problems across a wide range of scientific applications (Zhang et al., 2025; Farea et al., 2025a; Li, 2025; Kharazmi and Zayernouri, 2021; Murari et al., 2025; Rodrigues, 2024; Shaier et al., 2021; Farea et al., 2025b). Recent advances in fractional modeling and numerical solution techniques have demonstrated the benefits of fractional calculus for capturing memory and nonlocal effects in complex systems. For example, analytical and numerical methods for nonlinear fractional reaction-diffusion equations, such as those arising in blood flow modeling via Laplace-Residual Power Series methods have been developed to efficiently approximate solutions that would be challenging for classical approaches (Ali et al., 2025b). Fractional differential equations have also been applied to brain metabolite dynamics in circadian rhythm models where Caputo-Fabrizio derivatives and series solution techniques provide existence, uniqueness and convergence results beyond integer-order formulations (Ziada and Botros, 2025). In addition, fractional calculus has been used to model nonlinear, multi-dimensional DNA systems, highlighting how fractional models can effectively capture long-range interactions and memory effects that are absent in integer order descriptions (Ali et al., 2025a). The theoretical and numerical developments in these works support the use of fractional models and motivate the incorporation of classical fractional solvers as benchmarks in data-driven methods such as PINNs.

PINNs have been effectively applied to epidemic models such as SIR and SIRD, accurately inferring both state trajectories and time varying transmission rates from noisy outbreak data (Millevoi et al., 2023). Fractional extensions of PINNs further advance this framework by embedding fractional derivatives directly into the loss function, enabling representation of memory effects inherent in fractional PDEs. For example, PINN formulations have been used to solve time fractional Black-Scholes and related fractional diffusion equations by integrating non-integer operators into network residuals (Nuugulu et al., 2025). Other recent work explores enhanced PINN variants for β-conformable fractional differential equations, showing that specialized architecture variants like NRPINN can improve solution quality without domain discretization (Bulut and Yigider, 2025). Moreover, expanded PINN capabilities across biological and epidemiological dynamical systems illustrate the versatility of physics informed approaches when applied to ODEs and coupled systems characterized by known governing laws (Farea et al., 2025b). Nevertheless, critical limitations remain in existing PINN and fractional PINN approaches:

Automatic differentiation cannot directly compute nonlocal fractional operators, which will requires numerical discretization, auxiliary grids or transform techniques that increase complexity and cost, particularly in time fractional problems (Nuugulu et al., 2025; Bulut and Yigider, 2025).

Conventional training of PINNs can exhibit uneven optimization among the different loss components (physics residuals, observational data and boundary/initial conditions), which may cause training to stall or converge to suboptimal solutions. While recent studies introduce adaptive loss-weighting and sampling schemes to address this, more comprehensive methodological advances are still an unresolved research problem (Farea et al., 2025a) and some applications related problem is given in Lawal et al. (2022), Sopasakis and Kalliadasis (2019), Srivastava (2024), and Farea et al. (2025c).

Techniques using Monte Carlo sampling or structured sampling modules can reduce grid dependence, but often at the cost of estimator variance and sensitivity to hyperparameters, which limits robustness in large scale, high dimensional systems (Ahmad et al., 2025).

Many existing fractional PINN studies focus on specialized PDEs or prototype systems, rather than unified, scalable frameworks capable of handling multi scale biological systems with intertwined memory effects and noisy observational data, such as those found in Alzheimer's progression.

These limitations underscore the need for robust and scalable fractional PINN frameworks that can efficiently balance physics and data objectives, improve training stability and enable reliable parameter inference even under data scarcity. To bridge these gaps, foundational and emerging studies collectively underscore the adaptability and strength of PINNs in modeling complex dynamical systems, ranging from epidemic scenarios to tumor progression and systems biology. PINNs unite mechanistic knowledge and data driven learning with fractional, time varying and multi phase dynamics making them an ideal tool for realistic disease modeling across domains.

The structure of this paper is organized as follows. In Section 2, we introduce the mathematical preliminaries of fractional calculus and present the physics informed neural network (PINN) framework used to approximate fractional-order dynamical systems. Section 3 formulates the proposed fractional-order Alzheimer's disease model and provides a detailed theoretical analysis, including positivity, boundedness, existence, uniqueness, equilibrium points and stability properties. Section 4 presents sensitivity analysis and numerical simulations that illustrate the influence of key biological parameters and fractional order on disease progression. In Section 5, a fractional PINN-based optimal control framework is developed and its performance is analyzed. In addition, the proposed fractional PINN is benchmarked against a classical Grünwald Letnikov numerical solver to verify its ability to reproduce controlled fractional dynamics. Finally, the concluding section summarizes the findings and gives potential directions for future research.

2 Mathematical framework: Physics Informed Neural Networks

We consider a system of nonlinear ordinary or fractional differential equations governing the evolution of biological or physical quantities. Let D=[0,T] denote the temporal domain of interest and let u(t) ∈ ℝ^*n* represent the vector of state variables (e.g., neuronal population, amyloid-β concentration, etc.). The general form of the governing system is given by:

F(t,u(t),CDtσu(t);θ)=0, t∈D,(1)

with initial condition:

u(0)=u0∈ℝn.(2)

Here:

• F is a (possibly nonlinear) differential operator parameterized by θ,
• DtσCu(t) denotes the Caputo fractional derivative of order σ ∈ (0, 1],
• θ is a set of known or learnable parameters.

Definition 2.1. A Physics Informed Neural Network (PINN) is a neural network uϕ:Darrowℝn with parameters ϕ, trained to approximate the true solution u(t) such that:

• It minimizes the discrepancy with available data,
• It satisfies the governing differential equation(s).

• Td={ti(d)}i=1Nd⊂D be the set of data points,
• Tc={tj(c)}j=1Nc⊂D be the set of collocation points.

The total loss function used to train u_ϕ is:

L(ϕ)=Ldata(ϕ)+λphysLphys(ϕ),(3)

where:

Ldata(ϕ)=1Nd∑i=1Nd‖uϕ(ti(d))-uobs(ti(d))‖2,(4)

Lphys(ϕ)=1Nc∑j=1Nc‖F(tj(c),uϕ(tj(c)),CDtσuϕ(tj(c));θ)‖2,(5)

and λ_phys > 0 is a regularization parameter.

Remark 1. In practice, the Caputo derivative is approximated numerically using a discrete convolution formula such as the Grünwald-Letnikov scheme:

DtσCu(tn)≈1Δtσ∑k=0nwk(σ)u(tn−k),(6)

where:

• Δt is the time step size,
• wk(σ)=(−1)k(σk) are fractional weights,
• (σk)=Γ(σ+1)Γ(k+1)Γ(σ−k+1).

This allows the PINN to model nonlocal memory effects present in biological systems.

Definition 2.2. The residual function evaluated at the collocation point tj(c) is given by:

Rϕ(tj(c)):=F(tj(c),uϕ(tj(c)),CDtσuϕ(tj(c));θ).(7)

The physics loss term is then defined as:

Lphys(ϕ)=1Nc∑j=1Nc‖Rϕ(tj(c))‖2.(8)

Remark 2. Under the activation functions (e.g., tanh, ReLU), neural networks are universal function approximators. Thus, given sufficient, u_ϕ(t) can approximate the true solution u(t) arbitrarily well, provided that the optimization landscape is well conditioned.

Definition 2.3. (Kilbas et al., 2006) Assume that [a, b] ⊂ ℝ. Then the fractional integral of order σ for g ∈ L¹([a, b], ℝ) can be written as:

Itσg(t)=1Γ(σ)∫0t(t-v)(1-σ)g(v)dv

where t > 0, σ > 0 and integral on the right side is point-wise defined on ℝ⁺, ℝ⁺ = [0, ∞).

Definition 2.4. (Podlubny, 1999) Let g be a continuous function on [0, T]. The derivative of Caputo can be written as

Dtσ0Cg(t)=1Γ(n−σ)[∫0t(t−v)n−σ−1dndvng(t)(v)dv]

where n = ⌊σ⌋+1 and ⌊σ⌋ be the integer of β. Where 0 < σ < 1 then the Caputo derivatives will be:

Dtσ0Cg(t)=1Γ(n−σ)[∫0t(t−v)−σg′(t)(v)dv]

Lemma 1. (Bansal et al., 2023) Let t≥t₀ and let F:[t0,∞)→ℝ+ be a continuous function. For any σ ∈ (0, 1) and any constant F^* ∈ ℝ⁺, the following inequality holds:

DtσC​(F(t)−F−Fln​F(t)F) ≤ (1−FF(t))DtσCF(t),

where DtσC denotes the Caputo fractional derivative of order σ.

3 Fractional order model

Amyloid and tau are the two primary proteins that are hypothesized to obstruct brain cell to cell communication. The amyloid beta peptide accumulates in Alzheimer's disease and deposits as plaques around brain vasculature and neuronal cells. This deposition is linked to a decrease in neural function, which impairs memory and cognition and affects daily functioning including speaking, writing and thinking. Once amyloid beta has accumulated to a certain degree, aberrant tau begins to surge, following the initial appearance of amyloid beta clusters. Then, a positive feedback loop takes place leading to an increase in aberrant tau and amyloid beta formation. Numerous bacteria and viruses, as well as tau tangles and amyloid beta plaques can cause neuro-inflammation. The resident innate immune cells of the central nervous system known as microglia, have the ability to trigger the activation of inflammatory pathways and change their physiological function, which can accelerate the progression of disease. Microglia are believed to become activated when harmful amyloid beta and tau proteins are present. Toxic proteins and other detritus are removed from dead and dying cells by microglia. Neuronal dysfunction, damage and loss may arise from chronic inflammation caused by microglia's inability to keep up with everything that needs to be cleaned.

The model is given by:

D0σCFN(t)=ΠN+ρTμ−αFNAβ−ϕ1FND0σCIN(t)=αFNAβ−β1INMρ−(γ+ϕ2)IND0σCAβ(t)=γIN−β2AβMρ−(dβ+κ)AβD0σCTμ(t)=κAβ−β3TμMρ−(dμ+ρ)TμD0σCMρ(t)=(β1IN+β2Aβ+β3Tμ)Mρ−ϕ3Mρ(9)

The model parameters in Table 1 represents the biologically meaningful processes governing neuronal survival, protein aggregation and immune response. Parameters such as α, κ and γ quantify the interaction strength between amyloid beta, tau protein and neuronal compartments, while β₁, β₂ and β₃ describe microglial clearance rates of infected neurons and toxic proteins. Some parameters are assumed or estimated due to their limited availability of precise experimental measurements, however, their values are chosen within biologically plausible ranges consistent with existing literature. The small initial values assigned to amyloid beta and tau protein concentrations reflects an early stage pathological conditions and ensure physiological realism of the simulated trajectories.

Table 1. Initial variables states and parameters values.

3.1 Positivity and boundedness of solutions

For numerous forms of differential operators, including non-integer and integer orders, we provide a comprehensive analysis that substantiates the requirements for maintaining the positiveness of the suggested model solutions. To accomplish this, Now, define the Norm as,

‖B‖∞=supt∈DB|B(t)|(10)

D_*B* represents the domain of B. The definition given above can be used to generate the following inequality for the function F_*N*:

FN(t)=ΠN+ρTμ-αFNAβ-ϕ1FN, ∀t≥0 ≥-(αAβ+ϕ1)FN, ∀t≥0 ≥-(αsupt∈DO|Aβ|+ϕ1)FN, ∀t≥0 ≥-(α‖Aβ‖∞+ϕ1)FN, ∀t≥0(11)

This yields

FN(t)≥FN0exp-(α‖Aβ‖∞+ϕ1)t, ∀t≥0(12)

For the other functions:

IN(t)=αFNAβ-β1INMρ-(γ+ϕ2)IN, ∀t≥0 ≥-(β1Mρ+γ+ϕ2)IN, ∀t≥0 ≥-(β1‖Mρ‖∞+ϕ2)IN, ∀t≥0(13)

This yields

IN(t)≥IN0exp-(β‖Mρ‖∞+ϕ2)t, ∀t≥0.(14)

3.2 Positive solutions with non local operators

Here, we show that solutions with non-local operators are positive for a fractional calculus model. If every initial condition is met for non-local operators, then every solution is positive. The Caputo derivative gives,

FN(t)≥FN(0)Eσ(−(α‖Aβ‖∞+ϕ1)tσ), ∀t≥0IN(t)≥IN(0)Eσ(−(β‖Mρ‖∞+ϕ2)tσ), ∀t≥0Aβ(t)≥Aβ(0)Eσ(−(β2‖Mρ‖∞+dβ+κ)tσ), ∀t≥0Tμ(t)≥Tμ(0)Eσ(−(β3‖Mρ‖∞+dμ+ρ)tσ), ∀t≥0Mρ(t)≥Mρ(0)Eσ((β1‖IN‖∞+β2‖Aβ‖∞+β3‖Tμ‖∞ −ϕ3)tσ), ∀t≥0.(15)

3.3 Existence and uniqueness analysis

Now we use fixed point theory to investigate the existence and uniqueness of our model. Banach's contraction theorem will guaranties the model's singularity, whereas Schauder's fixed point theorem will guaranties its existence. Theorems will show that the our model has a unique solution are important since they suggest that your problem can be solved in a unique way. We can estimate a solution using numerical methods once we are certain that it exists. By using a fractional derivative for 0 < η ≤ 1 in the Caputo sense, system (9) can be made more generic. Let

ρ1(t,FN)=ΠN+ρTμ-αFNAβ-ϕ1FNρ2(t,IN)=αFNAβ-β1INMρ-(γ+ϕ2)INρ3(t,Aβ)=γIN-β2AβMρ-(dβ+κ)Aβρ4(t,Tμ)=κAβ-β3TμMρ-(dμ+ρ)Tμρ5(t,Mρ)=(β1IN+β2Aβ+β3Tμ)Mρ-ϕ3Mρ(16)

With the fractional integral and initial condition, the situation is:

FN(t)=FN0+1Γ(σ)∫0t(t-v)σ-1ρ1(t,FN)dvIN(t)=IN0+1Γ(σ)∫0t(t-v)σ-1ρ2(t,IN)dvAβ(t)=Aβ0+1Γ(σ)∫0t(t-v)σ-1ρ3(t,Aβ)dvTμ(t)=Tμ0+1Γ(σ)∫0t(t-v)σ-1ρ3(t,Tμ)dvMρ(t)=Mρ0+1Γ(σ)∫0t(t-v)σ-1ρ4(t,Mρ)dv(17)

Let

Θ(t)={FN(t)IN(t)Aβ(t)Tμ(t)Mρ(t),Θ0={FN0IN0Aβ0Tμ0Mρ0,ℵ(ξ,Θ(ξ))={ρ1(t,FN)ρ2(t,IN)ρ3(t,Aβ)ρ4(t,Tμ)ρ5(t,Mρ)(18)

Thus,

Θ(t)=Θ0+1Γ(σ)∫0t(t-v)σ-1ℵ(v,Θ(v))dv(19)

Now take a Banach space with a norm ℘[0, ℝ] = χ:

‖FN,IN,Aβ,Tμ,Mρ‖=maxt∈[0,ℝ][|FN,IN,Aβ,Tμ,Mρ|](20)

Let a mapping defined as ◇: χ → χ

◇Θ(t)=Θ0+1Γ(σ)∫0t(t-v)σ-1ℵ(v,Θ(v))dv(21)

Furthermore, we subject a nonlinear function to the next proposition:

• (P1) Constants ψ_*m*, ψ_*m* > 0 exist such that

• (P2) Every Θ,Θ̄∈χ has a constant L_*m* > 0 according to which

Theorem 1. If the hypotheses (P1) is correct, then the system (9) has at least one solution.

Proof. Let χ = C([0, τ], ℝ^*n*) with the sup-norm

‖Θ‖=supt∈[0,τ]‖Θ(t)‖.(24)

Assume that the nonlinear function ℵ:[0, τ] × ℝ^*n* → ℝ^*n* satisfies the growth and Lipschitz conditions

‖ℵ(t,Θ)‖≤ψm(‖Θ‖+1), ‖ℵ(t,Θ1)−ℵ(t,Θ2)‖≤L‖Θ1−Θ2‖,(25)

for all t ∈ [0, τ] and Θ, Θ₁, Θ₂ ∈ χ, where ψ_*m*, L ≥ 0 are constants.

Define

ψ≥maxt∈[0,τ]‖Θ0‖+ψmτσΓ(σ+1)1−ψmτσΓ(σ+1),(26)

and let

B:={Θ∈χ:‖Θ‖≤ψ}.(27)

Clearly, B is closed, convex and bounded in χ.

Now define the operator T:χ→χ by

(TΘ)(t)=Θ0+1Γ(σ)∫0t(t-v)σ-1ℵ(v,Θ(v))dv, t∈[0,τ].(28)

T maps B into B. For Θ ∈ B and t ∈ [0, τ],

‖(TΘ)(t)‖≤‖Θ0‖+1Γ(σ)∫0t(t−v)σ−1‖ℵ(v,Θ(v))‖dv ≤‖Θ0‖+1Γ(σ)∫0t(t−v)σ−1ψm(‖Θ‖+1)dv ≤‖Θ0‖+(ψm‖Θ‖+ψm)τσΓ(σ+1) ≤‖Θ0‖+(ψmψ+ψm)τσΓ(σ+1)≤ψ,(29)

by the choice of ψ. Thus, T(B)⊆B.

T is continuous and relatively compact. Let t₁<t₂ ∈ [0, τ]. For Θ ∈ B,

‖(TΘ)(t2)−(TΘ)(t1)‖ =1Γ(σ)‖∫0t2(t2−v)σ−1ℵ(v,Θ(v))dv −∫0t1(t1−v)σ−1ℵ(v,Θ(v))dv‖≤1Γ(σ)∫0t1[(t2−v)σ−1−(t1−v)σ−1]‖ℵ(v,Θ(v))‖dv +1Γ(σ)∫t1t2(t2−v)σ−1‖ℵ(v,Θ(v))‖dv ≤(ψmψ+ψm)t2σ−t1σΓ(σ+1).(30)

As t₂→t₁, the right-hand side tends to 0, so TΘ is equicontinuous on [0, τ]. Since ‖(TΘ)(t)‖≤ψ for all t, the set T(B) is uniformly bounded. By the Arzelà-Ascoli theorem T(B) is relatively compact in χ. Continuity of T follows from dominated convergence.

3.3.1 Existence of a solution

The operator T is continuous, maps B into a relatively compact subset of B and B is closed, bounded and convex. Hence, by Schauder's fixed point theorem T has at least one fixed point Θ^* ∈ B. Thus, the system has at least one solution.

3.3.2 Uniqueness under a contraction condition

If the Lipschitz constant L satisfies

LτσΓ(σ+1)<1,(31)

then for Θ₁, Θ₂ ∈ B,

‖(TΘ1)(t)-(TΘ2)(t)‖≤1Γ(σ)∫0t(t-v)σ-1‖ ℵ(v,Θ1(v))-ℵ(v,Θ2(v))‖dv ≤L‖Θ1-Θ2‖τσΓ(σ+1).(32)

Taking the supremum over t ∈ [0, τ] gives

‖TΘ1-TΘ2‖≤LτσΓ(σ+1)‖Θ1-Θ2‖.(33)

Thus, T is a contraction. By the Banach fixed point theorem the fixed point is unique.

Therefore, the system has at least one solution (by Schauder's theorem) and under the contraction condition this solution is unique (by Banach's theorem).

Theorem 2. Assume that the system (9) has a unique solution if the conditions (P2) are satisfied.

Proof. Assume that Θ̄,Θ̄1∈χ. Then,

‖◇(Θ̄)-◇(Θ̄1)‖=maxt∈[0,τ]|1Γ(σ)∫0t(t2-v)σ-1ℵ(v,Θ̄(v))dv. .-1Γ(σ)∫0t(t1-v)σ-1ℵ(v,Θ̄1(v))dv|(34)

≤τσΓ(σ+1)Lm|Θ̄-Θ̄1|.(35)

◇ is consequently a contraction. The system (9) has a unique solution according to the Banach fixed point theorem.

3.4 Equilibrium point analysis

In this part, the equilibrium points will be examined. The equilibrium points of the system (9) are: When IN0=0, the disease free equilibrium points (DEF) will be,

E0=(FN0,IN0,Aβ0,Tμ0,Mρ0)(36)

If IN0=0, then FN0=ΠNϕ1,Aβ0=0 and Tμ0=0, represent the disease free equilibrium points (DEF), then,

E0=(ΠNϕ1,0,0,0)(37)

Now, the endemic equilibrium points are:

E=(FN,IN,Aβ,Tμ)FN=γINκρ+ΠNdμdβ+ΠNdβ+ΠNdμκ+ΠNρκαdμγIN+αγINρ+ϕ1dμdβ+ϕ1dβρ+ϕ1dμκ+ϕ1κρ IN=IN Aβ=γINdβ+κ Tμ=κγIN*(dμ+ρ)(dβ+κ)(38)

The proposed system is absolutely stable.

3.5 The reproduction number

We now compute the basic reproduction number R₀ for the Alzheimer's model using the next-generation matrix method.

The infected states are

X(t)=(IN(t),Aβ(t),Tμ(t))⊤.

At the disease-free equilibrium (DFE), we have

IN=Aβ=Tμ=0, FN=ΠNϕ1.

Linearizing the equations around the DFE gives:

DtαIN=αinfFN*Aβ-(γab+ϕ2)IN,DtαAβ=γabIN-(db+κ)Aβ,DtαTμ=κAβ-(dm+ρ)Tμ.

We write the system as

DtαX=(F-V)X,

with

F=(0αinfFN*0γab000κ0), V=(γab+ϕ2000db+κ000dm+ρ).

Since V is diagonal, its inverse is

V−1=diag​((γab+ϕ2)−1, (db+κ)−1, (dm+ρ)−1).

Hence the next-generation matrix is

K=FV-1=(0αinfFN*db+κ0γabγab+ϕ2000κdb+κ0).

The characteristic polynomial is

det(K-λI)=λ(-λ2+αinfFN*db+κ·γabγab+ϕ2).

Thus the eigenvalues are

λ1=0, λ2,3=±αinfFN*db+κ·γabγab+ϕ2.

The spectral radius of K (the largest eigenvalue in absolute value) is

R0=αinfγab(db+κ)(γab+ϕ2)·ΠNϕ1.

Figure 1 shows how variations in key parameters influence the basic reproduction number R₀, reflecting the system's sensitivity.

3.6 Iterative and stability analysis via Caputo operator

Theorem 3. Let (𝔉, |·|) be a Banach space and 𝔈: 𝔉 → 𝔉 be a mapping that satisfies:

‖Em-En‖≤𝔍‖X-Em‖+π‖m-n‖,(39)

for all m, n ∈ 𝔉, where 0 ≤ 𝔍 and 0 ≤ π < 1. Then 𝔈 is Picard 𝔈-stable.

Theorem 4. Assume that 𝔈 is a self map with the following definition:

E[FNα(t)]=FNα+1(t)=FNα(0)+L-1[1EσL{ΠN+ρTμα-αFNαAβα-ϕ1FNα}],E[INα(t)]=INα+1(t)=INα(0)+L-1[1EσL{αFNαAβα-β1INαMρα-(γ+ϕ2)INα}],E[Aβα(t)]=Aβα+1(t)=Aβα(0)+L-1[1EσL{γINα-β2AβαMρα-(dβ+κ)Aβα}],E[Tμα(t)]=Tμα+1(t)=Tμα(0)+L-1[1EσL{κAβα-β3TμαMρα-(dμ+ρ)Tμα}],E[Mρα(t)]=Mρα+1(t)=Mρα(0)+L-1[1EσL{(β1INα+β2Aβα+β3Tμα)Mρα-ϕ3Mρα}].(40)

The iteration is 𝔄-stable in L¹(α, β) if the following criteria are met:

(1−ρw1(φ)−α(K1+K2)w2(φ)−ϕ1w3(φ))<1,(1+α(K1+K2)w2(φ)−β1(K3+K4)w4(φ)−(γ+ϕ2)w5(φ))<1,(1+γw6(φ)−β2(K5+K6)w7−(dβ+κ)w8(φ))<1,(1+κw8(φ)−β3(K7+K8)w9−(dμ+ρ)w10(φ))<1,(1+β1(K3+K4)(w4(φ)+β2(K5+K6)w7)+β3(K7+K8)w9−ϕ3w11(φ))<1.(41)

Proof. We examine the following for (α, β) ∈ ℕ × ℕ in order to show that 𝔈 has a fixed point:

E[FNβ(t)]-E[FNα(t)]=FNβ(t)-FNα(t)+L-1[1EσL{ΠN+ρTμβ-αFNβAββ-ϕ1FNβ}]-L-1[1EσL{ΠN+ρFNαAβα-ϕ1FNα}].(42)

E[INβ(t)]-E[INα(t)]=INβ(t)-INα(t)+L-1[1EσL{αFNβAββ-β1INβMρβ-(γ+ϕ2)INβ}]-L-1[1EσL{αFNαAβα-β1INαMρα-(γ+ϕ2)INα}].E[Aββ(t)]-E[Aβα(t)]=Aββ(t)-Aβα(t)+L-1[1EσL{γINβ-β2AββMρβ-(dβ+κ)Aββ}]-L-1[1EσL{γINα-β2AβαMρα-(dβ+κ)Aβα}].(43)

E[Tμβ(t)]-E[Tμα(t)]=Tμβ(t)-Tμα(t)+L-1[1EσL{κAββ-β3TμβMρβ-(dμ+ρ)Tμβ}]-L-1[1EσL{κAβα-β3TμαMρα-(dμ+ρ)Tμα}].(44)

E[Mρβ(t)]-E[Mρα(t)]=Mρβ(t)-Mρα(t)+L-1[1EσL{(β1INβ+β2Aββ+β3Tμβ)Mρβ-ϕ3Mρβ}]-L-1[1EσL{(β1INα+β2Aβα+β3Tμα)Mρα-ϕ3Mρα}].(45)

By calculating the norm of both sides of the above equations, we obtain:

‖E[FNβ]−E[FNα]‖≤‖FNβ−FNα‖+ℒ−1[1Eσℒ{‖ρ(Tμβ−Tμα)‖−‖αFNβ(Aββ−Aβα)‖+‖αAβα(FNβ−FNα)‖+‖−ϕ1(FNβ−FNα)‖+⋯.(46)

‖E[INβ(t)]−E[INα(t)]‖≤‖INβ(t)−INα(t)‖+ℒ−1[1Eσℒ{‖−(γ+ϕ2)(INβ−INα)‖+‖αFNβ(Aββ−Aβα)‖+‖αAβα(FNβ−FNα)‖+⋯.(47)

‖E[Aββ(t)]−E[Aβα(t)]‖≤‖Aββ(t)−Aβα(t)‖+ℒ−1[1Eσℒ{‖γ(INβ−INα)‖−‖β2Aββ(Mρβ−Mρα)‖+‖β2Mρα(Aββ−Aβα)‖+⋯.(48)

‖E[Tμβ(t)]−E[Tμα(t)]‖≤‖Tμβ(t)−Tμα(t)‖+ℒ−1[1Eσℒ{‖κ(Aββ−Aβα)‖−‖β3Tμβ(Mρβ−Mρα)‖+‖β3Mρα(Tμβ−Tμα)‖+⋯.(49)

‖E[Mρβ(t)]−E[Mρα(t)]‖≤‖Mρβ(t)−Mρα(t)‖+ℒ−1[1Eσℒ{‖(β1(INβ−INα)‖+‖β2(Aββ−Aβα)‖ +‖β3(Tμβ−Tμα)(Mρβ−Mρα)‖+⋯.(50)

Let,

‖FNβ(t)-FNα(t)‖≅‖INβ(t)-INα(t)‖≅‖Aββ(t)-Aβα(t)‖≅‖Tμβ(t)-Tμα(t)‖≅‖Mρβ(t)-Mρα(t)‖.(51)

To get the following relation:

‖E[FNβ(t)]-E[FNα(t)]‖≤(1-ρw1(φ)-α(K1+K2)w2(φ)-ϕ1w3(φ))‖FNβ-FNα‖,(52)

‖E[INβ(t)]−E[INα(t)]‖≤(1+α(K1+K2)w2(φ)−β1(K3+K4)w4(φ)−(γ+ϕ2)w5(φ))‖INβ−INα‖,(53)

‖E[Aββ(t)]-E[Aβα(t)]‖≤(1+γw6(φ)-β2(K5+K6)w7-(dβ+κ)w8(φ))‖Aββ-Aβα‖,(54)

‖E[Tμβ(t)]−E[Tμα(t)]‖≤(1+κw8(φ)−β3(K7+K8)w9−(dμ+ρ)w10(φ))‖Tμβ−Tμα‖,(55)

‖E​[Mρβ]−E​[Mρα]‖≤(1+β1​(K3+K4))(w4(φ)+β2​(K5+K6)w7+β3​(K7+K8)w9−ϕ3w11(φ))‖Mρβ−Mρα‖.(56)

Also, the convergent sequences F_{N_β}, I_{N_β}, A_{β_β}, T_{μ_β} and M_{ρ_β} are bounded. Then we obtain positive constants K₁, K₂, K₃, K₄ and K₅ for all t such that

‖FNα‖<K1, ‖INα‖<K2, ‖Aβα‖<K3, ‖Tμα‖<K4, ‖Mρα‖<K5, (β,α)∈ℕ×ℕ.(57)

Then from Equations 56–60, 61, we get

‖E[FNβ(t)]−E[FNα(t)]‖≤(1−ρw1(φ)−α(K1+K2)w2(φ)−ϕ1w3(φ))‖FNβ−FNα‖.(58)

Where 𝔴₁, 𝔴₂, 𝔴₃, 𝔴₄ and 𝔴₅ are functions of L-1{1EσL}.

Similarly,

‖E[INβ]−E[INα]‖≤(1+α(K1+K2)w2(φ)−β1(K3+K4)w4(φ)−(γ+ϕ2)w5(φ))‖INβ−INα‖,‖E[Aββ(t)]−E[Aβα(t)]‖≤(1+γw6(φ)−β2(K5+K6)w7−(dβ+κ)w8(φ))‖Aββ−Aβα‖,‖E[Tμβ(t)]−E[Tμα(t)]‖≤(1+κw8(φ)−β3(K7+K8)w9−(dμ+ρ)w10(φ))‖Tμβ−Tμα‖,‖E​[Mρβ]−E​[Mρα]‖≤(1+β1​(K3+K4))(w4(φ)+β2​(K5+K6)w7+β3​(K7+K8)w9−ϕ3w11(φ))‖Mρβ−Mρα‖.(59)

Where

(1−ρw1(φ)−α(K1+K2)w2(φ)−ϕ1w3(φ))<1,(1+α(K1+K2)w2(φ)−β1(K3+K4)w4(φ)−(γ+ϕ2)w5(φ))<1,(1+γw6(φ)−β2(K5+K6)w7−(dβ+κ)w8(φ))<1,(1+κw8(φ)−β3(K7+K8)w9−(dμ+ρ)w10(φ))<1,(1+β1(K3+K4)(w4(φ)+β2(K5+K6)w7+β3(K7+K8)w9−ϕ3w11(φ)))<1.(60)

As 𝔈 has a fixed point and using Equations 62, 63, let

π=(0,0,0,0,0),(61)

ℑ={(1−ρw1(φ)−α(K1+K2)w2(φ)−ϕ1w3(φ)),(1+α(K1+K2)w2(φ)−β1(K3+K4)w4(φ)−(γ+ϕ2)w5(φ)),(1+γw6(φ)−β2(K5+K6)w7−(dβ+κ)w8(φ)),(1+κw8(φ)−β3(K7+K8)w9−(dμ+ρ)w10(φ)),(1+β1(K3+K4)(w4(φ)+β2(K5+K6)w7+β3(K7+K8)w9−ϕ3w11)).(62)

This completes the proof.

3.7 Local stability analysis

Definition 3.1. The Hartman-Grobman theorem states that if the linearization of the equations produces no zero or imaginary eigenvalues, then there is a continuous function with a continuous inverse in the region of this point into ℝ^*n*.

The Jacobian matrix of the model (F_*N*, I_*N*, A_β, T_μ) is given by:

[-ϕ10-αΠNϕ1ρ0-γ-ϕ2αΠNϕ100γ-dβ-κ000κ-dμ-ρ]

All of the eigenvalues, which are calculated using Maple software, were found to have negative real values after verification. Thus, the equilibrium point can be considered locally stable.

3.8 Global stability analysis using Lyapunov function

Theorem 5. Consider the Alzheimer's fractional order system with the endemic equilibrium

E=(FN,IN,Aβ,Tμ,Mρ).

If the condition

Υ1<Υ2,

holds, then the equilibrium E^* is globally asymptotically stable in the invariant region.

Proof. Assume the Volterra-type Lyapunov function:

L=L1(FN-FN-FNlogFNFN)+L2(IN-IN-INlogININ) +L3(Aβ-Aβ-AβlogAβAβ)+L4(Tμ-Tμ-TμlogTμTμ) +L5(Mρ-Mρ-MρlogMρMρ*)

Assuming the positive constants L_*i*, i = 1, 2, 3, 4, 5 and replacing Equation 53 in the system, Lemma 1 can be applied.

DtσCL≤L1(FN−FNFN)cDtσFN+L2(IN−ININ)cDtσIN +L3(Aβ−AβAβ)cDtσAβ+L4(Tμ−TμTμ)cDtσTμ +L5(Mρ−Mρ*Mρ)cDtσMρ.

Now, modify the preceding formula's derivative to obtain:

DtσCL≤L1(FN−FNFN)(ΠN+ρTμ−αFNAβ−ϕ1FN) +L2(IN−ININ)(αFNAβ−β1INMρ−(γ+ϕ2)IN) +L3(Aβ−AβAβ)(γIN−β2AβMρ−(dβ+κ)Aβ) +L4(Tμ−TμTμ)(κAβ−β3TμMρ−(dμ+ρ)Tμ) +L5(Mρ−Mρ*Mρ)((β1IN+β2Aβ+β3Tμ) Mρ−ϕ3Mρ).

Now, by replacing FN=FN-FN,IN=IN-IN,Aβ=Aβ-Aβ,Tμ=Tμ-Tμ,Mρ=Mρ-Mρ*, we obtain:

DtσCL≤L1(FN−FNFN)(ΠN+ρ(Tμ−Tμ)−α(FN−FN) (Aβ−Aβ)−ϕ1(FN−FN)) +L2(IN−ININ)(α(FN−FN)(Aβ−Aβ)−β1(IN−IN) (Mρ−Mρ)−(γ+ϕ2)(IN−IN)) +L3(Aβ−AβAβ)(γ(IN−IN)−β2(Aβ−Aβ) (Mρ−Mρ)−(dβ+κ)(Aβ−Aβ)) +L4(Tμ−TμTμ)(κ(Aβ−Aβ)−β3(Tμ−Tμ) (Mρ−Mρ)−(dμ+ρ)(Tμ−Tμ)) +L5(Mρ−MρMρ)(β1(IN−IN)+β2(Aβ−Aβ) +β3(Tμ−Tμ))(Mρ−Mρ)−ϕ3(Mρ−Mρ*).

Now, assume L₁ = L₂ = L₃ = L₄ = L₅ = 1 and by rearranging the equation above, we obtain:

DtσCL≤ΠN−ΠNFNFN+ρTμ−ρTμ−ρTμFNFN +ρTμFNFN−αAβ(FN−FN)2FN+αAβ(FN−FN)2FN −ϕ1(FN−FN)2FN+αFNAβ−αFNAβ−αFNAβ +αFNAβ−αFNAβININ+αFNAβININ +αFNAβININ−αFNAβININ−β1Mρ(IN−IN)2IN +β1Mρ(IN−IN)2IN −(γ+ϕ2)(IN−IN)2IN+γIN−γIN−γINAβAβ +γINAβAβ−β2Mρ(Aβ−Aβ)2Aβ +β2Mρ(Aβ−Aβ)2Aβ−(dβ+κ)(Aβ−Aβ)2Aβ +κAβ−κAβ−κAβTμTμ +κAβTμTμ−β3Mρ(Tμ−Tμ)2Tμ+β3Mρ(Tμ−Tμ)2Tμ −(dμ+ρ)(Tμ−Tμ)2Tμ +β1IN(Mρ−Mρ)2Mρ−β1IN(Mρ−Mρ)2Mρ +β2Aβ(Mρ−Mρ)2Mρ−β2Aβ(Mρ−Mρ)2Mρ +β3Tμ(Mρ−Mρ)2Mρ−β3Tμ(Mρ−Mρ)2Mρ−ϕ3 (Mρ−Mρ*)2Mρ.

By using the following assumption,

DtσCL≤Υ1−Υ2

where

Υ1=ΠN+ρTμ+ρTμFNFN+αAβ(FN-FN)2FN+αFNAβ +αFNAβ+αFNAβININ +αFNAβININ+β1Mρ(IN-IN)2IN+ϕ2(IN-IN)2IN +γIN+γINAβAβ+β2Mρ(Aβ-Aβ)2Aβ+κAβ +κAβTμTμ+β3Mρ(Tμ-Tμ)2Tμ+β1IN(Mρ-Mρ)2Mρ +β2Aβ(Mρ-Mρ)2Mρ+β3Tμ(Mρ-Mρ*)2Mρ,

and

Υ2=ΠNFNFN+ρTμ+ρTμFNFN+αAβ(FN-FN)2FN +ϕ1(FN-FN)2FN+αFNAβ+αFNAβ +αFNAβININ+αFNAβININ+β1Mρ(IN-IN)2IN+(γ+ϕ2) (IN-IN)2IN+γIN +γINAβAβ+β2Mρ(Aβ-Aβ)2Aβ+(dβ+κ)(Aβ-Aβ)2Aβ +κAβ+κAβTμTμ +β3Mρ(Tμ-Tμ)2Tμ+(dμ+ρ)(Tμ-Tμ)2Tμ +β1IN(Mρ-Mρ)2Mρ+β2Aβ(Mρ-Mρ)2Mρ +β3Tμ(Mρ-Mρ)2Mρ+ϕ3(Mρ-Mρ)2Mρ.

If

Υ1<Υ2arrowcDtσL<0,

now, let FN=FN,IN=IN,Aβ=Aβ,Tμ=Tμ,Mρ=Mρ, then it can be concluded that E=(FN,IN,Aβ,Tμ,Mρ*) is the compact invariant set in

{(FN,IN,Aβ,Tμ,Mρ*)∈Δ:cDtσL=0}.

Hence, if Υ₁ < Υ₂, then the equilibrium points E^* are globally asymptotically stable in the invariant region.

3.9 Chaos stabilizing system

Utilizing the points of equilibrium, the linear output technique can be used to stabilize the suggested system (Equation 9), which is considered as a controlled-design fractional-order system (Nisar et al., 2024a,b).

Let ω₁, ω₂, ω₃, ω₄ and ω₅ be the controlled parameters, while FN,IN,Aβ* and Tμ* are the proposed model equilibrium points and a Jacobian matrix can be constructed to be as:

-ϕ1-ω10-αΠNϕ1ρ0-γ-ϕ2-ω2αΠNϕ100γ-dβ-κ-ω3000κ-dμ-ρ-ω4

Assuming ω₁ = 1, ω₂ = 2, ω₃ = 3, ω₄ = 4 and ω₅ = 5, with the values of the parameters ϕ₁ = 0.00003, α = 0.08, Π_*N* = 700, ρ = 0.025, γ = 0.00017, ϕ₂ = 0.00019, d_β = 9.51, κ = 0.025, d_μ = 0.277.By using Maple software, the roots are:

λ1=-1.00003000000000, λ2=-4.30200000000000,λ3=-25.8118817160782, λ4=-11.3265217160782(64)

Since, all of of the eigenvalues in Equation are real and negative numbers, the equilibrium points are asymptotically stable.

As illustrated in Figure 2, the dynamics of the Alzheimer's disease model compartments over 200 days are shown. The temporal evolution of functional neurons (F_*N*), infected neurons (I_*N*), amyloid beta (A_β), tau protein (T_μ) and microglia (M_ρ) are depicted.

3.10 Logarithmic sensitivity analysis

We evaluate how small perturbations in key parameters affect model outcomes. The logarithmic sensitivity is defined as:

Slog(t)=∂ln y(t)∂ln p≈ln y+(t)-ln y-(t)2ln(1+δ)(65)

where:

• y₊(t): output with parameter p(1+δ)
• y₋(t): output with parameter p(1−δ)
• δ = 0.1 (10% perturbation)

This formulation captures both the direction and magnitude of sensitivity in a normalized way.

4 Results interpretation

The time-dependent log-sensitivities of the Alzheimer model variables to ±10% perturbations in parameters are illustrated in Figure 3. The analysis reveals the following observations:

• *Functional neurons (F_N)* (Figure 3a): High negative sensitivity to α (amyloid beta toxicity) indicating vulnerability to A_β buildup.
• *Infected neurons (I_N)* (Figure 3b): Sensitive to both α and β₁; infection rises with A_β and is controlled by microglia.
• *Amyloid beta (A_β)* (Figure 3c): Strongly impacted by α, β₁ and β₂ reflecting production and clearance balance.
• *Tau protein (T_μ)* (Figure 3d): Follows the dynamics of A_β and shows sensitivity to κ and β₃.
• *Microglia (M_ρ)* (Figure 3e): Fluctuates mildly, reflecting its stabilizing regulatory role.

Logarithmic sensitivity analysis using NSFD reveals nonlinear time varying dependence of disease progression on key biological parameters (see Table 2). This method guides strategic control interventions and enhances the understanding of neuro degenerative dynamics.

Table 2. Summary.

4.1 Mapping model states to measurable biomarkers

The state variables of the proposed fractional-order model shows direct correspondence with clinically measurable biomarkers used in Alzheimer's disease diagnosis and progression assessment. The functional neuron population F_*N* is associated with neuroimaging-] derived indicators such as brain volume, cortical thickness and cognitive performance scores including MMSE and ADAS-Cog. Amyloid beta concentration A_β corresponds to cerebrospinal fluid (CSF) Aβ₄₂ levels and amyloid PET imaging, while tau protein T_μ aligns with CSF phosphorylated tau biomarkers and tau PET tracers. The microglial compartment M_ρ reflects neuroinflammatory activity and can be linked to TSPO-PET imaging markers. This mapping enables the proposed framework to bridge mechanistic modeling with observable clinical data thus facilitating data-driven inference and model calibration.

4.2 Simulations

The simulation of the Alzheimer's disease model using the NSFD scheme with ±10% parameter variations in neuronal infection rate, initial Tau protein concentration, Neuron production rate and the amyloid to tau conversion rate reveals biologically consistent dynamics across all compartments are shown in Figures 4–8. Functional neurons show a gradual decline interspersed with fluctuations, where increased neuronal production mitigates losses while higher infection rates accelerate degeneration. Infected neurons (I_*N*) exhibit spike like growth patterns whose magnitude and frequency are strongly affected by neuronal infection rate and Neuron production rate suggesting that even small changes in infection or neurogenesis processes can significantly shift neuronal vulnerability. Amyloid beta (A_β) demonstrates recurrent aggregation peaks, with the timing and amplitude of these cycles markedly influenced by all four parameters highlighting their central role in driving pathological progression. Tau protein (T_μ) oscillations are particularly sensitive to the initial tau burden (T_*m*0) and the amyloid tau interaction rate (b₃), where larger values accelerate accumulation and intensify spikes. Microglia (M_ρ) display fluctuating activation, with bursts modulated by both neuronal production and infection processes, reflecting their downstream immune response to amyloid and neuronal damage. Overall, these results indicate that modest variations in critical parameters substantially reshape the long term trajectories of all biomarkers reinforcing the model's ability to capture the delicate balance between neuronal survival protein aggregation and immune activation in Alzheimer's disease.

Figures 9–13 represents the dynamics of the Alzheimer's disease model under different values of the fractional order α ∈ {0.7, 0.8, 0.9, 1.0}, simulated over a period of 500 days. The evolution of the biological components: functional neurons (F_*N*), infected neurons (I_*N*), amyloid beta (A_β), tau protein (T_μ) and microglia (M_ρ) is presented to demonstrate the influence of fractional memory on system behavior.

For the functional neurons (F_*N*), the population increases over time with lower fractional orders (e.g., α = 0.7) showing faster growth compared to the integer order case (α = 1.0). This indicates that systems with stronger memory effects (smaller α) promote enhanced neuronal recovery and resilience. In contrast, infected neurons (I_*N*) decay monotonically toward zero for all α values with smaller α leading to a more rapid decay. This suggests that fractional dynamics facilitate faster infection clearance potentially due to increased microglial responsiveness or reduced propagation delay.

The concentration of amyloid beta (A_β) remains nearly constant around zero with very small amplitude (on the order of 10⁻⁶), implying that the production and clearance of amyloid are balanced and largely insensitive to variations in α within the observed time frame. Similarly, the tau protein (T_μ) exhibits a rapid exponential decay to zero across all cases, with slightly faster clearance for smaller α, indicating that memory-based fractional effects contribute to accelerated tau removal. Finally, the microglia population (M_ρ) starts around 0.02g/ml and decays over time, showing faster stabilization for smaller α values. This behavior implies that fractional memory enhances early microglial activation and leads to quicker stabilization, whereas the integer order dynamics sustain activation longer.

Overall, decreasing the fractional order α increases the rate of stabilization across all biological compartments, reflecting stronger memory and nonlocality effects. Smaller α values correspond to subdiffusive, memory-driven dynamics that accelerate system responses and stabilization, while larger α values approach the behavior of classical ordinary differential equations with slower transitions. Thus, the fractional order model captures the intrinsic temporal memory and delayed responses characteristic of Alzheimer's disease progression more effectively than the integer order formulation.

5 Optimal control: physics informed neural network approach

5.1 Modeling framework and neural network methodology

The modeling framework is based on a system of nonlinear fractional-order differential equations describing the progression of Alzheimer's disease under therapeutic intervention. The model consists of five biologically relevant state variables: functional neurons (F_*N*), inflamed neurons (I_*N*), amyloid-β concentration (A_β), tau protein concentration (T_μ), and microglial activity (M_ρ). Two time-dependent control variables, z₁(t) and z₂(t), are incorporated to represent therapeutic interventions targeting amyloid accumulation and neurodegeneration, respectively.

The governing equations are formulated using the Caputo fractional derivative, which captures long-term memory effects and nonlocal temporal dependence–key features of chronic neurodegenerative processes. For numerical implementation, the Caputo derivative is approximated using the Grünwald-Letnikov (GL) discretization, which expresses the fractional derivative as a weighted convolution of historical states, enabling stable evaluation of memory-dependent dynamics.

Synthetic reference trajectories for all state variables were generated by numerically solving the governing system using fixed parameter values and initial conditions. These reference trajectories serve as ground truth for training and evaluation. To assess robustness, additive Gaussian noise of varying intensity (0%–20%) was injected into the reference data during evaluation. This strategy enables a systematic investigation of stability and generalization under noisy observations without altering the underlying dynamics.

Three learning frameworks were investigated:

• Fractional-order PINN (proposed method),
• Integer-order PINN obtained by setting the derivative order to unity,
• A purely data-driven neural network serving as a baseline model.

All models were trained and evaluated under identical network architectures, optimization settings and time discretizations to ensure a fair comparison.

The Physics-Informed Neural Network (PINN) incorporates the governing equations directly into the loss function through physics-based residuals. In contrast, the integer-order PINN enforces classical integer-order dynamics, while the baseline model minimizes only a mean squared error loss without any physics constraints. This unified experimental design allows the influence of fractional-order memory and physical supervision to be isolated and quantified.

5.2 Neural network architecture and governing equations

All three models employ the same fully connected feedforward neural network architecture, which maps scalar time inputs to vector-valued outputs:

[FN(t),IN(t),Aβ(t),Tμ(t),Mρ(t),z1(t),z2(t)].

The network consists of multiple hidden layers with nonlinear activation functions, enabling the approximation of complex temporal dynamics. Using an identical architecture across all models ensures that observed performance differences arise solely from the imposed physical constraints rather than architectural bias.

For the fractional PINN, each state output is constrained by the governing fractional-order equations. The Caputo derivative is approximated using the Grünwald-Letnikov scheme. Let u(t) denote a network-predicted state variable and wj(α) represent the GL weights corresponding to fractional order α. The fractional derivative at time t_*n* is approximated as

Dtαu(tn)≈1hα∑j=0nwj(α)u(tn-j),(66)

where h is the uniform time step size. These residuals are incorporated into the training loss to enforce the fractional dynamics.

For the integer-order PINN, the same formulation is applied with derivative order α = 1, thereby enforcing classical integer-order dynamics while retaining the physics-informed training structure. The baseline neural network is trained purely in a data-driven manner, minimizing the mean squared error between predicted and reference state trajectories without learning control variables or enforcing physical laws. The key computational aspects of the implementation include: Fractional derivatives computed via convolution based Grünwald-Letnikov kernels, Automatic differentiation for residual evaluation, Separate evaluation of fractional order, integer-order and baseline data-driven models. Monte Carlo noise injection at multiple noise levels (0%–20%),

6 Simulation results and comparative analysis

To assess the effectiveness of the Physics Informed Neural Network (PINN) relative to a purely data driven Baseline Neural Network (NN) simulations were carried out for all compartments of the Alzheimer's disease model. Both models were trained on synthetic data, where Gaussian noise was deliberately added only for the baseline model to simulate measurement error. The PINN model in contrast, was trained not only to fit the noisy data but also to respect the underlying fractional order governing equations thereby enforcing physical consistency.

The results (Figure 14) were visualized through overlaid plots for each compartment variable, including both neural and inflammatory dynamics as well as the control interventions. These plots offer a clear comparison of the predicted trajectories from both models over a simulation period of 1000 days.

The variable F_*N*(t), representing functional neurons displays similar dynamics in both models, showing a sigmoidal growth with saturation around day 500. Both the baseline and PINN models closely approximate the expected trajectory indicating robust prediction for this compartment.

In contrast, the inflamed neurons I_*N*(t) show significant model divergence. The baseline model starts from a higher initial value and undergoes a rapid decay within the first 100 days before stabilizing around 0.002. The PINN, however begins with a lower value and steadily increases, reflecting a biologically plausible growth trajectory constrained by the governing equations. This suggests that the PINN captures inflammation dynamics more accurately under partial observability.

The amyloid beta concentration A_β(t) behaves distinctly across the models. The PINN output shows a sharp decline and remains suppressed throughout the simulation, whereas the baseline model maintains an increasing or constant trend. This discrepancy highlights the value of physics informed constraints in enforcing expected therapeutic behavior.

For the tau protein T_μ(t), the baseline model maintains a nearly constant value around 0.7, while the PINN shows a rapid early drop to nearly zero, remaining at that level throughout the simulation. This indicates that the PINN enforces realistic decay under the modeled treatment scenario while the baseline lacks the capacity to do so based on noisy data alone.

Microglial activity M_ρ(t) also shows differing behavior. The baseline model maintains a high constant value (0.8), lacking any dynamic response. Conversely, the PINN begins at a value of 1 rapidly declines to 0.8 in the first 100 days and gradually stabilizes near 0.8. The decline aligns with the model's representation of therapeutic effects modulating inflammation.

The control interventions z₁(t) and z₂(t), responsible for targeting amyloid and inflammatory pathways show a consistent pattern. The baseline model outputs nearly constant values of 0.7 for both, indicating no learned adaptation. The PINN however predicts an immediate drop to near zero values and maintains that level implying the network has identified minimal intervention as sufficient under the given conditions, as governed by the physics based dynamics.

This plot illustrates the sensitivity of the three models to increasing observation noise. As the noise level increases from 0 to 20 percent, the fractionalorder PINN exhibits only a mild increase in final loss, maintaining values several orders of magnitude lower than the integer order PINN. In contrast, the integer order PINN shows a rapid degradation in performance, with final loss increasing sharply as noise grows, indicating poor robustness to perturbations.

The baseline model achieves very low loss at 0 percent noise but does not encode any physical structure; therefore, its apparent stability does not reflect correct dynamical behavior and deteriorates in predictive quality when evaluated beyond the training regime. This highlights that low loss alone is insufficient without physical consistency.

The average training time is shown in Figure 15 for all three models. Both the fractional and integer-order PINNs require comparable computational effort due to the inclusion of physics-based residuals in the loss function. The baseline model trains significantly faster because it optimizes only a data-driven objective. Importantly, the modest additional cost of the fractional PINN is justified by its substantial gains in robustness, stability and interpretability.

The MSE for the functional neuron population (Figure 15) clearly demonstrates the superiority of the fractional PINN. Across all noise levels, the fractional model consistently produces orders of magnitude lower error than the integer order PINN. The integer-order model suffers from severe error amplification under noise, indicating that classical integer dynamics are insufficient to capture the long-term memory effects intrinsic to neurodegenerative processes. The baseline model maintains low error numerically but lacks physiological reliability due to the absence of governing constraints.

The bar chart (Figure 16) quantifies robustness as the ratio

Robustness=Loss at 20% noiseLoss at 0% noise.

The fractional PINN exhibits a robustness ratio approximately 2–3 orders of magnitude smaller than the integer order PINN, confirming its superior resistance to noise. The integer order model shows extreme sensitivity, with loss exploding as noise increases. The baseline model appears stable numerically but fails to generalize physically meaningful dynamics.

The convergence curves Figure 16 show that the fractional PINN converges smoothly and stably, even at higher noise levels. Loss decreases monotonically with minimal oscillations, indicating well-conditioned optimization guided by fractional physics constraints. In contrast, the integer order PINN exhibits oscillatory behavior and early stagnation, particularly under noisy conditions. This reflects instability in enforcing integer order dynamics that do not align with the true memory-dependent system.

The fractional PINN closely reproduces the GL reference solution across the entire simulation horizon, exhibiting smooth trajectories and consistent control-induced behavior. In contrast, the integer-order PINN deviates noticeably from the GL solution and displays oscillatory dynamics, particularly during transient phases where memory effects dominate.

As shown in Table 3, the final loss and training time statistics for the baseline, integer-order PINN and fractional-order PINN models are presented. The mean squared error (MSE) metrics for these models are detailed in Table 4.

Table 3. Final loss and training time statistics for baseline, integerorder PINN and fractional order PINN models.

Table 4. Mean squared error (MSE) metrics for baseline, integer-order PINN and fractional-order PINN models.

These results collectively demonstrate that incorporating fractional-order memory into a physics-informed neural network significantly enhances robustness, stability and physiological interpretability compared to both integer order PINNs and purely data-driven models, particularly under noisy and realistic observation conditions.

6.1 Comparison with a classical fractional numerical solver

To further validate the proposed fractional-order PINN framework and address concerns regarding numerical reliability, we compare its predictions with a classical Grünwald Letnikov (GL) finite-difference scheme applied to the controlled fractional-order Alzheimer's disease model. The GL method serves as a well-established reference for simulating Caputo-type fractional dynamics, explicitly accounting for the long-term memory effects inherent in neurodegenerative processes.

Figure 17 illustrates the evolution of the functional neuron population under optimal control obtained using the GL solver. For Therapeutic interventions for Alzheimer's disease, Optimal Control emerges as the most effective strategy. This approach not only increases the concentration of functional neurons over time but also significantly reduces the levels of infected neurons, amyloid beta, tau protein and microglia activity. By dynamically adjusting control variables Optimal Control achieves the desired therapeutic outcomes more efficiently.

On the other hand, Constant Control demonstrates moderate effectiveness. While it maintains a steady state of functional neurons and reduces infected neurons to some extent, it does not match the level of improvement seen with Optimal Control. The static nature of Constant Control limits its ability to adapt and optimize therapeutic effects, making it less effective in addressing the complexities of Alzheimer's disease progression.

In summary, these results demonstrate the superiority of the PINN approach in capturing biologically meaningful dynamics across all compartments, especially under conditions of noise and limited data. The baseline model, trained solely on data, exhibits static or overly smooth output that fails to reflect system behavior accurately. In contrast, PINN leverages both data and domain knowledge, resulting in more robust, interpretable and clinically consistent predictions. Conceptual radar chart comparing common strengths and limitations of different modeling approaches across multiple research dimensions, with higher values denoting stronger relative performance shown in Figure 18.

7 Conclusion

We developed a fractional order compartmental model of Alzheimer's disease that integrates amyloid beta accumulation, tau pathology, microglial response and neuronal dynamics, using Caputo derivatives to capture long-term memory effects. Analytical results demonstrated positivity, boundedness and stability properties while the reproduction number R₀ provided a threshold criterion linking amyloid production, clearance, tau conversion and the pool of susceptible neurons. Sensitivity and elasticity analyzes highlighted amyloid infectivity and clearance as the dominant factors shaping both disease initiation and long term trajectories, with tau progression playing a secondary role in sustaining pathology. Numerical simulations confirmed biologically consistent outcomes and the integration of Physics Informed Neural Networks (PINNs) showed strong potential for reconstructing hidden states and predicting optimal interventions under sparse or noisy data conditions. Future work will focus on patient-specific parameter identification, multi-objective control strategies, incorporation of real clinical datasets, extension to other neurodegenerative and epidemiological models, parameter calibration using longitudinal biomarker datasets (PET imaging, CSF and plasma amyloid/tau measures) to validate fractional dynamics and reduce uncertainty.

Correction note

A correction has been made to this article. Details can be found at: 10.3389/fninf.2026.1821637.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Miodrag Zivkovic, Singidunum University, Serbia

Reviewed by: Zhandos Zulpykhar, L.N. Gumilyov Eurasian National University, Kazakhstan

Monica Wadie, Delta University for Science and Technology, Egypt