This paper proposes an efficient discretization method for a novel fractional-order adaptive exponential (FrAdEx) integrate-and-fire model to study fractional-order dynamics of neuronal activity. The discretization is based on an extension of the L1-type method and accurately handles the model's exponential growth and spiking mechanisms. The new method employs an implicit scheme with adaptive time-stepping to robustly handle the stiff systems caused by exponential terms. The implicit nonlinear system can be solved exactly without iterative methods, improving efficiency while maintaining accuracy. The paper provides a complete numerical error model for the discretization scheme, which can be extended to other integrate-and-fire models with minor modifications. The numerical method is rigorously validated and used to investigate various spiking oscillations of the model. The study reveals that fractional-order models can predict biophysical activity, with phase portraits describing transitions between different firing types. This simple model demonstrates tremendous potential, possessing sufficient dynamical expressiveness to qualitatively reproduce multiple features from a biophysical dynamics perspective.
This research aims to solve three core problems in neuronal dynamics modeling:
The authors are driven by:
The core task is to numerically solve the following fractional-order impulsive differential equation system:
FrAdEx Model (non-dimensionalized form):
\frac{d^{\alpha_1} V}{dt^{\alpha_1}} = I - (V - E_L) + \exp(V) - w \\ \tau_w \frac{d^{\alpha_2} w}{dt^{\alpha_2}} = a(V - E_L) - w \end{cases}$$ **Reset Conditions**: $$\text{if } V > V_{peak} \text{ then } \begin{cases} V \leftarrow V_r \\ w \leftarrow w + b \end{cases}$$ Where: - $V(t)$: membrane potential - $w(t)$: adaptation variable - $0 < \alpha_i < 1$: fractional orders - $I(t)$: external current - Parameters: $(g_L, E_L, \Delta_T, V_T, \tau_w, a, b, V_r, V_{peak})$ **Challenges**: 1. Caputo fractional derivative defined as an integral operator with weakly singular kernel 2. Exponential terms causing stiff systems 3. Unknown state-dependent spike times $\{\tau_m\}$ ### Model Architecture #### 1. Piecewise Definition of Caputo Fractional Derivative For piecewise absolutely continuous functions $y \in PAC([0,T];\mathbb{R}^2)$, the piecewise Caputo derivative is defined as: $${}^{PC}D^{\alpha}_{0+}[y](t) = \frac{1}{\Gamma(1-\alpha)} \left[\sum_{j=0}^{m-1} \int_{\tau_j}^{\tau_{j+1}} \frac{y'(s)}{(t-s)^{\alpha}} ds + \int_{\tau_m}^{t} \frac{y'(s)}{(t-s)^{\alpha}} ds\right]$$ where $t \in (\tau_m, \tau_{m+1}]$ and $\tau_m$ are spike times. #### 2. L1-Type Discretization Using linear interpolation on each interval $[t_n, t_{n+1}]$: $$y(s) \approx \frac{t_{n+1}-s}{t_{n+1}-t_n}y_n^+ + \frac{s-t_n}{t_{n+1}-t_n}y_{n+1}^-$$ Obtaining the discrete scheme: $$\sum_{k=0}^{n} d_{n+1,k} \odot \frac{y_{k+1}^- - y_k^+}{\Delta t_k} = f(t_{n+1}, y_{n+1}^-)$$ Weight coefficients: $$d_{n+1,k} = \frac{(t_{n+1}-t_k)^{1-\alpha} - (t_{n+1}-t_{k+1})^{1-\alpha}}{\Gamma(2-\alpha)}$$ #### 3. Exact Solution via Lambert W Function Reformulating the implicit equation as: $$\hat{V}^- + c_2 = c_3 \exp(\hat{V}^-)$$ where $c_2, c_3$ are known coefficients. The solution is: $$\hat{V}^- = -c_2 - W[-c_3 \exp(-c_2)]$$ $$\hat{w}^- = c_0 \hat{V}^- + c_1$$ where $W(\cdot)$ is the Lambert W function, computable directly without iteration. #### 4. Adaptive Time-Stepping Based on error indicator: $$\chi_{n+1} = \|\Gamma(1+\alpha)\| \frac{(t_{n+1}-t_n)^{\alpha}}{t_{n+1}^{\alpha} - t_n^{\alpha}} \frac{\|y_{n+1}-y_n\|}{\|y_n\|}$$ Normalized to $[0,1]$ interval: $$\chi_{n+1} = \frac{\hat{\chi}_{n+1} - \chi_{min}}{\chi_{max} - \chi_{min}}$$ **Adaptive Strategy**: - If $0 < \chi_{n+1} < 1$: maintain step size $\Delta t_{n+1} = \theta \Delta t_n$ - If $\chi_{n+1} < 0$: increase step size $\Delta t_{n+1} = \rho \Delta t_n$ - If $\chi_{n+1} > 1$: decrease step size $\Delta t_{n+1} = \sigma \Delta t_n$ #### 5. Spike Time Estimation When $\hat{V}_{n+1}^-$ becomes complex (Lambert W function argument outside real domain), solve: $$c_3(\Delta t^*) \exp(-c_2(\Delta t^*) + 1) = 1$$ to obtain maximum allowable time step $\Delta t_{Lambert}$, ensuring real-valued solutions. ### Technical Innovations 1. **Lambert W Exact Solution**: - Avoids Newton-Raphson iteration - Direct analytical solution computation - Improved computational efficiency 2. **Integration of Adaptive Time-Stepping with Spike Handling**: - Automatic mesh refinement in exponential growth regions - Lambert W function constraints ensure numerical stability - First-order spike time estimation 3. **Error Analysis for Piecewise L1 Method**: - First analysis for fractional-order systems with state-dependent impulses - Proved global error is $O(\Delta t_{max})$ - Jump terms exactly cancel in error analysis 4. **Universal Design**: - Method applicable to general integrate-and-fire models - Minimal assumptions on reset conditions - Easy extension to other models ## Experimental Setup ### Datasets This paper does not involve traditional datasets but validates through numerical simulations: 1. **PIF Model**: Known analytical solution for convergence verification - Parameters: $C=100$ pFms$^{\alpha-1}$, $I=160$ pA, $V_{peak}=0$ mV, $V_r=-48$ mV - Test $\alpha \in \{0.5, 0.75, 0.95\}$ 2. **LIF Model**: Validates adaptive algorithm - Additional parameters: $g_L=3$ nS, $E_L=-50$ mV - $\alpha = 0.85$ 3. **FrAdEx Model**: Complete model validation - Multiple parameter configurations simulating different neuronal firing patterns - $\alpha \in [0.9, 0.999]$ ### Evaluation Metrics 1. **Relative $\ell^2$ Error**: $$E(x, x_{ref}) = \frac{\|x - x_{ref}\|_2}{\|x_{ref}\|_2}$$ 2. **Spike Time Error**: Comparing numerical spike times $\{t_n\}$ with exact/reference spike times $\{\tau_m\}$ 3. **Convergence Order**: Verifying $O(\Delta t_{max})$ convergence through log-log plots 4. **Computational Efficiency**: - Relationship between iteration count and computation time - Adaptive vs. fixed step-size efficiency comparison ### Comparison Methods 1. **Fixed Step-Size L1 Method**: Baseline for comparing adaptive method efficiency 2. **Exact Solution** (PIF model): Validates numerical method correctness 3. **Self-Convergence** (FrAdEx model): Using very fine mesh as reference solution ### Implementation Details 1. **Programming Environment**: Python + numpy + scipy 2. **Open-Source Library**: pycaputo (fractional calculus library developed by authors) 3. **Adaptive Parameters**: - Safety factor: $\theta = 1.0$ - Reduction factor: $\sigma = 0.5$ - Increase factor: $\rho \in [1.5, 2.0]$ - Minimum step size: $\Delta t_{min} = 10^{-5}$ - Initial step size: $\Delta t_0 = 10^{-2}$ 4. **Error Bounds**: - $\chi_{min} = \{2^{-k} | k=0,...,7\}$ - $\chi_{max} = \{2^{1-k} | k=0,...,7\}$ ## Experimental Results ### Main Results #### 1. PIF Model Convergence (Figure 3) - **Test Configuration**: $\alpha \in \{0.5, 0.75, 0.95\}$, $\Delta t \in \{10^{-2}, 5\times10^{-3}, 10^{-3}, 5\times10^{-4}\}$ - **Result**: All fractional orders achieve **first-order convergence** $O(\Delta t_{max})$ - **Verification**: Compared with exact analytical solution, error decreases from $10^{-3}$ to $10^{-5}$ - **Spike Count**: 6 spike times, all spike positions accurately estimated #### 2. LIF Model Adaptivity (Figure 4) - **Configuration**: $\alpha=0.85$, $\chi_{max} \in \{2^2, 2^{-2}, 2^{-6}\}$ - **Observations**: - In smooth regions, time step increases stepwise (from $10^{-5}$ to $10^{-1}$) - Near spikes, step size dramatically decreases to capture exponential growth - Stricter error bounds ($\chi_{max}=2^{-6}$) result in smaller average step sizes - **Efficiency**: Adaptive method reduces computation while maintaining accuracy #### 3. FrAdEx Model Self-Convergence (Figure 5) - **Configuration**: $\alpha=0.9$, evolution to $T=50$ (non-dimensional), 5 spikes - **Convergence**: - Global error exhibits clear first-order convergence - Error at each spike time decreases linearly with $\Delta t_{max}$ - Subsequent spike errors do not accumulate deterioration - **Error Range**: From $10^{-1}$ (coarse mesh) to $10^{-3}$ (fine mesh) #### 4. Computational Efficiency Comparison (Figure 6) - **Asymptotic Complexity**: - Both adaptive and fixed step-size methods are $O(N^2)$ (inherent complexity of fractional-order memory terms) - Experiments verify theoretical quadratic scaling - **Efficiency Advantage**: - Achieving $10^{-3}$ relative error: - Adaptive method: ~0.1 seconds - Fixed step-size: ~1 second (**10-fold difference**) - Adaptive method significantly reduces computation time at same accuracy ### Neuronal Firing Pattern Reproduction (Section 8) #### Experiment Group 1 (Figure 7): Parameter Set 1 - **$\alpha=0.999$**: **Chattering** - Manifests as densely spaced spike clusters - Phase portrait shows tight limit cycles - **$\alpha=0.98$**: **Fast spiking with broad spike afterpotential (SAP)** - Small curvature afterpotential following spikes - Phase portrait shows looser trajectories - **$\alpha=0.93$**: **Tonic spiking with sharp SAP** - Membrane potential increases monotonically after rapid repolarization - No adaptation, regular firing #### Experiment Group 2 (Figure 8): Parameter Set 2 - **$\alpha=0.999$**: **Tonic firing with broad SAP** - Regular action potential discharge - **$\alpha=0.98$**: **Regular spiking with sharp SAP** - **$\alpha=0.93$**: **Spike frequency adaptation** - Initial spike intervals short, gradually increasing - Phase portrait shows spiral convergence pattern #### Experiment Group 3 (Figure 9): Parameter Set 3 - **$\alpha=0.999$**: **Intrinsic bursting** - Initial bursting followed by regular single spikes - **$\alpha=0.95$**: **Tonic spiking with sharp SAP** - **$\alpha=0.9$**: **Regular spiking with frequency adaptation** ### Key Findings 1. **Fractional Order Impact on Firing Patterns**: - $\alpha \to 1$: Closer to integer-order behavior, exhibiting complex bursting and chattering - Decreasing $\alpha$: Firing patterns become more regular, showing stronger memory effects - Fractional order serves as control parameter, adjusting neuronal firing type 2. **Numerical Method Robustness**: - Maintains stability across all parameter configurations - Accurately captures phase portrait transitions - Handles multi-timescale dynamics from fast spiking to slow adaptation 3. **Biophysical Significance**: - FrAdEx model qualitatively reproduces multiple known neuronal firing patterns - Phase portraits clearly show transitions between firing types - Simple model possesses sufficient expressiveness ## Related Work ### Development of Integrate-and-Fire Models 1. **Classical Models**: - Lapicque (1907): Leaky integrate-and-fire (LIF) model - Izhikevich (2003): Quadratic integrate-and-fire model - Fourcaud-Trocmé et al. (2003): Exponential integrate-and-fire model 2. **AdEx Model**: - Brette & Gerstner (2005): First proposed AdEx model - Naud et al. (2008): Systematic study of AdEx firing patterns - Touboul & Brette (2008): AdEx dynamics and bifurcation analysis ### Fractional-Order Neuronal Models 1. **Fractional-Order LIF**: - Teka et al. (2014): First proposed fractional-order LIF model, studying spike-timing adaptation - Weinberg & Santamaria (2017): History-dependent neuronal activity 2. **Fractal Derivative Extensions**: - Souza et al. (2024): Extended AdEx model using local fractal derivatives - This paper's distinction: Uses Caputo fractional derivative, providing more rigorous mathematical framework ### Numerical Methods for Fractional-Order Differential Equations 1. **L1 Method**: - Li & Zeng (2015): Classical L1 method monograph - Li & Cai (2019): Theory and numerical approximation of fractional integrals and derivatives 2. **Non-Uniform Grid Methods**: - Li et al. (2017): Higher-order numerical methods - Yang & Zeng (2023): Modified L1 method 3. **Adaptive Methods**: - Jannelli (2020): Adaptive procedures for fractional differential equations - This paper's extension: Combines adaptive methods with state-dependent impulses ### Impulsive Differential Equations 1. **Fixed-Time Impulses**: - Wang et al. (2016): Survey of impulsive fractional differential equations - Extensive existing research 2. **State-Dependent Impulses**: - Lakshmikantham et al. (1994): Comparison principles for variable-time impulsive differential equations - **This paper's contribution**: First theory and numerical analysis for fractional-order systems with state-dependent impulses ### Advantages of This Work 1. **Mathematical Rigor**: Uses Caputo derivative rather than local fractal derivative 2. **Numerical Efficiency**: Lambert W exact solution without iteration 3. **Complete Error Analysis**: First for fractional-order systems with state-dependent impulses 4. **Generality**: Method extensible to other integrate-and-fire models ## Conclusions and Discussion ### Main Conclusions 1. **Model Capability**: - FrAdEx model successfully reproduces multiple neuronal firing patterns (chattering, tonic spiking, frequency adaptation, bursting, etc.) - Fractional order $\alpha$ serves as control parameter adjusting firing behavior - Smaller $\alpha$ enhances memory effects, leading to more regular firing patterns 2. **Numerical Method**: - Proposed L1-type method achieves first-order convergence in all test cases - Lambert W solving strategy significantly improves computational efficiency - Adaptive time-stepping successfully handles multi-timescale dynamics 3. **Theoretical Contributions**: - First complete error analysis for fractional-order systems with state-dependent impulses - Proved global error is $O(\Delta t_{max})$ - Method framework extensible to other integrate-and-fire models ### Limitations 1. **Convergence Order Restriction**: - Current method achieves first-order accuracy - Extending to higher-order methods for fractional differential equations requires significant effort - Higher-order spike time estimation for exponential growth models remains unclear 2. **Spike Accumulation Error**: - Theoretical analysis assumes spike count $m$ is not too large - Estimates may fail when $m \max(y(\tau_j^+) - y(\tau_j^-)) = O(\Delta t_{max}^{-1})$ - Long-time evolution may require additional considerations 3. **Parameter Selection**: - Adaptive algorithm parameters $(\chi_{min}, \chi_{max}, \theta, \sigma, \rho)$ require empirical adjustment - Lack robust methods for automatic parameter selection 4. **Model Validation**: - Qualitatively reproduces firing patterns but lacks quantitative comparison with experimental data - Requires more biophysical validation 5. **Computational Complexity**: - $O(N^2)$ complexity remains expensive for long-time simulations - While Fourier methods could improve to $O(N\log N)$, applicability to non-uniform discontinuous systems remains unclear ### Future Directions 1. **Higher-Order Methods**: - Develop second-order or higher-order numerical methods for fractional integrate-and-fire models - Explore higher-order spike time estimation in exponential growth regions 2. **Neuronal Networks**: - Extend to coupled FrAdEx neuronal networks - Study synchronization and collective behavior in networks - Address Lambert W solving challenges in coupled cases 3. **Theoretical Analysis**: - Complete dynamics and bifurcation analysis - Classification of firing patterns across parameter regions - Theoretical characterization of fractional order's impact on neuronal dynamics 4. **Experimental Validation**: - Quantitative comparison with real neuronal recording data - Parameter fitting and model selection - Determine biologically reasonable $\alpha$ value ranges 5. **Algorithm Improvements**: - Develop automatic parameter selection strategies - Explore fast algorithms (e.g., fast Fourier transform) on non-uniform grids - Parallel implementation ## In-Depth Evaluation ### Strengths 1. **Strong Innovation**: - First systematic study of numerical methods for fractional-order AdEx models - Novel and efficient Lambert W exact solving strategy - First complete error analysis for fractional-order systems with state-dependent impulses 2. **Theoretical Rigor**: - Complete mathematical framework (Theorems 6.1 and 6.2) - Rigorous error estimates and convergence proofs - Clear assumptions and applicability conditions 3. **Sufficient Experiments**: - Systematic validation on three models (PIF, LIF, FrAdEx) - Comprehensive testing of convergence, adaptivity, and efficiency - Multiple parameter configurations demonstrating neuronal firing patterns 4. **Practical Value**: - Open-source implementation (pycaputo library) - Method extensible to other models - Detailed algorithm pseudocode (Algorithm 1) 5. **Clear Writing**: - Well-structured, logical flow - Accurate mathematical notation - Rich, informative figures and tables ### Weaknesses 1. **Missing Higher-Order Methods**: - Only first-order method implemented, potentially insufficient accuracy for some applications - Feasibility of higher-order methods not explored 2. **Insufficient Biological Validation**: - Only qualitatively reproduces firing patterns - Lacks quantitative comparison with experimental data - No discussion of biological significance of $\alpha$ 3. **Insufficient Parameter Sensitivity Analysis**: - Adaptive algorithm parameter selection lacks systematic guidance - Parameter space not thoroughly explored 4. **Computational Efficiency**: - $O(N^2)$ complexity remains challenging for large-scale network simulations - Fast algorithm possibilities not explored 5. **Limited Theoretical Analysis**: - Spike accumulation error conditions are restrictive - No dynamics and bifurcation theoretical analysis provided ### Impact 1. **Field Contributions**: - Provides important tools for fractional-order neuronal modeling - Advances research on fractional-order systems with state-dependent impulses - Provides methodological framework for other integrate-and-fire models 2. **Practical Value**: - Open-source code promotes reproducibility and application - Method directly applicable to neuroscience research - Provides new modeling tools for computational neuroscience 3. **Reproducibility**: - Detailed algorithm descriptions - Open-source implementation (pycaputo) - Clear parameter specifications 4. **Potential Applications**: - Neuronal network simulation - Brain disease modeling (e.g., epilepsy) - Neuromorphic computing ### Applicable Scenarios 1. **Ideal Scenarios**: - Research requiring neuronal memory effect modeling - Multi-timescale neuronal dynamics studies - Small to medium-scale neuronal network simulations 2. **Inapplicable Scenarios**: - Applications requiring extremely high precision (current first-order method) - Large-scale network real-time simulation ($O(N^2)$ complexity) - Quantitative research requiring precise experimental data fitting 3. **Potential Extensions**: - Other integrate-and-fire model types (e.g., Izhikevich model) - Stochastic fractional-order neuronal models - Time-varying parameter systems ## Key References 1. **Brette & Gerstner (2005)**: Adaptive exponential integrate-and-fire model - Original AdEx model paper 2. **Teka et al. (2014)**: Fractional leaky integrate-and-fire model - Fractional-order LIF model 3. **Li & Zeng (2015)**: Numerical methods for fractional calculus - Classical textbook on fractional calculus numerical methods 4. **Jannelli (2020)**: Adaptive procedure for fractional differential equations - Adaptive time-stepping methods 5. **Fečkan et al. (2012)**: Impulsive fractional differential equations - Theoretical foundations of impulsive fractional differential equations --- **Overall Assessment**: This is a high-quality computational neuroscience paper making important contributions to fractional-order neuronal modeling and numerical methods. It features rigorous theory, sufficient experiments, and strong practical value. Main limitations include current implementation of only first-order methods and lack of quantitative biological validation. This work establishes important foundations for fractional-order integrate-and-fire model research, possessing significant academic value and application potential.