Dynamical behavior of reaction–diffusion neural networks and their synchronization arising in modeling epileptic seizure: A numerical simulation study

Abstract

Excessive synchronizations of neurons in the brain networks can be a reason for some episodic disorders such as epilepsy. In this paper, we simulate neural dynamic models and their synchronizations numerically by a new numerical algorithm. In order to overcome the complexity of the problem, a suitable combination of the Legendre spectral element method and operator splitting technique is proposed to convert a complex system of partial differential equations into sparse linear algebraic systems. The simulations are explored in different aspects such as accuracy, computing the wavefront and angular speeds, computing time, the condition number of obtained algebraic systems to demonstrate the proposed scheme effectiveness. The numerical simulations and comparisons confirm the method is suitable in terms of accuracy and speed.

Introduction

Epilepsy is a paroxysmal or episodic disorder which is studied at different scales, molecular to system levels in spatial scale or modeling of interictal spikes to disease evolution in temporal scale [1]. Studies at each scales have significant information for scientists investigating epilepsy. Dynamic models are the powerful approaches which give useful data about a complex system such as epilepsy. On the other hand, the (de)synchronization has a key role in some brain activities and study this phenomenon in the related dynamical systems can help to understand some cognitive processes [2]. Researches reveal that the abnormal synchronization of neurons in different parts of the brain might lead to epilepsy [3]. Thus, scientists explore the computational models of the dynamics of single neurons such as FitzHugh–Nagumo (FHN) model [4], [5], [6]; then, they study various kind of neural networks derived from the single neurons and explore the (de)synchronization in these networks [7], [8] to obtain some virtual data to complete the experimental data in the study of epilepsy.

In this paper, we focus on single and network of reaction–diffusion neural dynamics models. In other words, first, the dynamical modeling and behavior analysis of an FHN neuron is studied; then, synchronization of a network of FHN neurons which promising the understanding of cognitive processing has been presented by considering the unidirectional connectivity between two distant neurons. In order to explore the behavior of coupled FHN neurons, the influences of changing neuronal parameters for unlike neurons on the dynamical properties such as synchronization error have been studied.

The derived mathematical model is often complex and analytical solution is not available; therefore, a powerful computational algorithm is needed for effective simulations. Recently, a numerical method based on compactly supported radial basis functions (CSRBF) method has been proposed by Hemami et al. to simulate the model discussed in this work [6] and a numerical method based on the multi-domain pseudospectral method has been employed only for solving the dynamic model of a neuron by Olmos and Shizgal [5]. Moreover, Moghaderi and Dehghan used a mixed two-grid finite difference method to obtain the dynamical behavior of the one and two dimensional FHN model of a neuron [4]. There are also some other works which have investigated a simplified case of the FHN model, there is no recovery variable, by the various numerical approaches such as homotopy analysis [9], spectral element method [10], a semi-analytical method [11]. An FHN neural network is obtained by connecting a number of FHN reaction–diffusion systems to each other by a specific connection type. Since the FHN model is a complex dynamical system, its neural network is also a complex network. If the network components (FHN neuron) which are connected to each other act as the same, a synchrony is accrued in the system where control of the synchronization is necessary to prevent the malicious events in the network. So, first, we should detect the synchronization phenomena in the desired network; then, control that. Moreover, we face a complicated problem; thus, studying issues such as the synchronization in this system needs advanced numerical approaches. To the best of our knowledge, the combination of the operator splitting (OS) method with the spectral element method has not been used to simulate reaction–diffusion neural dynamics models; thus, we propose such an effective algorithm for numerical simulation of FHN model and studying the synchronization in the networks of FHN neurons.

Spectral methods (SM) are one of the useful numerical methods for applying in dynamical models because of their exponential convergence [12]. On the other hand, the finite element method is the other approach which is appropriate for simulation of the complex systems. But each of these methods has some drawbacks. In spectral methods, the coefficient matrices in the obtained system of algebraic equations are full and ill-conditioned. Moreover, there are more weaknesses such as limitations for applying the method in large or complex geometries, parallelization difficulties due to their global approximation [13]. On the other side, the low-order finite element is not helpful where an accurate simulation is important. Therefore, in 1984, Patera combined these methods for the first time to overcome the disadvantages of the methods and utilized the benefits of finite and spectral methods at the same time and called it the spectral element method (SEM) [14].

SEM is considered as a high-order finite element method which has a spectral convergence for the smooth functions ; besides, this method can decrease the complexity of computations in comparison with the classical spectral methods because of use of low-order polynomials in each element and sparse global operation matrices [15]. SEM is based on the Galerkin projection; in fact, first, the domain is divided into non-overlapping sub-domains; then, the Galerkin projection is implemented on each sub-domain [16]. Scientists have used SEM to simulate various models; for instance, Chen et al. introduced an SEM based on Laguerre functions to simulate Black–Scholes and Metron jump–diffusion models in an unbounded domain without any domain truncation [17]. In 2012, Yildrim and Karniadakis proposed a new hybrid spectral/DG approach for simulation of ocean waves [18]. In 2015, Fakhar-Izadi and Dehghan proposed SEM with a hierarchical basis to approximate second-order nonlinear two-dimensional partial differential equations [15]. They utilized fast Fourier transform to decrease the aliasing error. Dehghan and Sabouri investigated the predator–prey system by Legendre spectral element method and compared the obtained results with Legendre pseudospectral approach and finite element method to show the effectiveness of their proposed method [16]. In addition, Zayernouri and his collaborators developed SEM for solving the time–space fractional and fractional delay equations [19], [20], [21]. In 2014, Cantwell et al. used a high-order spectral/hp element method to approximate a kind reaction–diffusion problem and applied it to cardiac electrophysiology simulations [22]. Dehghan and Abbaszadeh examined fractional partial integral equations (fractional evolution equation) using a combination of difference scheme and SEM [23]. Sabouri and Dehghan used an implicit nodal SEM in complex geometries with triangular and deformed quadrilateral and hexahedral elements [24]. They implemented their proposed method on time-dependent nonlinear diffusion equations. In addition, in 2018, they improved a $hk$ mortar spectral element to solve the $p$-Laplacian problem and they study both of the $h$ and $k$ convergence for this method [25]. Kharazmi et al. proposed a new Petrov–Galerkin modal spectral element method to solve one-dimensional fractional elliptic problems [26]. In 2018, Fakhari-Izadi and Dehghan developed modal spectral element and finite difference methods for space and time discretization, respectively in curvilinear domains [27]. Recently, Abbaszadeh et al. developed a Legendre spectral element method to simulate the stochastic nonlinear system of advection–reaction–diffusion models effectively [28]. In addition, Abbaszadeh et al. employed spectral element method in conjunction with alternating direction implicit algorithm (ADI-SEM) for solving the multi-dimensional generalized modified anomalous sub-diffusion equation to reduce the computational complexity [29]. The interested reader of SEM can also see [30], [31], [32].

In the current paper, the Legendre spectral element method (LSEM) is proposed jointly with an appropriate operator splitting approach, which the main problem is divided to the three differential equations, two partial differential equations and a system of nonlinear ordinary differential equations. The proposed method is applied in the obtained dynamical systems. The Crank–Nicolson approach and LSEM are applied in PDEs for time and space discretization, respectively; then, the sub-problems reduce to the sparse and small dimension linear algebraic equations where can be solved by QR factorization technique efficiently. Finally, the fourth-order Runge–Kutta algorithm is implemented in the third obtained dynamical system (system of nonlinear ODEs) to discretize the time. Moreover, the LSEM combined with OS approach is employed to study the synchronization in FHN neural networks to find an approximate interval for the synchronized coupling strength parameter. Actually, an FHN network with $N$ neurons is modeled by a coupled system of $N$ FHN dynamical systems, and by using OS, this system is divided to $5N$ differential equations. Then, these subproblems are solved by a combination of the Crank–Nicolson, LSEM and RK4 approaches. Since there is no exact solution for the model, we compare our numerical results with the state-of-the-art methods. Moreover, as [6], a very accurate approximation with the proposed method are considered as the exact solutions to evaluate the approximation errors. The numerical results and computing times shed light on the efficiency of the proposed method. For instance, in one of the test cases the approximation error is about 10⁻⁵ and the computing time is 3 m. On the other hand, we explore the computational complexity of the proposed method and compare it with the pure spectral method and spectral method in conjunction with OS. It can be deduced that the spectral element method combined with OS algorithm is more appropriate than the pure spectral method. Although the SEM and spectral methods in conjunction with OS approach have the same computational complexity, the first method is more convenient due to its well-condition matrices. It is worth to mention that the main purpose of this manuscript is to show that the multi-domain spectral method especially spectral element method, which has never been applied to numerical simulation in mathematical neuroscience and cognitive neuroscience, can yield accurate and fast approximations of FHN single neuron and neural networks dynamical systems. Generally, our focus in this work is more devoted to developing an accurate, computationally efficient, stable, convergent technique to simulate the neuronal dynamical systems. So, we find a new numerical approach by combining Operator splitting (Trotter approach), Crank–Nicolson, Galerkin spectral element, fourth-order Runge–Kutta and QR factorization methods.

In sum, we propose an advanced and appropriate numerical method based on OS and LSEM schemes to simulate complex reaction–diffusion dynamical systems arisen in computational cognitive neuroscience. The model includes one and two time-dependent examples and reaction–diffusion dynamical networks. We examine the efficiency of the method analytically and numerically. Also, the stability and convergence of the proposed method are studied. It is shown that first, the OS scheme decreases the computational complexity dramatically, and second, the LSEM facilitates the calculations due to small number of the condition number of the obtained coefficient matrix compared to the spectral method. Moreover, we report the obtained numerical results with different test cases and compare them with the other methods which confirm the high accuracy and acceptable speed of the proposed algorithm. Finally, the proposed algorithm is utilized to simulate networks of the FHN neurons and examine the effects of the coupling strength parameter on the synchronization of the neurons.

The remainder of the paper is organized as follows: In Section 2, the models which are discussed in this work are described. Section 3, the details of the proposed method are explained and then the method is applied in the dynamical systems. The numerical simulations and results are discussed in Section 4. Finally, Section 5 concludes the paper.