A Fractional-Order Machine Learning Framework for Modeling Vertebral Column Pathology and Biomechanical Dynamics

1  Introduction

The human spinal column functions as the core structure of the skeleton while providing movement abilities and stability and safeguarding the spinal cord through its network of 33 vertebrae and intervertebral discs, and ligaments [1–5]. The natural balance between vertebral bones and discs becomes disrupted by spinal conditions like disk hernia and spondylolisthesis, which lead to persistent pain and functional decline that affects more than 80% of adults during their lifetime [6–8]. The conditions produce changes in biomechanical measurements, which include pelvic incidence and lumbar lordosis angle and sacral slope, pelvic tilt, and pelvic radius, because these measurements show how the spine is aligned and how weight is spread out [9–12]. The diagnosis process needs immediate treatment, but current methods fail to detect the intricate and memory-based spinal biomechanics. Conventional diagnostic methods depend on imaging techniques, which include X-rays and CT scans, and Magnetic Resonance Imaging (MRI), together with statistical models to detect structural defects [13–15]. Research shows that pelvic incidence determines spinal alignment through its direct relationship with lumbar lordosis and sacral slope measurements in people who have normal spinal structure [16–18]. Authors of [19] developed and evaluated an artificial intelligence algorithm to automate vertebral column segmentation in biplanar full-body radiographic images for degenerative scoliosis (DS) assessment. The research used 250 anonymous high-definition X-ray images from an institutional Picture Archiving and Communication System (PACS) to train and test a two-stage deep learning model based on U-Net (UNET) architecture, which included 200 images for training and 50 images for evaluation. The model first identified the spine region in full-body X-rays and then isolated spinal curvature, achieving high accuracy with Dice-Sørensen coefficients of 0.92 for anterior-posterior and 0.96 for lateral views, even in cases with complex spinal pathologies like lordosis, scoliosis, and spinal instrumentation. The researchers used a cross-sectional study [20] to study how spinopelvic alignment varies while creating a spinal balance classification system and developing formulas to predict lumbar lordosis (LL) based on pelvic incidence (PI). Sagittal parameters were evaluated using radiographic assessments, and participants were grouped into three distinct clusters with varying PI and LL averages using K-means clustering and a decision tree incorporating PI and sacral slope (SS). The prediction of LL and SS through cluster-specific linear regression models produced moderate correlation results. The authors of [21] studied how pelvic biomechanics, including pelvic incidence (PI) and pelvic tilt (PT) and lumbar lordosis (LL) and sacral slope (SS), and pelvic radius (PR), might predict spinal disorders like scoliosis and spondylolisthesis. Data from a centralized orthopedic database were used, with spine conditions classified as normal or abnormal based on imaging and clinical exams. Descriptive statistics, comparative analyses, decision trees, and logistic regression identified PI, LL, SS, and PR as significant predictors, with decision trees accurately classifying 69.5% of cases based on PI thresholds and models validated via ROC analysis and 10-fold cross-validation. The research identified abnormal sagittal spinopelvic alignment as the main cause of improper spinal load distribution while showing that higher Body Mass Index (BMI) increases the chance of developing degenerative central stenosis. The research shows how pelvic parameters lead to spinal problems, which creates new possibilities for better diagnostic methods and individualized treatment plans, yet additional forward-looking research needs to enhance current prediction systems and study movement-based evaluations. The methods described above treat biomechanical features as unchanging elements because they fail to consider the viscoelastic nature of spinal tissues, which results in previous stresses affecting present behavior. The first biomechanical studies used integer-order differential equation models, which failed to consider the delayed responses of intervertebral discs and ligaments because these models do not include memory effects, as fractional-order models do [22–24]. The method shows its greatest weakness when modeling spondylolisthesis because it fails to capture the body’s natural compensatory response, which results in increased lumbar lordosis.

Mathematical modeling has been instrumental in analyzing a wide range of phenomena with practical applications [25,26], where fractional-order calculus provides an effective way to solve these problems. The Caputo fractional derivative enables non-local memory-dependent system behavior, which scientists demonstrate through their work in biomechanical systems [27,28]. Fractional-order differential equations (FDEs) have been used to model viscoelastic tissues in other contexts, such as cartilage [29,30] and several other biological phenomena [31–33], but their application to spinal biomechanics remains underexplored. Medical diagnostics underwent a complete transformation through the use of machine learning (ML), which emerged as a leading innovation in healthcare [34,35]. Random Forest, Gradient Boost, XGBoost, and Deep Neural Networks have shown great promise for binary classification of abnormalities [36–38]. These models function well to detect patterns in large datasets, yet they struggle to provide explanations based on biomechanical principles. Studies about medical diagnostics show that Voting Ensembles enhance system stability through their use of ensemble methods [39–41]. Standalone ML methods do not include temporal information or tissue memory, which prevents them from effectively modeling disease progression.

Recent studies show that physics-informed neural networks, which combine ML with dynamic modeling through data-driven predictions and physical constraints, show promising results [42–44]. Yet, no prior work has fully bridged ML with fractional-order modeling for vertebral column pathology. The UCI vertebral column dataset contains 310 cases with six biomechanical features, which makes it an ideal resource for studying linked anatomical patterns between pelvic incidence and pelvic tilt and sacral slope as reported in clinical research.

The research presents an innovative fractional-order machine learning system that demonstrates the ability to model spinal conditions and biomechanical operations of the vertebral column. The system will use advanced ML algorithms to solve five linked Caputo FDEs, which will model tissue memory behavior and forecast abnormality probabilities. The idea of integrating ML algorithms with FDE models has recently been applied to some other biological phenomena and has proven quite positive [45–47]. The current framework processes the UCI dataset through Synthetic Minority Over-sampling Technique (SMOTE) oversampling and Z-score normalization, and Minimum Redundancy Maximum Relevance (MRMR) feature selection to identify pelvic incidence as a key feature. A diagnostic GUI system enables doctors to apply their findings in real-world clinical situations. The research uses an interdisciplinary approach that builds on existing studies to create a dynamic and interpretable system for spinal pathology diagnosis while showing potential for improved musculoskeletal health precision medicine. The rest of the paper is organized as follows: Section 2 presents the methodology that integrates the ML algorithm and the FDE modeling. Section 3 summarizes the data and data processing technique used. The ML algorithms are detailed in Section 4. The developed FDE model, informed by the ML best-performing model, is presented in Section 5, and its theoretical analysis is rigorously established in Section 6. The numerical analysis is detailed in Section 7, while the entire results and conclusion are presented in Sections 8 and 9, respectively.

2  Methodology

The research methodology combines machine learning techniques with fractional-order differential equations to analyze vertebral column abnormalities through the UCI vertebral column dataset, which contains 310 samples with six biomechanical measurements and two class labels that include 210 abnormal and 100 normal cases. The dataset underwent preprocessing through SMOTE oversampling, which created 1000 instances to address class imbalance, and then Z-score normalization was used for standardization. The MRMR feature selection process identified five important features, which include pelvic incidence and lumbar lordosis angle, pelvic radius and sacral slope, and pelvic tilt. The fractional-order model received input features that were min-max scaled between 0 and 1. The evaluation process involved training multiple machine learning models which included Random Forest with 500 trees minimizing Gini impurity and Gradient Boost with 300 trees at a learning rate of 0.02 minimizing exponential loss and XGBoost Approximate with 300 trees at the same learning rate minimizing logistic loss and DNN with cross-entropy loss and Rectified Linear Unit (ReLU) activation and 128/64 neuron layers and batch normalization and dropout rates of 0.3 and 0.2 optimized through Adam and Voting Ensemble that combined these models with logistic regression using majority voting. The models underwent evaluation through 10-fold cross-validation, which measured their performance by accuracy, precision, recall, F1 score, AUC, and computation time. The model consists of five interconnected Caputo fractional differential equations, which operate at an order of 0.8 to represent normalized biomechanical features, pelvic incidence, pelvic tilt, lumbar lordosis angle, sacral slope, and pelvic radius. Their correlations directly influence the fractional-order growth rates in the model, with strongly linked parameters exhibiting accelerated evolution under stress, reflecting how pathological changes in one propagate to others via viscoelastic tissue responses. The model uses logistic growth functions to simulate bounded systems and linear relationships to represent mechanical connections between body parts. The model uses logistic growth functions to represent bounded systems and linear relationships between body parts to simulate their mechanical interactions. The sum of pelvic tilt and sacral slope creates the approximation for pelvic incidence. The model uses pseudo-time, which orders patients based on their pelvic incidence measurements to show how their disease progresses. The Caputo fractional derivative models viscoelastic memory effects in spinal tissues based on the best-performing ML model shown in Fig. 1. The system was solved numerically using the Adams-Bashforth-Moulton predictor-corrector method on a uniform grid with a step size of 0.01 and a final time of 10, and the solution process used Fast Fourier Transform for discrete convolutions and an explicit method based on Grünwald-Letnikov for stability, with solutions restricted to the range of [0, 1] to maintain physical accuracy. The parameter optimization process took place through MATLAB’s fmincon using Sequential Quadratic Programming, which minimized a weighted mean squared error cost function based on MRMR-derived weights and regularization and coupling strength constraints while maintaining growth rates between 0.05 and 1.5 and coupling coefficients between −1 and 1, and fractional order between 0.2 and 0.8. The Gradient Boost machine learning model produces the abnormality probability, which gets interpolated through a piecewise cubic Hermite polynomial across the pseudo-time axis. The resulting values from this interpolation serve to determine how pathological factors impact pelvic tilt and lumbar lordosis, and pelvic radius within the fractional-order model. The model predicts how abnormal conditions affect the alignment of pelvic tilt and lumbar lordosis and pelvic radius, which results in compensatory biomechanical responses. The fractional model solutions were matched to normalized UCI dataset values through the use of pseudo-time, which organizes pelvic incidence to represent disease progression. Theoretical proofs of solution existence, uniqueness, boundedness, and stability were established using the Banach fixed-point theorem and Lyapunov functions, accommodating negative coupling coefficients for inhibitory biomechanical effects. These negative coefficients physiologically represent the body’s compensatory responses, such as flattening of the lumbar curvature to redistribute spinal load in response to pathology, thereby reducing stress on affected vertebrae and preventing further degeneration.

Figure 1: Schematic of the hybrid ML-FDE framework, illustrating bidirectional data flow between diagnostic features, growth trajectory, and predictive dynamics.

3  Data Summary

Table 1 presents a detailed overview of the dataset attributes, which serve as the foundation for this investigation. The table presents both the range values and mean ± standard deviation for all features, which are divided into Abnormal and Normal classes. The table also shows the normalized values through min-max scaling between 0 and 1 for the fractional-order model and machine learning analyses.

Data Processing

The dataset consists of data from column_2C.dat, which provides six features together with one class label. The class labels undergo a binary conversion process, which assigns Abnormal to 1 and Normal to 0. SMOTE operates on 67.7% of the Abnormal cases to create a balanced dataset, which results in 1000 instances. Z-score normalization standardizes features, and MRMR selects the top five features. The processed data serves as the basis for training and evaluation of ML models through 10-fold cross-validation.

4  Machine Learning Algorithms

Model Architectures

The machine learning models for the vertebral column diagnosis are defined by their mathematical formulations, primarily focusing on their objective functions and key components. Their core equations are summarized as follows:

{RF: ∑k=1Kpk(1−pk),pk = prob. of class k at tree node, 500 treesGB: ∑i=1nexp⁡(−yif(xi)),yi∈{−1,1} = class label, f(xi) = ensemble score, η=0.02, 300 treesXGB: ∑i=1nlog⁡(1+exp⁡(−yif(xi))),yi∈{−1,1} = class label, f(xi) = ensemble score, η=0.02, 300 treesDNN: −∑i=1n∑k=12yiklog⁡(y^ik),yik∈{0,1} = true class k,y^ik = predicted prob., ReLU σ(z)=max(0,z)VE: ∑i=1n[−yilog⁡(y^i)−(1−yi)log⁡(1−y^i)],yi∈{0,1} = true label, y^i=11+e−wTxi = predicted prob.(1)

Random Forest (RF) aggregates 500 decision trees, each minimizing Gini impurity, where pk is the probability of class k (Normal or Abnormal) at a tree node. Gradient Boost (GB) and XGBoost Approx (XGB) minimize exponential and logistic losses, respectively, with 300 trees and learning rate η=0.02, where yi is the true label (−1 = Normal, 1 = Abnormal) and f(xi) is the ensemble’s weighted prediction score for instance xi. The Deep Neural Network (DNN) uses cross-entropy loss, with yik as the true class indicator (1 if instance i is class k, else 0), y^ik as the predicted probability, ReLU activation σ(z)=max(0,z), softmax output, 128/64 neuron layers, batch normalization, and dropout (0.3, 0.2), optimized via Adam. The Voting Ensemble (VE) combines RF, GB, XGB, and logistic regression (log-reg), where the log-reg component minimizes log-loss with yi as the true label (0 = Normal, 1 = Abnormal) and y^i as the predicted probability via the sigmoid function, using majority voting for final predictions.

5  Fractional-Order Modeling

The fractional-order model emerged to understand spinal biomechanics through its application to vertebral column disorders like disk hernia and spondylolisthesis by using data from the UCI vertebral column dataset. The model operates through five interconnected fractional-order differential equations, which simulate how key biomechanical elements change over time while representing biological tissue memory effects through fractional derivatives. The Gradient Boost machine learning model was trained on these features to predict the abnormality probability P(t), which informs the coupling terms in the fractional-order system as presented in Section 4. The model employs the Caputo fractional derivative to simulate viscoelastic behavior and long-term memory effects of spinal tissue by creating pseudo-time through patient sorting based on pelvic incidence, which mirrors anatomical development and disease progression. The equations were developed to incorporate logistic growth terms, which control feature limits and linear coupling terms that simulate biomechanical relationships between pelvic incidence and pelvic tilt and sacral slope. The optimization process reached its best results by using MATLAB’s fmincon function to minimize weighted mean squared error, demonstrating an excellent match between the model and normalized feature data.

The system is given by:

{Dtα0cP1(t)=r1P1(t)(1−P1(t))+k1(P2(t)+P4(t)),Dtα0cP2(t)=r2P2(t)(1−P2(t))+k2P(t)+k6P3(t),Dtα0cP3(t)=r3P3(t)(1−P3(t))+k3P(t)+k7P4(t),Dtα0cP4(t)=r4P4(t)(1−P4(t))+k4P1(t),Dtα0cP5(t)=r5P5(t)(1−P5(t))+k5P(t),(2)

where Dtα0c(.) denotes the Caputo fractional derivative of order α∈(0,1).

The Caputo derivative is defined as [49,50]:

Dtα0cy(t)=1Γ(1−α)∫0t(t−s)−αy′(s)ds,(3)

where Γ(z) is the Gamma function defined as

Γ(z)=∫0∞uz−1e−udu,Re(z)>0.

The fractional-order model enables simulation of spinal biomechanical operations at a level that matches biological precision. The Caputo fractional derivative of order α shows how viscoelastic spinal tissues, including intervertebral discs and ligaments, maintain memory of previous stresses that affect their current alignment. The logistic terms riPi(1−Pi) function to maintain system stability by creating limits that duplicate the natural limits of biomechanical angles and distances. The mathematical expression k1(P2+P4) represents the connection between pelvic incidence and the sum of pelvic tilt and sacral slope, while k3 and k7 demonstrate how lumbar curvature changes because of abnormalities and sacral orientation, which patients with spondylolisthesis develop as compensatory mechanisms. The negative values of k3 and k7 represent inhibitory effects that decrease lumbar curvature because of pathological conditions and sacral alignment stabilization feedback, thus creating a realistic simulation of spinal movement.

6  Computational Analysis

In this section, we present the theoretical computations (existence, uniqueness, and stability analysis) to validate the fractional-order model.

Theorem 1: (Existence and uniqueness of solution): Consider the system of Caputo fractional-order differential equations of order α∈(0,1) given in system (2) with initial conditions Pi(0)=Pi0∈[0,1] for i=1,2,3,4,5, where ri>0, kj∈R (allowing negative values), and P(t)∈C([0,T],[0,1]) is a known continuous function. There exists a unique continuous solution P(t)=(P1(t),P2(t),P3(t),P4(t),P5(t))T on [0,T] for some T>0, and the solution remains bounded. The existence and uniqueness hold regardless of the signs of kj.

Proof: We prove the theorem using the Banach fixed-point theorem on the equivalent Volterra integral system. All computations are provided explicitly using inequalities and mathematical expressions.

Let P(t)=(P1(t),P2(t),P3(t),P4(t),P5(t))T∈R5, with P(0)=P0=(P10,P20,P30,P40,P50)T∈[0,1]5. The system is:

Dtα0cP(t)=f(t,P(t),P(t)),

where f=(f1,f2,f3,f4,f5)T with fi defined by the right hand side in system (2). Since P(t)∈C([0,T], [0,1]), we have:

‖P‖C([0,T])=supt∈[0,T]|P(t)|≤MP=1.

Applying the Riemann–Liouville integral operator Itα0RLy(t)=1Γ(α)∫0t(t−s)α−1y(s)ds to system (2) leads to an integral form:

P(t)=P0+1Γ(α)∫0t(t−s)α−1f(s,P(s),P(s))ds.

Consider the Banach space C([0,h],R5) for h>0 small, with norm:

‖P‖=supt∈[0,h]‖P(t)‖∞=supt∈[0,h]maxi=1,…,5|Pi(t)|.

Define the closed ball:

BR={P∈C([0,h],R5):‖P−P0‖≤R},

so for P∈BR, we have:

‖P(t)‖∞≤‖P(t)−P0‖+‖P0‖≤R+1.

Define the operator T:C([0,h],R5)→C([0,h],R5):

(TP)(t)=P0+1Γ(α)∫0t(t−s)α−1f(s,P(s),P(s))ds.

We show T maps BR to BR and is a contraction.

For P∈BR, compute:

‖(TP)(t)−P0‖∞=supt∈[0,h]‖1Γ(α)∫0t(t−s)α−1f(s,P(s),P(s))ds‖∞≤1Γ(α)supt∈[0,h]∫0t(t−s)α−1‖f(s,P(s),P(s))‖∞ds.

Bound ‖f(s,P,P)‖∞=maxi|fi|. For each fi on P∈[0,1+R]5, P∈[0,1]:

For f1:

|f1|=|r1P1(1−P1)+k1(P2+P4)|≤r1|P1(1−P1)|+|k1|(|P2|+|P4|).

The logistic term satisfies:

g(u)=u(1−u),g′(u)=1−2u,g′′(u)=−2<0,

so g(u) is maximized at u=0.5 with g(0.5)=0.25. For u∈[0,1+R], g(u)≤0.25, since g(1+R)=1+R−(1+R)2=−R−R2≤0. Thus:

|P1(1−P1)|≤14.

Since |P2|,|P4|≤1+R, we have:

|f1|≤r1⋅14+|k1|(1+R+1+R)=r14+2|k1|(1+R).

For f2:

|f2|≤r2|P2(1−P2)|+|k2||P|+|k6||P3|≤r24+|k2|⋅1+|k6|(1+R).

For f3:

|f3|≤r34+|k3|⋅1+|k7|(1+R).

For f4:

|f4|≤r44+|k4|(1+R).

For f5:

|f5|≤r54+|k5|⋅1.

Let rmax=maxiri, kmax=maxj|kj|. Then:

‖f‖∞≤max{r14+2|k1|(1+R),r24+|k2|+|k6|(1+R),r34+|k3|+|k7|(1+R),r44+|k4|(1+R),r54+|k5|}≤rmax4+3kmax(1+R)=:M.

Evaluate the integral:

∫0t(t−s)α−1ds=∫0tuα−1du=[uαα]0t=tαα,

so:

‖TP−P0‖≤MΓ(α)⋅hαα=MhαΓ(α)⋅α=MhαΓ(α+1).

Choose h>0 such that:

MhαΓ(α+1)≤R ⟹ h≤(RΓ(α+1)M)1/α.

Since M is finite (as ri,kj are fixed), for any R>0, there exists h>0 small enough, so T:BR→BR.

Step 2: T is a contraction.

For P,Q∈BR, compute:

‖(TP)(t)−(TQ)(t)‖∞≤1Γ(α)supt∈[0,h]∫0t(t−s)α−1‖f(s,P(s),P(s))−f(s,Q(s),P(s))‖∞ds.

Show f is Lipschitz in P. Compute for each component:

For f1:

f1(P)−f1(Q)=r1[P1(1−P1)−Q1(1−Q1)]+k1[(P2+P4)−(Q2+Q4)]=r1(P1−Q1−P12+Q12)+k1(P2−Q2+P4−Q4).

Since P12−Q12=(P1−Q1)(P1+Q1), we have:

P1−P12−Q1+Q12=(P1−Q1)−(P12−Q12)=(P1−Q1)(1−P1−Q1).

Thus:

|f1(P)−f1(Q)|≤r1|P1−Q1||1−P1−Q1|+|k1|(|P2−Q2|+|P4−Q4|).

Since |P1|,|Q1|≤1+R, |1−P1−Q1|≤1+2(1+R)=3+2R, so:

|f1(P)−f1(Q)|≤r1(3+2R)|P1−Q1|+|k1|(|P2−Q2|+|P4−Q4|)≤[r1(3+2R)+2|k1|]‖P−Q‖∞.

For f2:

f2(P)−f2(Q)=r2(P2−Q2)(1−P2−Q2)+k6(P3−Q3),|f2|≤r2(3+2R)|P2−Q2|+|k6||P3−Q3|≤[r2(3+2R)+|k6|]‖P−Q‖∞.

For f3:

|f3|≤[r3(3+2R)+|k7|]‖P−Q‖∞.

For f4:

|f4|≤[r4(3+2R)+|k4|]‖P−Q‖∞.

For f5:

|f5|≤r5(3+2R)‖P−Q‖∞.

Thus:

‖f(P,P)−f(Q,P)‖∞≤L‖P−Q‖∞,

where:

L=max{r1(3+2R)+2|k1|,r2(3+2R)+|k6|,r3(3+2R)+|k7|,r4(3+2R)+|k4|,r5(3+2R)}.

Note that L depends only on |kj|, so negative kj do not affect the Lipschitz condition. Then:

‖TP−TQ‖≤LΓ(α)supt∈[0,h]∫0t(t−s)α−1‖P(s)−Q(s)‖∞ds≤LΓ(α)⋅hαα⋅‖P−Q‖=LhαΓ(α+1)‖P−Q‖.

Choose h such that:

K=LhαΓ(α+1)<1 ⟹ h<(Γ(α+1)L)1/α.

Since L is finite, such h exists. Thus, T is a contraction on BR.

Now, applying the Banach Fixed-Point theorem, since T maps BR to BR and is a contraction, there exists a unique P∈BR such that TP=P, i.e., a unique solution on [0,h].

For global existence and boundedness, assume the solution exists on [0,Tmax) with Tmax≤T. If Tmax <T, then ‖P(t)‖∞→∞ as t→Tmax−. But:

|0cDtαPi|≤ri4+3kmax(1+R).

By the integral equation, for t∈[0,h]:

|Pi(t)|≤|Pi0|+1Γ(α)∫0t(t−s)α−1(ri4+3kmax(1+R))ds≤1+(rmax4+3kmax(1+R))hαΓ(α+1).

Choose R and h to be small so that the solution remains bounded. Negative kj may reduce |fi|, aiding boundedness. By continuation, the solution extends to [0,T], as the logistic terms prevent blow-up. Thus, the solution exists, is unique, and is bounded on [0,T], with negative kj permissible. ◻

Theorem 2: (Stability of the Fractional-Order Model): Consider the system of Caputo fractional-order differential equations of order α∈(0,1) given by (2), with initial conditions Pi(0)=Pi0∈[0,1] for i=1,2,3,4,5, where ri>0, kj∈R, and P(t)∈C([0,T],[0,1]) is a known continuous function representing the abnormality probability. For sufficiently small |kj|, the solution P(t)=(P1(t),P2(t),P3(t),P4(t),P5(t))T is uniformly ultimately bounded in a compact set containing [0,1]5, and any equilibrium points in [0,1]5 are locally asymptotically stable under perturbations, with explicit bounds ensuring stability regardless of the signs of kj.

Proof: We establish uniform ultimate boundedness and local asymptotic stability using a Lyapunov function, providing rigorous computations for the system in (2). Using the Caputo definition in (3), the system is:

Dtα0cP(t)=f(t,P(t),P(t)),

where f=(f1,f2,f3,f4,f5)T is defined by (2), and supt∈[0,T]|P(t)|≤1.

Uniform Ultimate Boundedness

Define the Lyapunov function:

V(P)=∑i=1512riPi2,

where ri>0 are the growth rates from (2). Since V(P)=∑i=1512riPi2≥0 and V(P)=0 only at P=0, it is positive definite. Compute the fractional derivative:

Dtα0cV(P(t))≤∑i=15∂V∂Pi⋅Dtα0cPi(t)=∑i=15Pirifi(t,P,P).

Substitute fi from (2):

∑i=15Pirifi=P1r1[r1P1(1−P1)+k1(P2+P4)]+P2r2[r2P2(1−P2)+k2P+k6P3]+P3r3[r3P3(1−P3)+k3P+k7P4]+P4r4[r4P4(1−P4)+k4P1]+P5r5[r5P5(1−P5)+k5P].

Simplify:

=∑i=15Pi2(1−Pi)+k1r1P1(P2+P4)+k2r2PP2+k6r2P2P3+k3r3PP3+k7r3P3P4+k4r4P1P4+k5r5PP5.

Bound the logistic terms:

Pi2(1−Pi)=Pi2−Pi3.

Since Pi∈[0,1], the function g(Pi)=Pi2−Pi3=Pi2(1−Pi) has derivative g′(Pi)=2Pi−3Pi2, with critical points at Pi=0 or Pi=23. Evaluate: g(0)=0, g(23)=49−827=427≈0.148, g(1)=0. Thus:

Pi2(1−Pi)≤427.

However, for large Pi, consider the behavior:

Pi2(1−Pi)=Pi2−Pi3≤−12Pi2,

for Pi≥1, since −Pi3≤−Pi2. Thus:

∑i=15Pi2(1−Pi)≤∑i=15(−12Pi2)=−12‖P‖22,forPi≥1.

For coupling terms, use |Pi|≤‖P‖2, |P|≤1:

|k1r1P1(P2+P4)|≤|k1|r1|P1|(|P2|+|P4|)≤|k1|r1‖P‖2⋅2‖P‖2=2|k1|r1‖P‖22,|k2r2PP2|≤|k2|r2⋅1⋅‖P‖2,|k6r2P2P3|≤|k6|r2‖P‖22,|k3r3PP3|≤|k3|r3‖P‖2,|k7r3P3P4|≤|k7|r3‖P‖22,|k4r4P1P4|≤|k4|r4‖P‖22,|k5r5PP5|≤|k5|r5‖P‖2.

Let rmin=miniri, kmax=maxj|kj|. Sum the quadratic terms:

|k1r1P1(P2+P4)+k6r2P2P3+k7r3P3P4+k4r4P1P4|≤(2|k1|r1+|k6|r2+|k7|r3+|k4|r4)‖P‖22.

Since ri≥rmin, bound:

2|k1|r1+|k6|r2+|k7|r3+|k4|r4≤2|k1|+|k6|+|k7|+|k4|rmin≤2kmax+3kmaxrmin=5kmaxrmin.

Sum the linear terms:

|k2r2PP2+k3r3PP3+k5r5PP5|≤(|k2|r2+|k3|r3+|k5|r5)‖P‖2≤|k2|+|k3|+|k5|rmin≤3kmaxrmin‖P‖2.

Thus:

Dtα0cV≤−12‖P‖22+5kmaxrmin‖P‖22+3kmaxrmin‖P‖2=(−12+5kmaxrmin)‖P‖22+3kmaxrmin‖P‖2.

Define a=12−5kmaxrmin, b=3kmaxrmin. For small kmax, ensure a>0:

5kmaxrmin<12 ⟹ kmax<rmin10.

Using optimized values rmin=0.05, kmax=0.1382, check:

5⋅0.13820.05=13.82>12,

so adjust by assuming smaller |kj|. Try kmax<0.0510=0.005:

a=12−5⋅0.0050.05=12−0.5=0,

so kmax<0.005 is too restrictive. So we use V=∑i=15Pi2:

Dtα0cV≤∑i=152Pifi=2∑i=15riPi2(1−Pi)+2∑coupling terms.

Bound:

2riPi2(1−Pi)≤2ri⋅427,for Pi∈[0,1],2riPi2(1−Pi)≤2ri(−Pi2)=−2riPi2,for Pi≥1.

So:

∑i=152riPi2(1−Pi)≤−2rmin‖P‖22,for ‖P‖2≥1.

Coupling terms:

2|k1P1(P2+P4)+⋯+k5PP5|≤2(2|k1|+|k6|+|k7|+|k4|)‖P‖22+2(|k2|+|k3|+|k5|)‖P‖2≤2⋅5kmax‖P‖22+2⋅3kmax‖P‖2.

Thus:

Dtα0cV≤−2rmin‖P‖22+10kmax‖P‖22+6kmax‖P‖2=(−2rmin+10kmax)‖P‖22+6kmax‖P‖2.

Let a=2rmin−10kmax, b=6kmax. Require a>0:

kmax<rmin5=0.055=0.01.

For kmax=0.01:

a=2⋅0.05−10⋅0.01=0.1−0.1=0,

so assume kmax<0.01. Then:

Dtα0cV≤−a‖P‖22+b‖P‖2.

For u=‖P‖2:

−au2+bu=−a(u2−bau)=−a((u−b2a)2−b24a2).

Choose R=ba=6kmax2rmin−10kmax. For u≥R:

Dtα0cV≤−a((u−b2a)2−b24a2)≤−a⋅b24a2=−b24a=−(6kmax)24(2rmin−10kmax)=−9kmax22rmin−10kmax.

For kmax=0.005, rmin=0.05:

9⋅0.00522⋅0.05−10⋅0.005=9⋅0.0000250.1−0.05=0.0002250.05=0.0045.

Thus:

Dtα0cV≤−0.0045,for ‖P‖2≥R≈6⋅0.0052⋅0.05−10⋅0.005=0.030.05=0.6.

Since V=‖P‖22, V≥0.36 implies V≥0.6. By the fractional comparison principle, V(t) is bounded, ensuring P(t) is uniformly ultimately bounded in {P:‖P‖2≤R+ϵ}, containing [0,1]5. Negative kj are permissible, as bounds use |kj|.

Local Asymptotic Stability

Equilibria satisfy f(t,P∗,P∗)=0 with P∗ constant. For small |kj|, approximate Pi∗≈0.5. The Jacobian is:

J=(−r1k10k100−r2k60000−r3k70k400−r400000−r5).

Characteristic polynomial:

det(λI−J)=(λ+r5)(λ+r2)det(λ+r1−k1−k10λ+r3−k7−k40λ+r4).

Compute:

det(λ+r1−k1−k10λ+r3−k7−k40λ+r4)=(λ+r1)(λ+r3)(λ+r4)+k1k7k4−k1k4(λ+r3).

For kj=0, eigenvalues are −ri. For kmax=0.005, perturbation |k1k7k4|≤0.0053=0.000000125. Gershgorin discs ensure Re(λi)<0 for small kmax, confirming local asymptotic stability by Matignon’s theorem [51,52]. Negative kj does not affect stability, as |kj| governs the perturbation.

Thus, solutions are uniformly ultimately bounded, and equilibria are locally asymptotically stable for small |kj|. ◻

7  Numerical Analysis

The Adams-Bashforth-Moulton predictor-corrector numerical method solves the fractional-order model for vertebral column data by applying the Caputo fractional derivative Dtα0c(.) with α between 0 and 1 at the optimized value of 0.8. The model tracks five normalized biomechanical features Y(t)=(P1(t),P2(t),P3(t),P4(t),P5(t))T which include pelvic incidence and pelvic tilt and lumbar lordosis angle and sacral slope and pelvic radius all measured within specific ranges and normalized to values between 0 and 1. The system operates under the following definition: