A Hybrid Machine Learning and Fractional-Order Dynamical Framework for Multi-Scale Prediction of Breast Cancer Progression

1  Introduction

Breast cancer continues to be a major global health challenge, with complex biological behaviors that make accurate prediction of disease progression particularly difficult [1–5]. The disease’s heterogeneity manifests not only in its molecular subtypes but also in its highly variable growth patterns and treatment responses [6,7]. Current diagnostic methods, while effective for detection, often struggle to provide comprehensive insights into the dynamic nature of tumor evolution [8,9]. This limitation has spurred significant interest in developing more sophisticated analytical approaches that can bridge the gap between immediate diagnostic needs and long-term prognostic predictions, particularly through the integration of computational methods with biological insights. The emergence of machine learning in medical diagnostics has revolutionized the approach to cancer classification [10–12]. These data-driven methods excel at identifying complex patterns in high-dimensional datasets, enabling highly accurate differentiation between malignant and benign lesions. Various algorithms, from traditional logistic regression to advanced deep neural networks, have demonstrated impressive performance in breast cancer diagnosis [13,14]. For instance, the study in [15] presented an AI-assisted deep learning model for early melanoma detection that achieved high accuracy while demonstrating potential for clinical implementation. However, while these models provide excellent snapshots of disease state at a given moment, they typically cannot model the temporal dynamics of tumor progression. This represents a significant limitation, as understanding how a tumor evolves is crucial for predicting its future behavior and selecting optimal treatment strategies. The static nature of conventional machine learning approaches contrasts with the dynamic reality of cancer biology, where past states continuously influence future progression. Mathematical modeling has been instrumental in analyzing a wide range of phenomena and practical applications [16–19]. They also offer complementary strengths to machine learning by providing frameworks to simulate biological processes over time [20–22]. Traditional differential equation models have been widely used to describe tumor growth dynamics [23]; however, they often rely on simplifying assumptions that may not fully capture the complexity of real-world cancer progression. This is where fractional calculus presents an intriguing alternative. Unlike classical calculus, fractional-order models incorporate memory effects and long-range dependencies, which align well with the observed behavior of biological systems [24–26], and have shown value across multiple disciplines [27–29]. The ability to account for historical influences on current growth patterns makes fractional calculus particularly suitable for modeling cancer progression, where a tumor’s past states significantly impact its future development. Recent applications of fractional-order models in biology and medicine have shown promising results, though their integration with data-driven approaches remains underexplored [30,31]. Recent advances have demonstrated the potential of combining these approaches. The authors in [32] developed a novel approach using a fractional recurrent neural network (RNN) to achieve observer-based synchronization in a fractional-order chaotic cancer cell model, while research [33] proposed a mathematical model integrated with machine learning techniques that achieved 97.65% accuracy in breast cancer prediction. Despite these advances, a critical gap remains in developing a unified framework that combines machine learning’s diagnostic precision with fractional calculus’s ability to model progression dynamics while maintaining clinical interpretability and utility.

Our work addresses this gap through a novel hybrid framework that begins with comprehensive machine learning analysis to identify the most predictive diagnostic features, then constructs a fractional-order dynamical system that captures the essential interactions governing tumor progression. This approach allows us to not only classify tumors with high accuracy but also to simulate their likely future behavior under different conditions. A key innovation is our use of the fractional order parameter as a potential biomarker for tumor aggressiveness, where different values correspond to distinct growth patterns observed clinically. The theoretical underpinnings of our model are rigorously developed, with careful attention to ensuring biological plausibility and clinical relevance, including proofs of existence and uniqueness of solutions. Practical implementation is achieved through a user-friendly graphical interface that integrates both classification capabilities and predictive modeling, representing a crucial step in translating complex computational methods into clinically actionable information. The broader implications of this work extend beyond breast cancer, offering a new paradigm for disease modeling that combines data-driven pattern recognition with mechanistic understanding. While advanced computational frameworks like ours aim to improve diagnostic precision, patient attitudes toward cancer management, including complementary and alternative medicine (CAM), remain critical to holistic care. Studies in low-resource settings, such as in [34], exploration of CAM perspectives among Pakistani cancer patients, reveal a gap between technological innovations and patient-centric adoption. Their findings underscore the need for culturally adaptable tools that bridge predictive analytics (like our hybrid ML-fractional model) with patient education, particularly in regions where breast cancer stigma or reliance on CAM may delay evidence-based treatment. By bridging these approaches, we move closer to truly personalized medicine, where treatment strategies can be tailored not just to a patient’s current disease state but to their predicted future progression. The rest of the paper is structured as follows: Section 2 introduces the methodology used in the study, providing the technical foundation for these advances. Section 3 details the various machine learning algorithms employed. In Section 4, the fractional-order modeling is formulated and extensively analyzed with analytic and numerical simulations. The results and implications are discussed in Section 5, while the conclusion, novelty, limitations, and future directions are given in Section 6, demonstrating the potential to transform our approach to breast cancer management.

2  Methodology

To investigate the dynamics of breast cancer progression, we developed an integrated analytical framework that combines machine learning classification with fractional-order dynamical modeling, implemented in MATLAB R2021a. This dual approach aims to leverage machine learning’s strength in data-driven pattern recognition alongside fractional calculus’s ability to capture memory-dependent biological processes, ensuring both predictive accuracy and biological interpretability. The methodology links machine learning outputs, such as feature importance, to the parameterization of mechanistic tumor growth models, creating a cohesive pipeline that bridges data science and mathematical biology. The dataset characteristics are summarized in Table 1, with its feature categories described in Table 2. The entire schematic flow is depicted in Fig. 1, with ultrasound image analysis illustrated in Fig. 2, showing a malignant breast tumor via diagnostic and breast ultrasound imaging, adapted from [35], highlighting irregular mass morphology, spiculated margins, and segmentation. The ML research design is illustrated in Fig. 3, with feature importance ranking and cross-validated performance metrics summarized in Tables 3 and 4, respectively. First, a comprehensive feature analysis using machine learning algorithms with the best performing algorithm, then the development of a fractional-order dynamical system informed by machine learning outputs, and finally, clinical implementation and prediction. The process begins with preprocessing the Wisconsin Diagnostic Breast Cancer dataset to ensure data quality and consistency. Subsequently, parallel workflows develop machine learning models for feature identification and fractional-order systems for tumor growth modeling. Machine learning algorithms identify key diagnostic features, which directly inform the parameterization of the fractional-order tumor growth model. This reciprocal relationship ensures that data-driven insights enhance the mechanistic understanding of cancer progression, while the fractional-order model provides interpretable dynamics grounded in biological principles. The resulting framework synthesizes empirical analysis and theoretical modeling, offering a novel paradigm for multi-scale cancer progression analysis with improved predictive performance and clinical applicability.

Figure 2: Plot detection of a malignant breast tumor in diagnostic imaging, using an ultrasound image of a malignant tumor from the real dataset [35]. Parameters for segmentation: threshold = 0.7, Sobel edge detection

Figure 3: Machine learning breast cancer diagnostic pipeline

2.1 Dataset Characteristics

•   All features are continuous, real-valued numbers

•   Features were computed from digitized images at 8-bit resolution

•   Data was collected at the University of Wisconsin Hospitals

•   Features represent characteristics of cell nuclei present in the images

•   The diagnosis (target) was confirmed by biopsy

2.2 Preprocessing Details

The MATLAB code performs the following preprocessing steps:

•   Conversion of diagnosis labels to binary (M = 1, B = 0)

•   Removal of patient ID column

•   Z-score normalization of all features

•   SMOTE resampling to address class imbalance

•   5-fold stratified cross-validation

3  Mathematical Formulations of the Machine Learning Models

The machine learning pipeline employs seven distinct algorithms to classify breast tumors, each with rigorous mathematical foundations. Logistic regression establishes a probabilistic framework via sigmoid transformations, while decision trees partition the feature space through recursive impurity minimization. Support vector machines construct optimal hyperplanes using kernel-based transformations, and ensemble methods (random forests, XGBoost) combine weak learners through bagging and gradient boosting. Neural networks implement hierarchical feature learning via layered nonlinear transformations. These formulations collectively enable robust prediction while maintaining interpretability through their mathematical structures.

Here are the enhanced subsections with detailed explanations for all components:

3.1 Logistic Regression

The logistic regression model estimates malignancy probability through the following formulation:

P(y=1|x)=σ(wTx+b)=11+e−(wTx+b),(1)

where

•   x∈R30 is the input feature vector containing:

–   10 core features (radius, texture, perimeter, etc.)

–   3 statistical measures per feature (mean, standard error, worst)

•   w∈R30 are the trainable weight parameters where:

–   wj represents feature j’s contribution to malignancy risk

–   Absolute magnitude indicates feature importance

–   Sign determines positive/negative correlation with malignancy

•   b∈R is the bias term encoding baseline log-odds

•   σ(z)=1/(1+e−z) is the sigmoid function with properties:

–   Bounds output between (0,1) for probability interpretation

–   Derivative σ′(z)=σ(z)(1−σ(z)) enables efficient gradient descent

The model is trained via maximum likelihood estimation with:

•   L2 regularization (ridge penalty) to prevent overfitting:

ℒreg=λ‖w‖22(2)

•   Optimized using limited-memory BFGS (L-BFGS) with:

–   Line search for step size determination

–   History size of 50 past updates for Hessian approximation

•   Convergence criteria:

‖∇ℒ‖2<10−5 or |ℒ(t+1)−ℒ(t)|<10−8(3)

The L-BFGS optimization ensures convergence by approximating the Hessian matrix with a limited-memory update, guaranteeing descent directions [38,39].

3.2 Decision Tree

The decision tree implements recursive partitioning through:

f(x)=∑m=1Mcm⋅I(x∈Rm)(4)

where

•   Rm⊂R30 are hyperrectangular regions formed by axis-aligned splits with:

–   Maximum depth of 30 (early stopped via validation)

–   Minimum 5 samples per leaf node

•   cm∈{0,1} is the class prediction for region Rm determined by:

cm=mode(yi|xi∈Rm)(5)

•   I(⋅) is the indicator function evaluating to 1 if x falls in Rm

Splits are determined by Gini impurity minimization:

G(Q)=∑k=12pk(1−pk),pk=1|Q|∑xi∈QI(yi=k)(6)

where

•   Q is the data subset at each node

•   pk is the empirical probability of class k in Q

•   Splits are chosen to maximize impurity reduction:

ΔG=G(Q)−|QL||Q|G(QL)−|QR||Q|G(QR)(7)

Advanced features include:

•   Cost-complexity pruning with complexity parameter α=0.01

•   Class-weighted splitting for imbalanced data:

Gw(Q)=∑k=12wkpk(1−pk)(8)

The tree provides interpretable decision paths:

Path(x)=⋂j=1d{xj≤θj or xj>θj}(9)

where θj are the learned split thresholds for each feature.

3.3 Support Vector Machine (SVM)

The Radial Basis Function (RBF) kernel SVM solves the dual optimization problem:

maxα∑i=1nαi−12∑i,jαiαjyiyjK(xi,xj)s.t.0≤αi≤C,∑iαiyi=0,(10)

where

•   αi∈R+ are Lagrange multipliers encoding support vector importance

•   C>0 is the regularization parameter controlling margin-violation tradeoff

•   yi∈{−1,1} are class labels (malignant/benign)

•   K(⋅,⋅) is the RBF kernel function

•   xi∈R30 are input feature vectors from WDBC dataset

with the RBF kernel:

K(xi,xj)=exp⁡(−γ‖xi−xj‖2),(11)

where

•   γ=12σ2 controls kernel width (inverse squared standard deviation)

•   ‖⋅‖ denotes Euclidean distance between feature vectors

The decision function is:

f(x)=sgn(∑i∈SVαiyiK(xi,x)+b),(12)

where

•   SV denotes the set of support vectors (αi>0)

•   b is the bias term computed from margin constraints

•   sgn(⋅) returns the predicted class (+1 or −1)

3.4 Random Forest

An ensemble of B=100 decision trees with bagging:

f^(x)=majority vote{fb(x)}b=1B,(13)

where

•   B is the number of trees in the ensemble (set to 100 via cross-validation)

•   fb:R30→{0,1} are individual decision tree classifiers

Each tree fb is trained on a bootstrap sample with:

•   Random feature selection at each split (typically 30≈5 features considered)

•   Gini impurity minimization for split criterion:

G(Q)=∑k=12pk(1−pk),pk=1|Q|∑(xi,yi)∈QI(yi=k)(14)

where Q is the data subset at each node and I(⋅) is the indicator function

•   Maximum tree depth determined by early stopping (typically 10–15 levels)

Key properties:

•   Out-of-bag (OOB) samples used for unbiased error estimation

•   Feature importance computed via mean decrease in Gini index

•   Inherent parallelization during training (trees independent)

Here’s the enhanced XGBoost subsection with detailed explanations of all terms and mechanisms:

3.5 XGBoost (Extreme Gradient Boosting)

The gradient boosted trees model minimizes the following regularized objective function:

ℒ(ϕ)=∑i=1nl(yi,y^i)+∑k=1KΩ(fk),(15)

where

l(yi,y^i) is the loss function measuring prediction error:

l(yi,y^i)=−[yilog⁡(σ(y^i))+(1−yi)log⁡(1−σ(y^i))](16)

with σ(⋅) as the logistic sigmoid function for binary classification

Ω(fk) is the regularization term controlling model complexity:

Ω(fk)=γTk+12λ‖wk‖2+α‖wk‖1(17)

where

–   Tk = number of leaves in tree fk

–   wk = vector of leaf weights (output values)

–   γ = minimum loss reduction required for split (controls tree growth)

–   λ = L2 regularization on leaf weights

–   α = L1 regularization on leaf weights (optional)

•   Each tree fk is learned additively to fit the pseudo-residuals:

fk(x)=−[∂l(yi,y^i(k−1))∂y^i(k−1)](18)

where y^i(k−1) is the prediction from previous k−1 trees

•   The prediction after K trees is:

y^i=∑k=1Kfk(xi)(19)

Key optimization features:

•   Tree Construction: Uses weighted quantile sketch for approximate greedy algorithm

•   Split Finding: Evaluates gain score for potential splits:

𝒢=12[(∑i∈ILgi)2∑i∈ILhi+λ+(∑i∈IRgi)2∑i∈IRhi+λ−(∑i∈Igi)2∑i∈Ihi+λ]−γ(20)

–   gi=∂y^(k−1)l(yi,y^(k−1)) (first-order gradient)

–   hi=∂y^(k−1)2l(yi,y^(k−1)) (second-order Hessian)

–   IL,IR = instance sets of left/right nodes

•   Shrinkage: Learning rate η applied to tree weights:

y^i(k)=y^i(k−1)+ηfk(xi)(21)

•   Column Subsampling: Randomly selects feature subsets to prevent overfitting

Implementation details for breast cancer prediction:

•   Early stopping with 50-round patience on validation AUC

•   Maximum tree depth of 6 (optimized via cross-validation)

•   Learning rate η=0.05 with 1000 boosting rounds

•   Minimum child weight (sum of instance weight Hessian) = 5

•   Column subsampling rate = 0.8 per tree

3.6 Neural Networks

The shallow (2 hidden layers) and deep (4 hidden layers) networks implement:

{h1=ReLU(BN(W1x+b1))h2=ReLU(BN(W2h1+b2))(shallow network)hl=ReLU(BN(Wlhl−1+bl))(for l=3 to 5 in deep network)P(y=1|x)=softmax(WLhL−1+bL)(22)

where

•   ReLU(z)=max(0,z) is the Rectified Linear Unit activation function

•   BN denotes batch normalization

•   Wl∈Rnl×nl−1 are weight matrices

•   bl∈Rnl are bias vectors

•   x∈Rd is the input feature vector (30 normalized features)

•   hl represents the hidden state at layer l

•   softmax converts final layer outputs to class probabilities

•   P(y=1|x) is the predicted probability of malignancy

with layer dimensions:

•   Shallow Network (2 hidden layers):

d→128→64→2(23)

–   Layer 1: 128 neurons with dropout (p = 0.2)

–   Layer 2: 64 neurons

–   Output: 2 neurons (softmax classification)

•   Deep Network (4 hidden layers):

d→256→128→64→32→2(24)

–   Layers 1–4: 256, 128, 64, 32 neurons

–   Dropout (p = 0.3) after first three layers, p = 0.2 after fourth

–   Output: 2 neurons (softmax classification)

All models use the categorical cross-entropy loss with L2 regularization:

ℒ=−1n∑i=1n∑c=12yi,clog⁡pi,c+λ∑l=1L‖Wl‖F2(25)

where

•   n is the number of training samples

•   λ=0.001 is the regularization strength

•   Training uses Adam optimizer with default parameters

Key implementation details:

•   Batch normalization after each hidden layer

•   Dropout rates: 20% for shallow, 30% (first three layers) and 20% (last layer) for deep

•   Training for 50 (shallow) and 80 (deep) epochs

•   Mini-batch size of 256 samples

•   Class-weighted sampling to address imbalance

4  Fractional-Order Tumor Dynamics Model

This section introduces a fractional-order differential system to model breast tumor progression. Unlike classical integer-order models, fractional calculus captures memory effects and long-range dependencies in biological systems, making it particularly suitable for cancer growth dynamics. Fractional derivatives capture memory effects via non-local operators, ideal for tumor growth, where past states influence future progression. Caputo’s definition ensures the physical interpretability of initial conditions [40,41].

4.1 Model Equations (Caputo Fractional Derivative)

The fractional-order model in (26) is constructed and informed by the machine learning algorithm and results in Section 3. The optimized parameters are obtained from fitting the model to the dataset. The system describes interactions between three key tumor features identified in the feature importance analysis (Table 3):

{Dtα0cR(t)=a1R(t)(1−R(t)K)−b1R(t)T(t)1+T(t)+c1C(t)1+C(t)(Radius growth)Dtα0cT(t)=a2T(t)(1−T(t)K/2)−b2R(t)T(t)1+R(t)+c2C(t)1+C(t)(Texture dynamics)Dtα0cC(t)=a3C(t)R(t)1+R(t)−b3C2(t)1+C(t)−c3R(t)T(t)C(t)(1+R(t))(1+T(t))(1+C(t))(Concavity evolution)(26)

where Dtα0c denotes the Caputo fractional derivative of order α∈(0,1], R(t) represents tumor radius (mm), T(t) describes tumor texture (gray-scale variation), and C(t) quantifies concavity score (contour irregularity). The Dtα0c of order α∈(0,1] is given by [42,43]:

Dtα0cg(t)=1Γ(1−α)∫0tg′(τ)(t−τ)αdτ,

where g(t) is an absolutely continuous function, and Γ(⋅) is the gamma function.

The Riemann-Liouville integral and Caputo derivative satisfy the following fundamental relation [44,45]:

Itα0RL(Dtα0cg(t))=g(t)−∑k=0n−1g(k)(0)tkk!,n−1<α≤n,n=⌈α⌉,

where Itα0RL is the Riemann-Liouville fractional integral, defined as:

Itα0RLg(t)=1Γ(α)∫0tg(τ)(t−τ)1−αdτ,α∈(0,1].(27)

The Caputo operator is chosen because it accommodates classical initial conditions and aligns with biological systems where memory effects depend on past states.

4.1.1 Biological Interpretation

The tumor radius R(t) follows fractional logistic growth with carrying capacity K, where the a1 term governs growth rate while texture exerts inhibitory effects (−b1RT) and concavity provides enhancement (+c1C). Texture dynamics T(t) grow proportionally to current values (a2T) but experience suppression from tumor size (−b2RT), with concavity contributing through a saturation term (c2C/(1+C)). Concavity evolution C(t) increases with tumor size (+a3CR) while being self-limiting (−b3C2), with the triple interaction term (−c3RTC) modeling complex biological feedback mechanisms.

Table 5 summarizes the model components and their relationships to measurable biological quantities. The fractional derivative Dtα0c provides key advantages by capturing long-term tumor history through memory effects, which integer-order models cannot represent. This formulation maintains direct interpretability through parameters that relate to feature importance rankings (Table 3) and correlation patterns. The model’s clinical relevance stems from its calibration compatibility with imaging data and potential integration with machine learning predictions for hybrid diagnostics. This mathematical framework provides a foundation for understanding tumor progression dynamics while complementing the machine learning pipeline developed in previous sections.

4.2 Functional Analytic Framework

Definition 1. The solution space is 𝒳=C([0,T],R+3) equipped with the weighted norm:

‖(R,T,C)‖𝒳=supt∈[0,T](e−λt[|R(t)|+|T(t)|+|C(t)|]),λ>0

where R+3={(R,T,C)∈R3:R,T,C≥0}.

This weighted norm construction ensures biological feasibility by enforcing non-negativity of tumor characteristics (radius, texture, concavity), a critical constraint since negative values would lack physical meaning in cancer progression. The exponential weighting accounts for clinical observation that recent tumor states disproportionately influence current growth rates.

Lemma 1. For α∈(0,1) and f∈AC([0,T]), the Caputo derivative satisfies:

1.   Semi-group property: Dtα0c(Itα0RLf)(t)=f(t)

2.   Fundamental theorem: Itα0RL(Dtα0cf)(t)=f(t)−f(0)

3.   Comparison principle: If Dtα0cu(t)≤Dtα0cv(t) and u(0)≤v(0), then u(t)≤v(t).

These operator properties mirror fundamental cancer behaviors: the semi-group property reflects self-consistent tumor evolution, the fundamental theorem ensures measurable changes from baseline pathology, and the comparison principle permits ranking of tumor aggression levels, analogous to histological grading systems.

Theorem 1 (Existence, Uniqueness, and Positivity). For the fractional tumor system in (26), with initial data (R0,T0,C0)∈R+3 and parameters

•   ai,bi,ci>0 for i=1,2,3,

•   K>max{2a2/b2,a3/b3},

•   α∈(0,1),

there exists a unique, non-negative solution (R,T,C)∈𝒳 for all t≥0.

Proof. Define the nonlinearities fi as the right-hand sides of (26). For (R,T,C),(R′,T′,C′)∈Bρ⊂𝒳 such that Bρ={x∈𝒳|‖x‖≤ρ,ρ>0}. Since 1/(1+R)(1+R′),1/(1+T)(1+T′),1/(1+C)(1+C′)<1, we have

|f1(R,T,C)−f1(R′,T′,C′)|≤(a1+2a1ρK+b1ρ(1+ρ))|R−R′|+b1ρ|T−T′|+c1|C−C′|≤η1(|R−R′|+|T−T′|+|C−C′|),

where η1=max{a1+2a1ρK+b1ρ(1+ρ),b1ρ,c1}. The other components satisfy similar bounds:

|f2(R,T,C)−f2(R′,T′,C′)|≤b2ρ|R−R′|+(a2+4a2ρK+b2ρ(1+ρ))|T−T|+c2|C−C′|≤η2(|R−R′|+|T−T′|+|C−C′|),

where η2=max{b2ρ,a2+4ρa2K+b2ρ(1+ρ),c2}.

|f3(R,T,C)−f3(R′,T′,C′)|≤(a3ρ+c3ρ2+2c3ρ3+c3ρ4)|R−R′|+(c3ρ2+3c3ρ3)|T−T′|+(a3ρ(1+ρ)+b3ρ(2+ρ)+c3ρ2+2c3ρ3+c3ρ4)|C−C′|≤η3(|R−R′|+|T−T′|+|C−C′|),

where η3=max{a3ρ+c3ρ2+2c3ρ3+c3ρ4,c3ρ2+3c3ρ3,a3ρ(1+ρ)+b3ρ(2+ρ)+c3ρ2+2c3ρ3+c3ρ4}. Using Lemma 1, convert the given system to the following Volterra form:

R(t)=R0+1Γ(α)∫0tf1(R,T,C)(t−s)1−αds,T(t)=T0+1Γ(α)∫0tf2(R,T,C)(t−s)1−αds,C(t)=C0+1Γ(α)∫0tf3(R,T,C)(t−s)1−αds.

Define the operator ℱ:𝒳→𝒳 componentwise via the integral formulation. For λ sufficiently large:

‖ℱ(R,T,C)−ℱ(R′,T′,C′)‖𝒳≤LλαΓ(α)‖(R,T,C)−(R′,T′,C′)‖𝒳

where L is the Lipschitz constant such that L=max{η1,η2,η3}. Choose λ>(L/Γ(α))1/α to ensure contraction.

For non-negative initial data, the right-hand sides satisfy:

•   f1(0,T,C)=c1C1+C≥0,

•   f2(R,0,C)=c2C1+C≥0,

•   f3(R,T,0)=0.

By the fractional comparison principle (Lemma 1), solutions remain non-negative.

Using the barrier method with:

Dtα0cR≤a1R(1−R/K) ⟹ R(t)≤max(R0,K),

and similar bounds for T(t), C(t), solutions cannot blow up in finite time. □

The contraction mapping guarantees clinically observed tumor trajectory stability, where small measurement errors don’t drastically alter predictions. The non-negativity preservation aligns with radiologically measurable tumor features, while the barrier method’s carrying capacity bound (K) models physiological space constraints in breast tissue

4.3 Global Stability Analysis

Theorem 2 (Global Asymptotic Stability of Tumor Progression Model). For the fractional-order system (26) with parameters satisfying:

1.   ai,bi,ci>0 for i=1,2,3,

2.   K>max(2a2b2,a3b3),

3.   Initial conditions (R(0),T(0),C(0))∈Ω:=R+3∩{R<K,T<K/2},

there exists a unique positive equilibrium E∗=(R∗,T∗,C∗) that is globally asymptotically stable in Ω for all α∈(0,1].

Proof. The equilibrium points satisfy the nonlinear algebraic system:

0=a1R∗(1−R∗K)−b1R∗T∗1+T∗+c1C∗1+C∗(28)

0=a2T∗(1−T∗K/2)−b2R∗T∗1+R∗+c2C∗1+C∗(29)

0=a3C∗R∗1+R∗−b3C∗21+C∗−c3R∗T∗C∗(1+R∗)(1+T∗)(1+C∗)(30)

Eq. (30) can be rewritten as:

C∗[a3R∗1+R∗−b3C∗1+C∗−c3R∗T∗(1+R∗)(1+T∗)(1+C∗)]=0

The nontrivial solution satisfies:

a3R∗1+R∗=b3C∗1+C∗+c3R∗T∗(1+R∗)(1+T∗)(1+C∗)

This implicitly defines a function C∗=ϕ(R∗,T∗) that is continuously differentiable by the Implicit Function Theorem. Substitute C∗=ϕ(R∗,T∗) into (28) and (29):

F1(R∗,T∗)=a1(1−R∗K)−b1T∗1+T∗+c1ϕ(R∗,T∗)R∗(1+ϕ(R∗,T∗))=0F2(R∗,T∗)=a2(1−T∗K/2)−b2R∗1+R∗+c2ϕ(R∗,T∗)T∗(1+ϕ(R∗,T∗))=0

The Jacobian determinant of (F1,F2) is non-zero in Ω when K>max(2a2/b2,a3/b3), guaranteeing a unique solution by the Inverse Function Theorem. Define the Volterra-type Lyapunov functional:

V(t)=Itα0RL[w1(R−R∗−R∗ln⁡RR∗)+w2(T−T∗−T∗ln⁡TT∗)+w3(C−C∗−C∗ln⁡CC∗)](31)

where wi>0 are weighting parameters to be determined. Using the composition rule for Caputo derivatives:

Dtα0cV(t)≤w1(1−R∗R)Dtα0cR+w2(1−T∗T)Dtα0cT+w3(1−C∗C)Dtα0cC=w1(1−R∗R)[a1R(1−RK)−b1RT1+T+c1C1+C]+w2(1−T∗T)[a2T(1−TK/2)−b2RT1+R+c2C1+C]+w3(1−C∗C)[a3CR1+R−b3C21+C−c3RTC(1+R)(1+T)(1+C)]

After extensive algebraic manipulation and application of Young’s inequalities, we obtain:

Dtα0cV(t)≤−∑i=13λi(Xi−Xi∗)2−Ψ(R,T,C)

where:

λ1=w1(a1K−b1ϵ12−b2w2ϵ22w1)λ2=w2(a2K/2−b1w12w2ϵ1−c3w3ϵ32w2(1+C∗))λ3=w3(b3−c1w12w3ϵ4−c2w22w3ϵ5−c3R∗2w3ϵ3(1+R∗)(1+T∗))

Here ϵi>0 are small constants chosen such that all λi>0 under the theorem’s parameter conditions, and Ψ contains higher-order terms that are positive definite.

Applying the fractional comparison principle and the properties of Mittag-Leffler functions:

V(t)≤V(0)Eα(−κtα)+μκ(1−Eα(−κtα))

where κ=min(λi)>0 and μ bounds the residual terms. As t→∞, Eα(−κtα)→0 exponentially, proving global convergence to E∗. □

The equilibrium conditions (28)–(30) encode balanced tumor microenvironments. The exponential convergence rate (κ) via Mittag-Leffler decay provides a theoretical basis for clinical timelines, where higher κ values predict faster stabilization post-treatment, and lower value characterizes aggressive tumors requiring prolonged therapy, matching chemotherapy response patterns.

4.4 Numerical Method

The numerical implementation of our fractional-order tumor model employs a predictor-corrector scheme adapted for Caputo fractional derivatives, implemented in MATLAB R2021a. The complete computational framework consists of the following components: Discretization Scheme:

Predictor Step:

{yk+1P=y0+1Γ(α)∑j=0kbj,k+1f(tj,yj)withbj,k+1=hαα[(k+1−j)α−(k−j)α](32)

Corrector Step:

{tk=t0+kh,k=0,1,…,Nwithh=tfinal−t0Nyk+1=y0+hαΓ(α+2)[f(tk+1,yk+1P)+∑j=0kaj,k+1f(tj,yj)]whereaj,k+1={kα+1−(k−α)(k+1)α,j=0(k−j+2)α+1+(k−j)α+1−2(k−j+1)α+1,1≤j≤k1,j=k+1(33)

4.4.1 Numerical Implementation of the Fractional Tumor Model

The discretized numerical scheme for the Caputo fractional-order system is given by:

{Predictor Step:Rk+1P=R0+hαΓ(α)∑j=0kbj,k+1[a1Rj(1−RjK)−b1RjTj1+Tj+c1Cj1+Cj]Tk+1P=T0+hαΓ(α)∑j=0kbj,k+1[a2Tj(1−TjK/2)−b2RjTj1+Rj+c2Cj1+Cj]Ck+1P=C0+hαΓ(α)∑j=0kbj,k+1[a3CjRj1+Rj−b3Cj21+Cj−c3RjTjCj(1+Rj)(1+Tj)(1+Cj)](34)

{Corrector Step:Rk+1=R0+hαΓ(α+2)[a1Rk+1P(1−Rk+1PK)−b1Rk+1PTk+1P1+Tk+1P+c1Ck+1P1+Ck+1P+∑j=0kaj,k+1(a1Rj(1−RjK)−b1RjTj1+Tj+c1Cj1+Cj)]Tk+1=T0+hαΓ(α+2)[a2Tk+1P(1−Tk+1PK/2)−b2Rk+1PTk+1P1+Rk+1P+c2Ck+1P1+Ck+1P]+∑j=0kaj,k+1(a2Tj(1−TjK/2)−b2RjTj1+Rj+c2Cj1+Cj)Ck+1=C0+hαΓ(α+2)a3Ck+1PRk+1P1+Rk+1P−b3Ck+1P21+Ck+1P−c3Rk+1PTk+1PCk+1P(1+Rk+1P)(1+Tk+1P)(1+Ck+1P)+∑j=0kaj,k+1(a3CjRj1+Rj−b3Cj21+Cj−c3RjTjCj(1+Rj)(1+Tj)(1+Cj))(35)

Weights Calculation:

aj,k+1={(k+1)α−kα,j=0(k−j+2)α+1−3(k−j+1)α+1+3(k−j)α+1−(k−j−1)α+1,1≤j≤k1,j=k+1bj,k+1=hαα[(k−j+1)α−(k−j)α]

Stability Control:

Δk=max(|Rk+1−Rk+1P|1+|Rk+1|,|Tk+1−Tk+1P|1+|Tk+1|,|Ck+1−Ck+1P|1+|Ck+1|)

If Δk>ϵ, then recompute with hnew=0.5h

Memory Management:

Mk=⌊k2⌋(Memory truncation index)

For j<Mk, weights are set to aj,k+1=0 and bj,k+1=0.

The numerical scheme maintains the following properties:

•   Conservation of non-negativity: Rk,Tk,Ck≥0 ∀k

•   Consistency error of order O(hmin(1+α,2))

•   Adaptive memory management for computational efficiency

•   Built-in stability verification through predictor-corrector difference

The implementation requires the following computational resources per time step:

{Storage: O(k) for history termsOperations: O(k2) for full history summationMatrix Inversions: 3×3 Jacobian matrices for Newton iterationsSpecial Functions: Γ(α) precomputed to machine precision(36)

Fractional Derivative Approximation:

Dtα0cy(tk+1)≈1hαΓ(2−α)∑j=0kwj[yk+1−j−yk−j]

where wj=j1−α−(j−1)1−α

Stability Control:

Δyk+1=min(ϵ‖yk‖‖f(tk,yk)‖,hα)⋅f(tk,yk)

with ϵ=10−6(relative tolerance)

Adaptive Step Refinement:

hnew=0.9⋅h⋅(ϵ‖yk+1−yk+1P‖)1/α

The algorithm implements the following computational workflow:

Initialization:

•   Set model parameters (a1,…,c3,K) from biological data

•   Define initial conditions y0=[R(0),T(0),C(0)]T

•   Initialize time grid t∈[t0,tfinal] with N=28 steps

Main Loop:

1.   For each fractional order α∈{0.7,0.75,…,0.95}:

•   Compute predictor yk+1P using historical values

•   Solve corrector equation via Newton-Raphson iteration:

{J=I−hαΓ(α+2)∂f∂y(tk+1,yk+1P)Δy=J−1[y0+hαΓ(α+2)f(tk+1,yk+1P)−yk+1P]yk+1=yk+1P+Δy(37)

•   Apply stability correction if ‖Δy‖>ϵ

Visualization:

•   Generate 2D temporal plots for each state variable

•   Construct 3D surface plots showing α-dependence

•   Compute phase portraits for variable pairs (R,T),(R,C),(T,C)

The numerical implementation provides the following key features:

•   O(h1+α) local truncation error

•   A-stability for α∈(0,1]

•   Memory optimization through variable step-size adaptation

•   Parallel computation across α values

The method preserves the non-negativity of solutions through the modified corrector step:

yk+1←max(yk+1,0)componentwise(38)

ensuring biologically realistic simulations.

5  Result and Discussion

The schematic diagram in Fig. 1 illustrates the directional flow of information in the hybrid framework, showing how the best-performing machine learning model informs the fractional-order model; from diagnostics, prognostics, to growth trajectory predictions. The comprehensive diagnostic pipeline shown in Fig. 3 demonstrates a systematic approach for breast cancer classification, beginning with data preprocessing and feature engineering before evaluating multiple machine learning models. The feature importance analysis in Fig. 4, and Table 3 reveals that worst-area characteristics dominate predictive power, with Area_worst (0.8659) and Radius_worst (0.8561) emerging as the most significant features, consistent with clinical understanding of tumor morphology. Fig. 5 further visualizes the distribution of these top-ranked features across diagnostic classes, showing a clear separation between malignant and benign cases for the most important predictors. Model performance evaluation, as shown in Table 4, yields important insights, with the deep neural network achieving superior accuracy (0.97716) and F1 score (0.96923) as shown in the cross-validated metrics, while maintaining reasonable computation time (5.3755 s). The confusion matrices in Fig. 6 provide granular detail on classification errors, revealing that most models achieve >95% accuracy with minimal false negatives. The bootstrap error analysis in Fig. 7 provides confidence in model stability, with tight confidence intervals around performance metrics. Correlation analysis of top features in Fig. 8 reveals expected relationships between morphological characteristics while maintaining sufficient independence for effective modeling. The SVM model demonstrates exceptional AUC (0.99572) despite longer training times, suggesting particular strength in probabilistic classification. These results are corroborated by the ROC analysis in Fig. 9a and precision-recall curves in Fig. 9b, which show consistently strong performance across all models, with the neural network and SVM maintaining high true positive rates across various thresholds.

Figure 4: Feature rank analysis

Figure 5: Top 6 ranked feature distribution

Figure 6: Confusion matrix for each model

Figure 7: Bootstrap error analysis (100 iterations)

Figure 8: Top ten feature correlation

Figure 9: Metric analysis: (a) ROC Comparison for each model, and (b) Precision-Recall Comparison for each model

The neural network architecture comparison in Fig. 10 illustrates the fundamental differences between shallow and deep approaches, with the deep network’s hierarchical structure enabling more sophisticated feature learning. Fig. 11 compares the decision boundaries of our best-performing shallow and deep neural networks, visualized in the feature space of Area_worst and ConcavePoints_worst. The shallow network (Fig. 11a) demonstrates a relatively larger boundary with moderate curvature, suggesting effective but simpler separation of benign (blue) and malignant (red) cases. In contrast, the deep network (Fig. 11b) exhibits a more strict, non-linear boundary that adapts precisely to local data distributions, particularly in regions where malignant cases exhibit extreme values of both features. The deep network generates a more articulated boundary with localized inflection points that better capture transitional cases. Both models demonstrate strong consensus in high-risk regions. This architectural advantage manifests in the decision boundary plots (Fig. 12), where the deep network captures more complex separations between classes. Feature sensitivity analysis in Fig. 12 confirms that both network types focus primarily on the same key features identified in the importance ranking, though the deep network shows more nuanced interactions between predictors. The implemented graphical user interface, shown in Fig. 13, demonstrates practical clinical applicability, allowing for the integration of these analytical insights into diagnostic workflows with real-time prediction capabilities. The GUI (Fig. 13) integrates real-time risk assessment with growth simulations, enabling clinicians to compare static ML predictions with dynamic trajectories for treatment planning. The feature values are being input, and a diagnosis prediction is being made. These results establish a robust framework for breast cancer diagnosis that strikes a balance between computational efficiency and clinical accuracy. The model fitting results demonstrate robust performance in capturing breast cancer growth dynamics, as evidenced by several quantitative metrics.

Figure 10: Neural network architectures for breast cancer classification: (a) Shallow network with 2 hidden layers (128 and 64 neurons) and 20% dropout; (b) Deep network with 4 hidden layers (256, 128, 64, and 32 neurons) featuring 30% dropout in first three layers and 20% in the last. Node colors indicate layer types: input (blue), hidden (green), and output (red). The visualization shows representative subsets of neurons for clarity, with actual layer sizes noted in parentheses.

Figure 11: Decision boundary comparison plot of Area_worst and ConcavePoints_worst for the best performing models: (a) Shallow Neural Network, and (b) Deep Neural Network

Figure 12: Feature sensitivity comparison plot of best performing models: (a) Shallow Neural Network, and (b) Deep Neural Network

Figure 13: Breast cancer predictor graphics user interface, using the best performing model

The result in Table 6 shows that the optimization process successfully converged with a final error of 0.2877, achieving first-order optimality of 5.58×10−7 within 108 iterations and 14.45 s of computation time. This efficient convergence indicates that the fractional-order model provides a mathematically sound framework for representing tumor progression. The close agreement between model predictions and clinical measurements shown in Fig. 14 further validates the model’s ability to reproduce observed growth patterns across multiple tumor characteristics. The optimized parameters listed in Table 7 reveal biologically meaningful interactions between tumor features. The radius growth rate (a1 = 0.4987) dominates the system dynamics, while the strong texture-radius competition (b2 = 1.0192) explains the inverse relationship observed in the phase portraits. Notably, the carrying capacity (K = 12.3588) suggests tumors in this model can grow to approximately 12 times their initial size before reaching resource limitations. The relatively small triple interaction term (c3 = 0.0082) confirms our earlier observation that concavity remains largely independent of other features during later growth stages. These parameter values collectively produce the characteristic growth curves seen in Fig. 15, where radius stabilizes while texture declines and concavity maintains relative stability. The model’s performance metrics and parameter estimates provide confidence in its ability to represent real breast cancer progression. The low final error and rapid convergence time suggest the framework is both accurate and computationally efficient for clinical applications. The parameter values align with known biological constraints, particularly the competition between tumor size and structural organization. This quantitative validation supports using the model for predicting tumor development patterns and potentially guiding treatment decisions based on individual growth characteristics. The successful fitting to clinical data shown in Fig. 14 demonstrates practical utility for modeling patient-specific tumor dynamics. The analysis of tumor growth dynamics reveals several key patterns in the evolution of tumor characteristics. As shown in Figs. 15a and 16a, the radius (R) exhibits an initial rapid growth phase followed by stabilization, with lower fractional orders (α = 0.70–0.80) demonstrating more gradual stabilization compared to higher orders (α = 0.90–0.95). This prolonged growth phase in low-α tumors may reflect the more aggressive behavior seen in certain breast cancer subtypes. The texture parameter (T) in Figs. 15b and 16b displays a consistent decreasing trend across all α values, with higher α values showing faster reduction to approximately T = 0.1 by 50 time units (5 weeks). This pattern suggests that tumors with classical growth dynamics develop more homogeneous internal structures more rapidly, while those with memory effects maintain structural irregularities for longer durations. The concavity (C) measurements in Figs. 15c and 16c show relatively stable values between 0.06 and 0.1 throughout the simulation period, with minor fluctuations being more pronounced for lower α values. This stability indicates that surface characteristics may be less sensitive to tumor growth dynamics compared to size and internal structure. The surface plots in Fig. 17 further demonstrate how these parameters evolve simultaneously with respect to both time and fractional order, revealing the complex interplay between tumor characteristics. The phase portraits provide additional insights into the relationships between these tumor features. Fig. 17a shows an inverse relationship between radius and texture, where increasing tumor size correlates with decreasing structural irregularity. This negative correlation aligns with clinical observations of larger tumors developing more homogeneous internal structures. The radius-concavity relationship in Fig. 17c demonstrates weak coupling, with concavity remaining relatively stable across different radius values. This suggests that surface deformations are governed by different biological mechanisms than those controlling tumor size. The texture-concavity relationship shown in Fig. 17b reveals nonlinear coupling, where texture dominates concavity modulation at early growth stages but becomes largely independent of structural changes beyond T≈0.2. These phase space trajectories provide valuable information about tumor development patterns, with low-α tumors maintaining higher concavity across all texture values, potentially serving as a biomarker for more aggressive subtypes. The differential responses based on α values support the utility of fractional-order modeling in capturing the spectrum of tumor behaviors observed clinically, particularly in representing the continuum from classical exponential growth to more complex, memory-dependent progression patterns seen in breast cancer. The prediction landscape in Fig. 18 illustrates our model’s ability to distinguish between malignant (red) and benign (blue) cases using fractional-order dynamics with α = 0.9, revealing a nonlinear decision boundary that reflects the complex interactions between tumor size and contour irregularity. The boundary’s curvature aligns with clinical knowledge, demonstrating a greater reliance on concavity for larger tumors while maintaining a robust separation of classes. Misclassified points near the transition zone correspond to diagnostically challenging borderline cases, highlighting the model’s potential for clinical decision support. The α parameter’s influence is evident in the boundary’s shape, where higher values produce more definitive classifications compared to the conservative predictions of lower α settings. This visualization validates how fractional calculus captures both morphological features and growth history, providing an interpretable framework that bridges machine learning predictions with biological understanding of tumor progression.

Figure 14: Model fitting to clinical data

Figure 15: Two dimensional plots of the equation dynamics in model (26): (a) Plot of R(t), (b) Plot of T(t), and (c) Plot of C(t)

Figure 16: Surface plots dynamics of model (26): (a) Plot of R(t), (b) Plot of T(t), and (c) Plot of C(t)

Figure 17: Phase plots of model (26): (a) Plot of T(t) over R(t), (b) Plot of C(t) over T(t), and (c) Plot of C(t) over R(t)