Space fractional derivatives have been used instead of integer operators with success in a variety of practical applications to study and describe transport processes in media characterized by spatial connectivity properties and high structural heterogeneity altering the classical laws of di usion. This study provides an implementation of space fractional FitzHugh-Nagumo model of excitable medium in two-dimensional space. This implementation is based on shifted Grunwald-Letnikov computational scheme which transforms the space fractional di usion equations to a system of algebraic equations. After that, the conjugated gradient method is used to solve the sparse matrix system. MPI and PETSc library is used for implementation. We investigate the in uence of some factors to our implementation performance and scalability such as core numbers, size of computation mesh, and fractional derivative orders. All tests are done on ccNUMA architecture with two CPUs.
Living system simulations make great impact on the understanding of biological processes. They help in checking the side e ects of drugs and allow to create a personalized model of human organs which help in planning of surgery and medication. Recently, some researchers have been looking for a way to apply fractional derivatives [1] to living system models. The main advantage of these approaches is to describe sub-di usion processes. The sub-di usion is a di usion process which gives a squared displacement of particles which are non-linearly dependent on time. Especially, this variant of di usion exists in tissue and subcellular space. The model of visco-elastic property of soft tissue which can explain oscillation damping more precisely is one of models which contains fractional derivatives. The review of these sorts of models are described in the introduction of Ph.D. thesis [2]. In addition, there are several other applications of models with fractional derivatives to life science: cancer growth [3], initial stage of HIV [4], and population growth [5].
Individual attention is required for the excitable medium models. These models describe neurons and muscle tissue (myocardium). The special problem is a more realistic description of di usion in consideration of subdi usion transport in subcellural space and tissue structure. In [6], the authors proposed a replacment of integer Laplacian in excitable medium model by fractional Laplacian. This replacement gives a precise reproduction of several physiological experiments in computer simulations. An important feature of fractional derivative in a new model behavior, which can not be repeated with changing di usion coe cient or introduction of nonlinear parts.
Unfortunately, the replacement of integer derivatives by fractional ones leads to a signi cant increase in the computational time. This fact motivates us to introduce a high performance implementation for fractional FitzHugh-Nagumo model [7] based on some computational methods [8]. A fractional FitzHughNagumo model is a generalizing of the standard model that describes the propagation of the electrical potential in heterogeneous cardiac tissue. The fractional model consists of a couple fractional Riesz space nonlinear reaction-di usion model and a system of ordinary di erential equations. This implementation is required for physiological experiments in the future.
We use PETSc library [9{11], which provides iterative solvers for sparse matrix. This library is used in more than 14 software packages for scienti c computing [9] and support multicore architectures, computation accelerators and systems with distributed memory. The implementation based on this popular and well-known library is more suitable for our goals.
Besides that, we investigate the in uence of some factors in our implementation performance and scalability such as core numbers, size of computation mesh, and fractional derivative orders. All tests were conducted on Cache coherent non-uniform memory access (ccNUMA) architecture with two six core central processing unit (CPU). The rest of this paper is as follows: the next section is devoted to describe FitzHugh-Nagumo model and fractional Laplalacian operator, then we proceed to describe conducted experiments. 2
The FitzHugh-Nagumo model is a reaction-di usion equation which describes waves propagation in an excitable medium. Initially, it was created as a mathematical description of the action potential propagation on an axon membrane. @u
= r (Kru) + u(1 u)(u a) v; (
= "( u v ); (
During a long period of time, the FitzHugh-Nagumo model was being used for research of auto-oscillatory process in the excitable media. For example, the spiral waves in a two dimensional domain and scroll waves in a three dimensional domain. These structures are models of arrhythmic activity in mammalian hearts.
The model is suitable for this study, because it is simple, contains only two phase variables and converges on big time steps. In addition, spiral waves provide test example for a performance analysis which is close to realistic problems. 2.1
Usually, researchers and engineers use derivatives and integrals with integer order of degree in applications. However, there are several de nitions of derivatives and integrals with real and complex numbers of degree. They are named fractional derivatives and fractional integrals. Some of these de nitions are given by Grunwald-Letnikov, Riemann-Liouville, Caputo, Marchaud, Hadamard, and more [12]. In some cases, many of these de nitions are equivalent to each other.
In this paper we use a fractional derivative with order 1 (1 < 1 2) which is de ned by Riesz [12].
2 cos( ;
= where v( ; y; t) ) 1 1 d ; d : x
Similar expressions are de ned for the space Riesz fractional derivative of order 2 (1 < 2 2) with respect to y.
After that, we can introduce a fractional Laplacian operator.
r = ; @jxj @jyj @jzj In [1, 13, 12], the authors describe the theory of fractional derivatives in detail.
After the replacement of integer Laplacian operator with the fractional one,
the fractional FitzHugh-Nagumo model appears:
@u @ 1 u @ 2 u
@t = Kx @jxj 1 + Ky @jyj 2 + u(1 u)(u a) v; (
= "( u v ); (
In our research, we use the numerical method which is described in [14, 8]. This
numerical method transforms (
The numerical simulation is formulated by considering positive integers m1, m2 and N such that hx = Dx=m1; hy = Dy=m2; = T =N be characteristics of space and time grids. The orthogonal uniform computation mesh includes m m nodes with step by space h. is a time step and n is a number of conducted time steps. uin;j is the value of unknown function u(xi; yj ; tn).
The GrunwaldLetnikov approximations for l.h.s and r.h.s fractional derivatives are given by: = =
1 (hx)a
1 (hx)a i+1 X g(l)u(xi l+1; yj ; tn) + O(hx); l=0 m i+1
X
g(l)u(xi+l 1; yj ; tn) + O(hx);
l=0
where coe cients gl are de ned by recurrent formula:
g(0) = 1; g(l) = (
Similar expressions are de ned for y axis direction. The Order of fractional derivative for y direction is denoted as 2.
Substitution of expressions (
i+1
uin;j + r(
l=0 vin;j = vi;j n 1
n 1 + si;j ;
r(
Kx
2 cos(
X l=0 m j+1
X l=0
Ky 2 cos( 2 2 )h 2
; g(l2)uin;j+l 1 = ui;j n 1 n 1 + fi;j ; where fin;j = uin;j (1 uin;j )(uin;j a) vin;j is nonlinear part in the rst model equation, and sin;j = "( uin;j vin;j ) is right hand part of second equations.
Boundary conditions in Dirichlet form are described in the following form: n n n n u0;j = um;j = ui;0 = ui;m = 0:
The linear system Au(n) = b(n 1) is constructed based on these expressions.
A is matrix of the system with size (m 1)(m 1) (m 1)(m 1) in left hand
(
(b) 1 = 2 = 2:0
side of equations (
Figure 1 presents matrix A with di erent values of 1 = 2. The problem became identical to usual FitzHugh-Nagumo model if 1 = 2 = 2:0.
In work [14] authors use Gauss-Seidel method to solve the equations. In our
research we use conjugated gradient method without precondition.
At the rst stage of the algorithm, arrays and the matrix are initialized by PETSc
library. PETSc library does a partition of arrays and the matrix across di erent
processes. The matrix and arrays are constructed in a parallel way. Each process
keeps in memory a single part of every structure. The formulas of matrix values
2 (curvature of the waves is due
are described in numerical scheme (
We consider three series of experiments to study performance of the simulation. In the rst series, we establish the e ect of order of fractional derivative to the performance. We take di erent values for 1; 2 and x the number of cores to one. The second series is concerned with the establishment of scalability for di erent numbers of threads (from 1 to 12). The fractional orders 1 = 2 are xed. The third series discusses the relation between number of mesh nodes and scalability. The order of fractional derivatives 1 = 2 chosen to be 1.9 and 2.0. All experiments were conducted on one node of Ural Federal University's cluster with two CPU Intel Xeon E5-2620 v2 @ 2.10GHz, which are joined in ccNUMA architecture.
In [6], the authors compare features of excitation wave behavior when the order of fractional derivative is equal to 1.9, 1.75, and 1.5. In this research, we Algorithm 1 The algorithm of initialization and main computational cycle 1: procedure Compute(t; s; f; T ) 2: Initiate(Ai;j ; uin;j ; vin;j) 3: ConstructTheMatrix(A) 4: SetInitialValues(u0; v0) 5: while t n < T do 6: b un 1 + t f (un 1; vn 1) 7: Solve(Aun = b) 8: vn vn 1 + t b 9: Print(t; un; vn) 10: n n + 1 11: end while 12: return un; vn 13: end procedure use the same parameters. Other variants and combinations are used to complete the study. If 1 = 2 = 2:0, fractional FitsHugh-Nagumo model is reduced to standard FitsHugh-Nagumo model with integer order of Laplacian. If 1 = 1:9 or 2 = 1:9, the simulation behavior is similar to the case of ordinary Laplacian, but performance drops. This happens because the matrix of the linear system is more complex for fractional derivative, see Fig. 1. This leads to a high number of steps in iterative CG solver.
The rst series of experiments are done on a square area with size 2:5 2:5
and 256 256 mesh nodes, step-by-time = 0:1. The parameters 1; 2 are
chosen in range from 1.1 to 2.0 with step 0.1. The simulation takes 1400 time
steps for each parameter combination. A spiral wave initiation protocol (
When 1 = 2; 1 2 [1:1; 2:0], the results are shown in Figure 3(a). For all possible combinations of parameters 1 2 [1:1; 2:0]; 2 2 [1:1; 2:0], the results are shown in Figure 3(b). The computational time grows with increasing 1; 2 and dramatically decreases when 1 = 2:0 or 2 = 2:0. The computational time for 1 = 2 = 2:0 is 35 times lower than in case of 1 = 2 = 1:9.
We can explain this result by reduction of the fractional problem to the problem with integer Laplacian when 1 = 2 = 2:0. This reduction leads to a simpli cation of the matrix (Fig. 1). The multi-diagonal matrix of fractional problem is reduced to banded matrix with ve diagonals. This simple matrix structure requires less iterations of CG solvers. The di erence in computational time for di erent derivative orders can be explained by features of the matrix or features of spiral wave behavior.
In the second series, we xed values of 1 = 2; 1; 2 2 [1:1; 2:0] and studied the scalability for number of cores from 1 to 12. Fractional derivative orders are set to be 1.1, 1.25, 1.5, 1.75, 1.9, 2.0. Simulation takes 1000 time steps on a 5:0 5:0 mesh with 512 512 mesh nodes.
The time of computations is shown on a Diagrams 4(a), 4(b) and scalability is shown on a Plot 4(c). The scalability is close to linear for 1 = 2 = 2:0, but in case 1 = 2 < 2:0 it increases up to a sixth core and stays on plate. Similar result is shown on two and three nodes of cluster. In our opinion, these results may be explained by counts of slow non-local memory readings from di erent cores in ccNUMA architecture.
In the third series, we study a dependency between scalability and mesh size. The computation is done on a di erent square meshes with sizes 1:25 1:25, 2:5 2:5, 5:0 5:0 and 128 128 , 256 256 , 512 512. This formulation of the problem is important, because the computational scheme 12 establishes dependency of each value to the values from adjoined row and column way to the mesh border. Expansion of the computational mesh leads to increasing of the number of nonzero elements in the matrix. Parameters of experiment are taken as: h = 0:009767 is step by space, = 0:1 is step by time, 1000 time steps were computed. The order of fractional derivatives 1 = 2 are set to be 1.9 and 2.0.
In this section we compare our implementation with some existed works in literature. In [15], the authors proposed some numerical methods for e cient computation of a class of matrix functions on GPU. This method allows to improve performance of exponential integrator and computation of fractional Laplacian. The authors describe computationally fractional-by-space reactiondi usion equations and fractional-by-space FitzHugh-Nagumo model on CPU and GPU. Libraries cuBLAS and cuSPARCE are used as parallelization technology. MATLAB is used for implementation.
In their study, 20000 time steps are computed during 175 369 seconds with CPU and 14 579 seconds with GPU. Approximately, it is about 8.768 seconds for one time step on CPU and 0.729 seconds for GPU. Our implementation computes 1000 time steps in 772 sec. with 12 CPU cores. It is about 0.771 steps in one second. Thus, performance of our implementation is faster than CPU implementation from the paper [15] and close in performance to GPU implementation.
In [16], the authors described parallel algorithm for solution fractional-byspace reaction-di usion equations in one dimensional case. Finite di erence method is used for fractional partial derivative discretization. Performance study was achieved with one core, one CPU (Intel Xeon X5540) and 64 cores. In the last case, MPI technology is used. Programming language was Fortran 90. Remarkably, computational scheme is based on Grunwald approximations. This is similar to our study, but we solve the problem in two dimensions.
In their study, performance increased by 3.38-3.58 times with 4 cores and 4.55{5.23 times with 8 cores. We computed minimal and maximal speedup achieved in our simulation for all variant of parameters excluding case of Laplacian with integer order. In our implementation, performance increased by 2.29 { 3.85 times on 4 cores, 3.5 { 5.0 times on 8 cores and 4.06 { 5.28 times on 12 cores. Based on these metrics, we can say that our implementation scalability is close to scalability of the implementation [16].
Besides that, the authors in [17] described program to solve Maxwell's fractional equation which uses PETSc, but performance and scalability tests are not described. In [18] and [19], a study of performance for several numerical methods without parallelization was performed. In [2] and [20, 21], a study of e ective parallelization techniques was done. The authors used PETSc. However, these works are not compatible with our work directly, because the authors solved the problem with another conditions and other computational methods. In addition, our study is more detailed. We study performance from hardware and model parameters together in one research. 5
In this paper, we describe a high performance implementation of the fractional FitzHugh-Nagumo model. It is based on numerical method [8] to approximate fractional reaction di usion equation and conjugate gradient method to solve the sparse linear system. Our solution is done with the aid of PETSc library which is used a lot for scienti c computing [9]. The simulation with fractional derivatives shows 17%-65% less performance than ones with integer order derivatives. We study the dependencies of simulation performance such as size of mesh and order of fractional derivative. In addition, we establish simulation scalability in ccNUMA architecture with two six-core CPU. The proposed implementation shows a linear scalability when a fractional problem is reduced to a problem with an ordinary Laplacian ( 1 = 2 = 2:0). However, the implementation speedup reaches plateau from six cores and more with other values 1; 2.
The proposed implementation can be used for research of physiological features of excitable medium model with fractional derivative in the future. Acknowledgements. The study was supported by the Russian Science Foundation (no. 14-35-00005). We used the computational cluster of Ural Federal University for this research.