Stochastically forced cardiac bidomain model

The bidomain system of degenerate reaction-diffusion equations is a well-established spatial model of electrical activity in cardiac tissue, with"reaction"linked to the cellular action potential and"diffusion"representing current flow between cells. The purpose of this paper is to introduce a"stochastically forced"version of the bidomain model that accounts for various random effects. We establish the existence of martingale (probabilistic weak) solutions to the stochastic bidomain model. The result is proved by means of an auxiliary nondegenerate system and the Faedo-Galerkin method. To prove convergence of the approximate solutions, we use the stochastic compactness method and Skorokhod-Jakubowski a.s.~representations. Finally, via a pathwise uniqueness result, we conclude that the martingale solutions are pathwise (i.e., probabilistic strong) solutions.

introduced the first mathematical model for the propagation of electrical signals along nerve fibers. This model was later tweaked to describe assorted phenomena in biology. Similar to nerve cells, conduction of electrical signals in cardiac tissue rely on the flow of ions through so-called ion channels in the cell membrane. This similarity has led to a number of cardiac models based on the Hodgkin-Huxley formalism [11,13,34,45,48,54]. Among these is the bidomain model [56], which is regarded as an apt spatial model of the electrical properties of cardiac tissue [13,54]. The bidomain equations result from the principle of conservation of current between the intra-and extracellular domains, followed by a homogenization process of the cellular model defined on a periodic structure of cardiac tissue (see, e.g., [13]). The bidomain model can be viewed as a PDE system, consisting of a degenerate parabolic (reactiondiffusion) PDE for the transmembrane potential and an elliptic PDE for the extracellular potential. These PDEs are supplemented by a nonlinear ODE system for the conduction dynamics of the ion channels. There are many membrane models of cardiac cells, differing in their complexity and in the level of detail with which they represent the biology (see [11] for a review). Herein we will utilize a simple model for voltage-gated ion channels [40].
The idiom "bidomain" reflects that the intra-and extracellular tissues are viewed as two superimposed anisotropic continuous media, with different longitudinal and transversal conductivities. If these conductivities are equal, then we have the so-called monodomain model (elliptic PDE reduces to an algebraic equation). The degenerate structure of the bidomain PDE system is due to the anisotropy of cardiac tissue [2,15]. Solutions exhibit discontinuous-like propagating excitation fronts. This, together with strongly varying time scales, makes the system difficult to solve by numerical methods.
The bidomain model is a deterministic system. This means that at each moment in time, the solution can be inferred from the prescribed data. This is at variance with several phenomena happening at the microscopic (cellular) and macroscopic (heart/torso) scales, where respectively channel noise and external random perturbations acting in the torso can play important roles. At the macroscopic level, the ECG signal, a coarse-grained representation of the electrical activity in the heart, is often contaminated by noise. One source for this noise is the fluctuating environment of the heart. In [39], the authors argue that such randomness cannot always be suppressed. Occasionally deterministic equations give qualitatively incorrect results, and it is important to quantify the nature of the noise and choose an appropriate model incorporating randomness.
At the cellular level, the membrane potential is due to disparities in ion concentrations (e.g., sodium, calcium, potassium) across the cell membrane. The ions move through the cell membrane due to random transitions between open and close states of the ion channels. The dynamics of the voltage potential reflect the aggregated behaviour of the individual ion channels, whose conformational changes control the conductance of each ionic current. The profound role of channel noise in excitable cells is summarised and discussed in [26]. Faithful modeling of channel noise gives raise to continuous-time Markov chains with voltage-dependent transition probabilities. In the limit of infinitely many ion channels, these models lead to deterministic Hodgkin-Huxley type equations. To capture channel noise, an alternative (and computationally much simpler) approach is to add well-placed stochastic terms to equations of the Hodgkin-Huxley type [26,38]. Indeed, recent studies (see [26] for a synthesis) indicate that this approach can give an accurate reproduction of channel fluctuations. For work specifically devoted to cardiac cells, see [19,39,45].
where χ m is the ratio of membrane surface area to tissue volume and c m > 0 is the (surface) capacitance of the membrane per unit area. The (nonlinear) function I ion (v, w) represents the ionic current per unit surface area, which depends on the transmembrane potential v and a vector w of ionic (recovery, gating, concentrations, etc.) variables. A simplified model, frequently used for analysis, assumes that the functional form of I ion is a cubic polynomial in v. The ionic variables w are governed by an ODE system, (1.7) Consisting of two (degenerate) parabolic PDEs, the system (1.5), (1.7) is occasionally referred to as the parabolic-parabolic form of the bidomain model. On the subject of well-posedness, i.e., existence, uniqueness, and stability of properly defined solutions, we remark that standard theory for parabolic-elliptic systems does not apply naturally. The main reason is that the anisotropies of the intra-and extracellular domains differ, entailing the degenerate structure of the system. Moreover, a maximum principle is not available. That being the case, a number of works [1,2,5,6,13,15,23,36,57] have recently provided well-posedness results for the bidomain model, applying differing solution concepts and technical frameworks.
1.3. Stochastic model & main results. The purpose of the present paper is to introduce and analyze a bidomain model that accounts for random effects (noise), by way of a few well-placed stochastic terms. The simplest way to insert randomness is to add Gaussian white noise to one or more of the ionic ODEs (1.5), leading to a system of (Itô) stochastic differential equations (SDEs): where W w is a cylindrical Wiener process, with noise amplitude α. Formally, we can think of α dW w as k≥1 α k dW w k (t), where {W w k } k≥1 is a sequence of independent 1D Brownian motions and {α k } k≥1 is a sequence of noise coefficients. Interpreting w as gating variables representing the fraction of open channel subunits of varying types, in [26] this type of noise is referred to as subunit noise. We will allow for subunit noise in our model, assuming for simplicity that the ionic variable w is a scalar and that the noise amplitude depends on the transmembrane potential v, α = α(v) (multiplicative noise). We will also introduce fluctuations into the bidomain system by replacing the PDEs (1.7) with the (Itô) stochastic partial differential equations (SPDEs) χ m c m dv − ∇ · (M i ∇u i ) dt + χ m I ion (v, w) dt = I i dt + β dW v χ m c m dv + ∇ · (M e ∇u e ) dt + χ m I ion (v, w) dt = −I e dt + β dW v , (1.9) where W v is a cylindrical Wiener process (independent of W w ), with noise amplitude β. Adding a stochastic term to the equation for the membrane potential v is labeled current noise in [26]. Current noise represents the aggregated effect of the random activity of ion channels on the voltage dynamics. Allowing the noise amplitude in (1.9) to depend on the membrane voltage v, we arrive at equations with so-called conductance noise [26]. The nonlinear term I ion (v, w) accounts for the total conductances of various ionic currents, and conductance noise pertains to adding "white noise" to the deterministic values of the conductances, i.e., replacing I ion by I ion +β(v) dWv dt , for some functionβ. Herein we include this case by permitting β in (1.9) to depend on the voltage variable v, β = β(v).
Our main contribution is to establish the existence of properly defined solutions to the SDE-SPDE system (1.8), (1.9). From the PDE perspective, we are searching for weak solutions in a certain Sobolev space (H 1 ). From the probabilistic point of view, we are considering martingale solutions, sometimes also referred to as weak solutions. The notions of weak & strong probabilistic solutions have different meaning from weak & strong solutions in the PDE literature. If the stochastic elements are fixed in advance, we speak of a strong (or pathwise) solution. The stochastic elements are collected in a stochastic basis Ω, F , {F t } t∈[0,T ] , P, W , where W = (W w , W v ) are cylindrical Wiener processes adapted to the filtration {F t } t∈[0,T ] . Whenever these elements constitute a part of the unknown solution, the relevant notion is that of a martingale solution. The connection between weak and strong solutions to Itô equations is exposed in the famous Yamada-Watanabe theorem, see, e.g., [31,37,44]. We reserve the name weak martingale solution for solutions that are weak in the PDE sense as well as being probabilistic weak.
We will prove that there exists a weak martingale solution to the stochastic bidomain system. Motivated by the approach in [2] (see also [5]) for the deterministic system, we use the Faedo-Galerkin method to construct approximate solutions, based on an auxiliary nondegenerate system obtained by adding εdu i and −εdu e respectively to the first and second equations in (1.9) (ε is a small positive parameter). The stochastic compactness method is put to use to conclude subsequential convergence of the approximate solutions.
Indeed, we first apply the Itô chain rule to derive some basic a priori estimates. The combination of multiplicative noise and the specific structure of the system makes these estimates notably harder to obtain than in the deterministic case. The a priori estimates lead to strong compactness of the approximations in the t, x variables (in the deterministic context [2]). In the stochastic setting, there is an additional (probability) variable ω ∈ D in which strong compactness is not expected. Traditionally, one handles this issue by arguing for weak compactness of the probability laws of the approximate solutions, via tightness and Prokhorov's theorem [31]. The ensuing step is to construct a.s. convergent versions of the approximations using the Skorokhod representation theorem. This theorem supplies new random variables on a new probability space, with the same laws as the original variables, converging almost surely. Equipped with a.s. convergence, we are able to show that the limit variables constitute a weak martingale solution. Finally, thanks to a uniqueness result and the Gyöngy-Krylov characterization of convergence in probability [27], we pass a la Yamada-Watanabe [31] from martingale to pathwise (probabilistic strong) solutions.
Martingale solutions and the stochastic compactness method have been harnessed by many authors for different classes of SPDEs, see e.g. [3,4,16,18,21,22,24,28,30,41,46,49] for problems related to fluid mechanics. An important step in the compactness method is the construction of almost surely convergent versions of processes that converge weakly. This construction dates back to the work of Skorokhod, for processes taking values in a Polish (complete separable metric) space [16,31]. The classical Skorokhod theorem is befitting for the transmembrane variable v, but not the intracellular and extracellular variables u i , u e . This fact is a manifestation of the degenerate structure of the bidomain system, necessitating the use of a Bochner-Sobolev space equipped with the weak topology. We refer to Jakubowski [32] for a recent variant of the representation theorem that applies to so-called quasi-Polish spaces, specifically allowing for separable Banach spaces equipped with the weak topology, as well as spaces of weakly continuous functions with values in a separable Banach space. We refer to [7,8,9,10,43,53] for works making use of Skorokhod-Jakubowski a.s. representations.
The remaining part of this paper is organized as follows: The stochastic bidomain model is presented in Section 2. Section 3 outlines the underlying stochastic framework and list the conditions imposed on the "stochastic" data of the model. Solution concepts and the accompanying main results are collected in Section 4. The approximate (Faedo-Galerkin) solutions are constructed in Section 5. In Section 6 we establish several a priori estimates and prove convergence of the approximate solutions, thereby providing an existence result for weak martingale solutions. A pathwise uniqueness result is established in Section 7, which is then used in Section 8 to upgrade martingale solutions to pathwise solutions.

STOCHASTIC BIDOMAIN MODEL
The spatial domain of the heart is given by a bounded open set Ω ⊂ R 3 with piecewise smooth boundary ∂Ω. This three-dimensional slice of the cardiac muscle is viewed as two superimposed (anisotropic) continuous media, representing the intracellular (i) and extracellular (e) tissues. The tissues are connected at each point via the cell membrane. In our earlier outline of the (deterministic) bidomain model, we saw that the relevant quantities are the intracellular and extracellular potentials as well as the transmembrane potential v := u i − u e (defined in Ω T ). The conductivities of the intracellular and extracellular tissues are encoded in anisotropic matrices M i = M i (x) and M e = M e (x). Cardiac tissue is composed of fibers, with the conductivity being higher along the fibers than in the cross-fibre direction. The cardiac fibers are organized in sheets of varying surface orientation, giving raise to three principal directions for conduction: parallel to the fibers, perpendicular to the fibers but parallel to the sheet, and perpendicular to the sheet. When the fibers rotate from bottom to top, we have cardiac tissue with rotational anisotropy. Without rotation (axisymmetric anisotropy), the conductivity tensors take the form is a unit vector giving the fiber direction and σ l j = σ l j (x), σ t j = σ t j (x) are coefficients describing respectively the intra-and extracellular conductivities along and transversal to the fibre direction. Whenever the fibers are aligned with the axes, M i and M e are diagonal matrices. For more details, see, e.g., [12,13,54]. In this paper, we do not exploit structural properties of cardiac tissue, that is, we assume only that M i and M e are general (bounded, positive definitive) matrices, cf. (2.6) below. The stochastic bidomain model contains two nonlinearly coupled SPDEs involving the potentials u i , u e , v. These stochastic reaction-diffusion equations are further coupled to a nonlinear SDE for the gating (recovery) variable w. The dynamics of (u i , u e , v, w) is governed by the equations where c m > 0 is the surface capacitance of the membrane, χ m is the surface-to-volume ratio, and I i , I e are stimulation currents. In (2.2), randomness is represented by cylindrical Wiener processes W v , W w with nonlinear noise amplitudes β, α (cf. Section 3 for details).
We impose initial conditions on the transmembrane potential and the gating variable: The intra-and extracellular domains are often assumed to be electrically isolated, giving raise to zero flux (Neumann type) boundary conditions on the potentials u i , u e [13,54]. From a mathematical point of view, Dirichlet and mixed Dirichlet-Neumann type boundary conditions are utilized in [1] and [2], respectively. Herein we partition the boundary ∂Ω into regular parts Σ N and Σ D and impose the mixed boundary conditions (j = i, e) where ν denotes the exterior unit normal to the "Neumann part" Σ N of the boundary, which is defined a.e. with respect to the two-dimensional Hausdorff measure H 2 on ∂Ω.
Observe that the equations in (2.2) are invariant under the change of u i and u e into u i + k, u e + k, for any k ∈ R. Hence, unless Dirichlet conditions are imposed somewhere (Σ D = ∅), the bidomain system determines the electrical potentials only up to an additive constant. To ensure a unique solution in the case Σ D := ∅ (∂Ω = Σ N ), we may impose the normalization condition To avoid making this paper too long, we assume that Σ D = ∅. Moreover, we stick to homogenous boundary conditions, although we could have replaced the right-hand sides of (2.4) by sufficiently (Sobolev) regular functions. Regarding the "membrane" functions I ion and H, we have in mind the fairly uncluttered FitzHugh-Nagumo model [20,42]. This is a simple choice for the membrane kinetics that is often used to avoid difficulties arising from a large number of coupling variables. The model is specified by where the parameter a represents the threshold for excitation, ǫ represents excitability, and κ, γ, δ are parameters that influence the overall dynamics of the system. For background material on cardiac membrane models and their general mathematical structure, we refer to the books [13,34,54].
In an attempt to simplify the notation, we redefine and set We also assume I i , I e ≡ 0, as these source terms do not add new difficulties. The resulting stochastic bidomain system becomes along with the initial and boundary conditions (2.3) and (2.4). The cylindrical Wiener processes W v , W w in (2.5) are defined in Section 3. With regard to the conductivity matrices in (2.5), we assume the existence of positive constants m, M such that for j = i, e, Motivated by the discussion above on membrane models, we impose the following set of assumptions on the functions I, H in (2.5): • Generalized FitzHugh-Nagumo model (GFHN): where I 1 , I 2 , h ∈ C 1 (R) and for all v ∈ R, for some positive constants c I,1 , c I,2 , c I,3 , c I,4 , c H,1 , c H,2 and c I > 0.
We end this section with a remark about the so-called monodomain model.
Remark 2.1. The stochastic bidomain model simplifies if the anisotropy ratios σ l i /σ t i and σ l e /σ t e are equal, cf. (2.1). Indeed, suppose σ l i /σ t i = σ l e /σ t e and moreover σ t i = λσ t e for some constant λ > 0 (and thus σ l i = λσ l e ). Then, in view of (2.1), it follows that M i = λM e , and hence the first two equations in (2.5) can be combined into a single equation; thereby arriving at the stochastic monodomain system where M := λ 1+λ M i . The system (2.7) is a significant simplification of the bidomain model (2.5), and even though the assumption of equal anisotropy ratios is very strong, the monodomain model is adequate in certain situations [14].

STOCHASTIC FRAMEWORK
To define the cylindrical Wiener processes W = W v , W w and the stochastic integrals appearing in (2.5), we need to briefly recall some basic concepts and results from stochastic analysis. For a detailed account we refer to [16,44] (see also [31,33]).
We consider a complete probability space (D, F , P ), along with a complete rightcontinuous filtration {F t } t∈[0,T ] . Without loss of generality, we assume that the σ-algebra F is countably generated. Let B be a separable Banach space (equipped with the Borel σ-algebra B(B)). Then a B-valued random variable X is a measurable mapping from (D, F , P ) to (B, B(B)), D ∋ ω → X(ω) ∈ B. The expectation of a random variable X is denoted E[X] := D X dP . We use the abbreviation a.s. (almost surely) for P -almost every ω ∈ D. The collection of all (equivalence classes of) B-valued random variables is denoted by L 1 (D, F , P ). Equipped with the norm X L 1 (D,F ,P ) = E [ X B ] it becomes a Banach space. For p > 1 we define L p (D, F , P ) similarly, with norm In this paper we deal with time-dependent functions and processes, so B is typically an "evolutionary" (Bochner) space like B = L qt ((0, T ); L qx (Ω)), q t , q x ∈ [1, ∞].
A stochastic process X = {X(t)} t∈[0,T ] is a collection of B-valued random variables The paths t → X(ω, t) of a (jointly) measurable process X are automatically Borel measurable functions. Of course, it is the use of filtrations and adaptivity that differentiates the theory of stochastic processes from the one of functions on product spaces. A stochastic process X is adapted if be a sequence of independent one-dimensional Brownian motions adapted to the filtration {F t } t∈[0,T ] . We refer to as a (Brownian) stochastic basis. When a filtration is involved there are additional notions of measurability (predictable, optional, progressive) that occasionally are more convenient to work with. Herein we use the (stronger) notion of a predictable process, but we could also have used that of a progressive process. A predictable process is a P T × B([0, T ]) i.e., the σ-algebra generated by all left-continuous adapted processes; P T is also generated by the sets A predictable process is adapted and measurable. Although the converse implication is not true, adaptive processes with regular (e.g. continuous) paths are predictable. Fix a separable Hilbert space U, equipped with a complete orthonormal basis {ψ k } k≥1 . We utilize cylindrical Brownian motions W evolving over U. Informally, the are defined by W := k≥1 W k ψ k . We have to be a bit more precise, however, since the series does not converge in U: Let X be a separable Hilbert space with inner product (·, ·) X and norm · X . For the bidomain model (2.5), a natural choice is X = L 2 (Ω). The vector space of all bounded linear operators from U to X is denoted L(U, X). We denote by L 2 (U, X) the collection of Hilbert-Schmidt operators from U to X, i.e., R ∈ L 2 (U, X) ⇐⇒ R ∈ L(U, X) and We turn L 2 (U, X) into a Hilbert space by supplying it with the inner product Returning to the convergence of W = k≥1 W k ψ k , it is always possible to construct an auxiliary Hilbert space U 0 ⊃ U for which there exists a Hilbert-Schmidt embedding J : U → U 0 . Indeed, given a sequence {b k } k≥1 ∈ ℓ 2 with b k = 0 for all k (e.g., b k = 1/k), introduce the space and equip it with the norm u U0 := k≥1 a 2 k b 2 k 1 2 , u = k≥1 a k ψ k . Clearly, the embedding J : U → U 0 , u → u, is Hilbert-Schmidt. Now, on this larger Hilbert space U 0 , the infinite series k≥1 W k ψ k converges; indeed, as n → ∞, To summarize, the stochastic process is referred to as a cylindrical Brownian motion evolving over U. The right-hand side of (3.4) converges on the Hilbert space U 0 , with the embedding U ⊂ U 0 being Hilbert-Schmidt. Via standard martingale arguments, W is a.s. continuous with values in U 0 , i.e., W (ω, ·, ·) belongs to C([0, T ]; U 0 ) for P -a.e. ω ∈ D, and also L 2 (D, F , P ; C([0, T ]; U 0 )).
Without loss of generality, we assume that the filtration {F t } t∈[0,T ] is generated by W and the initial condition. See [16,44] for details. Given a cylindrical Brownian motion W , we can define the Itô stochastic integral GdW as follows [16,44]: provided the integrand G is a predictable X-valued process satisfying G ∈ L 2 D, F , P ; L 2 ((0, T ); L 2 (U, X)) .
The stochastic integral (3.5) is an X-valued square integrable martingale, satisfying the Burkholder-Davis-Gundy inequality where C is a constant depending on p ≥ 1. In terms of the basis {ψ k } k≥1 , (3.6) reads For the bidomain model (2.5), we take X = L 2 (Ω). With this choice, we can give meaning to the stochastic terms appearing in the weak formulation of (2.5), with ϕ ∈ L 2 (Ω). Since W = k≥1 W k ψ k is a cylindrical Brownian motion, we can write knowing that the series converges in L 2 (D, F , P ; C([0, T ])), where β k (v) := β(v)ψ k are real-valued functions. Sometimes we denote the right-hand side by t 0 Ω β(v)ϕ dx dW v . We need to impose conditions on the noise amplitudes β = η, σ. For each v ∈ L 2 (Ω), we assume that β(v) : U → L 2 (Ω) is defined by for a constant C β > 0. Thus, β becomes a mapping from L 2 (Ω) to L 2 U, L 2 (Ω) and satisfies, via (3.2), and similarly for β(v 1 ) − β(v 2 ) 2 L2(U,L 2 (Ω)) . To be sure, by way of (3.8), we have . Given a predictable process v ∈ L 2 D, F , P ; L 2 ((0, T ); L 2 (Ω)) , the stochastic integral t 0 β(v) dW w is well-defined, taking values in L 2 (Ω). By (3.9), where we have also used the Cauchy-Schwarz and Burkholder-Davis-Gundy inequalities. Hence, (3.7) makes sense. It is possible to allow β = η, σ to be time dependent, β = β(t, v). In this case, β(t) : L 2 (Ω) → L 2 U, L 2 (Ω) must satisfy (3.8) for a.e. t ∈ [0, T ], with a constant C β that is independent of t. This does not entail additional effort in the proofs, but for simplicity of presentation we suppress the time dependency throughout the paper.
We conclude this section by collecting a few relevant probability tools. Let S be a separable Banach (or Polish) space. We denote by B(S) the family of the Borel subsets of S and by P(S) the family of all Borel probability measures on S. Each random variable X : D → S induces a probability measure on S via the pushforward X # P : We will construct weak martingale solutions by applying the stochastic compactness method to a sequence of approximate solutions. In one step of the argument, we show tightness of the probability laws of the approximations. By the Prokhorov theorem, this is equivalent to exhibiting weak compactness of the laws. A sequence {µ n } n≥ ⊂ P(S) is weakly (narrowly) convergent to µ ∈ P(S) as n → ∞ if  Relating to convergence of the approximate solutions, it is essential that we secure strong compactness (a.s. convergence) in the ω variable. To that end, we are in need of a Skorokhod a.s. representation theorem [31], delivering a new probability space and new random variables, with the same laws as the original ones, converging almost surely.
As alluded to before, our path space is not a Polish space since weak topologies in Hilbert and Banach spaces are not metrizable. Therefore the original Skorokhod theorem is not applicable; instead we must rely on the more recent Jakubowski version [32] that applies to so-called quasi-Polish spaces. "Quasi-Polish" refers to spaces S for which there exists a countable family of continuous functionals that separate points (of S) [32]. Quasi-Polish spaces include separable Banach spaces equipped with the weak topology, and also spaces of weakly continuous functions taking values in some separable Banach space. The basic assumption (3.10) gives raise to a mapping between S and the Polish space which is one-to-one and continuous, but in generalf is not a homeomorphism of S onto a subspace of S. However, if we restrict to a σ-compact subspace of S, thenf becomes a measurable isomorphism [32]. In this paper we use the following form of the Skorokhod-Jakubowski theorem [32], taken from [8,43] (see also [9,10]). (2) every Borel subset of a σ-compact set in S belongs to Σ; (3) every probability measure supported by a σ-compact set in S has a unique Radon extension to the Borel σ-algebra B(S); (4) if {µ n } n≥1 is a tight sequence of probability measures on (S, Σ), then there exist a subsequence {n k } k≥1 , a probability space (D,F ,P ), and Borel measurable S-valued random variablesX k ,X, such that µ n k is the law ofX k and X k → XP -a.s. (in S). Moreover, the law µ ofX is a Radon measure.
We will need the Gyöngy-Krylov characterization of convergence in probability [27]. It will be used to upgrade weak martingale solutions to strong (pathwise) solutions, via a pathwise uniqueness result.

Lemma 3.4 (Gyöngy-Krylov characterization)
. Let S be a Polish space, and let {X n } n≥1 be a sequence of S-valued random variables on a probability space (D, F , P ). For each n, m ≥ 1, denote by µ n,m the joint law of (X n , X m ), that is, Then {X n } n≥1 converges in probability (and P -a.s. along a subsequence) ⇐⇒ for any subsequence {µ m k ,n k } k≥1 there exists a further subsequence that converges weakly to some µ ∈ P(S) that is supported on the diagonal:

Remark 3.5. As a matter of fact, we need access to the "⇐=" part of the Gyöngy-Krylov lemma for quasi-Polish spaces S. Suppose for any subsequence
that converges in distribution to (X, X) as j → ∞, for some X ∈ S, that is, the joint probability laws µ m k j ,n k j converge weakly to some µ ∈ P(S × S) that is supported on the diagonal. Recalling the mappingf between S and the Polish space [−1, 1] L , cf. (3.11), and the continuous mapping theorem, it follows that the sequence converges in probability and thus, it is not difficult to see that this implies that X nj j≥1 converges P -a.s. as well.
We are going to need the following convergence result for stochastic integrals due to Debussche, Glatt-Holtz, and Temam [18].
Finally, we recall the "Kolmogorov test" for the existence of continuous modifications of real-valued stochastic processes.
Then there exists a continuous modification of X. The paths of X are γ-Hölder continuous for every γ ∈ 0, δ κ .

NOTION OF SOLUTION AND MAIN RESULTS
Depending on the (probabilistic) notion of solution, the initial data (2.3) are imposed differently. For pathwise (probabilistic strong) solutions, we prescribe the initial data as random variables v 0 , w 0 ∈ L 2 (D, F , P ; L 2 (Ω)). For martingale (or probabilistic weak) solutions, of which the stochastic basis is an unknown component, we prescribe the initial data in terms of probability measures µ v0 , µ w0 on L 2 (Ω). The measures µ v0 and µ w0 should be viewed as "initial laws" in the sense that the laws of v(0), w(0) are required to coincide with µ v0 , µ w0 , respectively.
Sometimes we need to assume the existence of a number q 0 > 9 2 such that (4.1) As a matter of fact, we mostly need (4.1) with q 0 > 2. One exception occurs in Subsection 6.5, where we use q 0 > 9 2 to conclude that the transmembrane potential v is a.s. weakly time continuous, cf. part (5) in the definition below (for w this holds with just q 0 > 2).
Let us define precisely what is meant by a solution to the stochastic bidomain model. For this, we use the space The laws of v 0 := v(0) and w 0 := w(0) are respectively µ v0 and µ w0 : The following identities hold P -almost surely, for any t ∈ [0, T ]: for all ϕ i , ϕ e ∈ H 1 D (Ω) and ϕ ∈ L 2 (Ω).
Our main existence result is contained in The proof of Theorem 4.4 is divided into a series of steps. We construct approximate solutions in Section 5, which are shown to converge in Section 6. The convergence proof relies on several uniform a priori estimates that are established in Subsections 6.1 and 6.2. We use these estimates in Subsection 6.3 to conclude that the laws of the approximate solutions are tight and that the approximations (along a subsequence) converge to a limit. The limit is shown to be a weak martingale solution in Subsections 6.4 and 6.5.
If the stochastic basis S in Definition 4.1 is fixed in advance (not part of the solution), we speak of a weak solution or weak pathwise solution. A weak solution is thus weak in the PDE sense and strong in the probabilistic sense. In this case, we prescribe the initial data v 0 , w 0 as random variables relative to S.
Weak solutions are said to be unique if, given any pair of such solutionsÛ ,Ũ for whicĥ U andŨ coincide a.s. at t = 0, We establish pathwise uniqueness by demonstrating that v(t), w(t) depend continuously on the initial data v 0 , w 0 in L 2 (D, F , P ; L 2 (Ω)). Moreover, using the Poincaré inequality, we conclude as well the pathwise uniqueness of u i , u e . As alluded to earlier, we use this to "upgrade" martingale solutions to weak (pathwise) solutions, thereby delivering Regarding the proof of Theorem 4.6, we divide it into two steps. A pathwise uniqueness result is established in Section 7 by exhibiting an L 2 stability estimate for the difference between two solutions. We use this result in Section 8 to upgrade martingale solutions to pathwise solutions.

CONSTRUCTION OF APPROXIMATE SOLUTIONS
In this section we define the Faedo-Galerkin approximations. They are based on a nondegenerate system introduced below. In upcoming sections we use these approximations to construct weak martingale solutions to the stochastic bidomain model.
We begin by fixing a stochastic basis and F 0 -measurable initial data v 0 , w 0 ∈ L 2 (D; L 2 (Ω)) with respective laws µ v0 , µ w0 on L 2 (Ω). For each fixed ε > 0, the nondegenerate system reads with boundary conditions (2.4). Regarding (5.2), we must provide initial data for u i , u e (not v = u i − u e as in the original problem). For that reason, we decompose (arbitrarily) the initial condition v 0 in (2.3) as v 0 = u i,0 − u e,0 , for some F 0 -measurable random variables u i,0 and u e,0 , such that the law of u i,0 − u e,0 coincides with µ v0 . We replace (2.3) by In some situations, we make use of the strengthened assumption To construct and justify the validity of the Faedo-Galerkin approximations, we employ a classical Hilbert basis, which is orthonormal in L 2 and orthogonal in H 1 D . We refer for example to [50,Thm. 7.7,p. 87] (see also [47]) for the standard construction of such bases. We operate with the same basis {e l } n l=1 for all the unknowns u i , u e , v, w. We look for a solution to the problem arising as the projection of (5.2), (2.3), (2.4) onto the finite dimensional subspace X n := Span{e l } n l=1 . The (finite dimensional) approximate solutions take the form We pick the coefficients , which are finite dimensional stochastic processes relative to (5.1), such that (ℓ = 1, . . . , n) where ε in (5.2) is taken as We need to comment on the finite dimensional approximations of the stochastic terms utilized in (5.8) we can write In (5.8), we utilize the finite dimensional approximation with β n and W n then defined by w n (0) = w n 0 := n l=1 a n l (0)e l , a n l (0) := (w 0 , e l ) L 2 (Ω) .

(5.11)
In (5.11), consider for example u n j,0 . Since u j,0 ∈ L 2 D, F , P ; L 2 (Ω) , we have (by standard properties of finite-dimensional projections, cf. (5.14), (5.16) below) u n j,0 → u j,0 in L 2 (Ω), P -a.s., as n → ∞, and u n j,0 . On this account, the dominated convergence theorem implies , P -a.s., and thus in L 2 ω (L 2 x ). For the basis {e l } ∞ l=1 , we introduce the projection operators (see e.g. [8, page 1636 The restriction of Π n to L 2 (Ω) is also denoted by Π n : i.e., Π n is the orthogonal projection from L 2 (Ω) to Span{e l } ∞ j=1 . We have Note that we have the following equality for any u * ∈ (H 1 D (Ω)) * and u ∈ H 1 D (Ω): Furthermore, as n → ∞, Using the projection operator (5.13), we may write (5.8) in integrated form equivalently as equalities between (H 1 D (Ω)) * valued random variables: where v n 0 = u n i,0 − u n e,0 and u n i,0 = Π n u i,0 , u n e,0 = Π n u e,0 , w n 0 = Π n w 0 . In coming sections we investigate the convergence properties of the sequences u n j n≥1 (5.17). Meanwhile, we must verify the existence of a (pathwise) solution to the finite dimensional system (5.8). Proof. Using the orthonormality of the basis, (5.8) becomes the SDE system (ℓ = 1, . . . , n) d c n ℓ + ε n c n i,ℓ = A i,ℓ dt + Γ ℓ dW v,n , d c n ℓ − ε n c n e,ℓ = A e,ℓ dt + Γ ℓ dW v,n , da n ℓ = A H,ℓ dt + ζ ℓ dW w,n , for the coefficients c n j = c n j (t) (j = i, e) and a n = a n (t), cf. (5.7), where Adding the first and second equations in (5.18) yields (ℓ = 1, . . . , n) (5.20) , a n (t)} be the vector containing all the unknowns in (5.19) and (5.20). For technical reasons, related to (5.22) and (5.23) below, we write the left-hand sides of the first two equations in (5.20) in terms of the ε n scaled quantities √ ε n c n i , √ ε n c n e . Moreover, we view the righthand sides of all the equations as functions of C n (involving the ε n scaled quantities), which can always be done since ε n > 0 is a fixed number. As a result, the constants below may depend on 1/ε n . Let be the vector containing all the drift terms, and , be the collection of noise coefficients. The vector {W v,n , W v,n , W v,n , W w,n } is denoted by W n . Then (5.19) and (5.20) take the compact form (5.21) dC n (t) = F (C n (t)) dt + G(C n (t)) dW n (t), C n (0) = C n 0 , where C n 0 = c n (0), √ ε n c n i (0), √ ε n c n e (0), a n (0) , cf. (5.11). If F, G are globally Lipschitz continuous, classical SDE theory [33,44,52] provides the existence and uniqueness of a pathwise solution. However, due to the nonlinear nature of the ionic models, cf. (GFHN), the global Lipschitz condition does not hold for the SDE system (5.21). As a replacement, we consider the following two conditions: • (local weak monotonicity) ∀C 1 , C 2 ∈ R 4n , |C 1 | , |C 2 | ≤ r, for any r > 0, for some r-dependent positive constant K r . • (weak coercivity) ∀C ∈ R 4n , there exists a constant K > 0 such that Below we verify that the coefficients F, G in (5.21) satisfy these two conditions globally (i.e., (5.22) holds independent of r). Then, in view of Theorem 3.1.1 in [44], there exists a unique global adapted solution to (5.21).
Let us verify the weak monotonicity condition. To this end, set where u n i,1 , u n e,1 , w n 1 and u n i,2 , u n e,2 , w n 2 are arbitrary functions of the form of (5.6), with corresponding time coefficients c n i,1 , c n e,1 , a n 1 and c n i,2 , c n e,2 , a n 2 , respectively. Moreover, set c n 1 := c n i,1 −c n e,1 , c n 2 := c n i,1 −c n e,1 , and C n k := c n k , √ ε n c n i,k , √ ε n c n e,k , a n k for k = 1, 2. We wish to show that i.e., that F is globally one-sided Lipschitz. This requires comparing the "dt-terms" in (5.19) and (5.20) corresponding to the vectors C n 1 and C n 2 , resulting in three different types of terms, linked to the M j (diffusion) part, the I (ionic) part, and the H (gating) part of the equations, that is, I F = I M F + I I F + I H F . First, Similarly, and therefore I I Adding the integrands gives and thus, cf. (2.6), I M F = − j=i,e M j ∇U n j · ∇U n j ≤ 0. Hence, F is globally one-sided Lipschitz. In view of (3.8), it follows easily that G is globally Lipschitz: for some constant K G (depending on n). Summarizing, condition (5.22) holds. In much the same way, again using assumptions (GFHN) and (3.8), we deduce that for some constants K F , K G ; that is to say, condition (5.23) holds.
6. CONVERGENCE OF APPROXIMATE SOLUTIONS 6.1. Basic apriori estimates. To establish convergence of the Faedo-Galerkin approximations, we must supply a series of apriori estimates that are independent of the parameter n (cf. Lemma 6.1 below). At an informal level, assuming that the relevant functions are sufficiently regular, these estimates are obtained by considering where ε n is defined in (5.9), multiplying the first equation by u i , the second equation by −u e , and summing the resulting equations. For the moment, let us assume that the noise W v is one-dimensional and η(v) is a scalar function. To proceed we use the stochastic (Itô) product rule. Hence, we need access to the equation for du i , which turns out to be ∇ · M e ∇u e − 1 2 + ε n I(v, w) dt Note that this equation "blows up" as ε n → 0 (the same is true for the du e equation below). The stochastic product rule gives where Similar computations, this time involving the equation where After some computations we find that Whence, adding (6.2) and (6.3), Adding to this the equation for dw 2 , resulting from (5.2) and Itô's formula, the estimates in Lemma 6.1 below appear once we integrate in x and t, make use of spatial integration by parts, the boundary conditions (2.4), and properties of the nonlinear functions I, H implying (6.14) below. Arguing at the level of finite dimensional approximations, we now convert the computations outlined above into a rigorous proof.
In view of the n-independent estimates in Lemma 6.1, passing if necessary to a proper subsequence, we can assume that the following (weak) convergences hold as n → ∞: w n ⋆ ⇀ w in L 2 D, F , P ; L ∞ ((0, T ); L 2 (Ω)) .
The next result, a consequence of Lemma 6.1 and a martingale inequality, supplies highorder moment estimates, useful when converting a.s. convergence into L 2 convergence. Corollary 6.2. In addition to the assumptions in Lemma 6.1, suppose (5.5) holds with q 0 defined in (4.1). There exists a constant C > 0, independent of n, such that Proof. In view of (6.16), we have the following estimate for any (ω, t) ∈ D × [0, T ]: for some constant C 1 independent of n.
We raise both sides of this inequality to the power q 0 /2, take the expectation, and apply several elementary inequalities, eventually arriving at where Arguing as in (6.18), using a martingale inequality (and some elementary inequalities), for any δ > 0, where we have used that Similarly, relying on (3.8), With δ chosen small, combining (6.24) and (6.25) in (6.23) gives for some constant C 9 > 0 independent of n. Set and note that (6.26) reads , v n (t), w n (t), t ∈ [0, T ], satisfy (5.8), (5.9), (5.10), (5.11). With u n = v n or w n , there is a constant C > 0, independent of n, such that for any sufficiently small δ > 0.
Proof. We assume that v n , u n i , u n e , w n and η n , σ n have been extended by zero outside the time interval [0, T ]. Recalling (5.6) (i.e., v n = u n i − u n e ), it follows that In view of (5.18), see also (5.17), Similarly, using the equation for w n , cf. (5.18) and also (5.2), Integrating over t ∈ (0, T − τ ) and summing the resulting equations gives (6.28) We examine these six terms separately. For the Γ 1 term, noting that thanks to (2.6), we obtain using Cauchy-Schwarz's inequality. Hence, by Young's inequality and (6.5), for some constant C 1 > 0 independent of n.
Next, take notice of the bound t+τ t I (v n (s, x), w n (s, x)) ds where we have used the inequality resulting from (GFHN) and Young's inequality. Due to (6.30), (6.4) and (6.5), 1 4 , and for this reason, in view of Young's inequality and (6.5), Finally, we treat the stochastic terms. By the Cauchy-Schwarz inequality, Applying E[·] along with the Cauchy-Schwarz inequality, we gather the estimate Similarly, Collecting the previous estimates (6.29), (6.32), (6.33), (6.34), and (6.35) we readily conclude from (6.28) that the time translation estimate (6.27) holds.
Remark 6.4. The proof of (6.27), cf. estimate (6.33), reveals that the amount of time continuity of w n is actually better than stated; it is of order δ 6.3. Tightness and a.s. representations. To justify passing to the limit in the nonlinear terms in (5.2), we must show that {v n } n≥1 converges strongly, thereby upgrading the weak L 2 convergence in (6.21). Strong (t, x) convergence is a result of the spatial H 1 D bound (6.5) and the time translation estimate (6.27).
On the other hand, to secure strong (a.s.) convergence in the probability variable ω ∈ D we must invoke some nontrivial results of Skorokhod, linked to tightness of probability measures and a.s. representations of random variables. Actually, there is a complicating factor at play here, namely that the sequences {u n i } n≥1 , {u n e } n≥1 only converge weakly in (t, x) because of the degenerate structure of the bidomain model. As a result, we must turn to the Skorokhod-Jakubowski representation theorem [32], which applies to separable Banach spaces equipped with the weak topology and other so-called quasi-Polish spaces. At variance with the original Skorokhod representations on Polish spaces, the flexibility of the Jakubowski version comes at the expense of having to pass to a subsequence (which may be satisfactory in many situations). We refer to [7,8,9,10,43,53] for works making use of Skorokhod-Jakubowski a.s. representations.
Following [3,41] (for example), the aim is to establish tightness of the probability measures (laws) generated by the Faedo-Galerkin solutions {(U n , W n , U n 0 )} n≥1 , where (6.36) U n = u n i , u n e , v n , w n , W n = W v,n , W w,n , U n 0 = u n i,0 , u n e,0 , v n 0 , w n 0 . Accordingly, we choose the following path space for these measures: where U 0 is defined in (3.3) and the tag "-weak" signifies that the space is equipped with the weak topology. The σ-algebra of Borel subsets of X is denoted by B(X ). We introduce the (X , B(X ))-valued measurable mapping Φ n defined on (D, F , P ) by Φ n (ω) = (U n (ω), W n , U n 0 (ω)) . On (X , B(X )), we define the probability measure (law of of Φ n ) We denote by L u n i , L u n e the respective laws of u n i , u n e on L 2 ((0, T ); H 1 D (Ω T ))-weak, with similar notations for the laws of v n on L 2 (Ω T ), w n on L 2 ((0, T ); (H 1 D (Ω)) * ), W v,n , W w,n on C([0, T ]; U 0 ), and u n i,0 , u n e,0 , v n 0 , w n 0 on L 2 (Ω). Hence, L n = L u n i × L u n e × L v n × L w n × L u n i,0 × L u n e,0 × L v n 0 × L w n 0 .
Again by the Chebyshev inequality (6.43), we infer that where we have used (6.5), (6.7), and (6.27). On the grounds of this and (6.44), we can choose R δ such that (6.39) holds.
Let us introduce the following stochastic basis linked toΦ n , cf. (6.48): S n = D ,F , F n Lemma 6.9 (Faedo-Galerkin equations). Relative to the stochastic basisS n in (6.53), the functionsŨ n ,W n ,Ũ n 0 defined in (6.45) satisfy the following equationsP -a.s.: for each t ∈ [0, T ], where ε n is specified in (5.9) andW v,(n) ,W w,(n) are defined in (6.54). Moreover, (6.56)ṽ n =ũ n i −ũ n e , dP × dt × dx a.e. inD × (0, T ) × Ω, and (by construction)Ũ n ,W n are continuous, adapted (and thus predictable) processes. Finally, the laws ofṽ n 0 andw n 0 coincide with the laws of Π n v 0 and Π n w 0 , respectively, Proof. We establish the first equation in (6.55), with the remaining ones following along the same lines. In accordance with Lemma 5.2 and (6.36), recall that (U n , W n , U n 0 ) is the continuous adapted solution to the Faedo-Galerkin equations (5.17) relative to S, cf. (5.1).
Note that I n = 0 P -a.s. and so E[I n ] = 0. If we could write I n = L n (Φ n ) for a (deterministic) bounded continuous functional L n (·) on X , cf. (6.48), then by equality of the laws, alsoẼ[Ĩ n ] = 0 and the result would follow. However, this is not directly achievable since the stochastic integral is not a deterministic function of W v,n . Hence, certain modifications are needed to produce a workable proof [3]. First of all, we do not consider I n but rather the bounded map I n /(1 + I n ). Noting that E[I n ] = 0 implies (6.57) E I n 1 + I n = 0, the goal is to show that (6.58)Ẽ Ĩ n 1 +Ĩ n = 0, from which the first equation in (6.55) follows.
We defineη n,ν k similarly (with v n replaced byṽ n ). An "integration by parts" reveals that i.e., thanks to the regularization of η n k (v n ) in the t variable, t 0 η n,ν k dW v k (s) can be viewed as a (deterministic) functional of W v k . Denote by I ν n ,Ĩ ν n the random variables corresponding to I n ,Ĩ n with η n k (v n ), η n k (ṽ n ) replaced by η n,ν k ,η n,ν k , respectively, and note that now Combining (6.60), (6.61), (6.62), (6.57) we arrive at (6.58).
Lemma 6.9 shows thatŨ n ,W n ,Ũ n 0 satisfy the Faedo-Galerkin equations (5.17); hence, they are worthy of being referred to as "approximations". The next two lemmas summarize the relevant convergence properties satisfied by these approximations. such that as n → ∞, passing to a subsequence if necessary, Proof. The claims in (6.65) follow from the estimates in (6.21) and the sequential Banach-Alaoglu theorem. The relatioñ v i =ũ i −ũ e , dP × dt × dx a.e. inD × (0, T ) × Ω, is a consequence of (6.56) and the weak convergences in L 2 ω,t,x ofṽ n ,ũ n i ,ũ n e . The limit functionsũ i ,ũ e ,ṽ,w are easily identified with the a.s. representations in Lemma 6.6.
Proof. The proof merges the a.s. convergences in (6.46), the high-order moment estimates in (6.51), and Vitali's convergence theorem. To justify the first claim in (6.66), for example, we consider the estimatẽ see (6.51). From this we infer the equi-integrability (w.r.t.P ) of .
Accordingly, the first claim in (6.66) follows from theP -a.s. convergence in (6.46) and Vitali's convergence theorem, with the remaining claims following along similar lines.
Regarding the third claim, note also that forW n =W v,n orW w,n , which follows from equality of the laws and a martingale inequality.
For each n ≥ 1,W v,n andW w,n are (independent) cylindrical Wiener processes with respect to the stochastic basisS n , see (6.53). SinceW v,n →W v ,W w,n →W w in the sense of (6.46) or (6.66), it is more or less obvious that also the limit processesW v ,W w are cylindrical Wiener processes. Indeed, we have Proof. The proof is standard, see e.g. [43,Lemma 9.9] Consequently, to conclude the proof of the lemma, we need to verify that WithΦ defined in (6.48), it is sufficient to show that for all bounded continuous functionals L s (Φ) on X depending only on the values of Φ restricted to [0, s]. Since the laws of Φ n andΦ n coincide, cf. (6.48), where the last equality is a result of the {F n t }-martingale property of W n = W v,n , W w,n . By (6.46), (6.67), and Vitali's convergence theorem, we can pass to the limit in (6.68) as n → ∞. This concludes the proof of the lemma.
By cause of the first part of (6.66), passing to a subsequence if necessary, we may assume that as n → ∞, v n →ṽ for dP × dt × dx almost every (ω, t, x) ∈D × [0, T ] × Ω.

UNIQUENESS OF WEAK (PATHWISE) SOLUTIONS
In this section we prove an L 2 stability estimate and consequently a pathwise uniqueness result. This result is used in the next section to conclude the existence of a unique weak solution to the stochastic bidomain model. Let S, u i , u e , v, w be a weak solution according to Definition 4.1. We need a special case of the infinite dimensional version of Itô's formula [16,35,44]: To compute the first term on the right-hand side, multiply the first equation in (2.5) by u i , the second equation by −u e , and sum the resulting equations. The outcome is v dv − j=i,e ∇ · M j ∇u j u j dt + vI(v, w) dt = vη(v) dW v . k≥1 Ω w σ k (v) dx dW w k .

(7.2)
Remark 7.1. Krylov [35] provides a rather simple proof of the Itô formula for the squared norm u t 2 H , in the context of Hilbert space valued processes u t relative to a Gelfand triple V ⊂ H ⊂ V * , see also [44,Theorem 4.2.5]. As part of the proof he also establishes the continuity of the process t → u t 2 H (more precisely, for a modification of u t ). The idea of his proof is to use an appropriate approximation procedure to lift the differential du t into H, apply the Itô formula for Hilbert space valued processes, and then pass to the limit in the approximations. Note that Definition 4.1 asks that the paths of v(t) are weakly time continuous but not that they belong to C([0, T ]; L 2 (Ω)). However, as alluded to, this would actually be a consequence of the results in [35].
We are now in a position to prove the stability result.  Given the equations in (7.4), we apply the Itô formula, cf. (7.1) and (7.2). Using (2.6) and adding the resulting equations, we obtain in the end the following inequality:  The stochastic integrals in (7.5) are square-integrable, zero mean martingales. Moreover, using the Poincaré inequality, we have t 0 Ω |u e | 2 dx ds ≤C t 0 Ω |∇u e | 2 dx ds, for some constantC > 0. Since u i = v + u e , we control u i as well. As a result of all this, there is a constant C > 0 such that for any t ∈ [0, T ]. The Grönwall inequality delivers (7.3) "without sup". The refinement (7.3) ("with sup") comes from the application of a martingale inequality (3.6), see (6.17) for a similar situation.

EXISTENCE OF WEAK (PATHWISE) SOLUTION
In this section we prove that the stochastic bidomain model possesses a unique weak (pathwise) solution in the sense of Definition 4.5, thereby proving Theorem 4.6. The proof follows the traditional Yamada-Watanabe approach (see for example [18,24,31,37,44]), combining the existence of at least one weak martingale solution (Theorem 4.4) with a pathwise uniqueness result (Theorem 7.2), relying on the Gyöngy-Krylov characterization of convergence in probability (Lemma 3.4).
Referring to Section 5 for details, recall that U n = u n i , u n e , v n , w n , W n = W v,n , W w,n , U n 0 = u n i,0 , u n e,0 , v n 0 , w n where W v = {W v k } k≥1 , W w = {W w k } k≥1 are cylindrical Wiener processes. Moreover, recall that L n is the probability law of the random variable Φ n : D → X = X U × X W × X U0 , Φ n (ω) = (U n (ω), W n (ω), U n 0 (ω)) . We intend to show that the approximate solutions U n converge in probability (in X U ) to a random variable U = (u i , u e , v, w) (defined on S). Passing to a subsequence if necessary, we may as well replace convergence in probability by a.s. convergence. We then argue as in Subsection 6.4 to arrive at the conclusion that the limit U of {U n } n≥1 is a weak (pathwise) solution of the stochastic bidomain model.
It remains to prove that {U n } n≥1 converges in probability. To this end, we will use the Gyöngy-Krylov lemma along with pathwise uniqueness. By Lemma 6.5, the sequence {L n } n≥1 is tight on X . For n, m ≥ 1, denote by L n,m the law of the random variable Φ n,m (ω) = (U n (ω), U m (ω), W n (ω), U n 0 (ω), U m 0 (ω)) , defined on the extended path space Clearly, we have L n,m = L U n × L U m × L W n × L U n 0 × L U m 0 (see Subsection 6.3 for the notation). With obvious modifications to the proof of Lemma 6.5, we conclude that {L n,m } n,m≥1 is tight on X E . Let us now fix an arbitrary subsequence {L n k ,m k } k≥1 of {L n,m } n,m≥1 , which obviously is also tight on X E .
Passing to a further subsequence if needed (without relabelling as usual), the Skorokhod-Jakubowski representation theorem provides a new probability space (D,F ,P ) and new X E -valued random variables (8.1) Ū n k ,Û m k ,W n k ,Ū n k 0 ,Û m k 0 , Ū ,Û ,W ,Ū 0 ,Û 0 on (D,F ,P ), such that the law of Ū n k ,Û m k ,W n k ,Ū n k 0 ,Û m k 0 is L n k ,m k and as k → ∞, Ū n k ,Û m k ,W n k ,Ū n k 0 ,Û m k 0 → Ū ,Û ,W ,Ū 0 ,Û 0 P -almost surely (in X E ).
Indeed, U n k 0 = Π n k U 0 and U m k 0 = Π m k U 0 , and so for any ℓ ≤ min(n k , m k ), Denote by µ n k ,m k and µ the joint laws of Ū n k ,Û m k and Ū ,Û , respectively. Then, in view of (8.1), µ n k ,m k ⇀ µ as k → ∞. SinceŪ 0 =Û 0P -a.s., Theorem 7.2 ensures that U =ÛP -a.s. (in X U ). Hence, since this implies we can appeal to Lemma 3.4, cf. Remark 3.5, to complete the proof.