Simulation of SPDEs for Excitable Media Using Finite Elements

Abstract

In this paper, we address the question of the discretization of stochastic partial differential equations (SPDEs) for excitable media. Working with SPDEs driven by colored 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 or Mitchell–Schaeffer models.

Appendix: Proof of Theorem 1

Recall that the domain $D$ is polyhedral such that

$$\begin{aligned} \overline{D}=\bigcup _{T\in {\fancyscript{T}}_h}T. \end{aligned}$$

Let $i\in \{0,0_a,1\}$. The process $(D_h(t),~t\in [0,\tau ])$ defined by

$$\begin{aligned} D_h(t)=W^Q_t-W^{Q,h,i}_t \end{aligned}$$

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

$$\begin{aligned} \mathbb {E}\left( \sup _{t\in [0,\tau ]}\Vert D_h(t)\Vert ^2\right) \le c_2\mathbb {E}\left( \Vert D_h(\tau )\Vert ^2\right) \end{aligned}$$

with $c_2$ a constant which does not depend on $h$ or $\tau$. We begin with the case $i=1$. Since the processes $W^Q$ and $W^{Q,h,1}$ are regular in space, we write

$$\begin{aligned} \mathbb {E}\left( \Vert D_h(\tau )\Vert ^2\right) =\mathbb {E}\left( \int _D\left( W^Q_\tau (x)-W^{Q,h,1}_\tau (x)\right) ^2\mathrm{d}x\right) . \end{aligned}$$

We use the definition of $W^{Q,h,1}$ in Definition 3 and the fact that $\sum _{i=1}^{N_h}\psi _i=1$ to obtain

$$\begin{aligned} \mathbb {E}\left( \Vert D_h(\tau )\Vert ^2\right)&= \mathbb {E}\left( \int _D\left( W^Q_\tau (x)-\sum _{i=1}^{N_h}W^Q_\tau (P_i)\psi _i(x)\right) ^2\mathrm{d}x\right) \\&= \mathbb {E}\left( \int _D\left( \sum _{i=1}^{N_h}\left( W^Q_\tau (x)-W^Q_\tau (P_i)\right) \psi _i(x)\right) ^2\mathrm{d}x\right) \\&= \mathbb {E}\left( \int _D\sum _{i,j=1}^{N_h}\left( W^Q_\tau (x){-}W^Q_\tau (P_i)\right) \left( W^Q_\tau (x){-}W^Q_\tau (P_j)\right) \psi _i(x)\psi _{j}(x)\mathrm{d}x\right) . \end{aligned}$$

By an application of Fubini’s theorem, exchanging over the expectation, integral and summation, we get

$$\begin{aligned} \mathbb {E}\left( \Vert D_h(\tau )\Vert ^2\right)&= \sum _{i,j=1}^{N_h}\int _D\mathbb {E}\left( \left( W^Q_\tau (x)-W^Q_\tau (P_i)\right) \left( W^Q_\tau (x)-W^Q_\tau (P_j)\right) \right) \psi _i(x)\psi _{j}(x)\mathrm{d}x\\&= \tau \sum _{i,j=1}^{N_h}\int _D\left( C(0){-}C(P_i-x){-}C(P_j-x){+}C(P_i-P_j)\right) \psi _i(x)\psi _{j}(x)\mathrm{d}x. \end{aligned}$$

For all $1\le i,j\le N_h$, if the intersection of the supports of $\psi _i$ and $\psi _j$ is not empty, then

$$\begin{aligned} \forall x\in \text {supp} \psi _i,\quad \forall y\in \text {supp}\psi _j, |x-y|\le K h. \end{aligned}$$

Thus, there exists $K>0$ such that, for all $i,j$, if $\text {supp} \psi _i\cap \text {supp} \psi _{j}\ne \emptyset$ and $x\in \text {supp} \psi _i\cap \text {supp} \psi _{j}$, a Taylor’s expansion yields

$$\begin{aligned} |C(0)-C(P_i-x)-C(P_j-x)+C(P_i-P_j)|\le K\max _{x\in \overline{D}}\Vert \mathrm{Hess}~C(x)\Vert h^2, \end{aligned}$$

where we have used the fact that $\nabla C(0)=0$. Then,

$$\begin{aligned} \mathbb {E}\left( \Vert D_h(\tau )\Vert ^2\right) \le K\tau \max _{x\in \overline{D}}\Vert \mathrm{Hess}~C(x)\Vert h^2. \end{aligned}$$

This ends the proof for the case $i=1$. The case $i=0$ can be treated similarly.

For the case $i=0_a$, we proceed as follows. The process $W^{Q,h,0_a}$ is the orthonormal projection of $W^Q$ on the space P0, thus, we have, using the Pythagorean theorem,

$$\begin{aligned} \mathbb {E}\left( \Vert D_h(\tau )\Vert ^2\right) =\mathbb {E}\left( \left\| W^Q_\tau -W^{Q,h,0_a}_\tau \right\| ^2\right) =\mathbb {E}\left( \left\| W^Q_\tau \right\| ^2-\left\| W^{Q,h,0_a}_\tau \right\| ^2\right) . \end{aligned}$$

Then, recalling that the processes $W^Q$ and $W^{Q,h,0_a}$ are regular in space and using the fact that the triangles $T\in {\fancyscript{T}}_h$ do not intersect, we obtain,

$$\begin{aligned} \mathbb {E}\left( \Vert D_h(\tau )\Vert ^2\right)&=\mathbb {E}\left( \int _DW^Q_\tau (x)^2\mathrm{d} x-\sum _{T\in {\fancyscript{T}}_h}\frac{1}{|T|}(W^Q_\tau ,1_T)^2\right) . \end{aligned}$$

By an application of Fubini’s theorem, exchanging over the expectation and summation, we get

$$\begin{aligned} \mathbb {E}\left( \Vert D_h(\tau )\Vert ^2\right) =\tau \left( C(0)|D|-\sum _{T\in {\fancyscript{T}}_h}\frac{1}{|T|}(Q1_T,1_T)\right) . \end{aligned}$$

(28)

Since $\overline{D}=\bigcup _{T\in {\fancyscript{T}}_h}T$ we have

$$\begin{aligned} C(0)|D|=\sum _{T\in {\fancyscript{T}}_h}\frac{1}{|T|}\int _T\int _TC(0)\mathrm{d}z_1\mathrm{d}z_2, \end{aligned}$$

hence, plugging in (28)

$$\begin{aligned} \mathbb {E}\left( \Vert D_h(\tau )\Vert ^2\right) =\tau \sum _{T\in {\fancyscript{T}}_h}\frac{1}{|T|}\int _T\int _T[C(0)-C(z_1-z_2)]\mathrm{d}z_1\mathrm{d}z_2. \end{aligned}$$

(29)

Thanks to the fact that $\nabla C=0$, a Taylor’s expansion yields

$$\begin{aligned} C(0)-C(z_1-z_2)=(z_1-z_2)\cdot \mathrm{Hess}~C(0)(z_1-z_2)+\mathrm{o}(|z_1-z_2|^2). \end{aligned}$$

(30)

Thus, thanks to (2), for all $z_1,z_2$ in the same triangle $T$

$$\begin{aligned} |C(0)-C(z_1-z_2)|\le K\max _{x\in \overline{D}}\Vert \mathrm{Hess}~C(x)\Vert h^2, \end{aligned}$$

where the constant $K$ is independent from $T\in {\fancyscript{T}}_h$. Plugging in (29) yields

$$\begin{aligned} \mathbb {E}\left( \Vert D_h(\tau )\Vert ^2\right) \le K\max _{x\in \overline{D}}\Vert \mathrm{Hess}~C(x)\Vert \tau h^2 \end{aligned}$$

for a deterministic constant $K$.