Thierry Dumont
email: tdumont@math.univ-lyon1.fr

Max Duarte
email: max.duarte@unice.fr

Stéphane Descombes
email: stephane.descombes@unice.fr

Marie-Aimée Dronne
email: marie-aimee.dronne@recherche.univ-lyon1.fr

Marc Massot
email: marc.massot@ecp.fr

Violaine Louvet
email: louvet@math.univ-lyon1.fr

Simulation of human ischemic stroke in realistic 3D geometry

Keywords: Ischemic stroke, reaction-diffusion equations, operator splitting, parallel computing 2010 MSC: 35A35, 35K57, 65L06, 65M08, 65M50, 65Y05, 92B05

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Stroke is a major public health problem since it represents the second leading cause of death worldwide and the first cause of acquired disability in adults. In the United States, this disease strikes once every 40 seconds and causes death every 4 minutes, with an estimated 41.6% death rate in 2007 [START_REF] Lloyd-Jones | Heart disease and stroke statistics-2009 update: A report from the American Heart Association Statistics Committee and Stroke Statistics Subcommittee[END_REF]. Most frequently (80%) strokes result from the occlusion of one or several brain vessels and are thus called ischemic strokes (in the other cases, strokes are hemorrhagic strokes). Ischemic stroke involves many pathophysiological mechanisms causing devastating neurological damage (see for review [START_REF] Dirnagl | Pathobiology of ischaemic stroke: an integrated view[END_REF][START_REF] Moskowitz | The science of stroke: mechanisms in search of treatments[END_REF]). Understanding these mechanisms is of the most importance to develop new therapeutic strategies since no treatments are currently available for most stroke patients. Currently, the only FDA-approved treatment for stroke patients is a thrombolytic agent (tPA) which can only be given to less than 10% of patients because of its narrow time-window and its hemorrhagic risks [START_REF] Graham | Tissue plasminogen activator for acute ischemic stroke[END_REF]. Many neuroprotective agents (aimed at blocking the ischemic cascade) have also been developed but, although they had given very promising results in preclinical studies in rodent models, they appeared ineffective or even noxious during the clinical trials in stroke patients (see for review [START_REF] Keyser | Clinical trials with neuroprotective drugs in acute ischaemic stroke: are we doing the right thing?[END_REF][START_REF] Barber | Understanding and managing ischemic stroke[END_REF][START_REF] Wahlgren | Neuroprotection in cerebral ischaemia: facts and fancies-the need for new approaches[END_REF][START_REF] Durukan | Acute ischemic stroke: Overview of major experimental rodent models, pathophysiology, and therapy of focal cerebral ischemia[END_REF]). This discrepancy between the results in rodents and in humans is partly due to the anatomic and histological differences between rodent and human brains. In this case, results in rodents are thus difficult to extrapolate to stroke patients. As a consequence, a mathematical model and its numerical simulations can help both to test some biological hypotheses concerning the involved mechanisms and to give new insights concerning the effects of these neuroprotective agents.

Previous works have been conducted on stroke modeling. One of these models [START_REF] Dronne | A mathematical model of ion movements in grey matter during a stroke[END_REF] is focused on the main mechanisms leading to cell death during the first hour of an ischemic stroke (such as ionic movements, glutamate excitotoxicity and cytotoxic edema). This model is based on a system of ordinary differential equations (ODEs) and is mainly an electrophysiological model. It describes the dynamics of membrane potentials, cell volumes and ionic concentrations (K + , Na + , Cl -, Ca 2+ and Glu -) in brain cells and in the extracellular space during a stroke. This model was used to study the role of various cell types during ischemia [START_REF] Dronne | Role of astrocytes in grey matter during a stroke: A modelling approach[END_REF] and to explore the effects of various neuroprotective agents in stroke patients [START_REF] Dronne | Examples of the influence of the geometry on the propagation of progressive waves[END_REF]. Other models have been developed to simulate and study spreading depressions during a stroke. This phenomenon is characterized by a slowly propagating depolarization of brain cells along with drastic disruption of ionic gradients [START_REF] Somjen | Mechanisms of spreading depression and hypoxic spreading depression-like depolarization[END_REF]. These spreading depressions have recently been observed in stroke patients [START_REF] Dohmen | Spreading depolarizations occur in human ischemic stroke with high incidence[END_REF] and are supposed to extent the ischemic damage [START_REF] Shin | Vasoconstrictive neurovascular coupling during focal ischemic depolarizations[END_REF]. Some models reproduce and study the behavior of spreading depressions in neuronal cells [START_REF] Revett | Spreading depression in focal ischemia: a computational study[END_REF][START_REF] Kager | Simulated seizures and spreading depression in a neuron model incorporating interstitial space and ion concentrations[END_REF]. Others describe these depolarization waves though neuronal and glial cells [START_REF] Dronne | A modelling approach to explore some hypotheses of the failure of neuroprotective trials in ischemic stroke patient[END_REF]. Other models study the influence of the human brain cortex geometry on the propagation of these spreading depressions [START_REF] Dronne | A modelling approach to explore some hypotheses of the failure of neuroprotective trials in ischemic stroke patient[END_REF][START_REF] Grenier | A numerical study of the blocking of migraine by Rolando sulcus[END_REF]. All these models are based on reaction-diffusion systems and in this paper we choose to use the mathematical model [START_REF] Dronne | A mathematical model of ion movements in grey matter during a stroke[END_REF].

The final goal of our work is to utterly describe and reproduce precocious mechanisms of stroke (i.e. ionic movements, glutamate excitotoxicity and cytotoxic edema) including the spreading depressions, for a realistic brain geometry. A first description of the algorithms used for the numerical solution of this stroke model on 1D and 2D geometries was presented in a previous article [START_REF] Descombes | Numerical simulation of a stroke: Computational problems and methodology[END_REF]. However, since we need to take into account the anatomic and histological specificities of human brain, this model must be simulated on a 3D realistic geometry, which implies to develop powerful numerical methods able to deal with a broad spectrum of spatial and temporal scales. This paper focuses on the methods developed for the numerical solution of this model, with much more insights on the mathematical and numerical methods than in [START_REF] Descombes | Numerical simulation of a stroke: Computational problems and methodology[END_REF]. The numerical method is based on operator splitting and explicit/implicit Runge-Kutta methods. A very important feature of this method is that no linear system (of large size) is solved. We then show, for the first time, numerical simulations in 3D obtained thanks to a particular implementation of parallelism in the framework of shared memory machines. Moreover, these 3D simulations are computed on realistic geometries, obtained from MRI of the human brain, on conventional computational resources, that is on nowadays reasonably powerful workstations; and they are shown to match the observed dynamics from MRI images in stroke patient. Since accuracy in 3D simulations is not yet optimal, the ability of extending the proposed numerical strategy to adaptive multiresolution is presented in the framework of preliminary computations in simple geometries, based on a strategy introduced in [START_REF] Duarte | New resolution strategy for multi-scale reaction waves using time operator splitting, space adaptive multiresolution and dedicated high order implicit/explicit time integrators[END_REF]. The idea is to increase the level of accuracy in order to match all the spatial scales, with a better computational efficiency; thanks to the fact that phenomenons in strokes are spatially localized, a local mesh adaptation (like multiresolution techniques) is the most suitable.

The paper is organized as follows: in a first part, we present the reaction-diffusion model of the precocious mechanisms. We then focus on numerical methods: we first mention the different approaches which can be used to discretize the system in time and explain why in the context of such a stiff and large system only very few are relevant. We then present our numerical strategy based on splitting methods; a grid adaptation technique is also proposed as a possible improvement of the numerical strategy, considering particular features of the phenomena. We present the parallel implementation on shared memory machines of the numerical strategy, and discuss the numerical validation of the results. In the next section, 2D and 3D numerical results of simulations with complex geometry are presented. Biological results obtained are compared with real observations and discussed in the penultimate section. Biomarkers are used in order to validate these computations. A brief and prospective study based on coupling the proposed strategy with adaptive multiresolution in space is conducted, whereas conclusion and future works are presented in the last section.

Stroke modeling through stiff Reaction-Diffusion systems

In this section, we describe the model on which our study is based. This model includes ionic movements, glutamate excitotoxicity, cytotoxic edema and spreading depressions [START_REF] Dronne | A mathematical model of ion movements in grey matter during a stroke[END_REF][START_REF] Dronne | Role of astrocytes in grey matter during a stroke: A modelling approach[END_REF]. It thus focuses on the first hour of a stroke, when the ionic exchanges are the main mechanisms leading to cell death. This model is based on a reaction-diffusion system (equations are given in what follows in Table 1).

In this model, brain tissue is composed of two cell types, namely neurons and glial cells, and of extracellular space. Two domains are considered: the white and the gray matter which differ in their glial cell composition (astrocytes in gray matter and oligodendrocytes in white matter) and in their "neuronal area" composition (neuronal somas in gray matter and neuronal axons in white matter). Human brain cortex is exclusively composed of gray matter whereas human brain space is mainly composed of white matter (except the gray kernels). For simplicity reasons, we consider in the model that brain cortex contains only gray matter and brain space contains only white matter. The ionic species considered in this model are K + , Na + , Cl -, Ca 2+ and the Glutamate (glu). They pass through neuronal and glial membranes via ionic channels (such as voltage-gated channels, receptor-channels, stretch-channels) and via ionic pumps and transporters (which are energy-dependent) (see figure 1). The ionic exchanges through voltagegated channels have been first modeled by Hodgkin and Huxley [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF].

Neuron (soma) Astrocyte

Extracellular space

3Na + 2K + Ca 2+ Cl - pump Ca 2+ pump Cl - Cl - Ca 2+ pump Ca 2+ pump Cl - Ca 2+

ATP ATP

Grey matter

gapjunctions

Figure 1: Ionic exchanges in gray matter between neurons, astrocytes and the extracellular space through voltage-gated channels, ionic transporters, receptor-channels and ionic pumps (from [START_REF] Dronne | A mathematical model of ion movements in grey matter during a stroke[END_REF]).

The main precocious mechanisms of ischemic stroke can be described as follows (see for review [START_REF] Dirnagl | Pathobiology of ischaemic stroke: an integrated view[END_REF][START_REF] Moskowitz | The science of stroke: mechanisms in search of treatments[END_REF]): after the stroke onset, the cells in the ischemic area do not receive enough oxygen to maintain their production of energy. As a consequence, the activity of the ionic pumps decreases, which results in variations of ionic concentrations in the cells and in the extracellular space. These ionic variations have several consequences:

• the alteration of membrane potentials, resulting in membrane depolarization and in the opening of the voltage-gated channels;

• the cell swelling due to water influx;

• the increase of the neuronal concentration of Ca 2+ , resulting in enzyme activation and leading cells towards necrosis;

• the increase of glutamate in the extracellular space, reinforcing the excitotoxic process;

• the increase of the concentration of K + propagating in the extracellular space and the increase of Ca 2+ in the astrocytic synticium, creating waves of cortical spreading depressions, opening further ionic channels and thus expanding the ischemic damage far from the ischemic core.

From this, we can understand the importance of studying these propagation phenomena and of exploring the potential effects of some neuroprotectors which modulate or block specific voltage-gated channels. Consequently, the model considers the following variables:

• the volume fractions f n and f a (by brain volume unit) of neurons and glial cells. The fraction of volume remaining for the extracellular space is thus 1f nf a . f n and f a ;

• the membrane potentials V n and V a of neurons and glial cells (taking zero as reference potential in the extracellular space);

• the concentrations of K + , Na + , Cl -, Ca 2+ and Glutamate in the 3 spaces (neurons, glial cells and extracellular space).

All the variables depend both from time and coordinates. Altogether, the mean field model has m = 19 unknowns written as a reaction-diffusion system of equations. However, there is no diffusion for 4 unknowns, namely f n , f a , V n and V a and not all ion concentrations diffuse in gray matter and in white matter. Since gray matter contains astrocytes (which are linked into an astrocytic syncytium thanks to gap-junctions), ions are able to diffuse in the astrocytic space as well as in the extracellular space in gray matter. On the other side, as the main glial cells in white matter are oligodendrocytes (which do not have the same properties as astrocytes), ions are considered to be only able to diffuse in extracellular space in white matter. As a consequence, the model contains 10 reaction-diffusion equations in gray matter (for the concentrations of K + , Na + , Cl -, Ca 2+ and Glu -in astrocytes and in the extracellular space) and 5 reaction-diffusion equations in white matter (for the concentrations of (K + , Na + , Cl -, Ca 2+ and Glu -in the extracellular space). To be more precise, in Table 1 we summarize the modeling equations of a human ischemic stroke, developed by Dronne et al. in [START_REF] Dronne | A mathematical model of ion movements in grey matter during a stroke[END_REF], using the following notations:

• S stands for the ionic species (Ca 2+ , Na + , K + , Cl -and Glu -respectively).

• The subscript e stands for extracellular, n for neuronal and a for glial medium (astrocytes in grey matter and oligodendrocytes in white matter).

• ε s,a and ε s,e are the diffusion coefficients for each ion "S " in glial cells and in the extracellular space respectively. In white matter ε s,a = 0.

• It s,n and It s,a are the global ionic currents for each ion S through neuronal membrane and through glial membrane respectively. For example in grey matter, It s,n is the sum of all the currents concerning the transport of ion S through neuronal membrane via the voltagegated channels, transporters and receptor-channels represented in figure 1. These current equations mainly rely on Hodgkin-Huxley equations and come from other electrophysiological models (neuronal and glial models: [START_REF] Walz | Role of Na/K/Cl cotransport in astrocytes, Can[END_REF][START_REF] Destexhe | Kinetic models of synaptic transmission[END_REF][START_REF] Yamada | Methods in Neuronal Modeling: From Ions to Networks[END_REF][START_REF] Rossi | Glutamate release in severe brain ischaemia is mainly by reversed uptake[END_REF][START_REF] Shapiro | Osmotic forces and gap junctions in spreading depression: a computational model[END_REF] and cardiac models [START_REF] Difrancesco | A model of cardiac electrical activity incorporating ionic pumps and concentration changes[END_REF][START_REF] Lemieux | Modeling the dynamic features of the electrogenic Na,K pump of cardiac cells[END_REF] for some currents.

Other parameters are given by: N imp,a : number of moles of impermeant anions in the glial cells (constant) N imp,n : number of moles of impermeant anions in the neuron (constant) n n : number of neurons in each volume unit n a : number of glial cells in each volume unit σ n : neuron surface σ a : glial cells surface z s : valence of ion S v : volume of each unit c n : neuron capacity c a : glial capacity F : Faraday's constant For α n and α a we follow [START_REF] Yi | A mathematical study of volume shifts and ionic concentration changes during ischemia and hypoxia[END_REF]: α n = α a = L i RT/v, with L i = 1.21 × 10 12 cm 3 /(Pa • min), R = 8.3145 J/(mol • K), and T = 310.15 K.

Variables Equations

Ionic concentrations

S n ∂S n ∂t = - n n σ n It s,n f n z s F v - S n f n ∂ f n ∂t S a ∂S a ∂t = div(ε s,a grad S a ) - n a σ a It s,a f a z s F v - S a f a ∂ f a ∂t S e ∂S e ∂t = div(ε s,e grad S e ) + n n σ n It s,n + n a σ a It s,a (1 -f n -f a ) z s F v - S e 1 -f n -f a ∂(1 -f n -f a ) ∂t
Proportions of intracellular volumes

f n ∂ f n ∂t = α n S n - S e + N imp,n v f n f a ∂ f a ∂t = α a S a - S e + N imp,a v f a Membrane potentials V n ∂V n ∂t = - σ n c n s It s,n V a ∂V a ∂t = - σ a c a s It s,a
Table 1: Model equations describing the dynamics of the ionic concentrations, cell volumes and membrane potentials with diffusion in neurons, in glial cells and in the extracellular space.

The PDEs in Table 1 can be written in the following general form:

∂u i ∂t (x, t) -div(ε i (x)grad u i (x, t)) = f i (u 1 (x, t), . . . , u m (x, t)), u i (x, 0) = u 0 i (x), 1 ≤ i ≤ m, x ∈ Ω.              (1)
The domain Ω corresponds to a human brain and is divided in gray and white matter. These two matters differ in several coefficients in the reaction term (corresponding to the cell composition) and in their diffusion coefficients, as previously described. There are no fluxes of ions in and out of the brain and thus, the boundary conditions are of Neumann homogeneous type. For the initial conditions u i (x, 0) = u 0 i (x), 1 ≤ i ≤ m, a classical medical hypothesis is that the system is in a stable equilibrium: thus we take, and must find, a stable constant solution of system [START_REF] Lloyd-Jones | Heart disease and stroke statistics-2009 update: A report from the American Heart Association Statistics Committee and Stroke Statistics Subcommittee[END_REF].

Let us mention some characteristics of the system which are very important in the choice of numerical schemes:

• The reaction term F = ( f 1 , ...., f m ) t is extremely stiff; that is to say that if we consider the system of differential equations du/dt = F(u), it is a stiff system according to the definition given in [START_REF] Hairer | Solving ordinary differential equations[END_REF]. To see this, we have performed, by numerical differentiation, a computation of the Jacobian matrix (∂ f i /∂u j ), 1 ≤ i, j ≤ m, near a stable stationary value F(u) = 0, and we found numerically negative eigenvalues with negligible imaginary parts but with real parts in the range from -10 8 to about -1. Moreover, it is impossible to separate fast and slow variables and even if this was possible, the voltage dependent gates would make this separation very local in time and space. We have to deal with the stiffness of the reactive term F, which is the core of the model and is a program of about 500 lines of C language.

• The diffusion coefficients ε i (x) are low: about 10 -3 given by a non-dimensional analysis.

The resulting splitting time step for a proper resolution of the propagating phenomenon resulting from the coupling with the reaction term will lead to the resolution of heat equation in a mildly stiff framework. Exploiting this fact turns out to be very important: as we will explain at paragraph 3.1.1 we can use stabilized explicit methods when solving the heat equation associated with the diffusion, with the advantage of good numerical performances, and an easy implementation of parallel computations.

The diffusion coefficients ε i (x) take two constant values in gray and white matter (respectively ε g i and ε w i ). The interface conditions between gray and white matter are classical:

ε g i grad u i (x, t) • n = ε w i grad u i (x, t) • n, ( 2
)
where n is a normal unit vector to the boundary between gray and white matter. These conditions become Neumann homogeneous boundary conditions whenever one of the diffusion coefficients is zero.

Numerical strategy: operator splitting and time integrators

One dimensional simulations are very useful to fit parameters such as the diffusion coefficients which are known in the literature only with limited accuracy; two dimensional ones are useful to validate numerical methods and programs, but only three dimensional simulations can be relevant from the medical point of view. From medical considerations, and also by some considerations on reaction-diffusion systems, we know that a precise description of the brain geometry is mandatory for the simulations, otherwise the plausible waves would be strongly perturbed, see for example [START_REF] Dronne | Examples of the influence of the geometry on the propagation of progressive waves[END_REF]. We then have to think of a strategy dedicated to three dimensional simulations with a very fine spatial discretization allowing to resolve the broad spectrum of spatial and temporal scales of the system (1). The method developed has to be fast, robust and must take into account the properties of the model.

We describe now the methods introduced in this work, based on a spatial discretization which will be applied in dimension 2 and 3.

Concerning the spatial discretization, we have chosen a finite volume approach with a 5 points stencil in 2D, and a 7 points stencil in 3D. Our experience is that, with uniform finite volumes, at least = 10 7 volumes are necessary for a realistic three dimensional simulation. The continuous unknown u is then replaced by a vector U belonging to R m× corresponding to the m unknowns at each point x i , 1 ≤ i ≤ . We use MRI pictures and we consider pixels as center of volumes of an uniform grid. When we apply this spatial discretization to the system (1), this yields a large system of ordinary differential equations. Let us write this system under the form

dU dt = A ε U + F(U), (3)
A ε being a matrix corresponding to the discretization of the diffusion operator; this is a classical 5 terms (resp. 7 terms) by line matrices in dimension 2 (resp. 3). We now present the different approaches which can be used to discretize this system in time and we explain why in the context of such a stiff and large system such as (1), only few are efficient. The first idea is to use directly a solver of systems of ODEs, the so called method of lines, but due to the stiffness of the nonlinear term, a large system of algebraic equations should be solved at each time step, which is too much time consuming. It is then better to use different discretizations in time for the linear and the nonlinear terms. A first method is to use an Implicit-Explicit method by treating the linear term implicitly and the nonlinear term explicitly. If we denote by δt the time step and U k the approximated solution at time kδt, the simplest method is the following:

U k+1 -U k δt + A ε U k+1 = F(U k ).
One must solve a linear system at each step since diffusion is taken implicitly but the nonlinear term is taken explicitly. This method is of order 1 in time. More accurate, but not really more expensive, methods of the same type and of order at most 6 are described and analyzed in [START_REF] Akrivis | Implicit-explicit multistep finite element methods for nonlinear parabolic problems[END_REF]. The main advantage of these methods is that only linear systems must be solved but the drawback is that, due to the explicit computation of the reaction terms, these methods are adapted only to systems with non stiff reaction terms. Let us recall that the system (1) is very stiff, and these methods can only work with time steps of the same order of the fastest time scale of the system which is about 10 -8 seconds. This would result in an prohibitive computing time, about 4 × 10 11 steps for simulating the first hour of the evolution of the stroke.

A better idea for the treatment of the linear and the nonlinear part in the context of a stiff nonlinear term is to "reverse" the numerical treatments: to solve explicitly the linear part and implicitly the nonlinear part. The discretization of the linear part is made using an explicit Runge-Kutta method with extended stability domain along the negative real axis. The papers [START_REF] Sommeijer | RKC: an explicit solver for parabolic PDEs[END_REF] and [START_REF] Verwer | An implicit-explicit Runge-Kutta-Chebyshev scheme for diffusion-reaction equations[END_REF] settled the foundation for these methods called IMEX methods and particular methods devoted to stiff non linear problems are presented in [START_REF] Verwer | RKC time-stepping for advection-diffusion-reaction problems[END_REF] and [START_REF] Shampine | IRKC: an IMEX solver for stiff diffusion-reaction PDEs[END_REF]. The main advantage of these methods is that they treat diffusion terms explicitly and the stiff reaction terms implicitly. Furthermore, the stiff reaction term is decoupled over space grids and yields small sized systems. These methods are usually very efficient; nevertheless, the computational requirements associated mainly with an implicit solver over the discretized domain with the same time step become soon critical when treating large computational domains.

Finally, the only possible methods which can solve system (1) seem to be the so called splitting methods that we describe in details now.

Splitting methods

The idea is as old as numerical analysis and was used and analyzed by the Soviet school in the 60's (see for example [START_REF] Marchuk | Splitting and alternating direction methods[END_REF]). At that time, the main interest was the economy of computer memory. The idea, applied to spatially discretized reaction-diffusion equations is to solve alternatively the reaction and the diffusion problems. For example, starting from some initial condition, we solve for a time step of δt:

dV dt = A ε (V),
with an initial condition

V(0) = V 0 .
Let us call D δt this procedure. Taking V(δt) as initial condition, we solve for the same δt:

dW dt = F(W), with W(0) = V(δt),
and we call R δt this second procedure. By taking the value W(δt), we repeat this procedure to obtain, for k > 0, W(kδt). The previous approximation is an approximation of order 1 in time of the solution of (3). Let us recall that a method is of order 1 (or more generally of order p) if the expansion in powers of δt of the numerical solution coincides with that of the true solution up to and including the order 1 (more generally the order p). The previous approximation is called a Lie method, but one can define different numerical splittings schemes: to obtain an approximation at time δt, one can apply successively R δt and D δt , or more generally apply successively R δt/2 , D δt and R δt/2 , or D δt/2 , R δt and D δt/2 . The last two approximations are called Strang methods [START_REF] Strang | On the construction and comparison of difference schemes[END_REF] and are of order 2.

Let us explain the main advantages of these methods: the reaction and diffusion are decoupled, the solution of the D δt problem is reduced to the solution of m independent diffusion equations, and thus the complexity is reduced. Concerning the R δt problems, one immediately see that they are decoupled in systems of ODEs of size m, as many systems as nodes in the finite volume mesh, and that all these systems are independent.

Assuming first that D δt and R δt can be solved exactly without time discretization, these splitting methods can be used to solve stiff systems of reaction-diffusion. Better performances are expected by ending the splitting scheme with the integration of the reaction part or more generally with the part involving the fastest time scales of the phenomenon (see [START_REF] Descombes | Operator splitting for nonlinear reaction-diffusion systems with an entropic structure: singular perturbation and order reduction[END_REF][START_REF] Descombes | Operator splitting for stiff nonlinear reaction-diffusion systems: order reduction and application to spiral waves[END_REF] and references therein).

Keeping in mind these theoretical studies and considering the various numerical alternatives previously discussed, Strang's splitting scheme ending with the reaction part remains as the most appropriate resolution scheme for general multiscale problems and so far, the best choice for our numerical study. All the numerical simulations in this article are then performed with this scheme.

Efficient choice of numerical methods for the sub-steps

The concept of splitting is very simple; we still have to describe the numerical methods used in each sub-step. The order of these methods must always be at least equal to 3, so that the dominant part of the error come from Strang's splitting, noticing that at each sub-step, we solve a Cauchy initial value problem: multistep methods like those based on backward differentiation formulae (see [START_REF] Hairer | Solving ordinary differential equations[END_REF] for details) are not adapted to splitting methods since they need more than one initial condition at each time step to perform the time integration. These initial conditions are often approximated by a less accurate procedure. Thus we chose Runge-Kutta methods for the two sub-steps:

1. For the R δt/2 sub-step, we have to solve a stiff but spatially decoupled system of ODEs. A dedicated method adapted to stiff systems of ODEs (and thus an implicit method) must be used. In this strategy, we chose an A-and L-stable scheme: the Radau5 solver [START_REF] Hairer | Solving ordinary differential equations[END_REF]. 2. For the diffusion sub-step, a stable method must be chosen. However, the value of the diffusion coefficient, the gradients of the solution as well as the value of the splitting time step will not lead to a strong stability constraint, which would require the use of an implicit method. In this context, we can use explicit methods with enlarged stability domains such as Rock4 [START_REF] Abdulle | Fourth order Chebyshev methods with recurrence relation[END_REF].

Let us emphasize that:

• For the R δt/2 sub-step, many other dedicated methods with both A-and L-stability [START_REF] Hairer | Solving ordinary differential equations[END_REF] can be used to handle the stiffness associated to the systems of ODEs: for instance, the Linearized Euler extrapolated method and Rosenbrock methods (see [START_REF] Hairer | Solving ordinary differential equations[END_REF]), but none of them is as fast and as robust as Radau5.

• One of the main advantages of the Rock4 method is that only matrix-vector products must be computed in opposition to implicit methods for which linear systems must be solved.

The number of matrix-vector products at each step in Rock4 is at least 6 and grows with the stiffness of the systems which is measured here by the products δtε i λ max (λ max is the dominant eigenvalue of the Laplacian operator). Thus, for a given diffusion operator (that is to say for given ε i λ max ) the efficiency of the Rock4 method is related to the splitting time step δt: with the time step we have chosen (as explained in section 6), we always perform the minimum number of matrix-vector products. Should δt or ε i be much larger, then the complexity of Rock4 (that is to say the number of matrix-vector products) could make it non competitive with a scheme involving the solution of linear systems. Let us recall that λ max is proportional to h -2 , h being the size of the smallest finite volume used, and thus does not depend of the spatial dimension, so that the 3 dimensional computations benefit largely from the use of the Rock4 method.

An adaptive mesh strategy: the adaptive multiresolution approach

Most software for the solution of partial differential equations use constant discretization in time and even in space. But stroke is a phenomenon which is localized in a small part of the brain, at least at the onset, and which exhibits propagating waves and large gradients. Thus, an adaptive mesh is certainly well fitted for such simulations, resulting in a better resolution of thin spatial structures and hopefully, in better numerical performances. The mesh adaptation must be automatically managed by the software. A numerical strategy has been proposed and evaluated in [START_REF] Duarte | New resolution strategy for multi-scale reaction waves using time operator splitting, space adaptive multiresolution and dedicated high order implicit/explicit time integrators[END_REF], where the combination of adaptive spatial multiresolution with the time integration solvers previously detailed, is described.

Historically, adaptive methods like Multi Level Adaptive Techniques (MLAT) [START_REF] Brandt | Multi-level adaptive solutions to boundary value problems[END_REF] or Adaptive Mesh Refinement (AMR) [START_REF] Bell | Three-dimensional adaptive mesh refinement for hyperbolic conservation laws[END_REF][START_REF] Berger | Adaptive mesh refinement for hyperbolic partial differential equations[END_REF][START_REF] Berger | Local adaptive mesh refinement for shock hydrodynamics[END_REF] were among the first to achieve this goal, using a set of locally refined grids where steep gradients are found. Furthermore, adaptive multiresolution methods, based on Harten's pioneering work [START_REF] Harten | Multiresolution algorithms for the numerical solution of hyperbolic conservation laws[END_REF], have been developed for 1D and 2D hyperbolic conservation laws [START_REF] Cohen | Fully Adaptive Multiresolution Finite Volume Schemes for Conservation Laws[END_REF][START_REF] Gottschlich-Müller | Adaptive finite volume schemes for conservation laws based on local multiresolution techniques[END_REF] and then extended to 3D parabolic problems [START_REF] Roussel | A conservative Fully Adaptive Multiresolution algorithm for parabolic PDEs[END_REF]. Consequently, high data compression might be achieved with all these methods. However, one major advantage of the adaptive multiresolution techniques is that the numerical analysis of the errors has already been conducted [START_REF] Harten | Multiresolution algorithms for the numerical solution of hyperbolic conservation laws[END_REF][START_REF] Cohen | Fully Adaptive Multiresolution Finite Volume Schemes for Conservation Laws[END_REF] and thus, a solid theoretical background has already been settled. Considering that we have a set of nested spatial grids, from the coarsest to the finest one, a multiresolution transformation allows to represent a set of function data as values on a coarser grid plus a series of differences (error estimators) at different levels of such nested grids. The information at consecutive levels are then related by inter-level transformations: which are the projection and prediction operators. The theoretical background of such configuration (see [START_REF] Cohen | Fully Adaptive Multiresolution Finite Volume Schemes for Conservation Laws[END_REF]) states that wavelet coefficients can be then defined as prediction errors, and they will retain the detail (local regularity estimator) information when going from a coarse to a finer grid (see figure 2 for an example of nested grids).

Finally, one of the main interests of carrying out such multiscale decomposition is that this new representation leads us to define a whole set of regularity estimators all over the spatial domain. Hence, an effective data compression might be performed in order to retain only a minimal quantity of nodes where it is strictly necessary. The evolution problem is then simulated on a dynamical adapted grid. The coupling between splitting techniques and multiresolution is very easy and yields an efficient algorithm even if the spatial discretization is locally very refined.

For an overview on adaptive multiresolution techniques, we refer to the books of Cohen [START_REF] Cohen | Wavelet methods in numerical analysis[END_REF] and Müller [START_REF] Müller | Adaptive multiscale schemes for conservation laws[END_REF].

Numerical software

We have developed two different softwares for the solution of the system (1). Both implement the time integration strategy defined above in 3.1: FM: (Fixed Mesh) a code using a fixed spatial discretization, with finite volumes of constant size [START_REF] Dumont | ZEBRE: Numerical Software for Reaction-Diffusion systems[END_REF]. This code takes into account complex geometries in the following way: starting from MRI images, we take each pixel as the center of a finite volume; it aims to be a framework for testing and exploiting numerical methods for 1, 2 and 3D reaction-diffusion systems. It will be used in order to obtain the main results of the present contribution, that is numerical simulations of the detailed and stiff stroke model in complex 3D geometries.

MR: (Multi Resolution) a code using an adaptive multiresolution method as defined above in 3.2. In the framework of multiresolution, an important amount of work is still required in order to optimally combine all the numerical methods described here, the most difficult aspects are related to programming features such as data and code structures, as indicated in [START_REF] Duarte | New resolution strategy for multi-scale reaction waves using time operator splitting, space adaptive multiresolution and dedicated high order implicit/explicit time integrators[END_REF]. Nowadays, this program can only solve problems in simple domains like squares and cubes; simulations with an adaptive multiresolution approach on a complex geometry are not yet available, and we will only present here 2D and 3D simulations in simplified geometries for the sake of assessing our results and perspectives in the field.

Let us remark that the (FM) code is a highly optimized and complete code for the simulation of reaction-diffusion equations. In particular, stroke simulations in complex geometry can be performed for the first time, with standard computing resources, and constitutes the major advance of our contribution. On the other hand, the second code (MR) allows to validate to some extents the previous numerical results, and it is meant to be a potential extension to (FM) in future developments.

Implementation and performances of the numerical methods on shared memory machines

We describe now the implementation of the 3D simulations which are performed on a uniform grid and a complex geometry with the code (FM). Let us emphasize the particular parallelism implementation that we have conceived in the framework of shared memory machines. All the computations have been performed on a 8 core (2x4) 64 bits machine (AMD Shanghai processors).

Since an implicit procedure is required to handle the associated stiffness, the reaction steps are by far the most time consuming parts of the computation. This step is naturally parallel, as we have to solve a large number of independent systems of differential equations, each of them corresponding to a single cell. But solvers like Radau5 use a time adaptive strategy, together with Newton method: that is, the computing time is varying from one point to another, and a fixed domain partition strategy, with an affectation of sub-domains to processors, is not optimal for load balancing. Therefore, our implementation uses threads, implemented in the C++ boostthread library [START_REF] Williams | The boost thread library. Portable C++ multi-threading[END_REF]. Let us describe it shortly in the next part.

We divide the set of finite volumes into small subsets S i , i = 1, ..., k where k should be much larger than the number of computing units. We build a stack of all the S i before each time step. By calling the procedure GetSubset, each thread gets one S i , while the stack is not empty. Threads join when the stack is empty.

Algorithm 5.1: Get a subset( ) procedure GetSubset( ) Lock l; if stack not empty f inished ← f alse; s ← stack.pop(); else f inished ← true; s ← null; return ( f inished, s)
The "Lock" object (a boost thread library object [START_REF] Williams | The boost thread library. Portable C++ multi-threading[END_REF]) is created when a thread enters the procedure, forbidding entry to other threads (this is a basic and mandatory feature of all threading libraries). The stack is popped, giving the thread a set of differential equations to be solved. When the thread exits the procedure, the lock object is automatically destroyed and another thread can enter. To minimize lock competition effects, the size of each S i is taken constant (about 100 nodes) plus a random size between zero and 10% of the constant.

Each thread is actually a procedure:

Algorithm 5.2: Thread( ) procedure ThreadProcedure( ) repeat ( f inished, s) = GetS ubset();
if not f inished for each x ∈ s do solve ODE system at x; until f inished;

And eventually, the solution of the reaction step with n threads is very simple, using the Boost library:

Algorithm 5.3: Reaction( ) procedure Reaction(n)
thread group thrds(n); for i ← 1 to n do thrds.create thread(T hreadProcedure); thrds. join all() thread group is a container, member of the boost-thread library, which manages a set of threads (n threads here). The threads are spawned with the method create thread and run in parallel, each thread launching the T hreadProcedure routine. The join all method acts as a meeting point for all the threads, and the execution waits until all threads have finished their computations.

Numerical results: implementation checkout and accuracy evaluation of the code

In order to check out the implementation of the method in the codes (FM) and (MR), which use two different spatial discretizations and data structures, we have conducted a detailed comparison on a numerical test-case. Nevertheless, since we use the same numerical methods for the sub-steps integration, this does not result in a full validation. We have considered a 2D case in a regular geometry of [0, 5] × [0, 5] (cm), using two resolutions, one on a fixed grid computed with the code (FM) and the other on the adapted grid obtained by multiresolution with the code (MR). Both methods, based on Strang's splitting, use Radau5 and Rock4 as time integrators for the reaction and for the diffusion problem. For the model parameters, we have considered only one domain, the gray matter. The time of integration was restricted to one hour, t ∈ [0, 3600] seconds. The splitting time step δt was chosen experimentally in order to obtain a good approximation of the velocity on a sufficiently 2D fine space discretization of 1024 2 . We thus compare both solutions with a L2 -norm1 at final time t = 3600, and define an error estimator e given by u(., t)u MR (., t) L 2 , where index MR denotes the use of multiresolution techniques. The spatial discretization consists of an uniform mesh of 256 2 points. In the case of the spatial adaptive method, there is a set of nested grids arranged in 8 different levels from the coarsest to the finest discretization. The latter corresponds to the uniform mesh previously considered of 256 2 points. In order to compare these results, we must consider the same spatial discretization for both solutions: this is easily achieved with the mentioned projection/prediction operations on the adapted grid. Tables 2, 3 and 4 show the minimum and maximum values of variables in the neurons, the astrocytes and in the extracellular space respectively, as well as the normalized L 2 difference of the numerical results e obtained by the adaptive multiresolution strategy and our proposed numerical strategy.

Variable

Figure 3 shows the evolution of the propagating phenomenon on an adapted grid for variable K + in the neurons. The refined regions clearly correspond to the wavefront area where the steep spatial gradients are present. Finally, figure 4 reveals in a qualitative way the different representations of the numerical solution on an uniform mesh and on the adapted one.

All these numerical results show a great accordance between the solutions of the two different codes: (MR) and (FM), in the 2D simulations. Let us recall that both codes rely on two well tested, robust and publicly available numerical routines: the Radau5 and Rock4 methods; therefore, one can consider that they only differ by the different spatial discretizations and data structures they use: the comparison can thus be considered as a (partial) cross validation.

Besides such a level of comparisons, (MR) will be shown to pave the way towards higher levels of refinement for a better resolution of the details of the dynamics in 3D at a reasonable cost. In fact, considering this 2D numerical test-case, finer spatial discretizations yield naturally better resolution of both the wave velocity and the dynamics of the wavefront, as seen into figure 5. See the corresponding (MR) adapted grids into figure 6. Even if it is clear that the wave is better resolved on the finest grid (1024 2 ), in particular it is somewhat faster, the qualitative value of the wave velocity is correctly captured even on the coarsest grid (256 2 ), which corresponds roughly to the 3D simulations with (FM). However, we can not yet simulate with (MR) all the complex geometries of the brain we are investigating in this paper; and thus, all the results in complex geometries will be performed with (FM).

Let us make now some comments on the performance of the numerical method on shared memory machines for both diffusion and reaction equation solvers in (FM) used in the next results section.

Performances of the diffusion equation solver

Only a poor man's parallelism is implemented for this step, each diffusion equation being solved by one thread. But actually, as we will see later, the computing time of this step is less than 10% of the total computing time (in dimension 2 or 3); each step, for one diffusion equation, needs only 6 matrix vector products (5 being the minimum for the Rock4 method, plus one for the error estimate).

Performances of the reaction solver

The main question concerns the efficiency of the multithreaded parallelism. Figure 7 shows the computing time with 1, 2, 4 and 8 threads, across 12 steps. The wall clock computing time is multiplied by the number of threads (unit is in number of CPU clock ticks); in case of perfect scalability, all the points with the same number of threads should be at the same ordinate. This is roughly the case, considering that the comparison of computations obtained with 1 and 8 threads reveals a loss of scalability of only 6%. We conclude that this implementation is very efficient. This is a consequence of the complexity of the right hand side (the reaction term): even one single numerical evaluation of F is much more time consuming than the overhead introduced by the thread mechanism. As a conclusion of this part, we can notice that our computing strategy combining splitting techniques with dedicated integration of each sub-step and multiresolution is compatible with parallelization.

Biological results

We present and discuss here some simulation results obtained with the code FM on the complex geometry of the human brain. We simulate an ischemic stroke beginning in the cortex (in gray matter) and study the propagation of the ischemic damage. The input of the model is the decrease of the ionic currents through the ionic pumps. Two variables have been chosen for the model validation: the potassium concentration in the extracellular space ([K + ] e ) and the ratio of apparent diffusion coefficient of water (rADCw).

• The potassium concentration cannot be measured in vivo in the brain of stroke patients but it can be measured ex vivo or in vitro on brain tissues. These concentration values give some insights on the severity of the damage. The physiological value of [K + ] e is about 5 mM. It was observed to be able to increase up to 35 mM in areas of moderate ischemia where depolarization waves can spread [START_REF] Kraig | Extracellular ionic variations during spreading depression[END_REF] and up to 75-90 mM in areas of severe ischemia where most cells are dead [START_REF] Hansen | The extracellular potassium concentration in brain cortex following ischemia in hypo-and hyperglycemic rats[END_REF]. The first step of the model validation is thus to compare the values of the [K + ] e obtained in the simulations with those values.

• The rADCw is a biomarker which can be estimated in the brain of stroke patients thanks to diffusion-weighted (DW) magnetic resonance (MRI) imaging. It reflects the severity of the cytotoxic edema and could thus be used to predict the ischemic damage and its extension [START_REF] Verheul | Comparison of diffusion-weighted MRI with changes in cell volume in a rat model of brain injury[END_REF][START_REF] Rosso | Prediction of infarct growth based on apparent diffusion coefficients: penumbral assessment without intravenous contrast material[END_REF]. The value of this ratio is supposed to be 1 in physiological conditions and is known to decrease in ischemic areas. In several studies, this value in stroke patients was shown to be between 0.75 and 0.9 in areas of moderate ischemia and between 0.5 and 0.75 in areas of severe ischemia [START_REF] Desmond | The value of apparent diffusion coefficient maps in early cerebral ischemia[END_REF][START_REF] Fiehler | Apparent diffusion coefficient decreases and magnetic resonance imaging perfusion parameters are associated in ischemic tissue of acute stroke patients[END_REF][START_REF] Røhl | Viability thresholds of ischemic penumbra of hyperacute stroke defined by perfusion-weighted MRI and apparent diffusion coefficient[END_REF][START_REF] Sakoh | Prediction of tissue survival after middle cerebral artery occlusion based on changes in the apparent diffusion of water[END_REF]. This biomarker can be related to the proportions of the intracellular volumes. It was shown to be proportional to the volume of the extracellular space [START_REF] Verheul | Comparison of diffusion-weighted MRI with changes in cell volume in a rat model of brain injury[END_REF]. Moreover, since the extracellular proportion was displayed to have a value of 0.2 in physiological conditions (i.e. when rADCw=1) [START_REF] Mcbain | Regional variation of extracellular space in the hippocampus[END_REF], rADCw can be expressed as follows: rADCw = 5 (1f nf a ). Since f n and f a are two variables of the model, this ratio can be calculated for each time and for each coordinate. Another step of the model validation is thus to compare the calculated values of rADCw obtained in the simulations to the experimental values.

We present in figures 8, 9, 10, 11, 12, and 13 some results of 2D and 3D simulations, showing the values of K + and rADCw biomarker in the extracellular space in different areas.

Let us make some biological comments about these results:

• First of all, we obtained depolarization waves after the simulation of a vessel occlusion in brain cortex and the depolarization waves spread in gray matter (i.e. in brain cortex) and not in white matter, which is consistent with MRI images obtained in human brain [START_REF] Dohmen | Spreading depolarizations occur in human ischemic stroke with high incidence[END_REF].

• Concerning the potassium concentration in the extracellular space, figures 8, 9 and 10 show that this concentration reaches values such as 77 mM in the areas where the vessel was occluded, which is in accordance with the results obtained in the infarcted core in some experimental studies [START_REF] Hansen | The extracellular potassium concentration in brain cortex following ischemia in hypo-and hyperglycemic rats[END_REF]. These figures also show that [K + ] e is about 20 mM in the areas where depolarization waves are spreading, which is consistent with the values obtained in penumbra (i.e. areas of moderate ischemia, able to recover) during spreading depressions in several studies [START_REF] Kraig | Extracellular ionic variations during spreading depression[END_REF]. We can also notice that, in the safe areas, [K + ] e remains at its physiological value (5 mM).

• Concerning rADCw, figures 11, 12 and 13 show that this ratio reaches values such as 0.6 in the areas where the vessel was occluded, which is in accordance with the values observed on MRI images in the infarcted core of the brain of the stroke patient [START_REF] Fiehler | Apparent diffusion coefficient decreases and magnetic resonance imaging perfusion parameters are associated in ischemic tissue of acute stroke patients[END_REF][START_REF] Røhl | Viability thresholds of ischemic penumbra of hyperacute stroke defined by perfusion-weighted MRI and apparent diffusion coefficient[END_REF][START_REF] Sakoh | Prediction of tissue survival after middle cerebral artery occlusion based on changes in the apparent diffusion of water[END_REF]. These figures also show that rADCw has values between 0.75 and 0.9 in the areas where depolarization waves are spreading, which is consistent with the values obtained in penumbra during spreading depressions in stroke patients [START_REF] Fiehler | Apparent diffusion coefficient decreases and magnetic resonance imaging perfusion parameters are associated in ischemic tissue of acute stroke patients[END_REF][START_REF] Sakoh | Prediction of tissue survival after middle cerebral artery occlusion based on changes in the apparent diffusion of water[END_REF]. We can also notice that, in the safe areas, rADCw remains at its physiological value of 1.

To conclude, the simulation results concerning the localization of spreading depressions and the values of [K + ] e and rADCw are consistent with those obtained in experimental studies or observed on MRI images in stroke patients. These results give thus a first step of validation for the model and for the numerical methods used in this study.

Toward better computational efficiency and improved accuracy: Adaptive Multiresolution

In the previous simulations, we notice that the simulated waves spread at a slightly slower speed. In several studies, spreading depressions were shown to spread at a rate of several millimeters per minute [START_REF] Martins-Ferreira | Perspectives on spreading depression[END_REF], which is not currently the case in our simulations. In fact, it is shown in [START_REF] Keener | Propagation and its failure in coupled systems of discrete excitable cells[END_REF] that traveling waves solutions of reaction-diffusion equations can disappear in the numerical solution if the spatial discretization is too coarse; the velocity of the traveling waves is a function of the mesh size, and coarse meshes might perturb the accuracy of the computed wave velocity. In particular, in the previous 3D simulations, the mesh we can use is not fine enough to obtain a correct level of accuracy for the wave velocities. In fact, coming back to the 2D numerical testcase of section 6, we have seen in figure 5, that a high number of volumes is needed to reproduce accurately the phenomenon, approximately 1000 per dimension.

We can also measure the computing time of the reaction at a typical step (see figure 14 (one clock tick is about 0.35 10 -9 second)). Clearly the most expensive nodes are about 37 times more expensive than the less ones! But on the other hand, the overwhelming part of the nodes are not expensive ones. Actually, 79% of the nodes (the less expensive ones, which cost less than 4.5 10 6 tics) takes 60% of the computing time.

Therefore, one way to improve both the performances and the accuracy of the resolution is to use an adaptive mesh: use a fine mesh in the ischemized zone, where the solution exhibits large gradients and wave propagation, and a coarser mesh far from this part of the domain: the multiresolution strategy, as implemented in the code (MR) is a step towards this goal.

In order to make this more concrete than just a statement, we eventually present here a first 3D multiresolution simulation in a cube, where the finest grid available has size 256 3 , since complex geometry is not yet at hand (figure 15). Therefore, considering the same computing resources, the computing time is reduced by a factor of about five with respect to the fixed mesh simulation with (FM) on the finest grid, even though trying to compare two so heterogeneous codes is a very delicate task. In fact, the code (MR) considers a very low percentage (≤ 10%) of 256 3 , which coupled with what has been presented before, allows to explain the gain in CPU time. Nevertheless, at the current state of development of the code (MR), a lot of work remains to be done like taking into account complex geometry, improving performances of multiresolution methods by using adapted data structures and optimized routines, and finally, a high performance distributed parallelism implementation. Once this is achieved, an adaptive mesh approach shall overcome the natural limitations of accuracy and performance of even very performing strategies such as the implemented in (FM), and will allow one to solve entirely the multiscale dynamics of this kind of phenomena.

Finally, increasing the accuracy of the discretization, of the numerical methods and even of the implementation, is not sufficient to generate more precise and predictive simulations. This should be carefully coupled with the development of a more precise modeling of coefficients and boundary conditions, and it is beyond the scope of the present paper.

Conclusions and future work

We have presented for the first time numerical 3D simulations of an ischemic stroke in a realistic brain geometry, based on the model of Dronne et al. [START_REF] Dronne | A mathematical model of ion movements in grey matter during a stroke[END_REF]. Results are encouraging from numerical and medical points of view. This is a first major step towards an usable tool for predicting the evolution of a stroke. The next steps are to improve both numerical performances and modeling. For this, a lot of work remains to be done from the model to practical implementations. Concerning the numerical methods, many parameters in the model are known only with a coarse approximation. Thus, numerical simulations must be conducted to explore the sensitivity of the model to these parameters. The ultimate way to improve the performances is to switch from multithreaded parallelism to distributed parallelism, on massive parallel computers.

From a medical point of view, this model is of the most importance since it could be used to simulate on a realistic human brain geometry several neuroprotective agents aimed at blocking the ischemic cascade and at reducing the ischemic damage. Since the model contains many pharmacological targets (such as ionic transporters, voltage-gated channels, channel-receptors and stretch channels), it could be used to assess and study the effects of various therapeutic agents or associations of therapeutic agents. Moreover, since the model includes both ionic movements through the cells and their diffusion, we will be able to study the effects of these neuroprotective agents both on the severity and on the extension of the damage in each brain area. Developing powerful numerical methods are thus of the most importance to be able to simulate the time and spatial evolutions of these phenomena on a realistic human brain geometry.

Figure 2 :

 2 Figure 2: Example of nested grids for multiresolution adaptive technique.

Figure 3 :

 3 Figure 3: Evolution of K + in the neurons at t = 100 (top left), t = 1000 (top right), t = 2000 (bottom left) and t = 3000 (bottom right).

Figure 4 :

 4 Figure 4: K + in the neurons at t = 3600 on a 256 2 uniform mesh (left) and the corresponding adapted grid (right).

Figure 5 :

 5 Figure 5: K + in the neurons at 3600 s. for a 2D mesh of 256 2 (top), 512 2 (center) and 1024 2 (bottom).

Figure 6 :

 6 Figure 6: 2D adapted meshes equivalent to 256 2 (left) and 1024 2 (right) spatial discretizations at the finest grid.

Figure 7 :

 7 Figure 7: Performances of the multithreaded reaction solver along 12 time steps. Abscissa: time step. Ordinate: computing time in CPU clock ticks.

Figure 8 :Figure 9 :

 89 Figure 8: Evolution of K + in the extracellular space over one hour, 2D simulation (in millimolar (mM)), from left to right, top to bottom.

Figure 10 :

 10 Figure 10: Evolution of K + in the extracellular space, 3D simulation. View in the plane P of figure 9 (in mM).

Figure 11 :Figure 12 :Figure 13 :

 111213 Figure 11: Evolution of rADCw over one hour, 2D simulation, from left to right, top to bottom, every 15 minutes.

Figure 14 :

 14 Figure 14: Local measurement of the computational cost of the reaction (in CPU clock tics).

Figure 15 :

 15 Figure 15: 3D simulations with (MR). K + in the neurons (left) and corresponding adapted grid (right) at 1000 (top) and 3600 (bottom) seconds.

Ca 2+ voltage-gated channel (CaHVA)

			pump						pump	2K +
			Na + /K +						Na + /K +
									3Na +
		Ca 2+	3Na + exchanger Na + /Ca 2+	exchanger Na + /Ca 2+	3Na +	Ca 2+
								Na + voltage-gated
				Na + Na + voltage-gated channel (NaP)	channel (NaP)	Na +
								Ca 2+ voltage-gated channel (CaHVA)	Ca 2+
		K +	channel (KDR, BK) K + voltage-gated	K + voltage-gated channel (KDR, BK, Kir)	K + K +
		K +	glu Na +	glutamate transporter	glutamate transporter	glu Na +
		K +	Ca 2+ Na +	receptor NMDA	receptor AMPA	Na +	K +
		K +	Na +	receptor AMPA	contransporter Na + /K + /Cl -	Na + 2Cl -K +
								exchanger	HCO 3	-
	glu	HCO 3	-	Cl -	exchanger -Cl -/HCO 3	Cl -/HCO 3 -	Cl -	glu
	Cl -	extra currents						extra currents	Cl -

Table 2 :

 2 Minimum and maximum values of variables in the neurons, computed with the (MR) code, and normalized L 2 difference e of numerical results between uniform mesh and (MR). t= 3600 s.

	Variable	Min value	Max value	e
	K +	124.2309	132.6962	1.2265 × 10 -3
	Na +	15.0751	24.3063	6.3289 × 10 -3
	Ca 2+	1.5921 × 10 -4	0.4149	3.1653 × 10 -2
	Cl -	6.7503	10.8147	1.0143 × 10 -2
	glu	2.5460	2.9870	1.0133 × 10 -3
	V a	-75.2476	-19.6358	1.0817 × 10 -2
	f a	0.7128	0.8088	1.1877 × 10 -3

Table 3 :

 3 Minimum and maximum values of variables in the astrocytes, computed with the (MR) code, and normalized L 2 difference e of numerical results between uniform mesh and (MR). t= 3600 s.

Table 4 :

 4 Minimum and maximum values of variables in the extracellular space, computed with the (MR) code, and L 2 normalized difference e of numerical results between uniform mesh and (MR). t= 3600 s.

The L

-norm f L 2 of f is given by Ω f (x) 2 dx 1/2 .

Acknowledgments

This research was supported by two ANR grants (Agence Nationale pour la Recherche -French National Research Agency): AVC in silico (ANR-06-BYOS-0002-02 -projects leaders: M.-A. Dronne and E. Grenier -2006Grenier - -2009) ) and Séchelles (ANR-09-BLAN-0075-01 -project leader: S. Descombes -2009Descombes - -2013)). We would like to thank Christian Tenaud (LIMSI -CNRS) for providing the basis of the multiresolution kernel of MR CHORUS, code that he developed at LIMSI for compressible Navier-Stokes equations.


Vincent Dubost

Tristan Cren

François Debontridder

Dimitri Roditchev

Cristian Vaju

Vincent Guiot

Laurent Cario

Benoît Corraze

Etienne Janod

Electric pulse induced electronic patchwork in the Mott insulator GaTa 4 Se 8

Keywords: numbers: 74.25.Gz, 74.72.Jt, 75.30.Fv, 75.40.-s

Following a recent discovery of the Insulator-to-Metal Transition induced by electric field in GaTa4Se8, we performed a detailed Scanning Tunneling Microscopy/Spectroscopy study of both pristine (insulating) and transited (conducting) crystals of this narrow gap Mott insulator. The spectroscopic maps show that pristine samples are spatially homogeneous insulators while the transited samples reveal at nanometer scale a complex electronic pattern that consists of metallic and superinsulating patches immersed in the pristine insulating matrix. Surprisingly, both kinds of patches are accompanied by a strong local topographic inflation, thus evidencing for a strong electron-lattice coupling involved in this metal-insulator transition. Finally, using a strong electric field generated across the STM tunneling junction, we demonstrate the possibility to trig the metal-insulator transition locally even at room temperature.

Many correlated materials exhibit exotic properties like metal-insulator transitions (MIT) in transitionmetal oxides or chalcogenides, high-Tc superconductivity in cuprates, colossal magnetoresistance in manganites. These phenomena are frequently associated with spatial electronic inhomogeneities or even electronic phase separation. The presence of patches with different magnetic, metallic or insulating nature was demonstrated in many real-space studies of systems such as manganites [1][2][3], cuprates superconductors [4][5][6][7][8], VO 2 [9][10][11][12], Ca substituted Sr 2 RuO 4 [13]. Generally, in these systems, the MIT is induced by tuning either the hopping amplitude via pressure (bandwidth-controlled transition) or bandfilling via chemical doping (doping-controlled transition) [16]. In addition to these ways of controlling the electronic state of correlated materials, we recently discovered that short electric field pulses could also drive the MIT in the Mott insulator GaTa 4 Se 8 [17][18][19]. Such a new way to electrically induce a MIT might be interesting for electronic devices, in particular if this MIT can be controlled at the nanometer scale at room temperature.

The Mott insulator GaTa 4 Se 8 adopts a deficient spinel structure (Fig. 1a) of rocksalt-type packing of Ta 4 Se 4 cubanes and GaSe 4 tetrahedra [20,21]. This compound belongs to a specific class of Mott insulators in which correlations do not take place at the atomic scale but rather at the scale of the Ta 4 clusters [22]. This makes these compounds easier to tune through the Metal-Insulator-Transition as it reduces the magnitude of Coulomb repulsion. GaTa 4 Se 8 is a rare example of stoichiometric Mott insulator that undergoes a pressure-induced MIT and even a superconductor transition [20]. Recently we found that an electric field pulse through GaTa 4 Se 8 induces a MIT (see Fig. 1b), accompanied at low temper-ature by granular superconductivity [17,23]. While the microscopic origin of this effect remains undecided, it was experimentally found in [23] that the temperature evolution of the resistance on the metallic side of the MIT in single crystals of GaTa 4 Se 8 can be accounted for by a two channel model of a granular metal in an insulating matrix. In this Letter report a detailed Scanning Tunneling Microscopy and Spectroscopy (STM/STS) study of this material on both sides of the MIT. On the insulating side we observed spatially homogeneous tunneling characteristics whereas on the metallic side of the transition the samples revealed nanometer scale electronic inhomogeneities accompanied by a strong local deformation. Moreover, we demonstrate the possibility to induce the MIT locally using the electric field generated across the STM junction.

Single crystals of GaTa 4 Se 8 of typical sizes around 300 µm were obtained by the selenium transport method [17]. The main steps for the preparation of the transited crystal used in this STM/STS study are summarized in Fig. 1b-d. First, two electric contacts were glued onto small GaTa 4 Se 8 single crystals of ∼300 microns (Fig. 1c). As expected for a Mott insulator, the resistance of the pristine sample displays a semiconducting/insulating behavior (see red curve in Fig. 1b). Application at 60 K of an electric field of about 30 kV/cm to this crystal induces a non-volatile MIT. After the electric pulse, the transited crystals present a low resistance state (see blue curve in Fig. 1b).

The crystals were then cleaved with a blade as sketched in Fig. 1d in order to obtain clean (100) surfaces. The cleavage plane was chosen in order to have the external electrode axis parallel to the cleaved surface. After the cleavage the samples were immediately introduced to the ultrahigh vacuum chamber (base pressure below 3.10 -11 mbar) holding the STM unit. The STM/STS experiments were provided at room temperature using electrochemically etched tungsten tips, thermally annealed under vacuum by direct current heating. The tunneling conductance dI/dV (V ) spectra were obtained by direct numerical derivation of the raw I(V ) spectra.

The large scale topographic STM images acquired on cleaved pristine crystals show atomic terraces which are several hundred of nanometers wide (Fig. 2a) [17]. The step heights are 0.5 nm, or its multiple, which corresponds to the thickness of slabs of Ta 4 Se 16 clusters and GaSe 4 tetrahedra (Fig. 1a). However, besides these steps, the topographic images of pristine GaTa 4 Se 8 crystal remains structureless as seen in figure 2c. The figures 2d-e present two conductance maps, taken at -600 mV and at 0 mV (Fermi level) respectively, which are extracted from the full dI/dV (V ) STS data set taken on the area of figure 2c. These maps reveal a spatially homogeneous insulating state. Typical conductance spectra, displayed in figure 2b, present a strong suppression of spectral weight around zero bias that is expected due to the gap at the Fermi level in the Mott insulating state. These data are consistent with the results of the optical conductivity measurements [24] that pointed out a gap of ∼ 100 meV. Note, that since the STS experiments were conducted at room temperature, the thermal broadening of 3.5k B T ≃ 90 meV affected the tunneling character- istics resulting in a smeared apparent gap and a finite conductance at zero bias, as shown in the inset of Fig. 2b.

If we now look at the STS data acquired on transited (metallic side) crystals, we find a very different picture (Fig. 3). First of all, the surface of transited GaTa 4 Se 8 crystals reveals specific topographic features at a scale of few tens of nanometers (see Figs. 3a and3f) which are absent in the pristine samples. We observed elongated blobs, with characteristic sizes 30-50 nanometers, organized in filaments and oriented along the direction of the electric field applied during the transition procedure via external electrodes. Second, the analysis of dI/dV (V ) conductance maps evidences that the MIT induced a complex electronic patchwork corresponding to the electronic phase separation at the nanoscale (Fig. 3bd). The first conductance map is taken at a voltage of -600 mV well above the gap (Fig. 3b), the second one is taken at -200 mV around gap edge (Fig. 3c) and the last one at the Fermi level (Fig. 3d). These conductance maps are strongly inhomogeneous. As underlined in figure 3c we encountered three typical zones: A (green), B (violet) and C (red) that correspond to three different electronic states which respective conductance spectra are displayed on Fig. 3e. Some conductance spectra corresponding to the zone A are featured in green. They are insulating-like and very similar to those observed in the pristine material, hence suggesting that these zones are barely affected by the electric pulse induced MIT. Spectra from the zone B, depicted in violet, have no conductance at zero bias (dI/dV (0)=0), and their corresponding gap of ∼ 400 meV is significantly larger than the gap in the pristine state. As they are more insulating than the pristine insulating medium we call these regions 'superinsulating'. Spectra from the zone C, represented in red, evokes those of a bad metal with a finite density of states at the Fermi level; hence these spectra will be qualified as 'metallic'. Interestingly, a contrast inversion appears between metallic and super-insulating patches when the maps acquired at voltages -200 mV (Fig. 3c) and -600 mV (Fig. 3b) are compared. This is due to the fact that the super-insulating patches are marked by a very low density of states below ±300 meV, while at -600 meV they represent the dominant contribution. On the contrary, the metallic patches, which density of states dominates at the Fermi level, have a lower density of states at -600 meV. It is also worth noting that both 'metallic' and 'super-insulating' spectra show an asymmetry between occupied and unoccupied states (i.e. between negative and positive bias voltage), whereas insulating spectra are almost symmetrical. Such an electronic patchwork is rep-resentative of the whole surface and has been encountered in numerous locations close or far from the macroscopic electrodes. This supports the conclusion drawn from four contact resistivity measurements [17] that the resistive switching in GaTa 4 Se 8 is globally a bulk effect: It is not related to any interfacial phenomenon to occur near the external electrodes used to transit the sample.

It is quite striking to observe that the topographic inhomogeneities (Fig. 3a) that appear through the pulse induce MIT transition, absent in the pristine material, are perfectly correlated to the induced metallic and superinsulating patches, as demonstrated on the 3D topography/spectroscopy map of Fig. 3f. In this image the 3D channel corresponds to the topography and the color channel corresponds to the conductance map at -200 mV. It is clear from this picture that the material inflates where new electronic patches appear, indicating that a very strong electron-lattice coupling is involved in the MIT. This strong electron-lattice coupling is reminiscent of the volume change encountered in the Mott metalinsulator transition such as in the canonical Mott insulator V 2 O 3 [25].

We have therefore strong indications that the electric field applied with macroscopic external electrodes (Fig. 1) creates an electronic phase separation concomitant to strong local sample deformations. At this step we may also expect that the strong electric field generated by the STM tip could locally modify the electronic state. Figure 4 shows successive topographic/spectroscopic maps measured on the same location of a transited crystal. The successive topography and conductance maps reveal that between two consecutive measurements, some zones switch back and forth between the super-insulating, insulating and metallic states (see the zooms on figure 4). Such switchings however, while being certainly due to an interaction between the STM tip and the sample, are not observed during purely topographic measurements, even after consecutive scans over 24 hours, if the bias voltage was kept lower than ∼ 700 mV. The fact that the switchings are observed only during spectroscopy measurements indicates that the electronic phase separation in GaTa 4 Se 8 is strongly affected by the voltage sweeps (±1V) applied trough the STM tip. This supports the existence of a threshold field, in agreement with the conclusions of macroscopic experiments [18,19]. Finally we note that the use of high voltage pulses (∼ 2 V) well above the threshold voltage to selectively induce new electronic phase appears to be impossible since the local inflation of the material during the pulse is so strong that the surface 'jumps' to the tip at every pulse and provokes tip crashes [26].

In conclusion, the comparison of our STM/STS data acquired on pristine and transited crystals gives a direct experimental evidence that the electric-field-induced metal-insulator transition recently observed in GaTa 4 Se 8 arises from the formation of a nanoscale electronic phase separation. The pristine crystals are homogeneous while the transited crystals display a peculiar electronic patchwork with a filamentary structure, that consists of a metallic/super-insulating network embedded in a barely pristine insulating matrix. The induced electronic inhomogeneities are associated with local topographic inflations indicating an important strain effect in the induced phase separation, probably originating from a strong electron-lattice coupling. The electronic phase separation can be reorganized by the electric field of the STM tip, suggesting that a local electric field might be used in this correlated systems to tune the Mott metal-insulator transition at the nanometer scale at room temperature.

FIG. 1 :

 1 FIG. 1: (Color online) (a)Representation of the crystal structure of GaTa4Se8 enhancing metal-metal bonding in the clusters (violet) Ta atoms : red, Se atoms : yellow, Ga-centered tetrahedra : green.(b) Temperature dependance of the resistance of the crystal used in this study. Red curve : before electric pulse, Blue curve : after electric pulse. (c) Photograph of the cleaved crystal used in this study. The arrows indicates the position of the external electrodes used to transit the crystal. (d) Schematic representation of the cleaving protocol.

FIG. 2 :

 2 FIG. 2: (Color online) (a) Topographic STM image (VT = -600 mV, IT = 100 pA) of the surface of a cleaved pristine crystal.(b) Representative dI/dV spectra acquired on the surface from panel (c-e). The image (c) is a small scale topographic image (VT = -600 mV, IT = 300 pA). The images (c,d) show the corresponding homogeneous conductance maps at -600 mV (d) and 0 mV (e).

FIG. 3 :

 3 FIG. 3: (Color online) (a) Topographic image of a cleaved transited GaTa4Se8 crystal (VT = -550 mV, IT = 250 pA). The figures b-d are the corresponding conductance map taken at -600 mV (b), at -200 mV and at the Fermi level 0 mV. These maps show clear electronic inhomogeneities. The figure (e) show the tunneling spectra corresponding to the zones A (green), B (violet), C (red) depicted on image c. The dI/dV spectra of the zone A (in green) are similar to the one of the insulating pristine samples, the spectra from zone B (violet) are more insulating and hence are called super-insulating while the spectra from zone C (in red) are 'metallic like'. In figure (f) the relief correspond to the topography while the color correspond to the conductance at -200 mV. This figures show clearly that the 'metallic' and super-insulating zones are topographically inflated.

FIG. 4 :

 4 FIG. 4: (Color online) The right part from a to f displays a series of consecutive Topographic/Spectroscopic measurements conducted on the same area (scale in nm). The topography was acquired with at VT = -630 mV, IT = 250 pA. The 3d canal correspond to topography while the color canal correspond to the conductance map at -200 mV. The left part of the figure, from g to l, show zooms on a same area, where switchings between insulating, metallic and super-insulating state are observed during consecutive frames. The topography is clearly correlated to the local electronic transitions.

The authors thank J. Martial at IMN for her skilled help in sample preparation. This work was supported by Grants ANR-05-JCJC-0123-01 (to L.C., B.C. and E.J.) and ANR-09-Blan-0154-01 from the French Agence Nationale de la Recherche.