Mathieu Bedez

Zakaria Belhachmi

Olivier Haeberlé

Renaud Greget

Saliha Moussaoui

Jean-Marie Bouteiller

Serge Bischoff

A fully parallel in time and space algorithm for simulating the electrical activity of a neural tissue

Keywords: GPU, neural activity, parareal algorithm, MPI, Partial differential equations 1 -Introduction

 Description of a method coupling parallelization in time and space on the GPU.  Resolution of a model describing the propagation of the electrical signal in a neural tissue using the parallel algorithm on GPUs.  Significant reduction in the calculation time to obtain an accurate solution.

A c c e p t e d M a n u s c r i p t 2

The main purpose of simulating EEG is to provide a better understanding of the connection between this global electrophysiological signal and the electrophysiological activity of subgroups of brain cells under physiological and/or pathological conditions. This electrical activity can be modeled using the bidomain model, now widely used in the simulation of the electrical signal during the heart activity, which constitutes the electrocardiogram. This model was introduced for the first time in the 1970s [START_REF] Schmitt | Biological Information Processing Using the Concept of Interpenetrating Domains[END_REF][START_REF] Tung | A bi-domain model for describing ischemic myocardial d-c potentials[END_REF]. Determining an analytical solution of this model, described via partial differential equations (PDEs), is not possible because the highly complex geometry. The digital resolution requests some simplifications in order to obtain a usable result within a reasonable computation time. These sacrifices are made, depending on the choice of the model, the precision of the mesh in time and space, or by fixing the problem within a portion of the definition of space. To avoid making concessions on the results accuracy, specific tools, which enable massively parallel computing, may be useful, distributing equitably the involved work, and thus reducing the computation time, while at the same time increasing the quality of results. In this study, we present a fully parallelization technique in time and space, and apply it to the monodomain model case [START_REF] Keener | Mathematical Physiology[END_REF], which constitutes a simplification of the bidomain model. This method will help us to determine the results quickly, and with an acceptable resolution in time and space. With this method, it is possible to run longer simulations within reasonable computation time, and also to perform simulations using very fine discretization grids, both in space and time. All this is not applicable in sequential methods, because calculation time would become excessively large We used GPU computation, an approach already used in neuroscience: in 2013, Hoang has developed a CPU/GPU simulation environment for large-scale biological networks: NeoCortical Simulator. The possibility of simulating spiking neural networks on one or several GPUs was investigated by developers [START_REF] Igarashi | Real-time simulation of a spiking neural network model of the basal ganglia circuitry using general purpose computing on graphics processing units[END_REF][START_REF] Wang | Simulation of large neuronal networks with biophysically accurate models on graphics processors[END_REF][START_REF] Brette | Simulating spiking neural networks on GPU[END_REF]. The use of GPUs yielded significant improvements compared to CPUs [START_REF] Baladron | Three Applications of GPU Computing in Neuroscience[END_REF]. Such implementations perform a parallelization at the neuron model level, whereas our study aims at performing a parallelization of the solving methods, at the global problem level. We performed several simulations to validate the different algorithms and the computational setups.

-Material and methods

-Bidomain model

As for the cardiac tissue, the neuronal tissue, defined by the domain B (see Fig. 1), can be modeled by decomposition into three distinct domains: the cells, called intracellular domain, the extracellular domain representing the outside of the cells, and the cellular membrane separating them. Each zone has an intracellular, an extracellular and a membrane potential ( , and respectively) and the voltage across the membrane is defined by their difference:

(1)

The intracellular, extracellular and membrane potentials are functions defined in by:

In the absence of any stimulus, all cells remain in resting state, ranging between -60 to -70 mV. When a cell is excited, a depolarization phase is initiated reaching a threshold value. At a given potential, sodium channels close spontaneously, terminating the depolarization phase. At the same time, potassium channels open, causing repolarization of the plasma membrane. Then the potassium channel closure period causes a transient hyperpolarization (or relative refractory period) of the membrane potential and reach its resting potential. Cell behavior has been described using models formulated from ordinary differential equations (see Section 2.3). Fig. 2 shows the variations of and to a position over time.

A c c e p t e d M a n u s c r i p t 3

The bidomain model might be used to describe the propagation of electrical signals along an axon, but also in the cells surrounding space. This model consists in two different types of equations: an evolution equation and an elliptic equation that describes the evolution of three potentials over time. Using a neuron model, and an appropriate conductivity tensor, these formulas can also be applied to neuronal tissue.

Throughout this research, the following assumption was made: in the intracellular and extracellular regions, current flow is resistive, resulting in a proportional relationship between the current density and the potential gradient (Fick s law), where the matrix is the conductivity tensor:

.

(2)

According to the current conservation law, the outflow from a domain must be equal to the inflow into the surrounding areas, which is described by the following equation:

.

(3)

The membrane current is defined with respect to the flow of the potential, by combining equations ( 2) and (3):

. ( 4
)
where and are conductivity tensors representing the intracellular and extracellular area, respectively, and is the unit outward normal vector of the intracellular domain. The behavior of the membrane is both resistive and capacitive. The perceived capacitive behavior of the membrane is due to an insulating bilipidic layer that separates the extra-cellular environment from the intra-cellular medium. The resistive part results from the transport of ions between the two environments by specific membrane proteins like voltage dependent channel. We consequently model this membrane with a resistor and a capacitor, connected in parallel, as schematically represented in Fig. 3. Associated currents are denoted and , respectively, and expressed in unit area per volume, while the membrane current is expressed in units of volume (Hodgkin and Huxley 1952). To maintain homogeneous units, one must multiply the ionic current and the membrane potential by area to volume ratio, which is denoted .

The macroscopic bidomain model is obtained from a homogenization process [START_REF] Neu | Homogenization of syncytial tissues[END_REF], which yields a coupled parabolic-elliptic system of equations. The first equation is:

in .

(5)

The second equation being:

in .

(6)

The ionic current is determined by the resolution of one or more ordinary differential equations (ODEs), e.g. as Hodgkin-Huxley (Hodgkin and Huxley 1952;[START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF] or Fitzhugh-Nagumo (FitzHugh 1955;[START_REF] Nagumo | An Active Pulse Transmission Line Simulating Nerve Axon[END_REF]) models, described in Section 2.3. We complete the model with Neumann boundary conditions for the extracellular potential and membrane potential allowing for the field to be perfectly isolated: 1, and the initial condition is given by: .

in (7) , (8)
(9)

-Monodomain model

The greatest difficulty in solving the bidomain problem is that the differential operator:

, is not explicitly defined from the expression of the membrane potential (equation 5), but only through the elliptic problem (equation 6). We therefore have two problems to solve:

(i) the dynamics problem according to (equation 5), (ii) the elliptic problem to update the extracellular potential (equation 6).

The monodomain approach was introduced to simplify this two-fold problem.

-Ionic current

To complete this model, we need to characterize the variable, which represents the ionic current. To this end, we propose to use well-established models that describe neuronal activity, and consequently define the dynamics of ionic currents. The selected model and its associated parameters are selected depending on the location within the brain. For our study, we use the Fitzhugh-Nagumo model, being a reduction of the Hodgkin-Huxley model (FitzHugh 1955;[START_REF] Nagumo | An Active Pulse Transmission Line Simulating Nerve Axon[END_REF].

The total ionic current described by the Hodgkin-Huxley model (Hodgkin and Huxley 1952;[START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF]) is given by the following expressions:

(13) , (14)
with: ,

with being the sodium current, being the potassium current, and representing the leakage current Variables and h (equations 11 and 12) are bounded between 0 and 1. Variables m and h describe the activation and deactivation of the sodium current, while n is the potassium A c c e p t e d M a n u s c r i p t current activation. The and functions correspond to the transitions between open and close states of the different gates involved in each ionic channel. Although this model precisely describes neuronal activity, it is complex and numerically costly to incorporate, as being composed of three ordinary differential equations (ODEs), as a whole to resolve. Introducing the Fitzhugh-Nagumo model, variables h and m vary over much slower time scales with respect to n and the membrane potential. Thus, one can combine the two variables in a single one, w, characterizing the degree of "refractoriness" of the system. The previous four-dimensional system (Hodgkin and Huxley 1952) is thus transformed into a two-dimensional one, having always the same dynamics characteristics.

It is composed of two variables:

the membrane potential and a recovery variable w:

(16) , ( 17
)
with the initial and boundary conditions (in the case of two spatial variables)

, where a, b and are the model parameters, and I is the external current being applied. Contrary to the Hodgkin-Huxley model, ionic currents do not explicitly appear, but the knowledge of the total ionic current is sufficient for the monodomain model.

Note that the values of the variables in the Fitzhugh-Nagumo model are bounded between -2 and 2, so it is necessary to rescale the obtained results in order to obtain physiological values. Thanks to the low computational complexity of this approach, it may be used in the modeling of very large neural networks.

-Parameters values

The parameter values, for the Fitzhugh-Nagumo model used here, are [START_REF] Fitzhugh | Impulses and Physiological States in Theoretical Models of Nerve Membrane[END_REF] :

a = 0.7 b = 0.8 = 13
We also define the parameter values for the monodomain model, assuming here that the conductivities , are the same in all directions (isotropic medium hypothesis):

-Spatial discretization

The spatial discretization of the PDE characterizing the monodomain is achieved using the finite difference method. We define the problem within where D is the dimension of space (1 or 2 in this study), and discretize space by introducing , the number of nodes along dimension and of step size . The differential operator then becomes: ,

A c c e p t e d M a n u s c r i p t 6 where is the approximation of and are the constant coefficients of the conductivity matrix . Other discretization formulas exist, but that one has the advantage of convenience.

-Time discretization

To solve the ordinary derivative equations, we introduce a time step , where N is the number of steps and T is the duration of the simulation. Time is approximated by discretized time is also associated with the discrete variables representing the unknowns of the equation system:

(membrane potential and recovery variable, respectively) for the Fitzhugh-Nagumo model.

Two schemes may be used to integrate these equations: Euler or 4 th oder Runge-Kutta. The time step to be chosen depends on the model, and also the model parameters values. For the complete resolution of monodomain model, one must sequentially solve each part, i.e:

1-Calculate the solution of the equation with the ionic current at time calculated above, using Euler and Runge-Kutta of order 4 integration: in .

(19)

2 -Determine the recovery variables to update the ionic current :

in . ( 20)

Each ionic current needs to be calculated at each node of the discretized space, resulting in a workload that varies depending on the location of the node, if one uses variable time step size.

In order to keep a high level of accuracy for the solution, while still maintaining acceptable computation duration, computations may be parallelized.

-Parallelization in time

To parallelize an ordinary differential equation or a partial differential equation, one can partition the time dimension. This technique appears at first look the least natural, since the resolution over time has the characteristics of being sequential. However, this technique has proven its interest in the case of very long temporal simulations. Lions et al. [START_REF] Lions | Real-time simulation of a spiking neuraComptes Rendus de l'Acadet al[END_REF] first proposed an algorithm to solve an evolution partial differential equation problem using evolutionary time multigrid. The convergence of this method has been proved in [START_REF] Bal | A "Parareal" Time Discretization for Non-Linear PDE's with Application to the Pricing of an American Put[END_REF][START_REF] Gander | On the Superlinear and Linear Convergence of the Parareal Algorithm[END_REF].

Consider the PDE problem:

, where A is the discretization of an operator. The time intervals are decomposed into N equal subintervals where . We introduce the operator that determines a coarse approximation of the solution of the problem (equation 16) with the initial condition . This propagator uses a step size . Let be an operator that determines a more accurate approximation of solution , using a time step The coarse propagator is used sequentially, then we solve in parallel, with the fine operator, the next equation over each subinterval:

, ( 22
)
with the initial condition: , where the functions (the coarse approximations) are defined for ; and the same initial boundary conditions . The standard algorithm begins with an initial coarse approximation with at time with the sequential calculation of , and then performs for each iteration the corrector step :

. ( 23
)
This parareal method stops when the next stop condition is satisfied: , (24) where is the tolerance.

Several implementations of the parareal algorithm were compared in [START_REF] Aubanel | Scheduling of tasks in the parareal algorithm[END_REF], including the distributed algorithm, presented as the most effective. Unlike the standard algorithm, the coarse propagation is shared across processes. At first iteration, process 0 performs the coarse resolution and sends the result to the next process, which in turn continues the coarse resolution. When the process has completed his M coarse resolution, it performs fine resolution using smaller time step than for the coarse resolution, followed by the corrector step from (equation 23).

The result is the sent to process M+1. If M+1 < N, then a convergence test is performed at each iteration. If process M converges, its work is finished. Otherwise, it starts a new iteration, taking into account the correction performed by the M-1 process. An example with 4 processes and 4 iterations processes is given in Fig. 4. The development of this algorithm is achieved through OpenMPI 1.4 (Message Passing Interface) that allows for exploiting the multi-processors characteristic of current CPUs. Another implementation of this algorithm was developed in CUDA, to handle the massively parallel architecture of GPUs. Unlike the MPI processes, which treat each a completely different piece of code, and therefore have to communicate their respective results, which can result in a blockage, GPU threads run the same piece of code, called the kernel, with different initial data.

-Spatial parallelization

One can also take benefit of the multi-threading capabilities of GPUs to parallelize the evolution equation on the monodomain model. The temporal resolution is carried out sequentially, whereas the spatial resolution is performed in parallel. To do so, we use the discretization of the gradient through the finite difference method. In its sequential resolution, we use the matrix formulation described in section 2.5, but with GPUs, we use the standard equation using the Taylor expansion to order 2: We define a CUDA grid having the same size as the mesh problem space, so that each CUDA thread can handle a single node. For example, in 1-D case, we restrict ourselves to a onedimensional grid containing 1-D blocks thread. The dimensions of grids and blocks are defined according to the size of the problem.

2.9 -Dynamic grids Version 5.0 of CUDA (Nvidia Corp. 2010), coupled with a GPU with a 3.5 or above microarchitecture version, has introduced dynamic parallelism. It is possible to launch kernels from threads running in the device i.e. threads can launch more threads. The interest of this technique is to launch an application in a coarse kernel, which in turn launches a finer kernel (see Fig. 5). Using this technique, the parallelism is coupled in time and space, thus enabling the use of a multithreaded version of the monodomain problem. The coarse kernel allows for parallel simulation time. Each thread of this coarse kernel is engaged to perform simulation in a time subinterval, and then all threads start the discrete calculation of the gradient (equation 20), and update the ionic currents (equation 13-14), on a fine kernel to parallelize computation space. By default, the grids launched within a thread block are executed sequentially: the next grid is launched only after the previous one has finished its computation. To implement a fully parallel algorithm, we must launch grid children competition. This is performed using CUDA streams. Kernels launched in different streams are executed concurrently. With the Kepler architecture, we can run up to 32 streams in competition (against 8 with the older Fermi architecture). In other words, the time interval can be split into 32 subintervals with the current GPU architecture.

-Results

The simulations are performed using the Python language and mpi4py, a full-featured Python bindings package for OpenMPI 1.4 that we deployed on the Rhenovia Pharma cluster composed of DL380 HP Servers, equipped with 256 processing cores clocked at 2.5 GHz and with 128 Go of RAM.

The GPU code was developed with the parallel computing platform and programming model provided by NVIDIA: CUDA 6.0. The compilation of GPU code was performed on a NVIDIA Kepler Tesla K20c board, with compute capability 3.5, running at a frequency of 706 MHz and equipped with 5120 MB of GDDR5 memory.

We define the acceleration of the algorithm as the ratio between the sequential resolution time and the parallel resolution one. This parameter is a simple measure of the gain between two methods of resolution. In this study, we compare the time resolution obtained with parallel policies from the sequential method. If the speed-up parameter is greater than 1, then the parallel method will be faster than the sequential method.

-Neurons Models

We initially validated the parareal algorithm on the system of ODEs modeling the activity of a neuron, following the Fitzhugh and Nagumo approach. We used the 4 th order Runge-Kutta scheme with a time step for the coarse solver, in order to ensure stability. For the fine solver, we choosed to use the explicit Euler scheme with a time step .

We determined the quadratic error between parareal algorithm results and those obtained with the sequential solution determined with explicit Euler scheme and a time step .

A c c e p t e d M a n u s c r i p t

The tests are based on the activity of a neuron for a period of 200 ms, to which is applied an external current over a period of 100 ms. The simulation time needed by the CPU to sequentially solve the Fitzhugh-Nagumo model is 92 seconds on one node of the Rhenovia Pharma cluster.

Then, we performed tests on 4, 8, 16, 20, and 32 different MPI processes.

Fig. 6(A-C) present the mean squared errors and corresponding time speed-up obtained by comparing the resolution time with sequential and parareal approaches, with MPI.

We note that for the Fitzhugh-Nagumo model, convergence takes place from the second or third iteration, according to the desired accuracy, i.e. when the mean square error is of the order of . In addition, a decrease in CPU time is observed when one increases the number of processors, without having to increase the number of iterations. Note that the error increases slightly when the number of processors increases, which can lead, in the case of a large number of processors, to a greater number of iterations needed to maintain convergence. These errors increase because there is a propagation effect. When the number of processors increases, the correction computed by one processor takes longer to spread towards the other processors.

We then tested the GPU version of this algorithm, applied to the same problem, and using the same simulation parameters. Fig. 6(B) presents the results obtained using the parareal algorithm on GPU. A grid is defined in one dimension containing a single block of threads. Using the same block, we can use the shared memory between threads and the synchronization between the threads within same block using the __syncthread() command. Shared and local memories are preferred, because the communication speed is significantly reduced when using global memory (up to 177 GB/s for global memory, while 1TB/s for shared memory is achievable). It also makes some blocking operations, so that the memory space dedicated to a specific thread is read too early, resulting in the recovery of a value, which has not yet been updated.

Our results indicate that the error between the reference (sequential resolution) solution and the GPU solution becomes negligible after no more than 3 iterations. The obtained error is notably larger than for the MPI parallelization (CPU) case with a small number of processes, but anyway remains below the desired accuracy. One important difference lies in the simulation time needed for additional iterations: with a GPU implementation, this time is very small, making computationally affordable the additional iterations needed to decrease the quadratic error. So, the GPU version of the parareal algorithm can drastically reduce the resolution time, while preserving a very high accuracy of the solution. Although some accuracy is to be sacrificed in order to quickly solve a problem, our results suggest that increasing the number of iterations reduces error, while still significantly decreasing the global simulation time. Fig. 6(C) shows the average speedup between the resolution time with MPI or CUDA and the resolution time in sequential. These speedups are determined by computing the ratio between the sequential resolution time and the parallel resolution time. We observe a significant time saving with the parallel resolution, compared to the sequential resolution. But the parareal CUDA implementation drastically decreases the computation time, by a factor 10 as compared to the MPI implementation.

-Monodomain 1-D model

In this section, we consider the parallelization of the monodomain model, applied on the interval representing the spatial domain, and for a simulated period of 100 ms.

The simulation consists in applying an external current during 90 ms, using a resolution , a coarse time step size and a fine time step size . The reference solution is calculated using the explicit Euler scheme with a time step size . This solution is obtained in 820 seconds using one node of the Rhenovia Pharma cluster.

The results are presented in Fig. 7(A) for different numbers of MPI processes and increasing iterations.

A c c e p t e d M a n u s c r i p t 10

We see that for an error in the order of , we need about 4 iterations to reach convergence. The mean square error appears to increase slightly when the number of processors increases, suggesting that a higher number of iterations may be needed for a large number of processes, depending on the accuracy and degree of smoothness the user wants.

The CUDA implementation of the monodomain resolution uses the dynamic grid (see section 3.5) to parallelize also the Laplacian calculation. To get the most from competition, we use the CUDA streams. Each stream will perform problem solving on a sub-interval, launching the Laplacian calculation on a children grid. The resolution of the problem is therefore completely parallelized. Fig. 7(B) presents the mean square errors for different sizes of sub-intervals (depending on the number of streams and hence the size of the parent grid). In all cases, the error is of the order of after the 3 rd iteration. Compared to the MPI resolution, we observe that the error after the first iteration is lower, being of the order of , while it is around for MPI. A more rapid convergence is observed for the CUDA implementation.

The speedup that depends on the number of MPI processes or CUDA streams, is shown in Fig. 7(C). The parareal simulation was completed when the error between two iterations respects a user-defined tolerance. Again, we see that there is indeed a reduction of computing time by a factor comprised between 1.6 and 9.5 when using the MPI implementation. When the CUDA implementation is used, we observe an even more dramatic reduction, by a factor of 80 to 138. The time parallelization combined with space parallelization therefore helps for solving this model almost immediately. As CUDA grids childrens depend on the size of the problem, the computation time of the monodomain model with a finer spatial resolution will not take much longer than this resolution.

-Monodomain 2-D model

We then apply the methods discussed earlier to parallelize calculations of 2-D simulation of the mono domain model. The spatial domain is defined by the square with a spatial resolution , the coarse and fine time step size being and , respectively.

The simulation represents a period of 100 ms, with application of an external current during 90 ms. The obtained results are summarized in Fig. 8, for varying numbers of MPI processes and CUDA streams.

Results from the MPI implementation are shown in Fig. 8(A ). It is observed that 6 iterations are needed to obtain an error of the order of between sequential and MPI parareal solution. Note that in this case, the speedup generally decreases by a factor of 2 compared to the 1-D case. It is therefore necessary to launch this integration on a large number of nodes in order to reduce the time of resolution, i.e. the time needed to reach a solution close to the reference solution.

Contrariwise, the computational time required to solve the monodomain model using dynamic grids is almost identical that the time required in 1-D, because the children grids has adapted to the size of problem. Thus, there is the same number of procedures to be performed by CUDA thread relative to the 1-D problem. Sequential solving time being logically longer, the speed-up is more important as shown Fig. 8(C). We observe that convergence seems faster than with MPI (see Fig. 8(B)), because the mean square error for the first iteration is small, while the MPI error is so significant that it has not been shown in Fig. 8(A). Indeed, it takes an average of four iterations to get to an error in the order of .

-Monodomain 3-D model

Finally, we tested our method by applying the monodomain model to a 3-D geometry. The space sector is defined by the cube with a spatial resolution , the steps of coarse and fine time are those defined in Section 3.3, i.e. and .

The simulation is performed over a period of 100 ms, by applying to a portion of the tissue an external current during 90 ms. Sparse matrices are used to reduce the memory space allocated in the sequential and MPI resolutions. The results obtained using the different methods are summarized in Fig. 9.

Fig. 9(A) shows the results obtained using the parareal method developed using the MPI technology. Again, there is a decrease of the error during the various iterations. Note also that the error is homogeneous with respect to the choice of the number of processors. Furthermore, it becomes acceptable between the third and fourth iteration, that is to say around . In the case of MPI, Fig. 9(C), we see that the speedup is of the same order of magnitude than those observed for 2-D geometry. The gain becomes non-negligible when the number of processors increases, which tends to confirm, once again, the usefulness of the parareal approach to solve the problem. Fig. 9(B) shows the results obtained when applying the method developed on GPU. The errors seem to again be homogeneous according to the number of streams used in the parallelization of the time scale. Also, one observes a continuous decrease of the error during the various iterations, even if they are obviously somewhat larger compared to the MPI application. An error of the order of is obtained after three or four iterations, which confirms the results obtained in the 2-D case. Fig. 9(C) shows the speed-up obtained using the fully parallel method. It is noted that the resolution time decreases greatly compared to sequential time and MPI, even if the gain is not as important as in 2-D. This is because the sequential resolution and MPI are performed using sparse matrices, thereby significantly reducing the calculation times, while the resolution in space GPU does not use the matrix formulation of the problem.

-Discussion

The size of the problem that can be set is limited by technical characteristics of the different technologies. In MPI, there is a limit to the size of the data, that commands Send and Broadcast can transfer. This limit depends on hardware characteristics of the machine (memory, processor), but also on the MPI version, which is used. For this reason, and especially in 3-D, the use of sparse matrices is required, thereby greatly reducing the memory to be used. Some data communications are however not compressible, limiting the maximum size of the problem. The same difficulty holds when using GPU architecture. The limits are obviously different depending on the graphics card and its architecture. The sizes of the different types of memory (local, shared, global) differ from an architecture to another. As we have seen in section 3.5, the number of streams that can be launched in competition is limited to 32 on the current Kepler architecture. This is reflected in our method by dividing by 32 the maximum time scale. The next architectures will push these limits further back, as was the case during the passage of the Fermi architecture to the Kepler one (with 8 versus 32 streams, respectively). Regarding the parallelization of space, the limit, if one attaches one thread per mesh node, depends on the maximum size of the GPU grid and the maximum size of a thread block. The graphics card used during our test phase (Tesla K20) defines a grid of maximum dimension (2147483647, 65535, 65535) and a maximum size of thread block (1024.1024, 64), already allowing a fine enough geometry for solving our present problem. All these specifications are to be considered on a power development, which will aim to solve a PDE problem in future, especially in 3-D, as the resolution of bidomain problem.

-Conclusion

We have applied parallelization methods to the monodomain model, which is a simplification of the more complete bidomain model. This parallelization allows for significant reduction of the simulation time by using fine resolutions in both time and space, using GPU technology. This A c c e p t e d M a n u s c r i p t method thus allows obtaining a better solution whether in space or time, but it can also perform much longer simulations on larger areas without catastrophic increase of the resolution time. With this technology, and using a fraction only of the resources of the graphics card, it was shown that the resolution of an ordinary system of differential equations is performed in less than one second of time. Then by applying the same process to 1-D, 2-D and 3-D monodomain models, we have shown that the resolution time also decreased by increasing the size of the resource used either in time or in space. Using these methods, we can simulate this model over much longer periods, but also in larger spaces. Furthermore, they allow for keeping a high level of resolution limiting the calculation time. GPU technology has also permitted to develop a resolution of the problem by parallelizing in space and also in time with dynamic grids. The ratio between the computation time in sequential and GPU implementation is very important compared to MPI (launched on multiple cluster nodes). This study showed that the use of a graphics card enables massive parallelization of a problem, thus requiring the full capabilities of the Tesla card. But in comparison, running MPI allows for a limited parallelization only. These methods have proven to be helpful for simulation of a monodomain model, and the logical extension will be to apply this method to a bidomain model for a better description of the propagation of electrical activity, but requiring considerable resources and especially time calculation.

-Acknowledgements

Mathieu Bedez thanks the entire team of Rhenovia Pharma for the biological expertise, which enabled us to do this study, as well as for the access to the Cluster computing, and Bertrand Simon, from MIPS laboratory, for giving access to its Tesla GPU K20c, to develop and demonstrate the concept of our method.

-Conflict of interest

Jean-Marie C. Bouteiller has a conflict of interest. The University of Southern California holds an equity interest in Rhenovia Pharma, and has also received licensing income from Rhenovia Pharma.

  unit normal vector of domain B. The different variables are shown in Fig.

  nothing else than equation 18 applied to monodomain model, and with the conductivity parameters given in Section 2.4.

FiguresFig. 1 -Fig. 3 -Fig. 4 -Fig. 5 -

 1345 Figures

Fig. 6

 6 Fig. 6 -(A) Mean square error for the Fitzhugh-Nagumo neuron model following the stop iteration of parareal algorithm, and for different numbers of MPI process. (B) Mean square error of the simulation of Fitzhugh-Nagumo model with parareal algorithm developed in CUDA. The error is determined for different stop iterations and size of CUDA grids. (C) Average speed-up of the CUDA resolution, compared to the parareal MPI resolution.

Fig. 7

 7 Fig. 7 -(A) Mean square error for the 1-D monodomain model between the sequential and the MPI parareal solutions at each iteration, up to convergence. (B) Mean square error for the 1-D monodomain model between the sequential and GPU solutions, at each iteration, up to convergence. (C) Speed-up between the sequential time calculation and the MPI/CUDA time calculation.

Fig. 8

 8 Fig. 8 -(A) Mean square error for the 2-D monodomain model between the sequential and the MPI parareal solutions at each iteration up to convergence. (B) Mean square error for the 2-D monodomain model between the sequential and GPU solutions, at each iteration, up to convergence. (C) Speed-up between the sequential time calculation and the MPI/CUDA time calculation.

Fig. 9

 9 Fig. 9 -(A) Mean square error for the 3-D monodomain model between the sequential and MPI solutions at each iteration of the parareal method. (B) Mean square error for the 3-D monodomain model between the sequential and GPU solutions, at each iteration of the parareal algorithm. (C) Speed-up between the sequential time calculation and MPI/CUDA time calculation.