Splitting physics-informed neural networks for inferring the dynamics of
  integer- and fractional-order neuron models
We introduce a new approach for solving forward systems of differential
equations using a combination of splitting methods and physics-informed neural
networks (PINNs). The proposed method, splitting PINN, effectively addresses
the challenge of applying PINNs to forward dynamical systems and demonstrates
improved accuracy through its application to neuron models. Specifically, we
apply operator splitting to decompose the original neuron model into
sub-problems that are then solved using PINNs. Moreover, we develop an $L^1$
scheme for discretizing fractional derivatives in fractional neuron models,
leading to improved accuracy and efficiency. The results of this study
highlight the potential of splitting PINNs in solving both integer- and
fractional-order neuron models, as well as other similar systems in
computational science and engineering.
Simin Shekarpaz, Fanhai Zeng, George Karniadakis
Splitting physics-informed neural networks for inferring the dynamics of integer- and fractional-order neuron models

###### Abstract

We introduce a new approach for solving forward systems of differential equations using a combination of splitting methods and physics-informed neural networks (PINNs). The proposed method, splitting PINN, effectively addresses the challenge of applying PINNs to forward dynamical systems and demonstrates improved accuracy through its application to neuron models. Specifically, we apply operator splitting to decompose the original neuron model into sub-problems that are then solved using PINNs. Moreover, we develop an \(L^{1}\) scheme for discretizing fractional derivatives in fractional neuron models, leading to improved accuracy and efficiency. The results of this study highlight the potential of splitting PINNs in solving both integer- and fractional-order neuron models, as well as other similar systems in computational science and engineering.

operator splitting, neuron models, fractional calculus. 1

simin_shekarpaz@brown.edu

fanhai_zeng@edu.cn

## 1 Introduction

The human brain is a complex system that involves the interactions of billions of neurons. Mathematical models can be used to simulate the neuronal activity in the brain as a system of differential equations, allowing researchers to better understand how the brain works. Studies related to spiking neurons are performed numerically or biophysically. In numerical studies, the main goal is to solve neural equations and investigate how the dynamic behavior changes for different inputs. Biophysical approaches focus on interpreting the dynamic behavior of spiking neurons according to available experimental observations [2, 58].
Another interesting aspect of spiking neuron models is that they can be formulated as fractional-order equations, which take into account long-term memory. The order of the derivative in these equations can affect the neuron's response [54, 59, 63], making this an important area of research. Recent works in both integer- and fractional-order neuron models are discussed in Section 3.3.

In this work we introduce a new approach for solving neuron models that combines operator splitting methods with physics-informed neural networks (PINNs). Operator splitting methods have been successfully applied in various fields of physics and engineering [6, 11, 15, 18, 23, 34, 51, 52], while PINNs provide a powerful tool for approximating the solution of differential equations. A general introduction to the splitting method can be found in [25, 43].

PINNs were first introduced by Raissi et al. [50]. In this method, the solution of a differential equation is approximated using a neural network, and the parameters of the network are determined by solving a minimization problem that includes residual functions at collocation points, as well as initial and boundary conditions.

PINNs have been applied successfully to a broad range of ordinary and partial differential equations, including fractional equations [48], integro-differential equations, stochastic partial differential equations [68], and inverse problems [44]. There have also been several extensions to the original PINN, such as Fractional PINN (FPINN) [48], physics-constrained neural networks (PCNN) [36, 69], variable hp-VPINN [30], conservative PINN (CPINN) [29], Bayesian PINN [66], parallel PINN [55], Self-Adaptive PINN [42], and Physics informed Adversarial training (PIAT) [53]. Innovations in activation functions, gradient optimization techniques, neural network structures, and loss function structures have driven recent advances in the field. Despite these advances, improvements are still possible, especially concerning unresolved theoretical and practical issues.

Our study makes two important contributions to the field of neural modeling. First, we propose a new method, called the splitting PINN, that employs the operator splitting technique to decompose the original spiking neuron model into sub-problems, which are then solved using PINNs. We demonstrate the effectiveness and accuracy of this method by applying it to integer- and fractional-order neuron models with oscillatory responses, for which vanilla PINN and FPINN formulations fail to predict the solutions. Second, we introduce a novel \(L^{1}\)-scheme for discretizing fractional derivatives in fractional neuron models, which leads to improved accuracy and efficiency in solving these complex models. Our results show that the combination of the splitting PINN method and the \(L^{1}\)-scheme accurately solves fractional neuron models and provides valuable insights into the underlying mechanisms of neural activity.

This paper is organized as follows: Section 2 provides an overview of the proposed method for a given system of differential equations. Section 3 introduces various neuron models and their properties, including the Leaky Integrate-and-Fire (LIF), Izhikevich, Hodgkin-Huxley (HH), and the fractional order Hodgkin-Huxley (FO-HH) models. The efficiency and accuracy of splitting PINNs and FPINNs are demonstrated in Section 4
by applying our algorithm to different neuron models. Finally, we present a discussion of the main results in Section 5. The results and conclusion of this paper will provide valuable insights into the potential of this new approach for solving neuron models and other similar systems in computational science and engineering.

## 2 Problem setup and solution methodology

Let us consider the general form of a nonlinear system of differential equations as follows,

\[\frac{dx_{i}}{dt}=f_{i}(t,x),\qquad i=1,\cdots,n,\qquad t\in[0,T]. \tag{2.1}\]

Before introducing the splitting PINN for solving the given systems of differential equations, we briefly review what splitting methods are, in general.

### Splitting method

For solving the given system (2.1), let us rewrite it as follows,

\[\frac{dx}{dt}=f(t,x(t)),\qquad x(0)=x_{0}\qquad x\in\mathbb{R}^{n},\]

then by using splitting method, \(x\) is decomposed into \((x^{*},x^{**})\), where \(x^{*},x^{**}\in\mathbb{R}^{d}\)\((d<n)\). So, we have

\[\frac{dx^{*}}{dt}=f(t,x^{*}(t)),\qquad x^{*}(0)=x_{0}^{*}, \tag{2.2}\]

\[\frac{dx^{**}}{dt}=f(t,x^{**}(t)),\qquad x^{**}(0)=x_{0}^{**}, \tag{2.3}\]

and \(x_{0}^{*}\) and \(x_{0}^{**}\) are the given vectors of initial conditions. After solving the above sub-systems by the given practical algorithms, we denote the solutions of (2.2) and (2.3) as

\[\phi_{\Delta t}^{(*)}x_{0}^{*}\qquad\text{and}\qquad\phi_{\Delta t}^{(**)}x_{ 0}^{**},\]

where \(\phi_{\Delta t}^{(*)}\) and \(\phi_{\Delta t}^{(**)}\) are \(x^{*}\)-flow and \(x^{**}\)-flow, respectively, and \(\Delta t\) is the step size. Then, by combining the solutions of sub-systems, the approximation operator in \(x\)-flow can be written as follows

\[\phi_{\Delta t}^{(*)}\circ\phi_{\Delta t}^{(**)}\qquad\text{or}\qquad\phi_{ \Delta t/2}^{(**)}\circ\phi_{\Delta t}^{(*)}\circ\phi_{\Delta t/2}^{(**)} \tag{2.4}\]

and \(\phi_{\Delta t}^{(*)}\) and \(\phi_{\Delta t}^{(**)}\) are interchangable. The first splitting is the Lie-splitting method [62], which is first-order, and the second is the Strang splitting method [56], which is a second-order method.
On interval \([0,T]\), we first split the original system into sub-systems (sub-problems) for each sub-interval \([t^{j},t^{j+1}]\)\((j=0,1,\cdots,J-1,\ t^{l}=T)\). Then the solution of the original system at time \(t^{j+1}\) can be approximated as follows

\[x(t^{j+1})=\phi^{(*)}_{\Delta t}\circ\phi^{(**)}_{\Delta t}x(t^{j})\]

where \(x(t^{j})\) is the accurate solution at \(t=t^{j}\).

### Physics informed neural network

Consider an initial value problem as follows,

\[\begin{split}\frac{dx}{dt}&=f(t,x(t)),\qquad x\in \Omega\subseteq\mathbb{R}^{n},\ t\in[0,T],\\ x(0)&=x_{0},\end{split} \tag{2.5}\]

where \(f\) is a nonlinear differential operator, and \(x\) is the unknown solution with known initial condition.

By using the PINN framework, the solution of the above equation is approximated by a fully connected neural network \(\mathcal{N}^{L}\) with \(L\) layers and \(N\) neurons, where the output of \(l\)-th hidden layer is defined as follows

\[\mathcal{N}^{l}(t)=W^{l}\sigma(\mathcal{N}^{l-1}(t))+b^{l},\qquad 2\leq l \leq L, \tag{2.6}\]

\(t\in\mathbb{R}\) is the input vector, \(\sigma(\cdot)\) is the activation function, \(\sigma(\mathcal{N}^{l-1}(t))\in\mathbb{R}^{N}\)\(\{W^{l}\in\mathbb{R}^{N\times N},b^{l}\in\mathbb{R}^{N}\}\) are the network parameters, \(\mathcal{N}^{1}(t)=W^{1}t+b^{1}\) and \(\mathcal{N}^{L}(t)\) is the output of the last layer which is used to approximate the solution. The unknown parameters can be learned by solving a minimization problem that consists of the residual error terms as follows

\[\min_{w,b}\qquad\frac{1}{N_{r}}\sum_{i=1}^{N_{r}}|\frac{dx}{dt}(t^{i})-f(t^{i },x(t^{i}))|^{2}+|x(0)-x_{0}|^{2}, \tag{2.7}\]

where \(N_{r}\) is the number of collocation points. The parameters are randomly initialized and optimized, and then the approximate solution is obtained. The framework uses automatic differentiation to calculate the derivatives of the solution, which eliminates the need for manual calculation or numerical discretization [5]. This capability is available in popular deep-learning frameworks such as TensorFlow and PyTorch.
Figure 1: **Overview of the Splitting PINN :** We first split the original system into subsystems (sub-problems) for each sub-interval \([t^{j},t^{j+1}]\). For each sub-interval, \(x(t^{j})\) is known, the sub-systems are solved using PINN and then the solutions are combined to obtain the approximate solution \(x(t^{j+1})\). To evaluate the error, we obtain the reference solution, \(x_{exact}(t^{j+1})\), by using a high-order numerical solver (for more details, see Appendix B). The algorithm proceeds until arriving at a given accuracy for each sub-interval.
## 3 Neuron models

Until now, various neuron models have been presented for the biological simulations of different parts of the brain. Among them, we can mention the leaky integrate-and-fire (LIF) model, the Izhikevich model, the Hodgkin-Huxley (HH) model, and FitzHugh-Nagumo (FHN) model, which model the membrane behavior [2]. These models can be classified into integer-and fractional-order models. Integer-order models can capture complex phenomena in the neuron system. However, they represent only one type of firing characteristic for constant parameters of the model. On the other hand, fractional-order models can exhibit different dynamic behavior of neurons for constant parameters [58]. This makes fractional-order models more versatile and capable of capturing a wider range of neuron behavior.

### Neuron models: integer order

#### 3.1.1 IF and LIF models

The neuron model of integrate-and-fire (IF model) is one of the best models due to the simplicity of calculations and closeness to human biological conditions. This model is a simplified version of the HH model, which is described by an equation and an assumption. Unlike the HH model, the IF model does not automatically generate an action potential. This model can be determined by

\[C_{m}\frac{dV}{dt}=I(t), \tag{3.1}\]

with the following spike condition: if \(V=V_{th}\), a spike at \(t_{spike}\) is generated and the membrane potential \(V(t)\) is set to \(V_{rest}\) for a refractory period \(\tau_{ref}\)[20, 21]. \(C_{m}\) is the membrane capacitance, \(V_{th}\) is the voltage threshold and \(V_{rest}\) is the resting membrane potential.

A generalized type of IF model is the leaky integrate-and-fire (LIF) model, which adds a leak to the membrane potential. This model is defined by the following equation,

\[\tau\frac{dV}{dt}=-(V-V_{rest})+RI(t), \tag{3.2}\]

where \(\tau=RC_{m}\) is the membrane time constant and \(R\) is the membrane resistance. Because of the important properties of this model, like its computational simplicity [27], accuracy in terms of the spiking behavior and spike times of neurons, and simulating speed [10, 39], this model has become one of the most popular and advantageous neuron models in neuromorphic computing [1, 7, 39, 45]. Also, the characteristic of membrane potential decay over time can be seen in the LIF model.

More complex types of IF models include exponential integrate-and-fire, quadratic integrate-and-fire, and adaptive exponential integrate-and-fire [8].
#### 3.1.2 Izhikevich model

Another neuron model for simulating the membrane behavior is the Izhikevich model. Two important features of this model are computational efficiency and biological plausibility. It reduces the more complex Hodgkin-Huxley model to a 2D system of ordinary differential equations of the form

\[\begin{split}\frac{dv}{dt}&=0.04v^{2}+5v+140-u+I(t), \\ \frac{du}{dt}&=a(bv-u),\end{split} \tag{3.3}\]

with the auxiliary condition

\[\text{if }v\geq v_{th},\text{ then}\left\{\begin{array}{l}v\gets c\\ u\gets u+d\end{array}\right\} \tag{3.4}\]

with \(u\) and \(v\) being dimensionless variables. The variable \(v\) is the membrane potential of the neuron, and \(u\) is the membrane recovery variable, which accounts for the activation of the \(K^{+}\) ion current and the inactivation of the \(Na^{+}\) ion current and provides negative feedback to the membrane potential. The dimensionless parameters \(a\), \(b\), \(c\), and \(d\) regulate the behavior of the neuron.

The auxiliary condition (3.4) triggers a reset of the neuron when the membrane potential surpasses the threshold, simulating a spike. This model is capable of reproducing the spiking and bursting behavior of neurons in real-time, making it a widely used model in simulations of large-scale neural networks.

#### 3.1.3 Hodgkin-Huxley model

From a biophysical perspective, the nerve cell's action potential is generated by the flow of ions through the cell membrane's ion channels. Hodgkin and Huxley described the dynamics of these membrane currents through a set of coupled differential equations based on their experiments on the giant squid axon [24].

The mechanism of the action potential can be understood with reference to Figure 2. A capacitor, resistor, and transistor were used to simulate the equivalent circuit. Changes in the action potential were observed by applying the current \(I(t)\) and adjusting the capacitance and leakage resistance of the sodium and potassium channels.

The circuit of the Hodgkin-Huxley (HH) model consists of four parallel branches: integrative branch, leaky branch, \(K^{+}\) channel, and \(Na^{+}\) channel. A system of four coupled differential equations was used to describe the membrane potential of a giant squid axon as follows

\[\begin{split}\frac{dV_{m}}{dt}=& F_{1}(t,V_{m},n,m,h), \qquad\frac{dn}{dt}=F_{2}(t,V_{m},n,m,h),\\ \frac{dm}{dt}=& F_{3}(t,V_{m},n,m,h),\qquad\frac{dh}{dt }=F_{4}(t,V_{m},n,m,h),\end{split} \tag{3.5}\]
where

\[F_{1} =\frac{1}{C_{m}}(-g_{L}(V_{m}-E_{L})-g_{K}n^{4}(V_{m}-E_{K})-g_{Na}m^ {3}h(V_{m}-E_{Na})+I(t)),\] \[F_{2} =\alpha_{n}(V_{m}(t))(1-n(t))-\beta_{n}(V_{m}(t))n, \tag{3.6}\] \[F_{3} =\alpha_{m}(V_{m}(t))(1-m(t))-\beta_{m}(V_{m}(t))m,\] \[F_{4} =\alpha_{h}(V_{m}(t))(1-h(t))-\beta_{h}(V_{m}(t))h,\]

and \(g_{Na}\), \(g_{K}\) and \(g_{L}\) are the maximum conductances of the \(Na^{+}\), \(K^{+}\) and leak currents. The variables \(\alpha_{x}\) and \(\beta_{x}\) of the current channel conductances are dependent functions of \(V_{m}(t)\) as follows:

\[\alpha_{n}(V_{m}) =\frac{0.1-0.01(V_{m}-V_{0})}{e^{1-0.1(V_{m}-V_{0})}-1}, \beta_{n}(V_{m}) =0.125e^{-(V_{m}-V_{0})/80},\] \[\alpha_{m}(V_{m}) =\frac{2.5-0.1(V_{m}-V_{0})}{e^{2.5-0.1(V_{m}-V_{0})}-1}, \beta_{m}(V_{m}) =4.0e^{-(V_{m}-V_{0})/80}, \tag{3.7}\] \[\alpha_{h}(V_{m}) =0.07e^{-(V_{m}-V_{0})/20}, \beta_{h}(V_{m}) =\frac{1}{1+e^{3-0.1(V_{m}-V_{0})}}.\]

### Neuron models: fractional order

In recent years, fractional differential equations have been developed to improve the modeling of many biological phenomena, including mechanical properties of viscoelastic tissue [40], the tissue-electrode interface [41], pharmacokinetics of drug delivery and absorption [13, 14, 49], and anomalous calcium sub-diffusion in micro-domains [57].

The main characteristic of using fractional derivatives is their non-locality. This means that the next state of the system depends on the current state and all historical states before it. This advantage makes the study of fractional order systems an active area of research.

Figure 2: Schematic diagram for the Hodgkin-Huxley model [24]. Left: ion channels on the membrane of the neuron. Right: simulated circuit with \(I\) denoting the input current.
#### 3.2.1 Fractional derivative definitions

In this Section, we present the definitions of fractional derivatives. There are different methods for defining fractional derivatives, among which we can mention the Grunwald-Letnikov derivative, Riemann-Liouville derivative, and Caputo derivative [33]. The models in this paper are defined using the Caputo fractional derivative.

**Definition 1**.: The Caputo fractional derivative of the function \(f(t)\) with order \(\alpha>0\) is defined as

\[{}_{C}D_{a,t}^{\alpha}f(t)=\frac{1}{\Gamma(n-\alpha)}\int_{a}^{t}(t-s)^{n- \alpha-1}f^{(n)}(s)ds, \tag{3.8}\]

where \(n-1<\alpha<n\) and \(n\) in a non-negative integer.

If \(a=0\), then we can use \(\frac{d^{\alpha}}{dt^{\alpha}}f(t)={}_{C}D_{0,t}^{\alpha}f(t)\).

#### 3.2.2 \(L^{1}\) scheme to approximate the fractional derivatives

An efficient method for approximating the Caputo derivative of order \(\alpha\)\((0<\alpha<1)\) is the \(L^{1}\) scheme, which was introduced by Oldham and Spanier [47]. Using this method, the Caputo derivative can be approximated by using the following formula

\[\frac{d^{\alpha}}{dt^{\alpha}}f(t^{j})\approx\delta_{t}^{\alpha}f^{j}=\sum_{k =0}^{j-1}b_{k}^{\alpha}\left[f(t^{j-k})-f(t^{j-1-k})\right], \tag{3.9}\]

where for the uniform time mesh \(t^{j}=j\Delta t,j\geq 0\), \(\Delta t\) is the step size. Then, \(b_{k}^{\alpha}\) is given by \(b_{k}^{\alpha}=\frac{(\Delta t)^{-\alpha}}{\Gamma(2-\alpha)}[(k+1)^{1-\alpha} -k^{1-\alpha}]\).

In [35], for a smooth \(f\), the error estimate of the above \(L^{1}\) scheme is

\[|\delta_{t}^{\alpha}f(t^{j})-\frac{d^{\alpha}}{dt^{\alpha}}f(t^{j})|\leq C \Delta t^{2-\alpha}, \tag{3.10}\]

where \(C=C(\alpha,f)\). This approximation has been used in many papers discussing fractional-order spiking neurons. The interested readers are referred to [46, 60, 2, 4]. Here, we will develop this method to discretize the fractional derivative, and then FPINN is used.

#### 3.2.3 Fractional Hodgkin-Huxley model

In fractional-order neural models, the neuron's dynamics depends on the order of the derivative, which can create different types of memory-dependent dynamics. The fractional order Hodgkin-Huxley (FO-HH) model is one of the neuron models that has attracted much attention.

The HH model has two basic problems; the first is that the Dielectric losses in the membrane are neglected. The second problem is that membrane capacity is considered
ideal. To overcome the above problems, a fractional model is proposed. The idea of fractional capacity is taken from Curie's empirical law [65], which can be written as follows

\[I_{c}(t)=C_{m}\frac{d^{q_{1}}V_{m}(t)}{dt^{q_{1}}}, \tag{3.11}\]

where \(V_{m}(t)\) is the excitation voltage, \(I_{c}(t)\) is the current in the capacitor, \(q_{1}\) is the order of differentiation, and \(C_{m}\) is the fractional capacitance [64]. The fractional order model provides a more accurate description based on long-term memory behavior. Motivated by the above discussion, we propose the following FO-HH model

\[\begin{split}\frac{d^{q_{1}}V_{m}}{dt^{q_{1}}}=& F_{1}(t,V_{m},n,m,h),\qquad\frac{d^{q_{2}}n}{dt^{q_{2}}}=F_{2}(t,V_{m},n,m,h), \\ \frac{d^{q_{3}}m}{dt^{q_{3}}}=& F_{3}(t,V_{m},n,m,h),\qquad\frac{d^{q_{4}}h}{dt^{q_{4}}}=F_{4}(t,V_{m},n,m,h), \end{split} \tag{3.12}\]

where \(q=(q_{1},q_{2},q_{3},q_{4})\) is the order of differentiation, and the other parameters are as in (3.5).

By defining the fractional model, we apply the dielectric loss in the membrane, and as a result, we will see the change in the refractory period with the same value of the given current in an integer order case. The refractory period is the time when the membrane is hyperpolarized and, hence, requires a stronger stimulus to produce a smaller action potential.

The modified HH model can be very effective in biomedical applications such as heart health analysis. For example, in an ECG waveform, the PR interval represents a refractory period. PR interval estimation is very important for cardiac diagnosis [12, 22].

### Prior works in neuron models

In [37], the authors compared the spiking rate patterns of five single neuron models, including LIF, Izhikevich, and Hodgkin-Huxley (HH) models, under different sustained current inputs. Numerical stability and accuracy were also considered. The multi-step methods for neuronal modeling, including the HH model, were proposed in [32]. In [31], the modified Khater (mK) method and B-spline scheme were proposed to find numerical solutions of the FitzHugh-Nagumo (FHN) equation, with a focus on finding different types of soliton wave solutions, studying their stability properties, and using them to obtain numerical solutions of the model.

In [19], the Hybrid Functions (HF) method was proposed as a solution for the HH model. The HF method was compared comprehensively with other algorithms, evaluating computational speed, absolute error, and integral time squared error. The finite difference scheme was used in [67] to solve the stochastic FHN model, including stability analysis and the calculation of explicit optimal a priori estimates for the existence of solutions. In [3], numerical solutions for the FHN and Izhikevich neuron models were
obtained using a non-standard finite difference scheme and GL discretization technique. The models were compared, and their behavior was analyzed in different fractional orders.

Fractional order modeling in neural systems is a relatively new area of research. The non-local definitions of fractional calculus used in these models provide a more realistic representation of neural systems and offer a deeper understanding of their behavior. In a study reported in [2], four numerical methods were applied to two fractional-order spiking neuron models, the FO-LIF and FO-HH models. The authors used a finite memory window version of the \(L^{1}\) approximation for comparison with well-known techniques such as the GL-based method, product integration approximation, and the Z-transform approach. In the four methods, the uniform mesh is used, and low-accurate solutions are obtained due to the singularity of the solution of fractional equations. In this paper, they have analyzed the spiking patterns, inter-spike interval adaptation, and steady-state spiking frequency for each numerical method under varying memory lengths. In a related study reported in [58], the authors have used a \(L^{1}\) scheme (linear interpolation based on uniform mesh) to discretize the Caputo fractional derivative. The first-order extrapolation is used to derive the linearized scheme, where the global error is \(O(\Delta t^{\alpha})\), and the error far from the origin is \(O(\Delta t)\). In this paper, they have investigated the effects of non-Markovian power-law voltage-dependent conductances on the generation of action potentials and spiking patterns in a Hodgkin-Huxley model. They used fractional derivatives to implement the slow-adapting power-law dynamics of the potassium and sodium conductance gating variables. The results showed that, with different input currents and derivative orders, a wide range of spiking patterns can be generated, such as square wave bursting, mixed mode oscillations, and pseudo-plateau potentials. These findings suggest that power-law conductances increase the number of spiking patterns a neuron can produce.

In [9], the dynamics and numerical simulations of a fractional-order coupled FHN neuronal model were discussed, and the stability properties of its equilibrium states were analyzed based on theoretical results. In [61], a non-standard finite difference scheme was used to solve the fractional Izhikevich neuron model, and a general formula for the synchronization of different Izhikevich neurons was proposed.

## 4 Results

In this Section, we use splitting PINN to solve the neuron models presented in Section 3, using the network architecture and hyperparameters specified in Table 1. The optimization algorithm used is Adam with a learning rate of 0.0001 and a scheduler. IN addition techniques such as adaptive activation function and feature expansion proposed in [28, 38] were also used to improve the computational efficiency of the proposed method.
The relative \(L^{2}\) norm of errors of the inferred solutions are also computed as follows,

\[\text{relative L}^{2}\text{ error}=\frac{\sqrt{\sum_{j=1}^{l}(V_{exact}(t^{j})-V_{app}(t^{j}))^{2}}}{ \sqrt{\sum_{j=1}^{l}(V_{exact}(t^{j}))^{2}}}.\]

Further details about the training procedure and hyperparameters can be found in the relevant Section.

**Data Validation:** In all of the examples, the reference solutions are obtained by using the adaptive time-stepping spectral collocation method, which is presented in Appendix B.

### LIF model

**PINN implementation:** The approximate solutions of the LIF model can be obtained by using PINNs. For solving this model, we discretize the time domain \([0,T]\) into subdomains \([t^{j},t^{j+1}]\) where \(j=0,1,\cdots,J-1\) and \(t^{j}=T\) and then PINNs are used to solve the model in each sub-interval. More detailed visual evaluations and the numerical results of this model are provided in Appendix A.

\begin{table}
\begin{tabular}{l|c c c c c} \hline Problems & Dep. & Wid. & Act. & Opt. & Iter. \\ \hline _LIF model_ & 5 & 40 & tanh & _Adam_ & 50000 \\ _LIF model with_ & 7 & 60 & tanh & _Adam_ & 10000 \\ _threshold voltage_ & & & & & \\ _Izhikevich model_ & 6,6 & 40,40 & tanh,tanh & _Adam_,_Adamax_ & 20000 \\ _HH model_ & 6,10 & 20,20 & tanh,sin & _Adam_,_Adamax_ & 20000 \\ _(step current function)_ & & & & & \\ _HH model_ & 6,10 & 20,20 & tanh,sin & _Adam_ & 20000 \\ _(constant current)_ & & & & & \\ _FO-HH model_ & 10,6 & 100,100 & tanh,sin & _Adam_ & 70000 \\ \((q_{i}=0.8)\) & & & & & \\ _FO-HH model_ & 10,6 & 100,100 & tanh,sin & _Adam_ & 50000 \\ \((q_{i}=0.6)\) & & & & & \\ _FO-HH model_ & 10,6 & 100,100 & tanh,sin & _Adam_ & 20000 \\ \((q_{i}=0.4)\) & & & & & \\ \hline \hline \end{tabular}
\end{table}
Table 1: PINN architectures and hyperparameters used for training. The first and second subcolumns under “Depth”, “Width”, and “Activation” correspond to the first and second sub-problems, respectively.
### Izhikevich model

**Splitting PINN implementation:** Consider \((u^{j},v^{j})\) to be the numerical solution at \(t=t^{j}=j\Delta t\) with step size, \(\Delta t\), then by using the splitting method, the following sub-problems on each sub-interval \([t^{j},t^{j+1}]\)\((j=0,1,\cdots,J-1,t^{j}=T)\) are obtained

\[\frac{du}{dt}=a(bv^{j}-u),\qquad u(t^{j})=u^{j}, \tag{4.1}\]

and

\[\frac{d\omega}{dt}=0.04v^{2}+5v+140-u^{j+1}+I(t),\qquad v(t^{j})=v^{j}. \tag{4.2}\]

This model is solved in the time interval \([0,100\ ms]\), with 2000 sub-intervals and 20 points per sub-interval. For solving sub-problems by using PINN, the fully connected neural networks are assumed to approximate the solutions, with the given architecture described in Table 1, \(v_{th}=30mV\),

\[I(t)=\begin{cases}0&\text{for }0<t<2,\\ 15&\text{for }t\geq 2,\end{cases}\]

and the choice of parameters given in [26], which can be described as follows:

* \(\mathbf{a}\sim 0.02\ \frac{1}{ms}\) : Time scale of the recovery variable \(u\). Smaller values result in slower recovery.
* \(\mathbf{b}\sim\ 0.2\ [dimensionless]\) : Sensitivity of the recovery variable \(u\) to the sub-threshold fluctuations of the membrane potential \(v\). Greater values couple \(v\) and \(u\) more strongly resulting in possible sub-threshold oscillations and low-threshold spiking dynamics. The case \(b<a\ (b>a)\) corresponds to saddle-node (Andronov-Hopf) bifurcation of the resting state.
* \(\mathbf{c}\sim-50\ mV\) : The after-spike reset value of the membrane potential \(v\) caused by the fast high-threshold \(K^{+}\) conductances.
* \(\mathbf{d}\sim 2\ mV\) : After-spike reset of the recovery variable \(u\) caused by slow high-threshold \(Na^{+}\) and \(K^{+}\) conductances.

\begin{table}
\begin{tabular}{l l l} \hline  & training error & testing error \\ \hline _membrane potential_ & \(0.115719\pm 0.008516\) & \(0.120458\pm 0.008221\) \\ \hline \(u\) & \(0.041201\pm 0.005980\) & \(0.043186\pm 0.005070\) \\ \hline \hline \end{tabular}
\end{table}
Table 2: mean relative \(L^{2}\) norm of errors for the Izhikevich model.
### Hodgkin-Huxley model

**Splitting PINN implementation:**

For solving the HH model by applying Splitting PINN on interval \([0,T]\), we have used the choice of parameters given in Table 3. Consider \((V^{j},n^{j},m^{j},h^{j})\) to be the solution at \(t^{j}=j\Delta t\) with step size \(\Delta t\). Then, the following sub-systems on each sub-interval \([t^{j},t^{j+1}]\) are obtained using the Lie splitting method. The first part is,

\[\frac{dn}{dt}=F_{2}(t,V^{j}_{m},n,m,h),\qquad n(t^{j})=n^{j},\qquad t \in(t^{j},t^{j+1}],\] \[\frac{dm}{dt}=F_{3}(t,V^{j}_{m},n,m,h),\qquad m(t^{j})=m^{j}, \qquad t\in(t^{j},t^{j+1}], \tag{4.3}\] \[\frac{dh}{dt}=F_{4}(t,V^{j}_{m},n,m,h),\qquad h(t^{j})=h^{j}\qquad t \in(t^{j},t^{j+1}].\]

where \(V^{j}_{m}\) is the known solution at time \(t=t^{j}\). The second part is given by

\[\frac{dV_{m}}{dt}=F_{1}(t,V_{m},n^{j+1},m^{j+1},h^{j+1}),\qquad V_{m}(t^{j})=V ^{j}_{m},\qquad t\in(t^{j},t^{j+1}], \tag{4.4}\]

Figure 4: Loss function of the Izhikevich model.

Figure 3: Izhikevich model: Left: membrane potential versus time. Right: recovery variable versus time.
and \((n^{j+1}\),\(m^{j+1}\),\(h^{j+1})\) is the solution obtained from the first part. The splitting approximation of the solution at \(t=t^{j+1}\) is \(\left(V_{m}^{j+1}\),\(n^{j+1},m^{j+1},h^{j+1}\right)\), and the solution for the whole domain is obtained by repeating this process.

Different kinds of input currents, including step current function and constant current, are considered for solving the HH model.

**Step current function:** For solving the HH model with the given step current in Figure 5(a) at the time interval \([0,20ms]\), we have used 800 sub-intervals in the splitting procedure with 40 training points for solving each sub-problem. The hyperparameters and architecture of PINNs for solving the sub-problems can be found in Table 1.

**Constant current:** This system is solved at time interval \([0,100ms]\), with 3000 sub-intervals and 30 points at each sub-interval. With the given hyperparameters and architectures in Table 1, PINN is used for solving each part, and the solutions are combined to obtain the approximate solutions.

The relative \(L^{2}\) norm of errors of approximate solutions are shown in Table 5.

\begin{table}
\begin{tabular}{l c c} \hline  & training error & testing error \\ \hline _membrane potential_ & \(0.001155\pm 0.000111\) & \(0.000886\pm 0.000149\) \\ \hline \(n\) & \(0.002720\pm 5.3\times 10^{-5}\) & \(0.002712\pm 5.3\times 10^{-5}\) \\ \hline \(m\) & \(0.014316\pm 0.000447\) & \(0.014071\pm 0.000451\) \\ \hline \(h\) & \(0.002685\pm 7.3\times 10^{-5}\) & \(0.002680\pm 7.2\times 10^{-5}\), \\ \hline \hline \end{tabular}
\end{table}
Table 4: mean relative \(L^{2}\) norm of errors for HH model with current step function.

\begin{table}
\begin{tabular}{l l l} \hline Parameter & Value & Description \\ \hline \(g_{Na}\) & \(120\)\(mS/cm^{2}\) & _Maximum Na\({}^{+}\) current conductance_ \\ \hline \(gK\) & \(36\)\(mS/cm^{2}\) & _Maximum K\({}^{+}\) current conductance_ \\ \hline \(g_{L}\) & \(0.3\)\(mS/cm^{2}\) & _Maximum leak current conductance_ \\ \hline \(E_{Na}\) & \(50\)\(mV\) & _Na\({}^{+}\) current reversal potential_ \\ \hline \(E_{K}\) & \(-77\)\(mV\) & _K\({}^{+}\) current reversal potential_ \\ \hline \(E_{L}\) & \(-54\)\(mV\) & _Leak current reversal potential_ \\ \hline \(C_{m}\) & \(1\)\(\mu F\cdot s^{q_{1}-1}/cm^{2}\) & _Membrane Capacitance_ \\ \hline \(V_{0}\) & \(-65\) & _Initial membrane potential_ \\ \hline \(m_{0}\) & \(0.0529\) & _Initial Na\({}^{+}\) current activation_ \\ \hline \(n_{0}\) & \(0.3177\) & _Initial Na\({}^{+}\) current inactivation_ \\ \hline \(h_{0}\) & \(0.5960\) & _Initial K\({}^{+}\) current activation_ \\ \hline \hline \end{tabular}
\end{table}
Table 3: Parameter values and their descriptions of the HH neuron model [58].
Figure 5: HH model: comparison of splitting PINN results with the reference solution. (a) membrane potential; (b) activation variable of potassium channel; (c) activation variable of the sodium channel; (d) deactivation variable of the sodium channel. The inset plot shows the step function input current.
Figure 6: HH model: absolute error of the normalized solutions using splitting PINN with a current step function. (a) membrane potential; (b) activation variable of potassium channel; (c) activation variable of the sodium channel; (d) deactivation variable of the sodium channel.

Figure 7: Loss function of the HH model with a current step function.
\begin{table}
\begin{tabular}{l c c} \hline \hline  & training error & testing error \\ \hline _membane potential_ & \(0.009925\pm 0.001957\) & \(0.009434\pm 0.001877\) \\ \hline \(n\) & \(0.005240\pm 0.000666\) & \(0.005229\pm 0.000661\) \\ \hline \(m\) & \(0.030731\pm 0.003931\) & \(0.030364\pm 0.003939\) \\ \hline \(h\) & \(0.007077\pm 0.000899\) & \(0.007066\pm 0.000898\) \\ \hline \hline \end{tabular}
\end{table}
Table 5: mean relative \(L^{2}\) norm of errors for HH model for constant current.

Figure 8: HH model: comparison of splitting PINN results with the reference solution for constant current input (\(I=10\)\(nA/cm^{2}\)). (a) membrane potential; (b) activation variable of potassium channel; (c) activation variable of the sodium channel; (d) deactivation variable of the sodium channel. Note that PINN fails to predict the solution.
### Fractional Hodgkin-Huxley model

**Splitting FPINN implementation:**

This problem is solved using the splitting fractional physics-informed neural network

Figure 10: Loss function of the HH model for constant current.

Figure 9: HH model: absolute error of the normalized solutions using splitting PINN with a constant current (\(I=10\)\(nA/cm^{2}\)). (a) membrane potential; (b) activation variable of potassium channel; (c) activation variable of the sodium channel; (d) deactivation variable of the sodium channel.
(FPINN) approach [48]. In this approach, the memory of the fractional derivative at each sub-interval encompasses the entire domain. So for each sub-interval \([t^{j},t^{j+1}]\)\((0\leq j\leq N_{1}-1)\), by using the developed \(L^{1}\) - scheme, we have the following sub-problems:

\[\delta_{t}^{q_{1}}V_{m}^{N_{2}j+l+1}=F_{1}(t^{N_{2}j+l+1},V_{m}^{N_{2}j+l+1},n, m,h), \tag{4.5}\]

where

\[\delta_{t}^{q_{1}}V_{m}^{N_{2}j+l+1}=\sum_{k=0}^{N_{2}j+l}b_{k}^{q_{1}}\left[V_ {m}(t^{N_{2}j+l+1-k})-V_{m}(t^{N_{2}j+l-k})\right], \tag{4.6}\]

and \(n\), \(m\) and \(h\) are the known solutions at \(t=t^{N_{2}j+l}\). The second subproblem is

\[\delta_{t}^{q_{2}}n^{N_{2}j+l+1} = F_{2}(t^{N_{2}j+l+1},V_{m},n^{N_{2}j+l+1},m,h),\] \[\delta_{t}^{q_{3}}m^{N_{2}j+l+1} = F_{3}(t^{N_{2}j+l+1},V_{m},n,m^{N_{2}j+l+1},h), \tag{4.7}\] \[\delta_{t}^{q_{4}}h^{N_{2}j+l+1} = F_{4}(t^{N_{2}j+l+1},V_{m},n,m,h^{N_{2}j+l+1}).\]

where \(V_{m}\) is the solution of the first sub-problem at \(t=t^{N_{2}j+l+1}\), \(l=0,1,\cdots,N_{2}-1\), and \(N_{2}\) is the number of residual points for each sub-interval.

Then, the solutions of sub-systems are combined to obtain the solution of the original system.

This model is solved with the proposed method for \(q_{i}=0.8,0.6,0.4\) at time interval \([0,100ms]\), with 2000 sub-intervals and 40 residual points in each sub-interval. Each part is solved via FPINN, with the given hyperparameters and network architectures in Table
Figure 11: FO-HH model: comparison of splitting FPINN results with the reference solution for constant current input (\(I=20\ nA/cm^{2}\)) and fractional order \(q_{1}=q_{2}=q_{3}=q_{4}=0.8\). (a) membrane potential; (b) activation variable of potassium channel; (c) activation variable of the sodium channel; (d) deactivation variable of the sodium channel. Note that FPINN fails to predict the solution.

Figure 12: Loss function of the FO-HH model for fractional order \(q_{1}=q_{2}=q_{3}=q_{4}=0.8\).
Figure 14: Loss function of the FO-HH model for fractional order \(q_{1}=q_{2}=q_{3}=q_{4}=0.6\).

Figure 13: FO-HH model: comparison of splitting FPINN results with the reference solution for constant current input (\(I=20\)\(nA/cm^{2}\)) and fractional order \(q_{1}=q_{2}=q_{3}=q_{4}=0.6\). (a) membrane potential; (b) activation variable of potassium channel; (c) activation variable of the sodium channel; (d) deactivation variable of the sodium channel. Note that FPINN fails to predict the solution.
Figure 16: Loss function of the FO-HH model for fractional order \(q_{1}=q_{2}=q_{3}=q_{4}=0.4\).

Figure 15: FO-HH model: comparison of splitting FPINN results with the reference solution for constant current input (\(I=20\)\(nA/cm^{2}\)) and fractional order \(q_{1}=q_{2}=q_{3}=q_{4}=0.4\). (a) membrane potential; (b) activation variable of potassium channel; (c) activation variable of the sodium channel; (d) deactivation variable of the sodium channel. Note that FPINN fails to predict the solution.
\begin{table}
\begin{tabular}{l c c c} \hline \hline  & \(q_{i}=0.8\) & \(q_{i}=0.6\) & \(q_{i}=0.4\) \\ \hline _membrane_ & \(0.084252\pm 0.011138\) & \(0.098086\pm 0.027659\) & \(0.043236\pm 0.001514\) \\ _potential_ & \(0.084460\pm 0.012576\) & \(0.096949\pm 0.028145\) & \(0.043241\pm 0.003771\) \\ \hline \(n\) & \(0.024353\pm 0.005469\) & \(0.023854\pm 0.008417\) & \(0.001909\pm 8.8\times 10^{-5}\) \\  & \(0.024360\pm 0.005527\) & \(0.023776\pm 0.008435\) & \(0.001933\pm 0.000177\) \\ \hline \(m\) & \(0.129660\pm 0.035021\) & \(0.093863\pm 0.031778\) & \(0.006009\pm 9.4\times 10^{-5}\) \\  & \(0.129083\pm 0.035903\) & \(0.093619\pm 0.032427\) & \(0.006063\pm 0.000214\) \\ \hline \(h\) & \(0.036292\pm 0.009040\) & \(0.040947\pm 0.013980\) & \(0.007536\pm 0.000374\) \\  & \(0.036336\pm 0.009121\) & \(0.040921\pm 0.013963\) & \(0.007592\pm 0.000649\) \\ \hline \hline \end{tabular}
\end{table}
Table 6: mean relative \(L^{2}\) norm of errors for FO-HH model for different orders of fractional derivatives. The first and second sub-rows correspond to the training and testing errors, respectively.

Figure 17: Plots of the membrane potential and the corresponding activation and deactivation variables of the FO-HH model with constant current \(I=20\ nA/cm^{2}\) for different orders of fractional derivatives, obtained using the splitting FPINN.
## 5 Discussion

We have presented a deep learning approach for solving nonlinear systems of differential equations. The performance and accuracy of splitting Physics-Informed Neural Networks (PINNs) are studied in the context of solving neuron models. In Section 4.2, the proposed method was used to solve the Izhikevich model. We showed the reference and learned solutions in Figure 3 for equi-distant training and random test datasets. It can be seen that the nonlinear solutions learned by the nonlinear splitting PINN match well with the reference solutions. The relative \(L^{2}\) norm of errors for the solutions are presented in Table 2, which shows the method's capability for this model. The convergence of the loss functions is shown in Figure 4. It can be seen that this function rapidly tends to values lower than \(10^{-5}\) and \(10^{-4}\) for the first and second sub-problems.

After solving the Izhikevich model, we focus on using the Splitting PINN method for the Hodgkin-Huxley model. In this model, applying input current and changes in voltage resulting from the opening and closing of ion channels lead to the generation of spikes. The amount of input current needed to generate a spike is at least \(2.7nA/cm^{2}\), and with increasing input current, the number of spikes increases. In Section 4.3, the voltage is obtained for two states of input current: step function current and constant input current. In the first case, when the current is applied over a given time interval, the voltage increases and produces an action potential (positive peak). After the spike, the potassium channel opens, and the sodium channel closes, causing the voltage to decrease and the neuron to enter a refractory period where the potential is below the resting potential of \(V_{rest}=-65mV\). The voltage then slowly returns to \(-65mV\). The splitting PINN method was used to calculate approximate solutions and compare them with reference solutions. The results are shown in Figure 5, and the absolute errors are given in Figure 6, which demonstrate the accuracy of the method. The relative \(L^{2}\) norm of errors for the current step function is presented in Table 4. The convergence of the loss functions is also shown in Figure 7 and tends to smaller values of \(10^{-6}\) and \(10^{-5}\) for the first and second sub-problems, respectively.

In the case of a constant input current, the approximate and reference solutions are shown in Figure 8. The proposed method's effectiveness is demonstrated by comparing the approximate results with the reference and numerical solutions obtained using PINN. The absolute errors of the solutions are shown in Figure 9, and the plots of the loss functions are displayed in Figure 10, which converge to values of \(10^{-6}\) for sub-problems. The relative \(L^{2}\) norm of errors of the solutions is presented in Table 5. These results demonstrate the capability and accuracy of the proposed algorithm in solving this model.

The fractional-order Hodgkin-Huxley model was also studied in Section 4.4. The system was modeled using Caputo's fractional derivative and simulated using the proposed method. The effect of memory, associated with fractional order, on firing activity was also examined. To address this, the \(L^{1}\) scheme was improved for solving the fractional-order model by proposing using full-domain memory, rather than sub-domain memory, to calculate the fractional derivative. The numerical results obtained using the improved
method for \(q_{i}=0.8,0.6,0.4\) are shown in Figures 11, 13 and 15, where the approximate results are compared with the reference solutions. The plots of the loss functions are also shown in Figures 12, 14, and 16, demonstrating the efficiency and convergence of the current method. The relative \(L^{2}\) norm of errors of the solutions is presented in Table 6, and as shown in the Section, the results are accurate and have converged to the reference solutions.

The voltage responses for the different orders of fractional derivatives under constant input current \(I=20\) are also displayed in Figure 17. The Figure demonstrates that various spike patterns can be produced. As \(q_{i}\) approaches 1, the firing frequency decreases and results in an increase in the number of spikes in the same time period. Furthermore, the first spike occurs at a later time. Following the cessation of the injected current, the memory-dependent spiking activity can also be observed by applying a step function current. Additionally, the regularity of solutions is observed to depend on the order of fractional derivatives in Figure 17. As \(q_{i}\) decreases, the regularity increases. Table 6 shows that, unregular solutions are more sensitive to parameter initialization, leading to a greater standard deviation. To address this issue, the network architecture can be improved.

The results in this section demonstrate high accuracy, and the solutions have converged to the reference values.

## 6 Conclusion

In this study, we introduced a novel method for solving neuron models represented as systems of differential equations. The Splitting PINN algorithm was demonstrated to be effective and accurate through comparisons with reference solutions and the mean \(L^{2}\) relative norm of errors. In addition, the results of the fractional-order Hodgkin-Huxley model highlighted the effect of memory on firing activity and voltage responses, showing that as \(q_{i}\) approaches 1, the firing frequency decreases while the number of spikes in the same time period increases.

This research provides valuable insights into the behavior of neuron membranes and the various spike patterns that can be generated. The performance of the proposed method demonstrates its superiority over the vanilla PINN algorithm in solving complex neuron models represented as systems of differential equations. This study contributes to the development of tools for investigating the behavior of neurons and the underlying mechanisms of neural activity.

In conclusion, the proposed Splitting PINN algorithm is a promising method for solving systems of differential equations and provides valuable insights into the behavior of neurons and their underlying mechanisms.
## Acknowledgments

We would like to thank Dr. Khemraj Shukla and Dr. Ehsan Kharazmi for helpful discussions. This work was supported by AFOSR MURI funding (FA9550-20-1-0358) and Graphs and Spikes for Earth and Embedded Systems (SEA-CROGS) project. Fanhai Zeng is supported by the National Natural Science Foundation of China (12171283), the National Key R&D Program of China (2021YFA1000202, 2021YFA1000200), the Science Foundation Program for Distinguished Young Scholars of Shandong (Overseas) (2022HWYQ-045).
