Simulation of SPDE's for Excitable Media using Finite Elements
mmanuscript No. (will be inserted by the editor)
Simulation of SPDE’s for Excitable Media using FiniteElements
Muriel Boulakia · Alexandre Genadot · Michèle Thieullen
Received: date / Accepted: date
Abstract
In this paper, we address the question of the discretization of Stochastic Par-tial Differential Equations (SPDE’s) for excitable media. Working with SPDE’s driven bycolored noise, we consider a numerical scheme based on finite differences in time (Euler-Maruyama) and finite elements in space. Motivated by biological considerations, we study numerically the emergence of reentrant patterns in excitable systems such as the Barkley orMitchell-Schaeffer models.
Keywords
Stochastic partial differential equation · Finite element method · Excitablemedia
Mathematics Subject Classification (2000) · · The present paper is concerned with the numerical simulation of Stochastic Partial Differ-ential Equations (SPDE’s) used to model excitable cells in order to analyze the effect ofnoise on such biological systems. Our aim is twofold. The first is to propose an efficient and
This work was supported by the Agence Nationale de la Recherche through the project MANDy, Mathemat-ical Analysis of Neuronal Dynamics, ANR-09-BLAN-0008-01M. BoulakiaSorbonne Universités, UPMC Univ Paris 06, UMR 7598, Laboratoire Jacques-Louis Lions, F-75005, Paris,France.CNRS, UMR 7598, Laboratoire Jacques-Louis Lions, F-75005, Paris, France.INRIA-Paris-Rocquencourt, EPC REO, Domaine de Voluceau, BP105, 78153 Le Chesnay Cedex.E-mail: [email protected]. GenadotUniversité Paris-Dauphine, UMR 7534, Centre De Recherche en Mathématiques de la Décision, 75775, Paris,France.E-mail: [email protected]. ThieullenSorbonne Universités, UPMC Univ Paris 06, UMR 7599, Laboratoire de Probabilités et Modèles Aléatoires,F-75005, Paris, France.E-mail: [email protected] a r X i v : . [ m a t h . P R ] N ov Muriel Boulakia et al. easy-to-implement method to simulate this kind of models. We focus our work on practicalnumerical implementation with software used for deterministic PDE’s such as FreeFem++[20] or equivalent. The second is to analyze the effect of noise on these systems thanks tonumerical experiments. Namely, in models for cardiac cells, we investigate the possibilityof purely noise induced reentrant patterns such as spiral or scroll-waves as these phenomenaare related to major troubles of the cardiac rhythm such as tachyarrhythmia. For numericalexperiments, we focus on the Barkley and Mitchell-Schaeffer models, both originally deter-ministic models to which we add a noise source.Mathematical models for excitable systems may describe a wide range of biological phe-nomena. Among these phenomena, the most known and studied are certainly the two fol-lowing ones: the generation and propagation of the nerve impulse along a nerve fiber andthe generation and propagation of a cardiac pulse in cardiac cells. For both, following theseminal work [22], very detailed models known as conductance based models have beendeveloped, describing the physio-biological mechanism leading to the generation and prop-agation of an action potential. These physiological models are quite difficult to handlemathematically and phenomenological models have been proposed. These models describequalitatively the generation and propagation of an action potential in excitable systems. Forinstance, the Morris-Lecar model for the nerve impulse and the Fitzhugh-Nagumo model forthe cardiac potential. In the present paper, we will consider two phenomenological models:the Barkley and the Mitchell-Schaeffer models.
Mathematically, the models that we will consider consist in a degenerate system of PartialDifferential Equations (PDE’s) driven by a stochastic term, often referred to as noise. Moreprecisely, the model may be written (cid:26) d u = [ ν∆ u + ε f ( u , v )] d t + σ d W , d v = g ( u , v ) d t , (1)on [ , T ] × D , where D is a regular bounded open set of R or R . This system is completedwith boundary and initial conditions. W is a colored Gaussian noise source which will bedefined more precisely later. System (1) is degenerate in two ways: there is no spatial op-erator such as the Laplacian neither noise source in the equation on v . All the consideredmodels have the features of classical stochastic PDEs for excitable systems. The generalstructure of f and g is also typical of excitable dynamics. In particular, in the models thatwe will consider, the neutral curve f ( u , v ) = v held fixed is cubic in shape.To achieve our first aim, that is to numerically compute a solution of system (1), we workwith a numerical scheme based on finite difference discretization in time and finite elementmethod in space. The choice of finite element discretization in space has been directed bytwo considerations. The first is that this method fits naturally to a general spatial domain:we want to investigate the behavior of solutions to (1) on domains with various geometry.The second is that it allows to numerically implement the method with popular softwareused to simulate deterministic PDE’s such as the finite elements software FreeFem++ orequivalent. The discretization of SPDE by finite differences in time and finite elements inspace has been considered by several authors in theoretical studies, see for example [16,11,28,29,32,40]. Other methods of discretization are considered for example in [1,23,24,27,31,42]. These methods are based on finite difference discretization in time coupled ei-ther to finite difference in space or to the Galerkin spectral method, or to the finite elementmethod on the integral formulation of the evolution equation. We emphasize that we do notconsider in this paper a Galerkin spectral method or exponential integrator, that is, roughlyspeaking, we neither use the spectral decomposition of the solution of (1) according to a imulation of SPDE’s for Excitable Media using Finite Elements 3 Hilbert basis of L ( D ) (or another Hilbert space related to D ) nor the semigroup attachedto the linear operator (the Laplacian in (1)), in order to build our scheme. We only usethe variational formulation of the problem in order to fit to commonly used finite elementsmethod for deterministic PDE’s, see [2,13,18]. Moreover, the present paper is more ori-ented toward numerical applications than the above cited papers, in the spirit of [37]. In[37], the author numerically analyzes the effect of noise on excitable systems thanks to aGalerkin spectral method of discretization on the square. In the present paper, we pursuethe same objective using the finite element method instead of the Galerkin spectral one. Webelieve that the finite element method is easier to adapt to various spatial domains. Let usnotice that a discretization scheme for SPDE’s driven by white noise for spatial domains ofdimension greater or equal to 2 may lead to non trivial phenomena, see [19]. Consideringcolored noises may also be seen as a way to circumvent these difficulties.As is well known, one can consider two types of errors related to a numerical scheme forstochastic evolution equations: the strong error and the weak error. The strong error forthe discretization that we consider has been analyzed for one dimensional spatial domains(line segments) in [40]. The weak error for more general spatial domains, of dimension 2or 3 for example, has been considered in [16]. In the present paper, we recall results aboutthe strong error of convergence of the scheme because we want to numerically investigatepathwise properties of the model. Working with spatial domains of dimension d , as obtainedin [28,29], the strong order of convergence of the considered method for a class of linear stochastic equations is twice less than the weak order obtained in [16]. This is what is ex-pected since this same duality between weak and strong order holds for the discretization offinite dimensional stochastic differential equations (SDE’s).Our motivation for considering systems such as (1) comes from biological considerations. Inthe cardiac muscle, tachyarrhythmia is a disturbance of the heart rhythm in which the heartrate is abnormally increased. This is a major trouble of the cardiac rhythm since it may leadto rapid loss of consciousness and to death. As explained in [21,25], the vast majority oftachyarrhythmia are perpetuated by a reentrant mechanism. In several studies, it has beenobserved that deterministic excitable systems of type (1) are able to generate sustained reen-trant patterns such as spiral or meander, see for example [26,6]. We show numerically thatreentrant patterns may be generated and perpetuated only by the presence of noise. We per-form the simulations on the Barkley model whose deterministic version has been intensivelystudied in [6,4,5] and the model of Mitchell-Schaeffer which allows to get more realisticshape for the action potential in cardiac cells [34,9]. For Barkley model, similar experi-ments are presented in [37] where Galerkin spectral method is used as simulation schemeon a square domain. In our simulations, done on a square with periodic conditions or on asmoothed cardioid, we observe two kinds of reentrant patterns due to noise: the first may beseen as a scroll wave phenomenon whereas the second corresponds to spiral phenomenon.Both phenomena may be regarded as sources of tachyarrhythmia since in both cases, areasof the spatial domain are successively activated by the same wave which re-enters in theregion.All the simulations in the present paper have been performed using the FreeFem++ finiteelement software, see [20]. This software offers the advantage to provide the mesh of thedomain, the corresponding finite element basis and to solve linear problems related to thefinite element discretization of the model on its own. One of the originalities of the presentwork is to use this software to simulate stochastic PDE’s.Let us emphasize that the generic model (1) is endowed with a timescale parameter ε . Thepresence of this parameter is fundamental for the observation of traveling waves in the sys-tem: ε enforces the system to be either quiescent or excited with a sharp transition between Muriel Boulakia et al. the two states. Moreover, the values of the timescale parameter ε and the strength of thenoise σ appear to be of first importance to obtain reentrant patterns. This fact is also pointedout by our numerical bifurcation analysis. Let us mention that noise induced phenomenahave been studied in [7] for finite dimensional system of stochastic differential equations.The theoretical study of slow-fast SPDEs, through averaging methods, has been consideredin [10,12,41] for SPDEs.In a forthcoming work, we plan to address the effect of noise on deterministic periodic forc-ing of the Barkley and Mitchell-Schaeffer models. We expect to observe as in [38] for theone dimensional case, the annihilation by weak noise of the generation of some waves ini-tiated by deterministic periodic forcing. We also want to investigate stochastic resonancephenomena in such a situation. On a theoretical point of view, we intend to derive the strongorder of convergence of the discretization method used in the present paper for non-linearequations and systems of equations such as the FitzHugh-Nagumo, Barkley or Mitchell-Schaeffer models but also on simplified conductance based models.The remainder of the paper is organized as follows. In Section 2, we begin with the def-inition of the noise source in system (1) and present its finite element discretization. InSection 3, we introduce a discretization scheme based on finite element method in spacefor a stochastic heat equation and we recall the estimates from [28,29] for the strong orderof convergence. Then we apply the method to Fitzhugh-Nagumo model. In Section 4, weinvestigate the influence of noise on Barkley and Mitchell-Schaeffer models. We show that noise may initiate reentrant patterns which are not observable in the deterministic case. Wealso provide numerical bifurcation diagrams between the noise intensity σ and the time-scale ε of the models. At last, some technical proofs are postponed to the Appendix. Q -Wiener processes. Q -Wiener processesLet D be an open bounded domain of R d , d = L ( D ) the set of square integrable measurable functions withrespect to the Lebesgue measure on R d . Writing H for L ( D ) , we recall that H is a realseparable Hilbert space and we denote its usual scalar product by ( · , · ) and the associatednorm by (cid:107) · (cid:107) . They are respectively given by ∀ ( φ , φ ) ∈ H × H , ( φ , φ ) = (cid:90) D φ ( x ) φ ( x ) d x , (cid:107) φ (cid:107) = (cid:18) (cid:90) D φ ( x ) d x (cid:19) . Let Q be a non-negative symmetric trace-class operator on H . Let us recall the definitionof a Q -Wiener process on H which can be found in [35], Section 4.4, as well as the basicproperties of such a process. Definition 1
Let Q be a non-negative symmetric trace-class operator on H . There exists aprobability space ( Ω , F , P ) on which we can define a stochastic process ( W Qt , t ∈ R + ) on H such that – For each t ∈ R + , W Qt is a H -valued random variable. – W Q starts at 0 at time 0: W Q = H , P -a.s. – ( W Qt , t ∈ R + ) is a Lévy process, that is, it is a process with independent and stationaryincrements: imulation of SPDE’s for Excitable Media using Finite Elements 5 – Independent increments: for a sequence t , . . . , t n of strictly increasing times, therandom variables W Qt − W Qt , . . . , W Qt n − W Qt n − are independent. – Stationary increments: for two times s < t , the random variable W Qt − W Qs has samelaw as W Qt − s . – ( W Qt , t ∈ R + ) is a Gaussian process: for any t ∈ R + and any φ ∈ H , ( W Qt , φ ) is a realcentered Gaussian random variable with variance t ( Q φ , φ ) . – ( W Qt , t ∈ R + ) is a H -valued pathwise continuous process, P -almost surely.We recall the definition of non-negative symmetric linear operator on H admitting a kernel. Definition 2
A non-negative symmetric linear operator Q : H → H is a linear operator de-fined on H such that ∀ ( φ , φ ) ∈ H × H , ( Q φ , φ ) = ( Q φ , φ ) , ( Q φ , φ ) ≥ . Let q be a real valued integrable function on D × D such that ∀ ( x , y ) ∈ D × D , q ( x , y ) = q ( y , x ) , ∀ M ∈ N , ∀ x i , y j ∈ D , ∀ a i ∈ R , i , j = , . . . M , M ∑ i , j = q ( x i , y j ) a i a j ≥ , that is q is symmetric and non-negative definite on D × D . We say that Q has the kernel q if ∀ φ ∈ H , ∀ x ∈ D , Q φ ( x ) = (cid:90) D φ ( y ) q ( x , y ) d y . Let Q : H → H be a non-negative symmetric operator with kernel q . Then Q is a trace classoperator whose trace is given by Tr ( Q ) = (cid:90) D q ( x , x ) d x . For examples of kernels and basic properties of symmetric non-negative linear operators onHilbert spaces, we refer to [35], Section 4.9.2 and Appendix A. Let us now state clearly ourassumptions on the operator Q . Assumption 1
The operator Q is a non-negative symmetric operator with kernel q given by ∀ ( x , y ) ∈ D × D , q ( x , y ) = C ( x − y ) , where C belongs to C ( D ) and is an even function on D satisfying: ∀ M ∈ N , ∀ x i , y j ∈ D , ∀ a i ∈ R , i , j = , . . . M , M ∑ i , j = C ( x i − y j ) a i a j ≥ . Particularly, ∇ C ( ) = and x (cid:55)→ C ( x ) − C ( ) | x | is bounded on a neighborhood of zero. For x ∈ D and t ∈ R + , one can show, see [35], Section 4.4, that we can define W Qt at thepoint x such that the process ( W Qt ( x ) , ( t , x ) ∈ R + × D ) is a centered Gaussian process with covariance between the points ( t , x ) and ( s , y ) given by E (cid:16) W Qt ( x ) W Qs ( y ) (cid:17) = t ∧ s q ( x , y ) . In this case, the correlations in time are said to be white whereas the correlations in spaceare colored by the kernel q . Muriel Boulakia et al.
Proposition 1
Under Assumption 1, the process ( W Qt ( x ) , ( t , x ) ∈ R + × D ) has a versionwith continuous paths in space and time.Proof This is an easy application of the Kolmogorov-Chentsov test, see [15] Chapter 3,Section 3.2. Note that the regularity of C is important to get the result. (cid:117)(cid:116) Remark 1
The more the kernel q smooth is, the more the Wiener process regular is. Forexample, let ι ∈ R d and define f ( x ) = ι · Hess C ( x ) ι for x ∈ D . Suppose that f is a twicedifferentiable function. Then, one can show that there exists a probability space on which ( W Qt , t ∈ R + ) is twice differentiable in the direction ι . Remark 2
Let us assume that there exists a constant α and a (small) positive real δ suchthat ∀ y ∈ D , | C ( ) − C ( y ) | ≤ α | y | + δ . Using the Kolmogorov-Chentsov continuity theorem, one can show that the process ( W Qt ( x ) , ( t , x ) ∈ R + × D ) has a modification which is γ -Hölder in time for all γ ∈ (cid:0) , (cid:1) and γ -Hölder in space forall γ ∈ (cid:16) , + δ (cid:17) . Thus if δ >
0, by Rademacher theorem, for γ =
1, this version is almosteverywhere differentiable on D . This is another way to obtain regularity in space for W Q without using another probability space. Q satisfies Assumption 1. Let us present our approximation ofthe Q -Wiener process W Q . We begin with the discretization of the domain D . Let T h be afamily of triangulations of the domain D by triangles ( d =
2) or tetrahedra ( d = T h is given by h = max T ∈ T h h ( T ) , where h ( T ) = max x , y ∈ T | x − y | is the diameter of the element T . We assume that there exista positive constant ρ such that ∀ h > , ∀ T ∈ T h , ∃ x ∈ T , T ⊂ B ( x , ρ h ) , (2)where B ( x , r ) stands for the euclidean ball centered at x with radius r . We assume furtherthat this triangulation is regular as in Figure 1 (see [36] p. 108 for a definition) wherea triangulation is displayed and the property (2) is illustrated. In the present work, weconsider two kinds of finite elements: the Lagrangian P0 and P1 finite elements. However,the method could be adapted to other finite elements. The basis associated to the P0 finiteelement method is B , T h = { T , T ∈ T h } , where the function 1 T denotes the indicator function of the element T . Let { P i , ≤ i ≤ N h } be the set of all the nodes associated to the triangulation T h . The basis for the P1 finiteelement method is given by B , T h = { ψ i , ≤ i ≤ N h } , where ψ i is the continuous piecewise affine function on D defined by ψ i ( P j ) = δ i j (Kro-necker symbol) for all 1 ≤ i , j ≤ N h . imulation of SPDE’s for Excitable Media using Finite Elements 7 ρ h Fig. 1
Meshing (triangulation) of a domain and illustration of (2)
Definition 3
The P0 approximation of the noise W Q is given for t ∈ R + by W Q , h , t = ∑ T ∈ T h W Qt ( g T ) T , (3)where g T is the center of gravity of T . The P1 approximation is W Q , h , t = N h ∑ i = W Qt ( P i ) ψ i . (4) We will also consider the following alternative choice for the P0 discretization W Q , h , a t = ∑ T ∈ T h | T | ( W Qt , T ) T . (5) W Q , h , a corresponds to an orthonormal projection on P0. These approximations are againWiener processes as stated in the following proposition. Proposition 2
For i ∈ { , a , } the stochastic processes ( W Q , h , it , t ∈ R + ) are centered Q h , i -Wiener processes where, for φ ∈ HQ h , φ = ∑ T , S ∈ T h ( T , φ ) q ( g T , g S ) S , Q h , a φ = ∑ T , S ∈ T h ( T , φ ) ( Q T , S ) | T || S | S and Q h , φ = N h ∑ i , j = ( ψ i , φ ) q ( P i , P j ) ψ j . Proof
The fact that for i ∈ { , a , } the stochastic processes ( W Q , h , it , t ∈ R + ) are Wienerprocesses is a direct consequence of their definition as linear functionals of the Wiener pro-cess ( W Qt , t ∈ R + ) , see Definition 3. The corresponding covariance operators are obtainedby computing the quantity E (( W Q , h , i , φ )( W Q , h , i , φ )) for i ∈ { , a , } and φ , φ ∈ H (the details are left to the reader). (cid:117)(cid:116) Muriel Boulakia et al.
The P0 approximation (5) of the noise has been considered for white noise in dimension 2in [11]. White noise corresponds to Q = Id H . In the white noise case, the associated Wienerprocess is not at all regular in space (the trace of Q is infinite in this case). In the presentpaper, we work with trace class operators and thus with noises which are regular in space.Notice that discretization schemes for SPDE driven by white noise for spatial domains ofdimension greater or equal to 2 may lead to non trivial phenomena. In particular, usualschemes may not converge to the desired SPDE, see [19]. Theorem 1 (A global error)
For any τ ∈ R + and i ∈ { , a , } we have E (cid:32) sup t ∈ [ , τ ] (cid:107) W Qt − W Q , h , it (cid:107) (cid:33) ≤ K τ h where K may be written K = ˜ K (cid:107) Hess C (cid:107) ∞ for some constant ˜ K which only depends on | D | .Proof The proof is postponed to Appendix A. (cid:117)(cid:116)
Let us comment the above result. Let us take, as it will be the case in the numerical experi-ments, the following special form for the kernel ∀ x ∈ D , C ξ ( x ) = a ξ e − b ξ | x | for three positive real numbers a , b , ξ , where D is a bounded domain in dimension 2. Thisis a so-called Gaussian kernel. To a particular ξ -dependent kernel C ξ , we associate thecorresponding ξ -dependent covariance operator Q ξ . Then, according to Theorem 1, in thisparticular case, we see that E (cid:32) sup t ∈ [ , τ ] (cid:107) W Q ξ t − W Q ξ , h , it (cid:107) (cid:33) = O (cid:18) τ h ξ (cid:19) for any τ ∈ R + . Thus, when ξ goes to zero, this estimation becomes useless since the righthand-side goes to infinity. In fact, when ξ goes to zero, C ξ converges in the distributionalsense to a Dirac mass, and W Q ξ tends to a white noise which is, as mentioned before, anirregular process. In particular, the white noise does not belong to H and this is why ourestimation is no longer useful in this case. The same phenomenon occurs for any boundeddomain of dimension d ≥
2. Let us mention that for white noise acting on steady PDEsand on particular domains (square and disc), the error considered in Theorem 1 have beenstudied in [11]: the regularity of the colored noised improved these estimates in our case.We also remark that the proof of Theorem 1 does not rely on the regularity of the functionsof the finite element basis of P0 or P1 here. The key points are that C is smooth enough,even and that ∑ i φ i =
1, where { φ i } corresponds to the finite element basis.To conclude this section, we display some simulations. In Figure 2 we show simulations ofthe noise W Q ξ with covariance kernel defined by ∀ ( x , y ) ∈ D × D , q ξ ( x , y ) = C ξ ( x − y ) = ξ e − π ξ | x − y | , (6)where ξ >
0. We use the same kernel as in [37] for comparison purposes. As alreadymentioned, when ξ goes to zero, the considered colored noise tends to a white noise. Onthe contrary, when ξ increases, the correlation between two distinct areas increases as well. imulation of SPDE’s for Excitable Media using Finite Elements 9 Fig. 2
Simulations of W Q ξ with ξ = , . , , This property is illustrated in Figure 2. In other words, ξ is a parameter which allows tocontrol the spatial correlation.In these simulations, we have discretized W Q ξ , h , with the P1 discretization which reads W Q ξ , h , = N h ∑ i = W Q ξ ( P i ) ψ i . We remark that the family { W Q ξ ( P i ) , ≤ i ≤ N h } is a centered Gaussian vector with co-variance matrix ( q ξ ( P i , P j )) ≤ i , j ≤ N h . Using some basic linear algebra, it is not difficult tosimulate a realization of this vector and to project it on the P1 finite element basis to obtainFigure 2.We now propose a log-log graph to illustrate the estimate of Theorem 1. In the case where D is the square [ , ] × [ , ] , let us consider the kernel q given by: ∀ ( x , y ) ∈ D × D , q ( x , y ) = f k p ( x ) f k p ( y ) , where for two given integers k , p ≥ f k p ( x ) = ( k π x ) sin ( p π x ) (if x = ( x , x ) ∈ D ). Then, the covariance operator is given by Q φ = ( φ , f k p ) f k p for any φ ∈ L ( D ) .Moreover, one can show that ∀ t ≥ , W Qt = β t f k p for some real-valued Brownian motion β . All the calculations are straightforward in thissetting. Suppose, in the finite element setting, that the square D is covered by 2 N triangles, N ∈ N . For any N ∈ N , we denote by W Q , N , the P0 approximation of W Q given by (3). Weshow in Figure 3 the log-log graph of the (discrete) function: N (cid:55)→ µ N = E ( (cid:107) W Q − W Q , N , (cid:107) ) . According to Theorem 1, we should have µ N = O (cid:16) N (cid:17) . This result is recovered numeri-cally in Figure 3. Fig. 3
Log-log graph of N (cid:55)→ µ N for N = , , ,
30. In green is the comparison with a line of slope − In this section, we first present our numerical scheme. The considered space-time discretiza-tion is based on the Euler scheme in time and on finite elements in space. In this section andin the next section, we will use the following notations. Let us fix a time horizon T . For N ≥ ∆ t = TN and denote by ( u n , v n ) ≤ n ≤ N a sequence of approxima-tions of the solution of (1) at times t n = n ∆ t , 0 ≤ n ≤ N . The scheme, semi-discretized intime, is based on the following variational formulation, for n ∈ { , . . . , N − } : (cid:16) u n + − u n ∆ t , ψ (cid:17) + κ ( ∇ u n + , ∇ ψ ) = ε ( f n , ψ ) + σ∆ t ( W Qn + − W Qn , ψ ) , (cid:16) v n + − v n ∆ t , φ (cid:17) = ( g n , φ ) (7)for ψ and φ in appropriate spaces of test functions. Here, f n and g n correspond to approxi-mations of the reaction terms f and g in (1). The way we compute f n and g n is detailed inthe sequel for each considered model. W Qn is an appropriate approximation of W Qt n based onone of the discretizations proposed in Definition 3.In Subsection 3.1, we consider the strong error in the case of a linear stochastic partialdifferential equation driven by a colored noise to study the accuracy of the finite elementdiscretization. We recall the estimates on the strong order of convergence for such a numeri-cal scheme obtained in [28] and we numerically illustrate this result for an explicit example.In Subsection 3.2, we present in more details the scheme for the Fitzhugh-Nagumo modelwith a colored noise source since this model is one of the most used phenomenologicalmodels in cardiac electro-physiology, see the seminal work [17] and the review [30].3.1 Linear parabolic equation with additive colored noiseLet us consider the following linear parabolic stochastic equation on ( , T ) × D (cid:26) d u t = Au t d t + σ d W Qt , u = ζ . (8) Remember that H = L ( D ) is a separable Hilbert space with the scalar product and the corresponding norm respectively denoted by ( · , · ) and (cid:107) · (cid:107) . We assume that W Q is a Q -Wiener process with an operator Q which satisfies Assumption 1. We impose the followingcondition on the operator A in (8). Assumption 2
The operator − A is a positive self-adjoint linear operator on H whose do-main is dense and compactly embedded in H. imulation of SPDE’s for Excitable Media using Finite Elements 11
It is well known that the spectrum of − A is made up of an increasing sequence of positiveeigenvalues ( λ i ) i ≥ . The corresponding eigenvectors { w i , i ≥ } form a Hilbert basis of H .The domain of ( − A ) is the set (cid:40) u = ∑ i ≥ ( u , w i ) w i , ∑ i ≥ λ i ( u , w i ) < ∞ (cid:41) , that we denote here by V . It is continuously and densely embedded in H . The V -norm isgiven by | u | = (cid:112) − ( Au , u ) for all u ∈ V . We define a coercive continuous bilinear form a on V × V by a ( u , v ) = − ( Au , v ) . Let ζ be a V -valued random variable. The following proposition states that problem (8) iswell posed. Proposition 3
Equation (8) has a unique mild solution:u t = e At ζ + σ (cid:90) t e A ( t − s ) d W Qs , Moreover u is continuous in time and u t ∈ V for all t ∈ [ , T ] , P -a.s. Proof
This result is a direct consequence of Theorem 5.4 of [15], Assumptions 1 and 2. (cid:117)(cid:116)
For h >
0, let V h be a finite dimensional subset of V with the property that for all v ∈ V , thereexists a sequence of elements v h ∈ V h such that lim h → (cid:107) v − v h (cid:107) =
0. For an element u of V , we introduce its orthogonal projection on V h and denote it Π h u . It is defined in a uniqueway by Π h u ∈ V h and ∀ v h ∈ V h , a ( Π h u − u , v h ) = . (9)Let I h be the dimension of V h . Notice that there exists a basis ( w i , h ) ≤ i ≤ I h of V h orthonormalin H with the following property: for each 1 ≤ i ≤ I h , there exists λ i , h such that ∀ v h ∈ V h , a ( v h , w i , h ) = λ i , h ( v h , w i , h ) , (see [36], Section 6.4). The family ( λ i , h ) ≤ i ≤ I h is an approximating sequence of the familyof eigenvalues ( λ i ) i ≥ so that λ i , h ≥ λ i , ∀ ≤ i ≤ I h . We study the following numerical scheme to approximate equation (8) defined recursivelyas follows. For u given in V h , find ( u hn ) ≤ n ≤ N in V h such that for all n ≤ N − (cid:26) ∆ t ( u hn + − u hn , v h ) + a ( u hn + , v h ) = σ∆ t ( W Q , hn + − W Q , hn , v h ) u h = u (10)for all v h ∈ V h where W Q , hn is an appropriate approximation of W Qn ∆ t in V h . The approximationerror of the scheme can be written as the sum of two errors. Definition 4
The discrete error introduced by the scheme (10) is defined by E hn = e hn + p hn where e hn = u hn − Π h u t n , p hn = Π h u t n − u t n (11)for 0 ≤ n ≤ N . For n ∈ { , . . . , N } , the error e hn is the difference between the approximated solution givenby the scheme and the elliptic projection on V h of the exact solution at time n ∆ t . The error p hn is the difference between the exact solution and its projection on V h at time n ∆ t .In order to give explicit bounds for the error defined above, let us choose our approximationin H specifically. This imposes to choose the space V explicitly. Assume that the operator A is such that V = H ( D ) and that V h is a space of P1 finite elements, see Section 2.2. Inthis P1 case, for n ∈ { , . . . , N } , we set W Q , hn = W Q , h , n ∆ t defined by Definition 3. The situationconsidered in the present section is also the one studied in example 3.4 and section 7 of [28].We have the following bound for the numerical error coming from such a numerical scheme. Theorem 2 (Example 3.4 and Corollary 7.2 in [28])
Let us assume that Assumptions 1and 2 are satisfied. Moreover, assume that we are in the P1 case: for n ∈ { , . . . , N } , W Q , hn = W Q , h , n ∆ t defined by Definition 3. Then, there exists ∆ t > such that for all n ∈ { , . . . , N } and ∆ t ∈ [ , ∆ t ] (cid:113) E ( (cid:107) E hn (cid:107) ) ≤ K ( h + √ ∆ t ) , (12) where K is a constant depending only on T and | D | . Let us recall the weak order of convergence of the considered scheme obtained in [16] but under weaker assumptions. Since C is a twice differentiable even function on D , ∆ C is abounded function on D and therefore, according to [16] Theorem 3.1, for any bounded realvalued twice differentiable function φ on L ( D ) , there exists a constant K depending onlyon T such that | E ( φ ( u hN )) − E ( φ ( u T )) | ≤ K ( h γ + ∆ t γ ) (13)for a given γ <
1. In our situation, it is more natural to consider the strong error since westudy pathwise behavior. For the method that we consider, estimates for the strong errorhave been obtained for one dimensional spatial domains and white noise in [40]. Many pa-pers exist for finite dimensional systems. The estimate of Theorem 2 lies in between thesetwo types of studies. The noise is colored but the spatial domain may be of any dimension.Notice that the order of weak convergence (13) is twice the order of strong convergence(12), as for finite dimensional stochastic differential equations.In the end of this subsection, we illustrate this error estimate in a simple situation. Weconsider the domain D = ( , l ) × ( , l ) for l >
0. We set A = ∆ and D ( A ) = H ( D ) ∩ H ( D ) ,thus V = H ( D ) . That is we consider the equation: d u t = ∆ u t d t + σ d W Qt , in D , u t = , on ∂ D , u = , in D (14)for t ∈ R + . W Q is a L ( D ) -valued Q -Wiener process defined by Definition 1 and Q satis-fies Assumption 1. The initial condition is zero, hence the solution of the correspondingdeterministic equation, without noise, is simply zero for all time. Following Proposition 3,equation (14) has a unique mild solution such that u t ∈ H ( D ) for all t ∈ [ , T ] , P -almostsurely. Moreover, u has a version with time continuous paths and such that, for any time T >
0: sup t ∈ [ , T ] E ( (cid:107) u t (cid:107) H ( D ) ) < ∞ . imulation of SPDE’s for Excitable Media using Finite Elements 13 We denote by ( e ∆ t , t ≥ ) the contraction semigroup associated to the operator ∆ . The mildsolution to equation (14) is defined as the following stochastic convolution u t = σ (cid:90) t e ∆ ( t − s ) d W Qs for t ∈ R + , P -almost-surely. In order to compute the expectation of the squared norm of u in L ( D ) analytically and also as precisely as possible numerically, we define the Hilbert basis ( e kp , k , p ≥ ) of L ( D ) which diagonalizes the operator ∆ defined on D ( A ) . For k , p ≥ ( x , y ) ∈ D e kp ( x , y ) = l sin (cid:18) k π l x (cid:19) sin (cid:16) p π l y (cid:17) . A direct computation shows that ∆ e kp = − λ kp e kp where λ kp = π l ( k + p ) . In the basis ( e kp , k , p ≥ ) of L ( D ) , the semigroup ( e ∆ t , t ≥ ) is given by e ∆ t φ = ∑ k , p ≥ e − λ kp t ( φ , e kp ) e kp for t ∈ R + and φ ∈ L ( D ) . Then for any t ∈ R + (c.f. Proposition 2.2.2 of [14]) E ( (cid:107) u t (cid:107) ) = σ (cid:90) t Tr (cid:16) e ∆ s Q (cid:17) d s = σ ∑ k , p ≥ − e − λ kp t λ kp ( Qe kp , e kp ) . In the sequel, we write Γ t = E ( (cid:107) u t (cid:107) ) . The above series expansion can then be implementedand we can compare this result with E ( (cid:107) u hn (cid:107) ) which is computed thanks to Monte-Carlosimulations. The Monte-Carlo simulation of E ( (cid:107) u hn (cid:107) ) consists in considering ( u h , pn ) ≤ p ≤ P , P ∈ N a sequence of independent realizations of the scheme (10) and define Γ ( P ) n ∆ t = P P ∑ p = (cid:107) u h , pn (cid:107) , (15)the approximation of Γ at time n ∆ t , n ∈ { , . . . , N } . We denote also by Γ ( P ) the continuouspiecewise linear version of Γ . Figure 4 displays numerical simulations of the processes ( Γ t , t ∈ R + ) and ( Γ ( P ) t , t ∈ R + ) . The simulations are done with l =
80. Moreover thedomain is triangulated with 5000 triangles giving a space step of about h = .
64 and anumber of vertices’s of about 2600. For this simulation, we choose P =
40 which is not bigbut Γ ( ) matches quite well with its corresponding theoretical version Γ , as expected by thelaw of large numbers. We remark also that for the same spatial discretization of the domain D , there is no particular statistical improvement to choose the P1 finite element basis insteadof the P0.3.2 Space-time discretization of the Fitzhugh-Nagumo modelWe write the scheme for the Fitzhugh-Nagumo which is a widely used model of excitablecells, see [17,30]. The stochastic Fitzhugh-Nagumo model, abbreviated by FHN model inthe sequel, consists in the following 2-dimensional system (cid:26) d u = [ κ∆ u + ε ( u ( − u )( u − a ) − v )] d t + σ d W Q , d v = [ u − v ] d t , (16) Fig. 4
Simulations of ( Γ t , t ∈ [ , ]) in green and its approximation Γ ( ) in red computed with: the P0 (onthe left) and P1 (on the right) approximations of the noise. For the simulation we choose the coefficient ofcorrelation ξ = q ξ defined by (6). The intensity of the noise is σ = .
15. The time step is0 .
05 whereas the space step is about 0 .
64. The two black curves are respectively Γ plus, respectively minus,the error introduced by the scheme which is expected to be of order √ ∆ t + h equals here to √ . + . on [ , T ] × D . In the above system, κ > diffusion coefficient, ε > time-scale coefficient, σ > a ∈ ( , ) a parameter. W Q is a Q -Wienerprocess satisfying Assumption 1. System (16) must be endowed with initial and boundaryconditions. We denote by u and v the initial conditions for u and v . Moreover we assume that u satisfies zero Neumann boundary conditions: ∀ t ∈ [ , T ] , ∂ u t ∂ n = , on ∂ D , (17)where ∂ D denotes the boundary of D and n is the external unit normal to this boundary.Noisy FHN model and especially, FHN with white noise, have been extensively studied. Werefer the reader to [8] where all the arguments needed to prove the following proposition aredeveloped. Proposition 4
Let W Q be a colored noise with Q satisfying Assumption 1. We assume thatu and v are in L ( D ) , P -almost surely. Then, for any time horizon T , the system (16)has a unique solution ( u , v ) defined on [ , T ] which is P -almost surely in C ([ , T ] , H ) × C ([ , T ] , H ) . The proof of this proposition relies on Itô Formula, see Chapter 1, Section 4.5 of [15], andthe fact that the functional defined by f ( x ) = x ( − x )( x − a ) , ∀ x ∈ R satisfies the inequality ( f ( u ) − f ( v ) , u − v ) ≤ + a − a (cid:107) u − v (cid:107) , ∀ ( u , v ) ∈ H × H , which implies that the map f − + a − a Id is dissipative. The local kinetics of system (16),that is the dynamics in the absence of spatial derivative, is illustrated in Figure 5. It describesthe dynamics of the system of ODEs (cid:26) d u = [ ε u ( − u )( u − a ) − v ] d t , d v = [ u − v ] d t , (18) imulation of SPDE’s for Excitable Media using Finite Elements 15 u ( t ) v ( t ) • • • Fig. 5
Phase portrait with nullclines of system (18) for a = . ε = .
1. The blue points correspond tothe three equilibrium points of the system. when the initial condition ( u , v ) is in [ , ] × [ , ] .We explicitly give the numerical scheme used to simulate system (16). Let us define thefunction k given by k ( x ) = ε (cid:0) − x + x ( + a ) (cid:1) , ∀ x ∈ R . This function corresponds to the non linear parts of the reaction term f . We use the followingsemi-implicit Euler-Maruyama scheme (cid:40) u n + − u n ∆ t = κ∆ u n + − a ε u n + + k ( u n ) − v n + + σ √ ∆ t W Q , n + , v n + − v n ∆ t = u n + − v n + , (19)where ( W Q , n ) ≤ n ≤ N + is a sequence of independent Q -Wiener processes evaluated at time1. Let ( H ∗ , B ( H ∗ ) , ˜ P ) be chosen so that the canonical process has the same law as W Q , n + under ˜ P . Then, for a given ( u n , v n ) ∈ H ( D ) × H , the equation ( ∆ t + a ε + ∆ t + ∆ t ) u n + − κ∆ u n + = k ( u n ) − + ∆ t v n + σ √ ∆ t W Q , n + has a unique weak solution u n + in H ( D ) , ˜ P -almost surely. This fact follows from Lax-Milgram Theorem and a measurable selection theorem, see Section 5 of the survey [39].Therefore, without loss of generality, we may assume in this section that the probabilityspace is ( C ([ , T ] , H ∗ ) , B ( C ([ , T ] , H ∗ )) , P ) such that under P , the canonical process hasthe same law as W Q . Remark 3
In the scheme (19), we could have chosen other ways to approximate the reac-tion term. For instance, it is also possible to work with a completely implicit scheme with k ( u n + ) instead of k ( u n ) in (19).Let us consider the weak form for the first equation of (19). We get, (cid:40) ( ∆ t + a ε + ∆ t + ∆ t )( u n + , ψ ) + κ ( ∇ u n + , ∇ ψ ) = ( k ( u n ) , ψ ) − + ∆ t ( v n , ψ ) + σ √ ∆ t ( W Q , n + , ψ ) , v n + − ∆ t + ∆ t u n + = + ∆ t v n (20) Fig. 6
Simulations of system (16) with ξ = σ = ε = . a = .
1. These figures must be read from theup-left to the down-right. The time step is 0 . ms and there is 0 . ms between each figure. for all ψ ∈ H ( D ) . Let h > ( ψ i , ≤ i ≤ N h ) be the P1 finite element basis defined inSection 2. For n ≥
0, we define the vectors u n = ( u n , i ) ≤ i ≤ N h , v n = ( v n , i ) ≤ i ≤ N h , W Qn + = ( W Q , n + ( P i )) ≤ i ≤ N h , which are respectively the coordinates of u n , v n and W Q , n + w.r.t. the basis ( ψ i , ≤ i ≤ N h ) .We also define the stiffness matrix A ∈ M N h ( R ) and the mass matrix M ∈ M N h ( R ) by A i j = ( ∇ ψ i , ∇ ψ j ) , M i j = ( ψ i , ψ j ) . System (20) can be rewritten as (cid:18) ( ∆ t + a ε + ∆ t + ∆ t ) M + κ A − ∆ t + ∆ t I I (cid:19) (cid:18) u n + v n + (cid:19) = (cid:18) − + ∆ t M + ∆ t I (cid:19) (cid:18) u n v n (cid:19) + (cid:18) G ( u n ) (cid:19) + (cid:18) σ √ ∆ t M
00 0 (cid:19) (cid:18) W Q , n + (cid:19) , where G ( u n ) = ( k ( u n ) , ψ i ) ≤ i ≤ N h ∈ R N h . As for the parabolic stochastic equation consideredin Section 3.1, one may expect a numerical strong error for this scheme of order E ( (cid:107) ( u t , v t ) − ( u t , n , v t , n ) (cid:107) ) = O ( h + √ ∆ t ) , (21)for ∆ t ≤ ∆ t . In (21) , ( u t , n , v t , n ) t ∈ [ , T ] is the interpolation of the discretized point which ispiecewise linear in time.We end up this section with Figure 6 which displays simulations of the stochastic Fitzhugh-Nagumo model (16) with zero Neumann boundary conditions on a cardioid domain and zeroinitial conditions. The kernel of the operator Q is given by equation (6) for some ξ >
0. Dueto a strong intensity of the noise source ( σ = imulation of SPDE’s for Excitable Media using Finite Elements 17 In this section, we focus on classical models for excitable cells, namely Barkley and Mitchell-Schaeffer models. We would like to observe cardiac arrhythmia, that is troubles that mayappear in the cardiac beats. Among the diversity of arrhythmia, the phenomena of tachycar-dia are certainly the most dangerous as they lead to rapid loss of consciousness and death.Tachycardia is described as follows in [25].The vast majority of tachyarrhythmias are perpetuated by reentrant mechanisms.Reentry occurs when previously activated tissue is repeatedly activated by the prop-agating action potential wave as it reenters the same anatomical region and reacti-vates it.In system (1), the equation on u gives the evolution of the cardiac action potential. Theequation on v takes into account the evolution of internal biological mechanisms leadingto the generation of this action potential. We will be more specifically interested by twosystems of this form: the Barkley and Mitchell-Schaeffer models.4.1 Numerical study of the Barkley model In the deterministic setting, a paradigm for excitable systems where reentrant phenomenasuch as spiral, meander or scroll waves have been observed and studied is the Barkley model,see [6,3,4,5]. This deterministic model is of the following form (cid:26) d u = [ κ∆ u + ε u ( − u )( u − v + ba )] d t , d v = [ u − v ] d t . (22)The parameter ε is typically small so that the time scale of u is much faster than that of v . Formore details on the dynamic of waves in excitable media, we refer the reader to [26]. TheBarkley model, like two-variables models of this type, faithfully captures the behavior ofmany excitable systems. The deterministic model (22) does not exhibit re-entrant patternsunless one imposes special conditions on the domain: for instance, one may impose thata portion of the spatial domain is a "dead zone". This means a region with impermeableboundaries where equations (22) do not apply: when a wave reaches this dead region, thetip of the wave may turn around and this induces a spiral behavior, see Section 2.2 of [26].One may also impose specific initial conditions such that some zones are intentionally hyper-polarized: the dead region is somehow transient in this case. As in [37] we add a colored noise with kernel of type (6) to equation (22) and so we consider (cid:26) d u = [ κ∆ u + ε u ( − u )( u − v + ba )] d t + σ d W Q ξ , d v = [ u − v ] d t , (23)where the kernel of Q ξ is given by (6) for ξ > D = [ , l ] × [ , l ] with periodic Fig. 7
Reentry is observed for system (23) with ξ = σ = . ε = . a = . b = .
01 and ν = ms , there is 0 . ms between each figures for a timestep of 0 . ms . boundary conditions: ∀ t ∈ R + , ∀ x ∈ [ , l ] u t ( x , ) = u t ( x , l ) , and ∂ u t ∂ n ( x , ) = ∂ u t ∂ n ( x , l ) , ∀ y ∈ [ , l ] u t ( , y ) = u t ( l , y ) , and ∂ u t ∂ n ( , y ) = ∂ u t ∂ n ( l , y ) , (24)where n is the external unit normal to the boundary. The numerical scheme is based on thefollowing variational formulation. Given u and v in H ( D ) , find ( u n , v n ) ≤ n ≤ N such thatfor all 0 ≤ n ≤ N − (cid:40) ( u n + − u n ∆ t , ψ ) + κ ( ∇ u n + , ∇ ψ ) = ε ( u n ( − u n )( u n − v n + ba ) , ψ ) + σ √ ∆ t ( W Q , n + , ψ ) , v n + − v n ∆ t = u n + − v n + (25)with boundary conditions u n ( x , ) = u n ( x , l ) , u n ( , x ) = u n ( l , y ) and for all ψ ∈ H ( D ) sat-isfying ψ ( x , ) = ψ ( x , l ) and ψ ( , y ) = ψ ( l , y ) for any ( x , y ) ∈ [ , l ] × [ , l ] . We have solvedthis problem using the P1 finite element methods, see Section 2.2.Our aim is to observe reentrant patterns generated by the presence of the noise source inthis system. Figure 7 displays simulations of (23) using the P1 finite element method. Weobserve the spontaneous generation of waves with a reentrant pattern. At some points inthe spatial domain, the system is excited and exhibits a reentrant evolution which is self-sustained: a previously activated zone is re-activated by the same wave periodically. Asexplained in [25] and quoted in Section 4.1.1, this phenomenon can be interpreted biolog-ically as tachycardia in the heart tissue. We observe that, as in [37], the constants a and b are chosen such that the deterministic version of system (23) may exhibit spiral pattern,see the bifurcation diagram between a and b in [5]. However, in our context, the generationof spiral is a phenomenon which is due solely to the presence of noise. In particular, thereis no need for a "dead region", as previously mentioned for the observation of spirals orreentrant patterns in a deterministic context. Figure 8 displays a simulation of system (23)on a cardioid domain with zero Neumann boundary conditions, see (17). We observe the imulation of SPDE’s for Excitable Media using Finite Elements 19a) b) c) d)e) f) g) h)i) j) k) l) Fig. 8
Simulations of system (23) with ξ = σ = . ε = . a = . b = .
01 and ν =
1. As forthe previous figure, the quiescent state is represented in green whereas the excited state is in violet. Anotherphenomena of re-entry is observed on this cardioid geometry with zero Neumann boundary conditions. Thereis 2 ms between each snapshot for a time step for the simulations equals to 0 . ms . spontaneous generation of a wave turning around itself like a spiral and thus reactivatingzones already activated by the same wave. To gain a better insight into these reentrant phenomena, a bifurcation diagram between ε and σ in system (23) is displayed in Figure 9. In this figure, the other parameters a , b , ν , ξ are held fixed. The domain and boundary conditions are the same as for Figure 7. Threedistinct areas emerge from repeated simulations: – the area NW (for No Wave) where no wave is observed. – the area W (for Wave) where at least one wave is generated on average. Such waves donot exhibit reentrant patterns. – the area RW (for Reentrant Wave) where waves with re-entry are observed. The wavehas the same pattern as in Figure 7.Let us mention that these three different area emerge from repeated simulations. The de-tection of re-entrant patterns is quite empirical here: we say that there is re-entrant patternsif we can actually see it on the figures. An automatic detection of re-entrant patterns maycertainly be derived from [33] even if rather difficult to implement in our setting.At transition between the areas W and RW, ring waves with the same pattern as reentrantwaves may be observed: two arms which join each other to form a ring. We also remark thatfor a fixed ε , when σ increases, the number of nucleated waves increases. On the contrary,for a fixed σ when ε increases the number of nucleated waves decreases. Let us notice thatfor small ε , that is when the transition between the quiescent and excited state is very sharp,small noise may powerfully initiate spike. However, we only observe reentrant patternswhen ε is large enough. Notice also that the separation curve between the zone NW and thetwo zones W and RW is exponentially shaped. This may be related to the large deviationtheory for slow-fast system of SPDE. .
01 0 . . . σ ε NWW RW
Fig. 9
Numerical bifurcation diagram between ε and σ of system (23) with ξ = a = . b = .
01 and ν = ε and σ -step are respectively 0 .
005 and 0 .
025 but has been refined around boundariesto draw the curve boundaries.
Fitzhugh-Nagumo model is the most popular phenomenological model for cardiac cells.However this model has some flaws, in particular the hyperpolarization and the stiff slopein the repolarization phase. The Mitchell-Schaeffer model [34] has been proposed to im-prove the shape of the action potential in cardiac cells. The spatial version of the Mitchell-Schaeffer model reads as follows d u = (cid:20) κ∆ u + v τ in u ( − u ) − u τ out (cid:21) d t + σ d W Q ξ , d v = (cid:20) τ open ( − v ) u < u gate − v τ close u ≥ u gate (cid:21) d t . (26)The numerical scheme is based on the following variational formulation. Given u and v in H ( D ) , find ( u n , v n ) such that for all 0 ≤ n ≤ N − ( u n + − u n ∆ t , ψ ) + κ ( ∇ u n + , ∇ ψ ) = ε ( v n τ in u n ( − u n ) − τ out u n , ψ ) + σ √ ∆ t ( W Qn + ( ) , ψ ) , v n + − v n ∆ t = τ open ( − v n ) u n < u gate − v n τ close u n ≥ u gate (27)for ψ ∈ H ( D ) . More precisely, we solve this problem with the P1 finite element method. Bifurcations have been investigated in Figure 10 for the same domain and boundary condi-tions as for the bifurcation diagram related to Barkley model (Figure 9). We choose to fixall the parameters except the intensity of the noise σ and τ close to investigate the influence ofthe strength of the noise and the characteristic time for the recovery variable v to get closed.From repeated simulations, five distinct areas emerge: imulation of SPDE’s for Excitable Media using Finite Elements 213 . . . . . σ τ close NWDW RW WT
Fig. 10
Numerical bifurcation diagram between τ close and σ of system (26) with ξ = τ in = . τ out = . τ open = . , u gate = .
13 and ν = .
03 held fixed. – the area NW (for No Wave) where no wave is observed. – the area W (for Wave) where at least one wave is generated on average. These waves do not exhibit reentrant patterns. However, these waves may be generated with the samepattern as reentrant waves: two arms which meet up and agree to form a ring. – the area RW (for Reentrant Wave) where waves with re-entry may be observed as inFigure 7. – the area DW (for Disorganized Wave) where reentrant waves are initiated but breakdown in numerous pieces resulting in a very disorganized evolution. In a sense, thisdisorganized evolution may be regarded as reentrant since previously activated zonemay be re-activated by one of these resulting pieces. – the area T (for Transition) is a transition area between reentrant waves and more disor-ganized patterns as observed in the area DW.We would like to end the paper with a short discussion about a problem we intend toaddress in future works. In [38], the authors present simulations of a stochastic spatially-extended Hodgkin-Huxley model. This is a celebrated model of propagation and generationof an action potential in a nerve fiber. They consider the case of a one dimensional nervefiber stimulated by a noisy signal. They show, using numerical experiments, that the pres-ence of weak noise in the model may powerfully annihilate the generation (but not the prop-agation) of waves. We reproduce these numerical experiments in our setting for the Barkleyand Fitzhugh Nagumo models in order to investigate if this phenomena may be observed intwo dimensions. Thus, we consider the Barkley model with a periodic deterministic inputand driven by a colored noise. We perform simulations of this model using growing noiseintensity. Unfortunately, we were not able to produce inhibition of a periodic deterministicsignal using weak noises. Of course, with stronger intensity of the noise, signals due tothe stochasticity of the system are generated and perturb the deterministic periodic signal.However, this phenomenon is not surprising. It seems to us that FitzHugh-Nagumo couldbe a better model to consider on the question of annihilate the generation of waves by weaknoise. We noted that a lot of care is required when choosing the mesh to use as well asthe parameters (intensity of the noise, intensity of the deterministic signal, duration of itsperiod), if we want to produce sound results. Acknowledgements
The authors would like to thank the anonymous reviewer for his/her valuable commentsand suggestions to improve the quality of the paper.
A Proof of Theorem 1
Recall that the domain D is polyhedral such that D = (cid:91) T ∈ T h T . Let i ∈ { , a , } . The process ( D h ( t ) , t ∈ [ , τ ]) defined by D h ( t ) = W Qt − W Q , h , it is a centered Wiener process. In particular, it is a continuous martingale and thus, by the Burkholder-Davis-Gundy inequality (see Theorem 3.4.9 of [35]) we have E (cid:32) sup t ∈ [ , τ ] (cid:107) D h ( t ) (cid:107) (cid:33) ≤ c E ( (cid:107) D h ( τ ) (cid:107) ) with c a constant which does not depend on h or τ . We begin with the case i =
1. Since the processes W Q and W Q , h , are regular in space, we write E ( (cid:107) D h ( τ ) (cid:107) ) = E (cid:18) (cid:90) D ( W Q τ ( x ) − W Q , h , τ ( x )) d x (cid:19) . We use the definition of W Q , h , in Definition 3 and the fact that ∑ N h i = ψ i = E ( (cid:107) D h ( τ ) (cid:107) ) = E (cid:32) (cid:90) D ( W Q τ ( x ) − N h ∑ i = W Q τ ( P i ) ψ i ( x )) d x (cid:33) = E (cid:32) (cid:90) D ( N h ∑ i = ( W Q τ ( x ) − W Q τ ( P i )) ψ i ( x )) d x (cid:33) = E (cid:32) (cid:90) D N h ∑ i , j = ( W Q τ ( x ) − W Q τ ( P i ))( W Q τ ( x ) − W Q τ ( P j )) ψ i ( x ) ψ j ( x ) d x (cid:33) . By an application of Fubini’s theorem, exchanging over the expectation, integral and summation, we get E ( (cid:107) D h ( τ ) (cid:107) ) = N h ∑ i , j = (cid:90) D E (cid:16) ( W Q τ ( x ) − W Q τ ( P i ))( W Q τ ( x ) − W Q τ ( P j )) (cid:17) ψ i ( x ) ψ j ( x ) d x = τ N h ∑ i , j = (cid:90) D (cid:0) C ( ) − C ( P i − x ) − C ( P j − x ) + C ( P i − P j ) (cid:1) ψ i ( x ) ψ j ( x ) d x . For all 1 ≤ i , j ≤ N h , if the intersection of the supports of ψ i and ψ j is not empty, then ∀ x ∈ supp ψ i , ∀ y ∈ supp ψ j , | x − y | ≤ Kh . Thus, there exists K > i , j , if supp ψ i ∩ supp ψ j (cid:54) = /0 and x ∈ supp ψ i ∩ supp ψ j , a Taylor’sexpansion yields | C ( ) − C ( P i − x ) − C ( P j − x ) + C ( P i − P j ) | ≤ K max x ∈ D (cid:107) Hess C ( x ) (cid:107) h , where we have used the fact that ∇ C ( ) =
0. Then, E ( (cid:107) D h ( τ ) (cid:107) ) ≤ K τ max x ∈ D (cid:107) Hess C ( x ) (cid:107) h . imulation of SPDE’s for Excitable Media using Finite Elements 23This ends the proof for the case i =
1. The case i = i = a , we proceed as follows. The process W Q , h , a is the orthonormal projection of W Q on thespace P0, thus, we have, using the Pythagorean theorem, E ( (cid:107) D h ( τ ) (cid:107) ) = E (cid:16) (cid:107) W Q τ − W Q , h , a τ (cid:107) (cid:17) = E (cid:16) (cid:107) W Q τ (cid:107) − (cid:107) W Q , h , a τ (cid:107) (cid:17) . Then, recalling that the processes W Q and W Q , h , a are regular in space and using the fact that the triangles T ∈ T h do not intersect, we obtain, E ( (cid:107) D h ( τ ) (cid:107) ) = E (cid:32) (cid:90) D W Q τ ( x ) d x − ∑ T ∈ T h | T | ( W Q τ , T ) (cid:33) . By an application of Fubini’s theorem, exchanging over the expectation and summation, we get E ( (cid:107) D h ( τ ) (cid:107) ) = τ (cid:32) C ( ) | D | − ∑ T ∈ T h | T | ( Q T , T ) (cid:33) . (28)Since D = (cid:83) T ∈ T h T we have C ( ) | D | = ∑ T ∈ T h | T | (cid:90) T (cid:90) T C ( ) d z d z , hence, plugging in (28) E ( (cid:107) D h ( τ ) (cid:107) ) = τ ∑ T ∈ T h | T | (cid:90) T (cid:90) T [ C ( ) − C ( z − z )] d z d z . (29)Thanks to the fact that ∇ C =
0, a Taylor’s expansion yields C ( ) − C ( z − z ) = ( z − z ) · Hess C ( )( z − z ) + o ( | z − z | ) . (30)Thus, thanks to (2), for all z , z in the same triangle T | C ( ) − C ( z − z ) | ≤ K max x ∈ D (cid:107) Hess C ( x ) (cid:107) h , where the constant K is independent from T ∈ T h . Plugging in (29) yields E ( (cid:107) D h ( τ ) (cid:107) ) ≤ K max x ∈ D (cid:107) Hess C ( x ) (cid:107) τ h for a deterministic constant K . References
1. Allen, E., Novosel, S., Zhang, Z.: Finite element and difference approximation of some linear stochasticpartial differential equations. Stochastics: Int. J. Probab. and Stoch. Proc. (1-2), 117–142 (1998)2. Babuska, I., Szabo, B., Katz, I.: The p-version of the finite element method. SIAM J. Numer. Anal. (3), 515–545 (1981)3. Barkley, D.: A model for fast computer simulation of waves in excitable media. Physica D: NonlinearPhenom. (1-2), 61–70 (1991)4. Barkley, D.: Linear stability analysis of rotating spiral waves in excitable media. Phys. Rev. Lett. (13),2090–2093 (1992)5. Barkley, D.: Euclidean symmetry and the dynamics of rotating spiral waves. Phys. Rev. Lett. (1),164–167 (1994)6. Barkley, D., Kness, M., Tuckerman, L.: Spiral-wave dynamics in a simple model of excitable media:The transition from simple to compound rotation. Physic. Rev. A (4), 2489–2492 (1990)7. Berglund, N., Gentz, B.: Noise-induced phenomena in slow-fast dynamical systems: a sample-pathsapproach, vol. 246. Springer Berlin (2006)4 Muriel Boulakia et al.8. Bonaccorsi, S., Mastrogiacomo, E.: Analysis of the stochastic Fitzhugh–Nagumo system. Infin. Dimens.Anal., Quantum Probab. Rel. Top. (3), 427–446 (2008)9. Boulakia, M., Cazeau, S., Fernández, M., Gerbeau, J.F., Zemzemi, N.: Mathematical modeling of elec-trocardiograms: a numerical study. Ann. Biomed. Eng. (3), 1071–1097 (2010)10. Bréhier, C.E.: Strong and weak order in averaging for SPDEs. Stoch. Proc. Appl. (7), 2553-2593(2012)11. Cao, Y., Yang, H., Yin, H.: Finite element methods for semilinear elliptic stochastic partial differentialequations. Numer. Math. , 181–198 (2007)12. Cerrai, S., Freidlin, M.: Averaging principle for a class of stochastic reaction–diffusion equations.Probab. Theory Rel. Fields (1-2), 137–177 (2009)13. Ciarlet, P., Lions, J.: Finite Element Methods, Handbook of Numerical Analysis , vol. 2. Elsevier (1991)14. Da Prato, G.: Kolmogorov Equations for Stochastic PDEs. Birkhäuser Basel (2004)15. Da Prato, G., Zabczyk, J.: Stochastic equations in infinite dimensions. Cambridge University Press,Cambridge (1992)16. Debussche, A., Printems, J.: Weak order for the discretization of the stochastic heat equation. Math.Comput. (266), 845–863 (2009)17. Fitzhugh, R.: Mathematical models of excitation and propagation in nerve. Biological Engineering(1969)18. Goudenège, L., Martin, D., Vial, G.: High order finite element calculations for the Cahn-Hilliard equa-tion. J. Sci. Comput. (2), 294–321 (2012)19. Hairer, M., Ryser, M., Weber, H.: Triviality of the 2d stochastic Allen-Cahn equation. Electron. J.Probab. , no. 39, 1–14 (2012)20. Hecht, F., Le Hyaric, A., Ohtsuka, K., Pironneau, O.: Freefem++, finite elements software21. Hinch, R.: An analytical study of the physiology and pathology of the propagation of cardiac actionpotentials. Progress in Biophys. and Mol. Bio. (1), 45 – 81 (2002)22. Hodgkin, A., Huxley, A.: Propagation of electrical signals along giant nerve fibres. Proc. Roy. Soc.London. (899), 177–183 (1952)23. Jentzen, A.: Pathwise numerical approximations of SPDEs with additive noise under non-global lipschitzcoefficients. Potential Anal. (4), 375–404 (2009)24. Jentzen, A., Röckner, M.: A Milstein scheme for SPDEs. arXiv preprint arXiv:1001.2751 (2012)25. Jordan, P., Christini, D.: Cardiac Arrhythmia. John Wiley & Sons, Inc. (2006)26. Keener, J.: Waves in excitable media. J. Appl. Math. (3), 528–548 (1980)27. Kovács, M., Larsson, S., Lindgren, F.: Strong convergence of the finite element method with truncatednoise for semilinear parabolic stochastic equations with additive noise. Numer. Algorithms (2-3),309–320 (2010)28. Kruse, R.: Optimal error estimates of Galerkin finite element methods for stochastic partial differentialequations with multiplicative noise. IMA J. of Numer. Anal. (1), 217–251 (2014)29. Kruse, R., Larsson, S.: Optimal regularity for semilinear stochastic partial differential equations withmultiplicative noise. Electron. J. Probab., 17(65), 1–19 (2012).30. Lindner, B., Garcia-Ojalvo, J., Neiman, A., Schimansky-Geier, L.: Effects of noise in excitable systems.Phys. Rep. (6), 321–424 (2004)31. Lord, G., Tambue, A.: A modified semi-implict Euler-Maruyama scheme for finite element discretizationof SPDEs. arXiv preprint arXiv:1004.1998 (2010)32. Lord, G., Tambue, A.: Stochastic exponential integrators for finite element discretization of SPDEs formultiplicative and additive noise. IMA J. Numeric. Anal. (2), 515–543 (2013)33. Lord, G., Thümmler, V.: Computing stochastic traveling waves. SIAM J. Scientific Comput. (1),24–43 (2012)34. Mitchell, C., Schaeffer, D.: A two-current model for the dynamics of cardiac membrane. Bull. Math.Bio. (5), 767–793 (2003)35. Peszat, S., Zabczyk, J.: Stochastic Partial Differential Equations with Lévy Noise. Cambridge UniversityPress (2007)36. Raviart, P., Thomas, J.: Introduction à l’analyse numérique des équations aux dérivées partielles. Masson(1983)37. Shardlow, T.: Numerical simulation of stochastic PDEs for excitable media. J. Comput. Appl. Math. (2), 429–446 (2005)38. Tuckwell, H., Jost, J.: Weak noise in neurons may powerfully inhibit the generation of repetitive spikingbut not its propagation. PLoS Comput. Bio. (5) (2010)39. Wagner, D.: Survey of measurable selection theorems: An update, Lecture Notes in Mathematics , vol.794. Springer Berlin Heidelberg (1980)40. Walsh, J.: Finite element methods for parabolic stochastic PDEs. Potential Anal. (1), 1–43 (2005)41. Wang, W., Roberts, A.: Average and deviation for slow–fast stochastic partial differential equations. J.Differ. Equ. (5) 1265–1286 (2012)42. Yan, Y.: Galerkin finite element methods for stochastic parabolic partial differential equations. SIAM J.Numeric. Anal.43