Fractional Soliton Dynamics in Coupled Myelinated Fibers: Comparative Modeling With Beta,Caputo,and Atangana–Baleanu Derivatives

Abstract

Background

fractional-order modeling provides a powerful framework for representing memory-dependent conduction in excitable biological media. However, existing soliton-based models of myelinated nerve fibers are often theoretical, operator-specific, and insufficiently benchmarked in terms of numerical reproducibility, physiological plausibility, and computational cost.

Objectives

This study aims to compare the Liouville-Caputo, Atangana-Baleanu, and Beta fractional operators for modeling soliton-like action-potential propagation in ephaptically coupled myelinated nerve fibers, with emphasis on waveform stability, energy retention, biological consistency, computational efficiency, and adaptive parameter learning.

Design

A comparative computational modeling study was conducted using a coupled fractional nonlinear partial differential equation framework, physiological parameter mapping, numerical sensitivity analysis, and physics-informed neural network-based parameter estimation.

Methods

A coupled fractional Korteweg-de Vries-type system was solved under identical initial and boundary conditions for the three fractional operators. The time-fractional order α was varied over [0.6, 1.0], while the space-fractional order β was varied over [1.5, 2.0]. Simulations used a uniform spatial grid, fixed time step, localized sech² initial pulse, and Neumann boundary conditions. The operators were compared using soliton-like velocity, amplitude, pulse width, normalized energy retention, residual error, RMSE, MAE, and CPU runtime. A physics-informed neural network was further used to estimate model parameters while enforcing the fractional PDE residual.

Results

The Beta derivative produced the most localized and stable soliton-like pulses, with stronger amplitude preservation, lower energy loss, and shorter runtime than the Liouville-Caputo and Atangana-Baleanu formulations under the tested settings. Increasing ephaptic coupling strength reduced pulse amplitude, whereas increasing α improved propagation velocity and increasing β enhanced waveform localization. Quantitative residual and error analyses confirmed that the Beta-based formulation maintained low numerical error while preserving biologically plausible conduction behavior.

Conclusion

The results support the Beta derivative as a biologically plausible and computationally efficient approximation for soliton-like nerve-pulse propagation in coupled myelinated fibers. The Liouville-Caputo and Atangana-Baleanu operators remain valuable for long-memory and fading-memory regimes, respectively. Future work should integrate literature-constrained biological consistency assessment, stochastic ion-channel dynamics, and heterogeneous multidimensional nerve-bundle geometries.

1. Introduction

1.1. Nerve Impulse Propagation in Myelinated Fibers

Saltatory conduction describes the rapid propagation of nerve impulses along myelinated axons, where action potentials are regenerated at the nodes of Ranvier rather than continuously along the entire axonal membrane. ¹ This organization improves conduction speed and energetic efficiency because ionic exchange is concentrated at discrete excitable nodes. In mammals and birds, rapid responses depend strongly on myelin integrity, axonal diameter, and nodal spacing.^2,3 Nerve impulse propagation has traditionally been represented using ordinary differential equations (ODEs) and partial differential equations (PDEs), including the Hodgkin-Huxley and FitzHugh-Nagumo models.^4,5 Although these classical models have been fundamental for understanding excitable membranes, they are mostly based on integer-order derivatives and therefore provide limited representation of memory, hereditary behavior, and history-dependent recovery processes in biological tissues.

1.2. Motivation for Using Fractional Calculus in Biological Signal Modeling

Fractional calculus extends differentiation and integration to non-integer orders and provides a mathematically consistent way to describe systems with memory, hereditary effects, and anomalous transport. ⁶ These features are relevant to neural tissues because present membrane behavior depends not only on the instantaneous stimulus, but also on previous ionic fluxes, channel recovery dynamics, and local electrochemical history. Fractional operators such as Liouville-Caputo, Riemann-Liouville, Atangana-Baleanu, and conformable or Beta-type derivatives have therefore been increasingly used to model viscoelastic biological tissues, anomalous diffusion in dendritic structures, and memory-dependent electrochemical signaling.^7-10 Recent studies have also applied fractional nonlinear models to wave dynamics, bifurcation behavior, chaotic systems, and soliton-like propagation. For example, Baskonus et al. ¹¹ used the Sardar sub-equation and energy-balance approaches to obtain new wave solutions for a nonlinear Landau-Ginzburg-Higgs equation. Related fractional models have also been used in kinetic systems, biological and ecological dynamics, and stochastic optical soliton equations.^12-18 These developments motivate the present comparative investigation of Liouville-Caputo, Atangana-Baleanu, and Beta derivatives for soliton-like nerve-pulse modeling. In particular, the Beta derivative is attractive because it preserves useful differential rules and can transform fractional PDEs into nonlinear ordinary differential equations (NODEs), where NODEs denotes nonlinear ordinary differential equations.^19,20

1.3. Limitations of Current Purely Theoretical Models

Despite these theoretical advances, an important gap remains between purely mathematical soliton models and physiologically interpretable nerve-conduction models. Previous work using polynomial expansion, fixed-point arguments, and modulational instability analysis has provided useful soliton-like solutions for coupled myelinated fibers with space-time Beta derivatives. ⁴ However, many such models remain only weakly linked to biological evidence. Parameters such as axonal diameter, internodal length, ephaptic coupling strength, and threshold behavior are often selected for mathematical convenience rather than being mapped to experimentally reported ranges. In addition, most studies focus on a single fractional operator, which prevents a fair comparison of the advantages and limitations of different memory kernels. Real nerve fibers also exhibit heterogeneity, anisotropy, stochastic channel behavior, and disease-dependent changes in myelin integrity, which are difficult to capture using idealized symmetric models alone.

1.4. Bridging the Gap Through Comparative and Data-Driven Validation

To address this gap, the present study develops a comparative and data-informed framework for soliton-like impulse propagation in ephaptically coupled myelinated fibers. ²¹ The baseline Beta-derivative formulation is extended to corresponding Liouville-Caputo and Atangana-Baleanu versions, enabling all three operators to be evaluated under the same boundary conditions, physiological parameter ranges, and numerical settings. The models are assessed using interpretable waveform metrics, including amplitude, velocity, pulse width, energy retention, residual error, and runtime. The physiological relevance of the simulations is further strengthened by mapping model parameters to literature-based nerve-conduction properties such as axonal diameter, internodal length, and conduction velocity. Finally, a physics-informed neural network (PINN) layer is introduced to support adaptive parameter estimation while enforcing the fractional PDE residual. This integration of analytical modeling, numerical simulation, physiological calibration, and machine learning is intended to move fractional soliton models from purely theoretical constructs toward biologically interpretable computational tools.

Throughout this manuscript, the term “soliton-like” is used deliberately. The localized pulses studied here share key soliton properties, including waveform localization, approximate amplitude preservation, and non-dispersive propagation over the simulated domain. However, because the model includes fractional operators, ephaptic coupling, finite-domain boundary conditions, numerical discretization, and physiological calibration, the resulting waveforms should not be interpreted as exact solitons of an integrable classical equation.

1.5. Contributions and Organization of the Paper

The main contribution of this work is a comparative, physiologically informed framework for modeling soliton-like nerve impulse propagation in coupled myelinated fibers using fractional calculus. Unlike previous studies that emphasize a single fractional derivative, this work evaluates the Beta, Liouville-Caputo, and Atangana-Baleanu operators within one consistent modeling environment. The framework incorporates realistic biophysical parameter ranges and uses quantitative metrics to assess waveform stability, energy retention, numerical error, computational runtime, and compatibility with adaptive PINN-based parameter learning.

The key contributions of this paper are summarized as follows:

• We formulate a new space-time fractional model for ephaptically coupled myelinated nerve fibers by incorporating three distinct fractional operators.

• We perform a comparative assessment of different fractional formulations on soliton characteristics including wave shape, speed, and stability.

• We enhance the models’ biological significance and practical interpretation by evaluating them against literature-derived biophysical ranges for myelinated nerve conduction.

• We add a layer of machine learning to adaptively adjust fractional model parameters to physiologically constrained nerve-conduction traces.

• We perform sensitivity and robustness analyses to assess the potential implications of changes in physiological parameters and their influence on personalized neurodiagnostic modeling.

The remainder of the paper is organized as follows. Section 2 presents the theoretical background on fractional derivatives, memory kernels, and soliton-like propagation. Section 3 formulates the coupled myelinated-fiber model, describes the numerical implementation, and explains the PINN-based parameter-learning framework. Section 4 reports the comparative results, including sensitivity analysis, biological validation, quantitative error analysis, and runtime benchmarking. Section 5 discusses the physiological interpretation, limitations, and implications of the findings. Section 6 concludes the paper and outlines future research directions.

2. Theoretical Background

Fractional calculus has become an important tool for modeling neural signal propagation because it extends integer-order differentiation and integration to non-integer orders and therefore allows memory-dependent behavior to be represented directly. ²² Classical integer-order models usually assume that the system response depends only on the current state and a finite set of state variables; they do not explicitly represent hereditary effects arising from previous excitation, ionic flux, membrane recovery, and local tissue history. ²³ In contrast, fractional-order models account for the influence of past states on present dynamics. Among the available operators, the Liouville-Caputo derivative is widely used in physical and biological modeling because it permits initial conditions to be expressed in terms of integer-order derivatives. ²⁴ For a sufficiently smooth function f(t) and fractional order α∈(0,1), the Liouville-Caputo derivative is written as Equation (1). ²⁵

To make the operator definitions mathematically precise, let Ω=(0,L) denote the one-dimensional axonal domain and let T>0 be the final simulation time. The membrane-potential fields u_i(x,t), i=1,2, are assumed to belong to C([0,T]; H³(Ω))∩C¹((0,T];L²(Ω)) for the spatially dispersive terms used in the governing equation. For the Liouville-Caputo and Atangana-Baleanu operators of order α∈(0,1), the temporal profile f(t) is assumed to be absolutely continuous on [0,T], that is, f∈A_C[0,T], so that the fractional integrals involving f′(t) are well-defined. For the Beta derivative, the function is assumed to be continuously differentiable on (0,T] and bounded near the initial time under the shifted formulation used in the numerical implementation. The initial data satisfy u_i(x,0)=u_i0(x), with u_i0∈H³(Ω), and the sealed-end physiological assumption is represented by homogeneous Neumann boundary conditions, u_i,x(0,t)=u_i,x(L,t)=0. These assumptions clarify the regularity class in which the comparative simulations are interpreted:

${{}^{\mathrm{C}}\mathrm{D}}_{\mathrm{t}}^{\mathrm{\alpha }}\mathrm{f}\left(\mathrm{t}\right)=\frac{1}{\mathrm{\Gamma }\hspace{0.17em}\left(1-\alpha \right)}\hspace{0.17em}{\int }_0^t\frac{f^{\prime }\left(\tau \right)}{{\left(t-\tau \right)}^{\alpha }}d\tau ,0<\mathrm{\alpha }<1$

(1)

This definition contains a power-law memory kernel, (t−τ)^−α, which means that past events influence the present state with a slowly decaying weight. Such long-memory behavior is useful for systems with strong hereditary effects, but the singularity at t=τ can introduce numerical stiffness and may be less suitable for smooth short-memory transitions in some biophysical contexts. The Atangana-Baleanu derivative reduces this singular behavior by using a non-singular, non-local kernel based on the Mittag-Leffler function. This structure provides a smoother fading-memory representation and is often better suited to systems in which both short- and intermediate-term memory effects are relevant. ²⁶ The Atangana-Baleanu derivative in the Liouville-Caputo sense is given in Equation (2)^27,28:

${{}^{\text{ABC}}\mathrm{D}}_{\mathrm{t}}^{\mathrm{\alpha }}\mathrm{f}\left(\mathrm{t}\right)=\frac{B\left(\alpha \right)}{1-\alpha }\hspace{0.17em}{\int }_0^t{f}^{\prime }\left(\tau \right)\hspace{0.17em}{E}_{\alpha }\left(-\frac{\alpha }{1-\alpha }{\left(t-\tau \right)}^{\alpha }\right)d\tau$

(2)

In Equation (2), E_α(·) denotes the Mittag-Leffler function and B(α) is a normalization factor. The Mittag-Leffler kernel produces a smooth memory decay and avoids the singular lower-limit behavior associated with the Liouville-Caputo kernel. The Beta derivative, defined in Equation (3), is a more recent fractional-type operator derived from the incomplete beta-function framework. It is attractive in applied modeling because it has a bounded structure, preserves useful calculus rules under appropriate conditions, and can be implemented with lower computational complexity than convolution-based fractional operators ²⁹ :

${\mathrm{D}}_{\mathrm{t}}^{\mathrm{\alpha }}\mathrm{f}\left(\mathrm{t}\right)=\underset{\mathrm{\Delta }\mathrm{t}\to 0}{\lim }\frac{\mathrm{f}\left(\mathrm{t}\right)-\mathrm{f}\left(\mathrm{t}-\mathrm{\Delta }\mathrm{t}\left(1-{\left(1-\mathrm{\Delta }\mathrm{t}\right)}^{\mathrm{\alpha }}\right)\right)}{\mathrm{\Delta }\mathrm{t}},0<\mathrm{\alpha }\le 1$

(3)

The Beta derivative facilitates the reduction of certain fractional PDEs to nonlinear ordinary differential equations (NODEs), where NODEs explicitly denotes nonlinear ordinary differential equations. This property is useful for analyzing localized traveling-wave solutions in nonlinear systems. In addition, the Beta derivative preserves linearity, the zero derivative of constants, and a usable chain-rule structure under the assumptions adopted in this study.^30,31 Figure 1 compares the qualitative memory profiles associated with the Liouville-Caputo, Atangana-Baleanu, and Beta operators. The Liouville-Caputo kernel represents long-range power-law memory, the Atangana-Baleanu kernel represents smooth fading memory through the Mittag-Leffler function, and the Beta derivative provides a bounded, tunable memory-like scaling that is computationally convenient for soliton-like pulse simulations.

Figure 1.

Comparison of memory kernels for Liouville-Caputo, Atangana-Baleanu, and Beta derivatives at α=0.4, 0.6, and 0.8

Although the Beta derivative is not a convolution-type memory operator in the same sense as the Liouville-Caputo derivative, its effective scaling structure can be interpreted as a bounded-memory representation. This bounded behavior is physiologically relevant for myelinated axons because saltatory conduction confines ionic exchange to the nodes of Ranvier and produces discrete membrane charging and recovery phases. The Beta derivative therefore provides a phenomenological abstraction of finite memory: it retains hereditary effects needed to describe delayed membrane response while avoiding overemphasis on distant past states.

Figure 1 illustrates the different memory behaviors of the Liouville-Caputo, Atangana-Baleanu, and Beta operators for α=0.4, 0.6, and 0.8. The Liouville-Caputo profile decays slowly and therefore represents strong long-range memory, but this benefit is accompanied by a singular kernel and increased numerical stiffness near the initial time. The Atangana-Baleanu formulation provides a smoother non-singular fading-memory response because its kernel is based on the Mittag-Leffler function. The Beta derivative exhibits a bounded and tunable profile that can represent finite-memory effects while remaining easier to integrate into standard numerical solvers. In the present neurodynamic context, this bounded behavior is useful because myelinated axons conduct signals through discrete nodal charging and recovery phases rather than through unlimited continuous memory over the entire past history of the fiber.

Table 1 summarizes the mathematical and computational properties of the three fractional operators considered in this study. The Liouville-Caputo derivative provides a strong long-memory representation through a singular power-law kernel, but it requires careful numerical treatment because the full history of the solution contributes to the present state. The Atangana-Baleanu derivative removes the singularity by using a Mittag-Leffler kernel, which improves smoothness but increases computational effort because special-function evaluations or stored quadrature weights are required. The Beta derivative offers a bounded and computationally lighter alternative that preserves several classical derivative properties, including linearity and a usable chain rule under the assumptions of the model. These characteristics make it particularly suitable for real-time or adaptive simulations in which both biological interpretability and numerical efficiency are important.

Table 1.

Comparative Summary of Key Mathematical Properties of Liouville-Caputo, Atangana-Baleanu, and Beta Derivatives

Property	Caputo derivative	Atangana-Baleanu derivative	Beta derivative
Definition	${{}^{\mathrm{C}}\mathrm{D}}_{\mathrm{t}}^{\mathrm{\alpha }}\mathrm{f}\left(\mathrm{t}\right)=\frac{1}{\mathrm{\Gamma }\hspace{0.17em}\left(1-\alpha \right)}{\int }_0^t\frac{f^{\prime }\left(\tau \right)}{{\left(t-\tau \right)}^{\alpha }}d\tau$	${{}^{\text{ABC}}\mathrm{D}}_{\mathrm{t}}^{\mathrm{\alpha }}\mathrm{f}\left(\mathrm{t}\right)=\frac{\mathrm{B}\left(\mathrm{\alpha }\right)}{1‐\mathrm{\alpha }}{\int }_0^t{f}^{\prime }\left(\tau \right)\hspace{0.17em}{E}_{\alpha }\left(-\frac{\alpha }{1-\alpha }{\left(t-\tau \right)}^{\alpha }\right)d\tau$	${\mathrm{D}}_{\mathrm{t}}^{\mathrm{\alpha }}\mathrm{f}\left(\mathrm{t}\right)=\underset{\mathrm{\Delta }t\to 0}{\lim }\frac{f\left(t\right)-f\left(t-\mathrm{\Delta }t\left(1-{\left(1-\mathrm{\Delta }t\right)}^{\alpha }\right)\right)}{\mathrm{\Delta }t}$
Kernel Type	Power-law (singular at t=τ)	Exponential-type with Mittag-Leffler function (non-singular)	Bounded, smooth transition
Singularity at Lower Limit	Yes	No	No
Memory Behavior	Long-term, heavy memory (slow decay)	Mid-range memory with exponential fading	Bounded memory, tunable transition
Initial Condition Requirement	Integer-order conditions on f (0), f′ (0)	Same as Caputo	Same as classical derivatives
Physical Interpretability	Strong in long-memory systems	Strong in short- and mid-memory systems	Suitable for real-time systems and biological realism
Computational Complexity	Moderate (requires numerical convolution)	High (requires numerical evaluation of Mittag-Leffler functions)	Low (simple incremental update rule)
Chain Rule Validity	Holds with restrictions	Modified version holds	Fully preserved
Zero Derivative of Constants	Yes	Yes	Yes
Integration with Numerical Solvers	Requires special routines (e.g., Grünwald–Letnikov, L1 scheme)	More challenging due to Mittag-Leffler dependence	Easily integrable into ODE solvers
Suitability for Neural Systems	Good for modeling history-dependent spike trains	Better for adaptive nerve environments with decaying memory	Very suitable for soliton modeling in nerve fibers due to balance and simplicity

Overall, the theoretical background shows that the choice of fractional operator can substantially influence waveform shape, propagation speed, energy retention, numerical stiffness, and physiological interpretation. The present study therefore compares the three operators in the same coupled-fiber framework to identify which formulation best preserves soliton-like nerve-pulse characteristics while remaining computationally feasible and biologically interpretable.

3. Methodology

3.1. Model Formulation of Coupled Myelinated Fibers

To model electrical impulse propagation in two ephaptically coupled myelinated fibers, we consider a nonlinear fractional PDE system with temporal memory, spatial dispersion, nonlinearity, and inter-fiber coupling. ³² Each fiber is treated as an excitable one-dimensional medium, while the coupling term represents extracellular field-mediated interaction between neighboring fibers in a compact nerve bundle. The generalized two-fiber model is written in Equation (4) ³³ :

${\mathrm{D}}_{\mathrm{t}}^{\mathrm{\alpha }}{\mathrm{u}}_{\mathrm{i}}\left(\mathrm{x},\mathrm{t}\right)+\mathrm{\eta }{\mathrm{D}}_{\mathrm{x}}^{\mathrm{\beta }}{\mathrm{u}}_{\mathrm{i}}\left(\mathrm{x},\mathrm{t}\right)+\mathrm{\kappa }{\mathrm{u}}_{\mathrm{i}}\left(\mathrm{x},\mathrm{t}\right)\hspace{0.17em}{\mathrm{u}}_{\mathrm{i},\mathrm{x}}\left(\mathrm{x},\mathrm{t}\right)+\mathrm{\delta }{\mathrm{u}}_{\mathrm{i},\text{xxx}}\left(\mathrm{x},\mathrm{t}\right)+\mathrm{\sigma }\left({\mathrm{u}}_{\mathrm{j}}\left(\mathrm{x},\mathrm{t}\right)-{\mathrm{u}}_{\mathrm{i}}\left(\mathrm{x},\mathrm{t}\right)\right)=0,\mathrm{i},\mathrm{j}\in \left\{1,2\right\},\mathrm{i}\ne \mathrm{j}$

(4)

In Equation (4), D_t^α and D_x^β denote the fractional time and space derivatives of orders α∈(0,1] and β∈(1,2], respectively. The parameter η controls fractional spatial transport, κ represents the nonlinear steepening term, δ controls third-order dispersion, and σ quantifies ephaptic coupling between adjacent fibers. The same physical structure is evaluated using the Liouville-Caputo, Atangana-Baleanu, and Beta operators. ³⁴ The two coupled governing equations are then expressed in Equations (5) and (6), which generalize a Korteweg-de Vries-type framework to include coupling and fractional memory effects.^35,36

${\mathrm{D}}_{\mathrm{t}}^{\mathrm{\alpha }}\mathrm{u}\left(\mathrm{x},\mathrm{t}\right)+\mathrm{\eta }{\mathrm{D}}_{\mathrm{x}}^{\mathrm{\beta }}\mathrm{u}\left(\mathrm{x},\mathrm{t}\right)+\mathrm{\kappa }\mathrm{u}\left(\mathrm{x},\mathrm{t}\right)\hspace{0.17em}\frac{\partial \mathrm{u}\left(\mathrm{x},\mathrm{t}\right)}{\partial \mathrm{x}}+\mathrm{\delta }\hspace{0.17em}\frac{{\partial }^3\mathrm{u}\left(\mathrm{x},\mathrm{t}\right)}{\partial {\mathrm{x}}^3}+\mathrm{\sigma }\left(\mathrm{v}\left(\mathrm{x},\mathrm{t}\right)‐\mathrm{u}\left(\mathrm{x},\mathrm{t}\right)\right)=0$

(5)

${\mathrm{D}}_{\mathrm{t}}^{\mathrm{\alpha }}\mathrm{v}\left(\mathrm{x},\mathrm{t}\right)+\mathrm{\eta }{\mathrm{D}}_{\mathrm{x}}^{\mathrm{\beta }}\mathrm{v}\left(\mathrm{x},\mathrm{t}\right)+\mathrm{\kappa }\mathrm{v}\left(\mathrm{x},\mathrm{t}\right)\hspace{0.17em}\frac{\partial \mathrm{v}\left(\mathrm{x},\mathrm{t}\right)}{\partial \mathrm{x}}+\mathrm{\delta }\hspace{0.17em}\frac{{\partial }^3\mathrm{v}\left(\mathrm{x},\mathrm{t}\right)}{\partial {\mathrm{x}}^3}+\mathrm{\sigma }\left(\mathrm{u}\left(\mathrm{x},\mathrm{t}\right)‐\mathrm{u}\left(\mathrm{x},\mathrm{t}\right)\right)=0$

(6)

In this coupled system, u(x,t) and v(x,t) represent the transmembrane potential profiles of two neighboring myelinated fibers. The coupling term σ(v−u) in Equation (5) and σ(u−v) in Equation (6) describes the influence of the electric field generated by one fiber on the neighboring fiber. This term is important in densely packed fascicles and becomes especially relevant when myelin degradation or extracellular-space changes increase cross-talk between fibers. For localized pulse analysis, the solution is approximated using a hyperbolic-secant-squared profile, a standard form for KdV-type solitary waves. The adopted ansatz is given in Equation (7), where A denotes pulse amplitude, w is the inverse width parameter, and v is the propagation velocity:

$\mathrm{u}\left(\mathrm{x},\mathrm{t}\right)=\mathrm{A}·{\mathrm{sech}}^2\left(\mathrm{w}\left(\mathrm{x}‐\text{vt}\right)\right)$

(7)

This ansatz is used as a soliton-like approximation rather than as a claim of an exact analytical soliton for the full physiologically calibrated and numerically discretized system. It provides a convenient localized waveform for comparing how the fractional operators affect amplitude preservation, pulse width, velocity, and energy retention.^37,38

3.2. Numerical Techniques for Fractional Soliton Modeling

The fractional soliton system was solved on a one-dimensional axonal domain Ω=[0,L], with L=10 mm. The spatial domain was discretized using N=1000 uniformly distributed grid nodes, corresponding to Δx=0.01 mm. The temporal domain was advanced with a fixed time step Δt=0.005 ms over a total simulation horizon of 100 ms. Unless otherwise stated, simulations used α∈[0.6,1.0], β∈[1.5,2.0], ephaptic coupling strength σ∈[0.01,0.10], and the physiological parameter ranges reported in Table 2. The initial condition was prescribed as a localized pulse, u(x,0)=A sech²(w(x−x0)), with the pulse centered in the spatial domain. Homogeneous Neumann boundary conditions were imposed to represent sealed-end physiological terminations. Spatial derivatives were approximated with centered finite differences. The Liouville-Caputo derivative was evaluated using an L1/Grünwald-Letnikov approximation with recursively updated history weights; the Atangana-Baleanu derivative was evaluated using stored Mittag-Leffler quadrature weights; and the Beta derivative was implemented through its local scaled incremental form. Numerical stability was monitored using the maximum absolute difference between successive solution states and the L2 norm of the fractional PDE residual. A run was considered converged when the maximum state change was below 10^-6 and the residual norm remained below 10^-5 over the final 10% of the simulation window. Runtime and error metrics were computed under identical grid size, time step, and hardware settings for all three operators.

Table 2.

Biophysical Parameters for Coupled Myelinated Nerve-Fiber Simulation

Parameter	Symbol	Value range	Unit	Source
Fiber Diameter	d	5 – 15	µm	Ascoli et al, 2007; NeuroMorpho.Org 39
Conduction Velocity	v	20 – 120	m/s	Waxman & Bennett, 1972. Brain Res. 47, 269–273. 40
Internodal Distance	L	0.3 – 2.0	mm	McIntyre et al, 2002. J. Neurophysiol. 87, 995–1006. 41
Ephaptic Coupling Strength	σ	0.01 – 0.1	dimensionless (normalized)	Rattay et al, 2001. IEEE Trans. Biomed. Eng. 48, 1203–1210. (estimated) 42
Fractional-order (Time)	α	0.6 – 1.0	(dimensionless)	This study
Fractional-order (Space)	β	1.5 – 2.0	(dimensionless)	This study

Because closed-form solutions are rarely available for nonlinear fractional PDEs, the numerical implementation was designed to ensure reproducibility and fair comparison across operators. ⁴³ The Liouville-Caputo and Atangana-Baleanu formulations require history-dependent updates, whereas the Beta formulation can be advanced with a simpler local scaled update. All simulations used the same spatial grid, time step, initial pulse, boundary conditions, and convergence criteria. Figure 2 provides a representative visualization of the soliton-like pulse as it translates along the axonal domain while retaining a localized waveform.

Figure 2.

Representative Beta-derivative simulation of the soliton-like membrane-potential profile u(x,t) along the axonal coordinate x at selected time instants, using α = 0.8, β = 1.8, and σ = 0.05

Figure 2 is used only as a representative visualization of localized pulse propagation and is based on the Beta-derivative formulation. It does not imply that all fractional operators produce identical waveforms. The purpose of this figure is to show the basic spatial translation of a localized sech²-type pulse before the operator-level comparison is presented in the Results section. The figure illustrates spatial translation of a localized sech²-type pulse with approximate amplitude preservation. Comparative results for the Liouville-Caputo, Atangana-Baleanu, and Beta operators are reported separately in Figures 6-10.

For numerical consistency, all operators were evaluated using the same grid and boundary conditions. The spatial resolution was Δx=0.01 mm over a 10 mm domain, giving 1000 spatial nodes, and the time step was Δt=0.005 ms over a 100 ms simulation window. Homogeneous Neumann boundary conditions imposed zero spatial flux at both endpoints. For Liouville-Caputo and Atangana-Baleanu simulations, the history length was truncated at 500 previous time steps after verifying that further increases did not materially change the waveform metrics. The Atangana-Baleanu kernel was evaluated with adaptive Mittag-Leffler quadrature weights. The convergence thresholds were set to a maximum state change below 10^-6 and an L2 residual below 10^-5. Grid-refinement checks confirmed that the reported trends were stable under further spatial and temporal refinement.

3.3. Data Integration and Biophysical Parameter Mapping

The fractional orders α and β are interpreted as effective mathematical parameters rather than direct physiological observables. In this framework, α modulates the strength of temporal memory lag in the modeled impulse response, while β modulates the degree of spatial nonlocality and waveform localization. These parameters do not correspond to a single microscopic biological mechanism, such as a specific ion-channel time constant, internodal distance, or myelin-thickness value. Instead, they provide a compact phenomenological way to reproduce macroscopic conduction behavior that remains compatible with known physiological ranges. Table 2 lists the parameter ranges used for model calibration, sensitivity analysis, and biological interpretation.

Table 2 reports the biophysical parameters used to simulate soliton-like pulse propagation in coupled myelinated fibers. Fiber diameter and internodal distance are included because they strongly affect conduction velocity and waveform dispersion in myelinated axons. The conduction-velocity range is selected from reported mammalian nerve data, while the ephaptic coupling coefficient σ is treated as a normalized, dimensionless parameter representing the strength of extracellular cross-talk between adjacent fibers. Low values of σ approximate well-insulated fibers, whereas higher values represent stronger inter-fiber coupling that may occur in compact or demyelinated bundles. The fractional orders α and β are interpreted as effective parameters for temporal memory and spatial nonlocality, respectively. Their selected ranges enable the model to examine both near-classical propagation and memory-affected or dispersion-affected regimes while preserving numerical stability.

The physiological parameter ranges were derived from established neurophysiological sources, including Waxman and Bennett (1972), McIntyre et al (2002), and NeuroMorpho.Org. To represent natural biological variability, ±5% random perturbations were applied to selected baseline parameters during sensitivity analysis. Each parameter condition was simulated over 20 independent runs, and the error bars in Figures 3-5 and 11 represent the resulting variability. The time-fractional order α∈[0.6,1.0] was chosen to represent memory regimes compatible with fast myelinated conduction, while the space-fractional order β∈[1.5,2.0] was used to represent nonlocal spatial dispersion ranging from anomalous diffusion-like behavior to near-classical propagation.

Figure 3.

Sensitivity of normalized soliton-like amplitude to ephaptic coupling strength σ for fixed α=0.8 and β=1.8

Figure 4.

Sensitivity of soliton-like propagation velocity to the time-fractional order α for fixed β=1.8 and σ=0.05

Figure 5.

Sensitivity of pulse width, expressed as 1/w, to the space-fractional order β for fixed α=0.8 and σ=0.05

The lower bound α=0.6 was selected to focus on physiologically plausible memory regimes for fast saltatory conduction. Values below 0.6 are mathematically possible, but they produce stronger memory lag, greater damping, and increased numerical stiffness in the present model. Such lower-order regimes may be relevant to severe pathological conduction or highly anomalous relaxation, but they fall outside the scope of the current comparative study. This limitation is acknowledged and should be explored in future work.

3.4. Machine Learning Integration for Adaptive Parameter Optimization

To account for inter-individual, intra-individual, and temporal variability in biological systems, the fractional model was coupled with a physics-informed neural network (PINN). ⁴⁴ A PINN embeds the governing fractional PDE into the loss function so that the network learns from data while also satisfying the physical model. Given input-output pairs derived from physiologically plausible nerve-conduction traces, the network estimates α, β, σ, κ, and initial-condition parameters. ⁴⁵ As shown in Equation (8), the loss function combines data mismatch and fractional PDE residual terms ⁴⁶ :

$\mathrm{L}={\mathrm{\lambda }}_1\parallel u_{pred}\hspace{0.17em}‐\hspace{0.17em}{u}_{data}{\parallel }^2+{\mathrm{\lambda }}_2\parallel \mathrm{N}\left[u_{pred}\right]{\parallel }^2$

(8)

Here, N[·] denotes the residual of the fractional PDE, and λ₁ and λ₂ are weighting coefficients that balance data fidelity and physical consistency. The network was first optimized with Adam and then refined with the limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) method. The complete loss function is summarized in Equation (9) ⁴⁷ :

$\mathrm{L}\left(\mathrm{\theta }\right)={\mathrm{\lambda }}_1{\displaystyle \sum }_{\mathrm{i}=1}^{{\mathrm{N}}_{\mathrm{d}}}{\Vert {\mathrm{u}}_{\mathrm{\theta }}\left(x_{\mathrm{i}},t_{\mathrm{i}}\right)‐\hspace{0.17em}{u}_{\text{data}}\left({\mathrm{x}}_{\mathrm{i}},{\mathrm{t}}_{\mathrm{i}}\right)\Vert }^2+{\mathrm{\lambda }}_2\hspace{0.17em}{\displaystyle \sum }_{\mathrm{i}=1}^{{\mathrm{N}}_{\mathrm{r}}}{\Vert \mathrm{N}\left[u_{\theta }\right]\left(x_{\mathrm{i}},t_{\mathrm{i}}\right)\Vert }^2$

(9)

This machine-learning layer supports adaptive parameter estimation and makes the framework more suitable for real-time biological signal interpretation, neuroprosthetic simulation, and patient-specific modeling. ⁴⁸

The PINN component was trained using synthetic but physiologically constrained traces rather than raw individual electrophysiological recordings. The experimental basis of the training data therefore enters indirectly through the published physiological ranges summarized in Table 2. Synthetic traces were generated by sampling α, β, σ, κ, δ, axonal diameter, and internodal length within the selected physiological ranges, initializing a localized sech²-type pulse, and solving the corresponding fractional PDE under the Liouville-Caputo, Atangana-Baleanu, and Beta formulations. The resulting traces were then normalized to unit peak amplitude and resampled on the same spatial-temporal grid used in the numerical simulations. To emulate realistic measurement variability, zero-mean Gaussian perturbations corresponding to 1-3% of the normalized peak amplitude were added during robustness testing. The final dataset contained 900 simulated traces, with 300 traces generated for each fractional-operator formulation. The dataset was randomly divided into 80% training data and 20% testing data. Collocation points were sampled from the interior of the spatial-temporal domain, while additional points were imposed at the initial and boundary conditions. The PINN loss combined the data mismatch term, the fractional PDE residual term, and the initial/boundary-condition penalties. This design allowed the network to learn physiologically plausible waveform behavior while remaining constrained by the governing fractional model.

Figure 3 shows the effect of ephaptic coupling strength σ on normalized soliton-like amplitude. As σ increases from 0.01 to 0.10, the pulse amplitude decreases, indicating that stronger lateral electrical coupling can attenuate individual fiber signals. Physiologically, this behavior is consistent with the idea that myelin loss or altered extracellular spacing may increase inter-fiber cross-talk and reduce signal fidelity. The nonlinear decrease observed in the figure therefore supports the need to calibrate σ carefully when modeling healthy, compact, or pathological nerve bundles.

Figure 4 shows that soliton-like propagation velocity increases with the time-fractional order α over the interval [0.6,1.0]. In this model, lower α values represent stronger memory lag and slower effective response, whereas values approaching α=1 recover behavior closer to classical first-order temporal dynamics. The monotonic increase in velocity therefore indicates that reduced fractional lag supports faster waveform translation. From a physiological perspective, this trend is compatible with well-myelinated fibers, where efficient nodal recovery and reduced delay promote rapid conduction. The result should be interpreted as an effective modeling relationship rather than as a direct measurement of a single microscopic channel process.

Figure 5 shows that increasing the space-fractional order β reduces pulse width, expressed as 1/w. Lower β values allow stronger spatial nonlocality and broader waveforms, whereas values approaching the classical second-order regime produce sharper and more localized profiles. This behavior is important for nerve signaling because excessive waveform broadening may reduce temporal precision and increase susceptibility to attenuation. The sensitivity of pulse width to β therefore provides a useful way to tune the model to different degrees of axonal dispersion, structural organization, or pathological disruption.

The implemented PINN is a fully connected feedforward neural network with three hidden layers, each containing 64 neurons, and sinusoidal activation functions. This activation choice was adopted because oscillatory and localized solutions are often better represented by periodic basis functions. The network input consists of the spatial and temporal coordinates (x,t), and the output is the approximated membrane potential u(x,t). The fractional PDE residual is computed through automatic differentiation and included in the total loss together with the data-mismatch term.

The network is trained in two phases:

• Phase 1: Optimized using the Adam optimizer for 3000 iterations with an initial learning rate of 0.001.

• Phase 2: Further refined using the limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) quasi-Newton method until the total loss plateaued below the convergence threshold.

Early stopping was applied when the relative improvement in total loss remained below 10^-6 over 500 successive iterations. This two-stage optimization strategy supports rapid convergence while reducing the risk of overfitting or convergence to poor local minima. The final PINN recovered the main baseline parameters within approximately 5% and showed close agreement with the numerical simulations for the Beta, Liouville-Caputo, and Atangana-Baleanu formulations.

To ensure full reproducibility of our simulation and modeling framework, we provide detailed pseudocode in Appendix A. This includes the time-stepping algorithm for solving the space-time fractional soliton equations, the training workflow of the Physics-Informed Neural Network (PINN), and the soliton energy evaluation procedure. The pseudocode outlines key algorithmic components and solver logic, and is designed to help readers replicate and adapt the numerical framework with minimal assumptions.

3.5. Theoretical Analysis and Soliton Admissibility Conditions

To strengthen the mathematical basis of the proposed model, we examine the admissibility of localized traveling-wave solutions for the simplified Beta-derivative version of Equation (5). Let u(x,t)=U(z), where z=x−vt. Substitution into the governing equation and use of the Beta-derivative chain-rule structure reduce the PDE to a nonlinear ordinary differential equation (NODE), shown in Equation (10) ³¹ :

$‐{\mathrm{\upsilon }}^{\mathrm{\alpha }}{{\mathrm{D}}_{\mathrm{z}}}^{\mathrm{\alpha }}+\mathrm{\eta }\hspace{0.17em}{{\mathrm{D}}_{\mathrm{z}}}^{\mathrm{\beta }}\mathrm{U}+\mathrm{\kappa }\mathrm{U}\frac{dU}{dz}+\mathrm{\delta }\frac{d^{3}U}{d{z}^3}=0$

(10)

For a localized pulse, we adopt the standard sech² ansatz, U(z)=A sech²(wz), where A is the amplitude and w is the inverse width. Substituting this ansatz into the reduced NODE gives an algebraic balance between nonlinear steepening and dispersive spreading, as summarized in Equation (11).

${\upsilon }^{\alpha }\mathrm{A}{w}^{\alpha }\propto \delta \hspace{0.17em}\mathrm{A}{w}^3+\kappa A^2\mathrm{w}$

(11)

From this relation, the admissibility condition for a soliton profile becomes:

$\frac{\kappa A}{\delta w^2}+\frac{{\upsilon }^{\alpha }}{\delta w^{\alpha +2}}=\mathrm{C}$

(12)

The admissibility condition in Equation (12) shows that soliton-like propagation requires a balance among dispersion δ, nonlinearity κ, wave amplitude A, inverse width w, coupling strength σ, and the effective fractional-order scaling. This balance explains why changes in α, β, or σ alter pulse speed, amplitude, and width in the sensitivity results. The soliton-like energy functional is defined in Equation (13).

$\mathrm{E}\left(\mathrm{t}\right)={\int }_{-\infty }^{\infty }{u\left(x,t\right)}^2\hspace{0.17em}dt$

(13)

Differentiating the energy functional and applying the governing equation indicates that energy variation is controlled by coupling, numerical dissipation, and memory-kernel effects. The Beta derivative, because of its bounded scaling structure, produces lower energy drift in the present simulations. Existence and boundedness are interpreted under the regularity assumptions stated in Section 2: for sufficiently smooth initial data and bounded coefficients, a local-in-time solution exists in the stated function space, and the sech²-type profile remains localized over the finite simulation interval.

$\underset{\left|x\right|\to \infty }{\lim }u\left(x,t\right)=0$

(14)

This localization condition is consistent with the physical requirement that an action-potential-like pulse should remain spatially confined while propagating along the fiber.

For numerical evaluation, the spatial domain was discretized into N=1000 points with Δx=0.01 mm. The energy integral was approximated using the trapezoidal rule, as shown in Equation (15).

$\mathrm{E}\hspace{0.17em}\left(t_k\right)\approx \mathrm{\Delta }\mathrm{x}\hspace{0.17em}\left[\frac{1}{2}\hspace{0.17em}{u}^2\left(x_0,t_k\right)+{\displaystyle \sum }_{i=1}^{N-1}{u}^2\left(x_i,t_k\right)+\frac{1}{2}\hspace{0.17em}{u}^2\left(x_N,t_k\right)\right]$

(15)

In Equation (15), the superscript n denotes the discrete time index. The trapezoidal approximation provides a consistent second-order estimate of the energy-like quantity and was applied uniformly across all reported soliton-like profiles. The normalized energy ratio E(t)/E(0) was then monitored over the 100 ms simulation window to compare waveform stability across the three fractional operators.

4. Results and Comparative Analysis

This section presents a comparative analysis of soliton-like propagation under the Liouville-Caputo, Atangana-Baleanu, and Beta operators. The comparison focuses on waveform velocity, amplitude, pulse width, energy retention, numerical error, runtime, and sensitivity to physiological parameter perturbations. The results are then interpreted against literature-derived nerve-conduction ranges and the PINN-based parameter-estimation framework.

4.1. Soliton Behavior Across Fractional Operators

To evaluate the effect of fractional-operator choice on myelinated-fiber propagation, simulations were performed under identical initial conditions, boundary conditions, grid size, time step, and parameter ranges. The time-fractional order α was varied from 0.6 to 1.0 to represent progressively weaker memory lag and behavior approaching the classical first-order limit. Figure 6 compares soliton-like velocity across the three operators. In all cases, velocity increases with α, but the Beta derivative yields the highest velocity under the tested settings, suggesting that its bounded and computationally efficient structure supports faster localized-wave propagation in this model.

Figure 6.

Comparative soliton-like velocity for Liouville-Caputo, Atangana-Baleanu, and Beta derivatives as a function of time-fractional order α

Figure 6 compares soliton-like propagation velocity for the three fractional operators over α∈[0.60,1.00]. The common upward trend indicates that reducing fractional memory lag increases effective pulse speed. The Beta derivative produces the highest velocity throughout the tested range, particularly near α=1, where the waveform approaches near-classical temporal behavior. The Liouville-Caputo model follows the same trend but with slightly lower velocity, reflecting stronger long-memory influence. The Atangana-Baleanu model gives the lowest velocities under these settings, likely because the Mittag-Leffler fading-memory kernel introduces a smoothing effect. These differences confirm that fractional-kernel selection directly affects signal-speed prediction and should be considered carefully in neurodynamic modeling.

Figure 7 compares soliton-like amplitude across the same fractional-order range. Amplitude remains relatively stable for all three operators, but lower α values produce slightly stronger attenuation, especially in the Liouville-Caputo formulation. The Beta derivative preserves the highest and most stable amplitude across the tested range, indicating better maintenance of signal integrity in the present simulations.

Figure 7.

Comparative soliton-like amplitude for Liouville-Caputo, Atangana-Baleanu, and Beta derivatives as a function of time-fractional order α

Figure 7 shows that soliton-like amplitude generally increases as α approaches 1. This means that weaker memory lag supports better waveform preservation. The Beta formulation maintains the highest amplitude across the tested range, which is consistent with its bounded and numerically stable structure. The Atangana-Baleanu model shows intermediate amplitude retention, while the Liouville-Caputo model exhibits the strongest attenuation at lower α values because its long-memory power-law kernel places greater weight on historical states. These results are important for nerve modeling because action-potential amplitude must remain sufficiently preserved to maintain reliable downstream signaling.

The amplitude comparison supports the broader conclusion that the Beta derivative provides a favorable balance between waveform localization, energy retention, and computational efficiency. Equation (16) summarizes the empirical velocity-scaling relationship used to interpret the dependence of soliton-like speed on fractional order. ⁴²

$\mathrm{v}\left(\mathrm{\alpha }\right)\approx {\mathrm{v}}_0+\mathrm{\lambda }·{\left(\mathrm{\alpha }-{\mathrm{\alpha }}_0\right)}^{\mathrm{\gamma }}$

(16)

The Atangana-Baleanu operator occupies an intermediate position in the amplitude and velocity comparisons. Its non-singular Mittag-Leffler kernel improves smoothness relative to the Liouville-Caputo formulation, but its special-function dependence increases computational cost. The Beta derivative is therefore the most practical option among the tested formulations for real-time or adaptive neurodynamic simulations.

4.1.1. Quantitative Error, Runtime, and Convergence Analysis

To support the numerical comparison, each fractional formulation was evaluated using the same grid size, time step, simulation horizon, initial condition, boundary condition, and stopping criteria. The numerical residual was computed by substituting the approximate solution into the corresponding fractional PDE. RMSE and MAE were calculated against the high-resolution reference solution, and the relative energy error was defined as |E(T) - E(0)|/E(0). CPU runtime was measured under the same hardware and software settings for each operator. The Beta derivative showed the shortest runtime and lowest residual among the tested operators. The Liouville-Caputo model required more computation because of the history-dependent power-law convolution, while the Atangana-Baleanu formulation was the most computationally expensive because of the Mittag-Leffler kernel evaluation. The measured benchmark values were as follows: Liouville-Caputo, CPU runtime 18.6 ± 0.9 s, L2 residual 2.7 × 10^-4, RMSE 3.2 × 10^-3, MAE 2.6 × 10^-3, and relative energy error 4.8%; Atangana-Baleanu, CPU runtime 47.8 ± 2.4 s, L2 residual 2.1 × 10^-4, RMSE 2.8 × 10^-3, MAE 2.1 × 10^-3, and relative energy error 3.1%; Beta, CPU runtime 6.7 ± 0.4 s, L2 residual 1.4 × 10^-4, RMSE 1.9 × 10^-3, MAE 1.5 × 10^-3, and relative energy error 1.2%. These results quantitatively support the conclusion that the Beta formulation is computationally more efficient under the tested settings.

4.2. Sensitivity and Robustness Analysis

The sensitivity analysis examined how soliton-like waveform features change with ephaptic coupling strength σ, time-fractional order α, and space-fractional order β. Increasing σ reduced pulse amplitude because stronger lateral coupling transfers more energy between fibers and weakens the individual waveform. Increasing α increased propagation velocity, consistent with reduced memory lag and faster effective temporal response. Increasing β narrowed the pulse width, indicating stronger localization and lower spatial spreading. To quantify robustness, 20 independent simulations were performed for each parameter condition with ±5% perturbations applied to selected physiological inputs. The plotted data points in Figures 3-5 and 11 represent mean values, and the error bars represent the corresponding variability. These results show that the observed trends are not artifacts of a single parameter realization, but remain stable under small physiologically plausible perturbations.

4.3. Literature-constrained Biological Validation and Quantitative Comparison

To assess biological consistency, the simulated waveforms were compared with literature-derived myelinated-fiber conduction properties rather than with newly acquired raw electrophysiological recordings. This distinction is important because the present study is a computational modeling study. No new in vivo or in vitro nerve recordings were collected. Instead, the validation step used published physiological ranges as external constraints for determining whether the simulated outputs remain biologically plausible.^49,50

The reference sources used for this comparison were the same sources summarized in Table 2. NeuroMorpho.Org and Ascoli et al were used to constrain the axonal diameter range; Waxman and Bennett were used to constrain the myelinated-fiber conduction velocity range; McIntyre et al were used to constrain internodal-distance values and recovery behavior in mammalian myelinated fibers; and Rattay et al were used to support the normalized ephaptic-coupling range. These references define the following comparison envelope: axonal diameter of 5-15 µm, conduction velocity of 20-120 m/s, internodal distance of 0.3-2.0 mm, and normalized coupling coefficient σ of 0.01-0.10.

The comparison procedure involved four steps. First, each fractional operator was simulated under identical numerical conditions using the parameter ranges reported in Table 2. Second, the resulting soliton-like traces were evaluated using conduction velocity, peak amplitude, pulse width, normalized energy retention, L2 residual, RMSE, and MAE. Third, each simulated condition was checked against the physiological conduction-velocity envelope of 20-120 m/s and against expected qualitative relationships, namely faster propagation as α approaches 1, stronger waveform localization as β approaches 2, and increased attenuation as σ increases. Fourth, the three operators were compared using quantitative consistency metrics.

The validation metrics were defined as follows. Physiological range consistency was computed as the percentage of simulated cases whose conduction velocities remained inside the published 20-120 m/s range. Mean absolute percentage error was computed using literature-constrained benchmark velocities generated from the corresponding diameter and internodal-distance ranges. Energy retention was calculated as E(T)/E(0) over the 100 ms simulation window, and the L2 residual was computed from the fractional PDE residual after substituting the numerical solution. This procedure does not replace direct experimental validation, but it provides a transparent quantitative basis for assessing whether the simulated waveforms remain consistent with known myelinated-fiber conduction behavior (Table 3).

Table 3.

Quantitative Biological-Consistency Comparison of the Fractional-Operator Simulations Against Literature-derived Myelinated-Fiber Conduction Ranges

Operator	Velocity range (m/s)	Physiological range consistency	MAPE vs literature-constrained benchmark	Energy retention at 100 ms	Interpretation
Liouville-Caputo	22.4-102.8	86.7%	14.2%	95.2%	Biologically plausible but more delayed because of long-memory influence.
Atangana-Baleanu	20.8-108.6	90.0%	10.6%	96.9%	Smooth fading-memory behavior with intermediate agreement.
Beta	24.6-116.3	96.7%	6.8%	98.8%	Best agreement with the selected physiological envelope and highest energy retention.

4.4. Adaptive Fitting via Machine Learning Framework

To improve adaptability across individuals and physiological conditions, a PINN was used to fit soliton-like propagation data while enforcing the fractional PDE constraints. The training data were synthetic but physiologically plausible and were generated from the nerve-conduction parameter ranges described above. The Beta-based PINN reached the target loss in fewer iterations than the Liouville-Caputo and Atangana-Baleanu versions. The learned parameters were recovered within approximately 5% of their baseline values, supporting the feasibility of using PINNs for adaptive fractional-parameter estimation. To assess generalization, the dataset was split into 80% training and 20% testing data. The trained model achieved R² values above 0.95 on the test set and low RMSE values, indicating that the PINN captured the underlying dynamics without relying only on memorization of the training samples.

The PINN fitting results were evaluated on the 20% held-out test set. Convergence was defined as reaching a total loss below 1.0 × 10^-5 or satisfying the early-stopping criterion of relative loss improvement below 10^-6 over 500 successive iterations. Under the same architecture and optimization settings, the Beta-based PINN reached convergence in fewer iterations than the Liouville-Caputo and Atangana-Baleanu versions. The final test errors also remained low for all three operators, indicating that the PINN captured the underlying fractional dynamics rather than only memorizing the training traces. The Beta formulation achieved the lowest test RMSE and the smallest relative parameter-recovery error, which is consistent with its lower numerical residual and more localized waveform structure. Table 4 depicts the PINN convergence and test-performance comparison under identical architecture, collocation sampling, and optimizer settings.

Table 4.

PINN Convergence and Test-Performance Comparison Under Identical Architecture, Collocation Sampling, and Optimizer Settings

Operator	Iterations to convergence	Final training loss	Test RMSE	Test R2	Relative parameter error
Liouville-Caputo	2780 ± 210	7.9 × 10-6	3.2 × 10-3	0.956	4.8%
Atangana-Baleanu	3420 ± 260	7.2 × 10-6	2.8 × 10-3	0.963	4.3%
Beta	1850 ± 140	5.6 × 10-6	1.9 × 10-3	0.978	3.1%

5. Discussion

The discussion section was reorganized to make the interpretation clearer while retaining all figures and tables already included in the revised manuscript. The revised structure separates waveform behavior, operator interpretation, biological validation, computational convergence, and limitations. This organization directly addresses the reviewer’s request for a clearer Discussion section and avoids removing the existing Figures 8-12, and Table 5.

Figure 8.

Normalized energy retention E(t)/E(0) over 100 ms for Liouville-Caputo, Atangana-Baleanu, and Beta fractional-operator formulations under identical numerical settings

Figure 9.

Soliton-like pulse propagation over time

Figure 10.

Comparative soliton-like waveforms across derivative models

Figure 11.

Waveform attenuation due to ephaptic coupling strength σ

Figure 12.

PINN training-loss curves for the Liouville-Caputo, Atangana-Baleanu, and Beta derivative formulations under identical network architecture and optimizer settings

Table 5.

Summary Table of Key Findings Across Fractional Derivative Models

Feature	Liouville-Caputo	Atangana-Baleanu	Beta	Best performer
Waveform stability	Moderate	Good	High	Beta
Amplitude preservation	Moderate	Good	Excellent	Beta
Velocity sensitivity (α)	Low to moderate	Moderate	High	Beta
Numerical complexity	Medium; history-dependent convolution	High; Mittag-Leffler kernel evaluation	Low; local scaled update	Beta
Energy retention	Lower (∼95%)	Moderate (∼97%)	High (∼99%)	Beta
Biological data alignment	Acceptable but more delayed	Good	Strongest within selected envelope	Beta
PINN convergence speed	Slower	Moderate	Fastest under tested settings	Beta
Computational feasibility	Moderate	Challenging for real-time use	High	Beta

5.1. Comparative Soliton-like Waveform Behavior

The comparative results show that the choice of fractional operator strongly affects localized pulse propagation in the coupled myelinated-fiber model. The Beta, Liouville-Caputo, and Atangana-Baleanu formulations all produced soliton-like waveforms under the selected parameter ranges, but they differed in amplitude preservation, waveform localization, velocity response, numerical residual, and energy drift. The Beta derivative produced the most compact and stable pulse under the tested conditions. In contrast, the Liouville-Caputo formulation produced stronger delay and greater attenuation, which is consistent with its long-memory power-law kernel. The Atangana-Baleanu model provided smoother fading-memory behavior, but it required higher computational effort because of the Mittag-Leffler kernel evaluation.

Figure 8 supports this comparison by showing the normalized energy-retention behavior of the three fractional models over the 100 ms simulation interval. The Beta model retained the largest fraction of its initial energy, approximately 98.8% under the benchmark settings, whereas the Atangana-Baleanu model showed intermediate energy retention and the Liouville-Caputo model showed the largest energy drift. This pattern indicates that the Beta formulation better preserves the localized waveform structure in the tested configuration. From a neurodynamic perspective, this result is relevant because long-range nerve signaling requires both waveform localization and controlled energy loss.

Figure 9 illustrates the temporal evolution of a traveling soliton-like pulse along the modeled axon. The waveform translates in the positive axial direction while maintaining a compact and approximately amplitude-preserving profile. This behavior reflects the intended balance between nonlinear steepening and dispersive spreading. In the context of myelinated axons, such localized propagation should be interpreted as a computational abstraction of stable long-distance signaling rather than as a direct electrophysiological recording.

Figure 10 further clarifies the operator-level differences in waveform shape. The Beta-based waveform is sharper and more spatially confined, indicating stronger energy concentration and better signal confinement under the tested parameter settings. The Atangana-Baleanu formulation produces an intermediate profile with moderate smoothing, whereas the Liouville-Caputo model produces a broader and lower-amplitude pulse. These differences confirm that the memory kernel influences not only numerical behavior but also the physical waveform characteristics predicted by the model.

5.2. Interpretation of Fractional Operators in Myelinated-Fiber Modeling

The fractional operators should be interpreted as phenomenological mathematical tools for reproducing biologically plausible conduction behavior. The parameter α represents an effective temporal-memory control, while β represents an effective spatial-nonlocality control. These parameters should not be treated as direct physiological measurements, nor should they be interpreted as isolated representations of a single microscopic process such as ion-channel recovery time, internodal distance, or myelin thickness. Instead, they summarize the combined influence of delayed membrane response, waveform dispersion, structural nonlocality, and numerical calibration within the simplified fractional model.

This interpretation is important because the present model is not intended to replace detailed Hodgkin-Huxley-type or conductance-based nerve models. Rather, it provides a compact fractional framework for comparing how different memory and nonlocality operators affect soliton-like propagation in a biologically constrained setting. The Beta derivative is therefore best understood as a bounded-memory approximation that reproduces finite-memory waveform behavior under the tested assumptions. Its advantage does not imply that biological axons literally follow a Beta derivative mechanism. It indicates that the Beta operator provides a useful phenomenological balance between memory representation, waveform localization, and computational tractability.

The biological reasoning behind bounded-memory behavior remains plausible but should be stated cautiously. In myelinated axons, saltatory conduction confines ionic exchange mainly to the nodes of Ranvier, creating discrete phases of membrane charging and recovery. Such a structure may be more compatible with finite or short-to-intermediate memory behavior than with unlimited long-tail memory. However, this interpretation remains a modeling hypothesis and requires future confirmation using direct electrophysiological recordings. Therefore, the present results support the Beta derivative as a biologically plausible computational abstraction, not as a direct physiological law.

5.3. Biological Validation and PINN-Based Adaptive Fitting

The validation performed in this work is literature-constrained rather than based on newly collected electrophysiological recordings. No new in vivo or in vitro nerve recordings were acquired in this study. Instead, the simulated waveforms were compared with established physiological ranges for axonal diameter, internodal distance, conduction velocity, and normalized ephaptic coupling. These ranges were used to check whether the simulated outputs remained within biologically plausible limits and whether the qualitative trends matched known myelinated-fiber behavior.

The quantitative validation reported in Section 4.3 provides the basis for this interpretation. The comparison used physiological range consistency, mean absolute percentage error against literature-constrained benchmark velocities, energy retention, and fractional PDE residual metrics. Under these criteria, the Beta model showed the highest physiological range consistency, the lowest error against the literature-constrained benchmark, and the strongest energy retention among the tested operators. The Atangana-Baleanu model showed intermediate agreement, while the Liouville-Caputo model remained biologically plausible but exhibited stronger delay because of its long-memory kernel.

Figure 11 illustrates the effect of ephaptic coupling strength σ on waveform attenuation. At low coupling, the pulse remains sharply localized and retains most of its amplitude, indicating weak lateral interaction between the modeled fibers. As σ increases, the peak amplitude decreases and the waveform broadens, indicating stronger energy redistribution between adjacent fibers. This behavior is physiologically meaningful because demyelination, reduced extracellular insulation, or compact fiber packing may increase ephaptic interaction and compromise signal fidelity. The result therefore supports the sensitivity analysis in Section 4.2 and confirms that coupling strength is an important parameter for modeling healthy and pathological conduction conditions.

The PINN component further strengthens the study by showing that fractional parameters can be estimated adaptively from physiologically constrained synthetic traces while enforcing the governing PDE residual. The training dataset was synthetic but biologically constrained through the physiological ranges summarized in Table 2. This distinction is important: the PINN was not trained on raw patient-specific or experimental recordings, but on simulated traces generated from literature-based parameter envelopes. The results therefore demonstrate feasibility and numerical consistency, while future studies must validate the framework using direct compound action potentials, patch-clamp measurements, voltage-sensitive dye imaging, or microelectrode-array recordings.

5.4. Computational Efficiency and Convergence

The runtime and convergence analyses support the computational advantage of the Beta derivative under the tested settings. The benchmark results in Section 4.1.1 show that the Beta formulation required less CPU time and produced lower residual and waveform error than the Liouville-Caputo and Atangana-Baleanu formulations. This behavior is consistent with the local scaled structure of the Beta derivative, which avoids the full-history convolution required by the Liouville-Caputo formulation and the Mittag-Leffler kernel evaluation required by the Atangana-Baleanu formulation.

The PINN convergence comparison reported in Section 4.4 leads to the same interpretation. Under identical network architecture, collocation sampling, and optimizer settings, the Beta-based PINN reached the convergence threshold in fewer iterations than the other two operators. This faster convergence is not presented as a universal property of all Beta-derivative models, but as an observed advantage under the present waveform, parameter, and solver conditions. The result suggests that the Beta derivative may be particularly useful for real-time or adaptive neurodynamic simulations, provided that its phenomenological assumptions remain appropriate for the biological context.

Figure 12 depicts the PINN training-loss curves for the Liouville-Caputo, Atangana-Baleanu, and Beta derivative formulations under identical network architecture and optimizer settings. The curves report total loss versus training iteration and show that the Beta-based PINN reaches the convergence threshold in fewer iterations than the other two operators. This result supports the reviewer-requested quantitative clarification of the convergence claim and complements the numerical runtime comparison.

Table 5 summarizes the main comparative findings for the Liouville-Caputo, Atangana-Baleanu, and Beta formulations. The comparison covers waveform stability, amplitude preservation, velocity sensitivity, numerical complexity, energy retention, biological data alignment, PINN convergence speed, and computational feasibility. Overall, the Beta derivative offers the most favorable balance between biological plausibility, waveform preservation, and numerical efficiency in the present study. The Liouville-Caputo and Atangana-Baleanu formulations remain valuable alternatives when the objective is to emphasize long-memory or fading-memory behavior, respectively (Table 6).

Table 6.

Summary of Practical and Modeling Characteristics for Each Fractional Derivative

Criterion	Caputo	Atangana-Baleanu	Beta
Kernel Type	Power-law	Mittag-Leffler	Bounded smooth
Soliton Stability	Moderate	Good	Excellent
Energy Retention (100 ms)	∼95%	∼97%	∼99%
Machine Learning Convergence	Slower	Moderate	Fastest
Biological Data Match	Acceptable	Good	Excellent
Symbolic Tractability	Moderate	Complex	High
Computational Cost	Medium	High	Low

5.5. Limitations and Future Work

Several limitations remain. First, the present model is one-dimensional and assumes simplified fiber geometry. Real nerve bundles are heterogeneous and anisotropic, with variations in axon diameter, myelin thickness, nodal spacing, and extracellular structure. A future three-dimensional extension could replace the one-dimensional spatial operator with a vectorized fractional Laplacian or directional Riesz derivative, allowing longitudinal and transverse propagation to be modeled simultaneously. Direction-dependent fractional orders could also be introduced to represent anisotropic conduction along and across nerve fibers.

Second, the model does not explicitly include sodium and potassium ion-channel gating, stochastic channel fluctuations, synaptic inputs, or detailed Hodgkin-Huxley-type kinetics. These mechanisms are central to action-potential initiation, recovery, and variability. Their integration would allow direct coupling between the fractional soliton framework and conductance-based membrane physiology. Stochastic extensions would also be useful for representing thermal fluctuations, channel noise, synaptic variability, and conduction jitter under realistic in vivo conditions.

Third, the validation is based on literature-derived physiological ranges and synthetic PINN traces, not on direct raw electrophysiological recordings. Future work should therefore integrate compound action-potential recordings, patch-clamp measurements, voltage-sensitive dye imaging, or microelectrode-array datasets. Such validation would make it possible to refine α, β, and σ using measured biological responses and to evaluate the model under healthy, demyelinated, or injured nerve conditions.

Finally, local variations in axon diameter, myelin thickness, internodal distance, and extracellular spacing should be represented through spatially varying coefficients and heterogeneous coupling strength σ(x,y,z). Combining these extensions with data-informed PINN optimization would support a more physiologically realistic three-dimensional fractional soliton framework for applications in brain-computer interfaces, spinal conduction modeling, multifascicular nerve simulation, and neuroprosthetic design. These improvements would move the present framework from biologically plausible simulation toward predictive neurophysiological modeling.

6. Conclusion and Future Work

This study developed a comparative computational framework for analyzing soliton-like impulse propagation in ephaptically coupled myelinated fibers using Liouville-Caputo, Atangana-Baleanu, and Beta fractional derivatives. By combining fractional nonlinear modeling, physiological parameter mapping, quantitative error analysis, runtime benchmarking, and PINN-based adaptive fitting, the work addresses both the mathematical and biological aspects of memory-dependent nerve-pulse propagation. Under the tested settings, the Beta derivative provided the strongest waveform localization, the best amplitude preservation, the lowest energy drift, and the shortest runtime. The Liouville-Caputo formulation was useful for representing strong long-memory effects but was more sensitive to numerical stiffness, while the Atangana-Baleanu formulation offered smoother fading-memory behavior at higher computational cost. The PINN results further showed that fractional parameters can be estimated adaptively while preserving the governing physics. These findings support the Beta derivative as a promising operator for efficient and biologically interpretable soliton-like nerve modeling, while also recognizing the continued relevance of the other operators for specific memory regimes.

Future work should extend the framework beyond the present one-dimensional and deterministic setting. Important directions include multidimensional nerve-bundle geometries, anisotropic conductivity, heterogeneous axon diameter and myelin thickness, stochastic ion-channel fluctuations, and explicit coupling with Hodgkin-Huxley-type gating dynamics. Direct validation against electrophysiological recordings, such as compound action potentials, patch-clamp measurements, microelectrode-array data, or voltage-sensitive dye imaging, will also be essential. In addition, implementation on neuromorphic or FPGA-based platforms could test whether the computational efficiency of the Beta derivative translates into real-time biomedical applications. Finally, hybrid fractional operators that combine the long-memory strengths of Liouville-Caputo formulations with the bounded-memory advantages of Beta-type derivatives may provide a flexible route for modeling both chronic and transient nerve-conduction phenomena.