Stochastic electromechanical bidomain model

We analyse a system of nonlinear stochastic partial differential equations (SPDEs) of mixed elliptic-parabolic type that models the propagation of electric signals and their effect on the deformation of cardiac tissue. The system governs the dynamics of ionic quantities, intra and extra-cellular potentials, and linearised elasticity equations. We introduce a framework called the active strain decomposition, which factors the material gradient of deformation into an active (electrophysiology-dependent) part and an elastic (passive) part, to capture the coupling between muscle contraction, biochemical reactions, and electric activity. Under the assumption of linearised elastic behaviour and a truncation of the nonlinear diffusivities, we propose a stochastic electromechanical bidomain model, and establish the existence of weak solutions for this model. To prove existence through the convergence of approximate solutions, we employ a stochastic compactness method in tandem with an auxiliary non-degenerate system and the Faedo–Galerkin method. We utilise a stochastic adaptation of de Rham’s theorem to deduce the weak convergence of the pressure approximations.

Our main interest lies in the mathematical investigation of the interplay between propagation of electrical potentials in cardiac tissue (under stochastic effects) and its elastic mechanical response.These distinct processes in cardiac function are intricately linked by several complex phenomena occurring at various spatio-temporal scales.At the macroscopic level, the bidomain equations describe the propagation of electrical potentials in the heart by conserving electrical fluxes between the extra-and intra-cellular domains, which are separated by a membrane functioning as a capacitor.The bidomain equations are influenced by differing conductivities within these domains, which also vary based on the orientation of the cardiac tissue fibers.A multi-continuum description of the heart can be obtained by homogenization arguments, wherein both constituents (intra-and extracellular potentials) coexist at each material point.Several references, such as [21,52], describe the bidomain equations, whereas [7,23,29,41,45] provide insight into homogenization arguments for the bidomain model.
At the macroscopic level, muscle deformation can be described using the equations of motion for a hyperelastic material in its reference configuration.The muscle tissue itself is active, meaning that it can contract without external loads, but instead through intrinsic mechanisms that occur primarily at the microscale.There are different approaches to incorporating these effects, but one common method is to decompose stresses into active and passive components, resulting in what is known as the active stress formulation.This approach has been widely applied in various fields, as seen in studies such as [27,40,54].
Another approach to incorporating the intrinsic mechanisms in muscle contraction is to use the active strain formulation [19,39].This involves factoring the deformation gradient into active and passive components, with the active deformation representing the fiber contraction driven by the depolarization of the cardiomyocytes.This approach directly incorporates micro-level information on fiber contraction and direction into the kinematics, without the need for intermediate transcription in terms of stress.In this formulation, the strain energy function depends on auxiliary internal state variables that represent the level of tissue activation across scales [49].We adopt the active strain approach in our work, but [2,50] provide comparisons between the two (active stress / active strain) approaches in terms of numerical implementation, constitutive issues, and stability.
The analysis of deterministic macroscopic cardiac models has been primarily focused on the study of solutions to the bidomain equations coupled with phenomenological or physiologically-based ionic models.A variational approach was first introduced in [22] and has since been extended in various directions.These extensions include the nonlinear, degenerate conductivity tensors [4], more complex ionic models [12,56], and the existence of global classical solutions [36].Additionally, further references are available in [3,11,21,26,34].The aforementioned works provide well-posedness results for the bidomain model, utilizing differing solution concepts and technical frameworks.
Existence theories for nonlinear elasticity can be found in [20], while applications of these theories to the particular case of hyperelastic materials and cardiac mechanics (and their discretizations) are explored in [40].Despite the extensive literature on numerical methods and models for cardiac electromechanics, including [27,54], rigorous studies on the solvability and stability of solutions are still lacking.
The aim of this paper is to explore stochastic electromechanical coupling by considering the dependence of the electrical potentials, governed by the stochastic bidomain equations of [5], on the deformation gradient.This dependence arises after transforming from Eulerian to Lagrangian coordinates and by virtue of the Piola identity.In turn, the active part of the deformation is assumed to incorporate the influence of calcium kinetics into the equations of linear structural mechanics.This modeling approach enables coupling in the opposite direction (feedback), achieving stochastic electromechanical coupling.
This work improves the stochastic bidomain model for cardiac electrophysiology [5] by incorporating the mechanical activity of the heart.By doing so, the resulting model gains complexity in its mathematical analysis and becomes more realistic and applicable from a practical standpoint.To achieve a realistic model, it is crucial to blend microscale effects into the model that capture the heart's electrical and mechanical behaviors.Our approach incorporates "current noise" into the bidomain equations to model the stochastic behavior of ion channels on cellular voltage dynamics.Furthermore, we include "subunit noise" into the equations governing the gating variables, capturing the stochastic transitions of ion channels between different states [28].The model considers the mechanical contractions of the heart, which has a feedback effect on the electrical activity by changing the electrical conductivity due to tissue strain [19].This necessitates coupling the stochastic bidomain equations with the elasticity equations, showcasing the complex interaction between the electrical signals of the heart and its mechanical actions.
Realistic models play a crucial role in bridging the gap between organ-level phenomena and their biophysical underpinnings.They also furnish clinicians with quantitative insights for refining diagnostic, prognostic, and therapeutic strategies [47,9,42].A case in point is electrocardiography, where the traditional analysis of ECG (electrocardiogram) and BSPM (body surface potential mapping) data, based on static cardiac models, ignores the heart's dynamic contractions and the influence of random effects.This limitation can lead to inaccuracies in interpreting electrical signals.Incorporating heart contractions and stochastic factors into our model, as [37] suggests, allows for a more dynamic and accurate mapping of electrical potentials, enhancing our understanding of cardiac functions beyond the constraints of static models.
We will revisit the stochastic bidomain model later in this introduction, but for now, let us focus on the elasticity problem.
1.2.Nonlinear elasticity problem.We model the heart as a homogeneous continuous material occupying a bounded domain O R ⊂ R d with Lipschitz continuous boundary ∂O R , where d = 3.The deformation of the heart is described by the equations of motion written in the reference configuration O R , and the current configuration is denoted by O.The goal is to determine the deformation field φ : O R → R d , which maps a material particle initially located at position X to its current position x := φ(X).The tensor gradient of deformation is given by F := D X φ, where D X denotes the gradient operator with respect to the material coordinates X, noting that det(F) > 0.
We assume that the cardiac tissue is a hyperelastic incompressible material, and as such, there exists a strain stored energy function W = W(X, F), which is differentiable with respect to F. This function allows us to obtain constitutive relations between strains and stresses.Notably, we assume the incompressibility of the material, which is enforced via a scalar Lagrange multiplier p associated with the incompressibility constraint, interpreted as "hydrostatic pressure".By assuming this constraint, we minimize the total elastic energy subject to the condition det(F) = 1.Furthermore, the first Piola stress tensor P, which represents force per unit undeformed surface, is given by: where Cof (•) is the cofactor matrix.
We seek to find the deformation field φ and Lagrange multiplier p that satisfy the balance equations for deformations and pressure given by (1.1) where g is a prescribed body force, and n represents the unit outward normal vector to ∂O R .We complete these equations with the Robin boundary condition where α > 0 is a constant parameter.Our choice of boundary condition in (1.2) can be adjusted to replicate the global motion of the cardiac muscle [49].
To obtain a precise form of the equation in (1.1), we require a particular constitutive relation defining W. Our study is limited to Neo-Hookean materials, where which gives ∂W ∂F = µF and P = µF − p Cof (F), where µ is an elastic modulus.While simplified, this description of the passive response of the muscle exhibits a nonlinear strain-stress relationship arising from the incompressibility constraint.Another form of strain-stress non-linearity results from the anisotropy inherited from the active strain incorporation (see, e.g., [40,49,50] for more information).
1.3.Stochastic bidomain model.The deterministic bidomain equations, first introduced in [55] and further discussed in the books [21,52], are derived from the application of Ohm's electrical conduction law and the continuity equation, which ensures the conservation of electrical charge, to both the intracellular and extracellular domains.Recently, a stochastic extension of the bidomain equations has been proposed in [5].In this work, we will refer to these equations as the stochastic bidomain model (or equations).
Let T > 0 be a fixed final time, and let O ⊂ R 3 be a bounded open subset that represents the heart.We make the assumption that the intracellular and extracellular current densities can be expressed in terms of potentials v i = v i (t, x) and v e = v e (t, x), respectively, at points (t, x) ∈ O T := (0, T ) × O.We define the transmembrane potential and the gating or recovery variable as v = v(t, x) := v i − v e and w = w(t, x), respectively.
In the global coordinate system, cardiac electrical conductivity is represented by the orthotropic tensors ) for k ∈ {e, i} and s ∈ {l, t, n}, denote the intra-and extracellular conductivities along the directions of the fibers.Here, d s = d s (x) for s ∈ {l, t, n}, represents the direction of the fibers, which is a local quantity.
The stochastic bidomain equations [5] take the following form: where c m and χ are model parameters, and the externally applied stimulation current is represented by the function I app .In (1.3), W κ is a cylindrical Wiener process with a noise amplitude β κ for κ = v, w.The multiplicative stochastic term β v (v) dW v added to the equations for the membrane potential v is commonly referred to as current noise, and it represents the aggregated effect of the random activity of ion channels on the voltage dynamics.On the other hand, the gating variable w represents the fraction of open channel subunits of varying types, and the multiplicative noise β w (w) dW w added to the w-equation is sometimes referred to as subunit noise.For a discussion of different ways of adding noise to the Hodgkin-Huxley equations, we refer readers to [28].
We complete the SPDE system (1.3) with homogeneous Neumann boundary conditions for all fields.It is worth noting that the choice of the membrane model to be used is reflected in the functions H(v, w) and I ion (v, w).For a phenomenological description of the action potential, the FitzHugh-Nagumo model [38,25], given by assumption (A.4) below, will be used this paper.
The available literature on the mathematical theory of the stochastic bidomain system (1.3) is currently limited.The work [5] confirmed the global well-posedness of weak solutions for multiplicative noise.For strong well-posedness results, assuming additive noise, we refer to the works [30] and [33].
1.4.Stochastic electromechanical equations.The active strain model [19] is used to couple the electrical and mechanical behaviors.This model factorizes the deformation (gradient) tensor F into a passive component F p at the tissue level and an active component F a at the cellular level such that F = F p F a .The tensor F p accounts for the deformation of the material to ensure compatibility and possible tension due to external loads.The active tensor F a represents the distortion that dictates deformation at the fiber level and is assumed to depend directly on the electrophysiology, as described in [1]: where γ κ , κ ∈ {l, t, n} are quantities that depend on the electrophysiology equations.
The active strain approach used for decomposing the deformation tensor F assumes the existence of a virtual intermediate configuration located between the reference and the current frames.In this configuration, the strain energy function depends only on the macroscale deformation F p , as expressed in [1]: Moreover, the strain energy can be rewritten such that the Piola stress tensor becomes , see [1].It should be noted that mechanical activation is primarily influenced by the release of intracellular calcium [43].Furthermore, local strain dynamics closely follow calcium release rather than transmembrane potential.Physiologically-based activation models require a dependence of γ κ not only on calcium, but also on local stretch, local stretch rate, sliding velocity of crossbridges, and other force-length experimental relationships [49].However, for simplicity, we will restrict ourselves to a phenomenological description of local activation in terms of gating variables, where the scalar fields γ l , γ t , and γ n are functions of a common scalar γ as expressed in Here, each γ κ : R → [−Γ κ , 0] is a Lipschitz continuous monotone function and Γ κ is small enough to ensure that det(F a ) stays uniformly bounded away from zero for γ ∈ R and κ ∈ {l, t, n}.More explicitly, the functions γ κ are assumed to have the form: where γ + := max(0, γ) and γ R is a reference value.The scalar field γ is governed by the ordinary differential equation where β, η 1 , η 2 ≥ 0 are physiological parameters.
In order to couple mechanical and electrical behavior, the stochastic bidomain equations are modified by changing variables from the current configuration (Eulerian coordinates) to the reference configuration (Lagrangian coordinates).This modification introduces a conduction term that depends on the deformation gradient F. As a result, the stochastic electromechanical activity in the heart is governed by the following system of coupled random elliptic equations and stochastic degenerate parabolic reaction-diffusion equations, along with two differential equations: (1.4) where C −1 a is the tensor obtained from the active factor of the deformation gradient.To complete the system of equations (1.4), we need to specify suitable initial conditions for v, w and γ, as well as boundary conditions for v i , v e and the elastic flux a(•, •, •, •).
1.5.Linearizing the mechanical equations.To facilitate the mathematical analysis of the stochastic electromechanical system (1.4), we employ a linearization technique for both the incompressibility condition det(F) = 1 and the flux in the equilibrium equation.
Specifically, we use the following approach for the determinant: Since det(F) = 1, we can use the approximations We can express the displacement u as the difference between the current position ϕ(X) and the reference position X.With this notation, the linearized incompressibility condition takes the form ∇ • u = 0. To linearize the flux in equation (1.5) with respect to F, we employ a Taylor expansion of Cof(F) around the identity matrix I: As a result, we arrive at a(x, γ, F, p) = µFC −1 a (x, γ) − pI.We can use the displacement gradient ∇u to write the first equation in (1.4) as ).Therefore, we can reformulate the last equation to obtain a Stokes-like equation of the form where along with the algebraic relation The system (1.7) is fully specified by the inclusion of boundary conditions, which includes the linearization of the condition (1.2): for some α > 0. By selecting appropriate boundary conditions (1.8), it becomes feasible to impose the compatibility constraint (1.11) mentioned below, among other benefits.
We also impose the following initial conditions: The following list of assumptions made in the model described by (1.7), (1.8), and (1.9) play a crucial role in the subsequent analysis: (A.1) σ(x, γ) is a symmetric, uniformly bounded, and positive definite tensor: Furthermore, the function σ(•, •) is in the class C 1 ( Ō × R).Meanwhile, the function g, see (1.6), is an element of (L 2 (O T )) 3 .
(A.4) The functions H, I ion are given by the generalized FitzHugh-Nagumo kinetics: where for some positive constants c I,1 , c (1.10) The following combatibility condition holds (1.11) O v e (t, x) dx = 0 for a.e.t ∈ (0, T ).
To present our main result, we need to establish some additional assumptions about the stochastic data of the system (1.7).Specifically, we will assume the existence of a complete probability space (Ω, F, P ), which is equipped with a complete right-continuous filtration {F t } t∈[0,T ] .Consider a sequence {W k } ∞ k=1 of independent finite-dimensional Brownian motions adapted to the filtration {F t } t∈[0,T ] .We will refer to the collection of all these objects as the (Brownian) stochastic basis, denoted by All relevant processes are defined on a stochastic basis S and they are assumed to be measurable (predictable) with respect to the filtration {F t } t∈[0,T ] .For those who need background information on stochastic analysis and SPDEs, including stochastic integrals, Itô's chain rule, and martingale inequalities, we recommend consulting [24,46].To learn about key concepts related to probability measures (on topological spaces), weak compactness and tightness, we suggest reading [10].For a primer on quasi-Polish spaces and the Skorokhod-Jakubowski representation theorem [32], readers can consult [16,17,44].For properties of Bochner spaces, such as L r (Ω; X) = L r Ω, F, P ; X , where X is a Banach space and r ∈ [1, ∞], we suggest consulting [31].
Fixing a Hilbert space U equipped with a complete orthonormal basis {ψ k } k≥1 , we always consider cylindrical Brownian motions W evolving over U. Specifically, we define W := k≥1 W k ψ k , where {W k } k≥1 are independent standard Brownian motions on R.
The vector space of all bounded linear operators from U to a separable Hilbert space X with inner product (•, •) X and norm |•| X is denoted by L(U, X).In particular, we are interested in the collection of Hilbert-Schmidt operators from U to X, denoted by L 2 (U, X).
For the stochastic electromechanical bidomain system (1.7), a natural choice for X is the separable Hilbert space L 2 (O): To ensure that the infinite series k≥1 W k ψ k converges, we can construct an auxiliary Hilbert space U 0 that contains U, and a Hilbert-Schmidt embedding J from U to U 0 .With this construction, the series converges in the space L 2 Ω; C([0, T ]; L 2 (U 0 , L 2 (O))) [24].
Considering a cylindrical Wiener process W with noise coefficient β, the stochastic integral β dW is defined in the sense of Itô as follows [24]: for any predictable integrand β satisfying The stochastic integral This explains how we interpret the stochastic integrals in (1.7), with (β, . Next, we gather the conditions that need to be imposed on the noise coefficients. for a constant C β > 0. As a result of (1.12), one has The following theorem encapsulates our main global existence result for the stochastic electromechanical bidomain model.However, we will postpone the exact definition of the solution concept until Section 2.
The organization of the remaining part of the paper is as follows: First, Section 2 introduces the concept of a weak martingale solution.Then, in Section 3, we present the construction of approximate (Faedo-Galerkin) solutions for an auxiliary non-degenerate system.Next, Section 4 presents a series of a priori estimates.These estimates play a crucial role in establishing the convergence of the approximate solutions for the auxiliary problem.This convergence is tied to demonstrating the tightness of the probability laws of the Faedo-Galerkin solutions, as discussed in Section 5. Moving forward, Section 6 establishes the existence of a solution for the non-degenerate problem.Subsequently, in Section 7, we prove the existence of a solution for the original problem.This section focuses on demonstrating the weak convergence of the pressure variable, which does not exhibit a uniform L 2 bound.By interconnecting and building upon one another, these sections form a proof of the main result (Theorem 1.1).Finally, some numerical examples are presented in Section 8.

NOTION OF SOLUTION
To specify the initial data for the problem under consideration, it is important to keep in mind that we are considering a probabilistic weak notion of solution.For probabilistic strong (pathwise) solutions, we prescribe the initial data as random variables.Specifically, we require that v 0 ∈ L 2 (Ω; L 2 (O)) and w 0 , γ 0 ∈ L 2 (Ω; H 1 (O)).On the other hand, if we are interested in martingale or probabilistic weak solutions, we cannot rely on the knowledge of the underlying stochastic basis.Therefore, we must specify the initial data in terms of probability measures.More precisely, we prescribe the initial data as probability measures µ v0 on L 2 (O) and µ w0 , µ γ0 on H 1 (O).In this context, the measures µ v0 , µ w0 and µ γ0 should be viewed as the "initial laws" because the laws of v(0), w(0), and γ(0) will be required to coincide with µ v0 , µ w0 , and µ γ0 , respectively.
We will utilize the following notion of solution for our stochastic system (1.7).

APPROXIMATE SOLUTIONS
In this section, we define the Faedo-Galerkin approximations for the stochastic system (1.7).We construct approximate solutions relative to a stochastic basis To define the Faedo-Galerkin solutions, we will use classical Hilbert bases that are orthonormal in L 2 and orthogonal in H 1 , see for example [48, p. 138-145].More precisely, consider bases {ψ l } l∈N and {ζ l } l∈N such that Span {ψ l } is dense in L 2 (O) 3 and H 1 (O) 3 , while Span {ζ l } is dense in L 2 (O) and H 1 (O).In order to impose the compatibility condition (1.11), we introduce the modified functions We observe that Span {µ l } is dense in the space Moreover, we can apply the Gram-Schmidt process to orthonormalize the basis {µ l }.The resulting L 2 (O)-orthonormal basis is still denoted by {µ l }.Furthermore, Span {µ l } is dense in H 1,0 (O).For n ∈ N, we introduce the finite dimensional spaces The basis functions µ l and ζ l are conventionally derived from the eigenfunctions of the Neumann-Laplace operator.This condition ensures that ∂e l /∂n = 0 on ∂O, where e l represents either µ l or ζ l for each l ∈ N. Drawing from elliptic regularity theory, we recognize that every such eigenfunction e l is a member of H 2 (O).Furthermore, the degree of smoothness of e l within Ō is determined by the smoothness of ∂O.For instance, if ∂O is C ∞ , then e l ∈ C ∞ .For more details, see the summary in [6, p. 9-10] (for example).
If we were to apply the Faedo-Galerkin method directly to the highly degenerate system (1.7), we would end up with a set of coupled stochastic differential equations (SDEs) and algebraic equations that would need to be solved at each time.Unfortunately, the classical SDE theory does not provide a straightforward guarantee of the existence and uniqueness of solutions for this type of system.To overcome this challenge, we will use a time regularization approach inspired by [4].This approach considers an auxiliary non-degenerate system obtained by adding temporal regularizations to the first four equations in (1.7).Specifically, we add ε∂ t u ε , ε∂ t p ε , εdv i,ε , and −εdv e,ε , respectively, where ε > 0 is a small parameter that we later send to zero.
Fixing a small number ε > 0, we first look for a solution of the following non-degenerate system where In the system (3.2), the elasticity equation has been augmented to incorporate the dynamic velocity term ε∂ t u ε .Meanwhile, the addition of the term ε∂ t p ε in the second equation represents an "artificial compressibility".The system (1.7) also includes the boundary conditions (1.8).However, when it comes to specifying the initial data for U ε , we need to exercise extra caution since (1.9)only provides initial data for v ε , w ε , and γ ε .Therefore, we must also supply auxiliary initial data for u ε , p ε , v i,ε , and v e,ε .We will explore this issue in greater detail later on.
We can now move on to the next step, which involves applying the Faedo-Galerkin method to the non-degenerate system (3.2).This method consists of projecting (3.2) onto the finite-dimensional spaces in (3.1).To accomplish this, we consider approximate solutions that take the following form: The time-dependent coefficients u n l (t), p n l (t), v n i,l (t), v n e,l (t), w n l (t), γ n l (t) are determined by requiring that the following (Galerkin) equations hold for l = 1, . . ., n: Let us comment on the (finite dimensional) approximations of the stochastic integrals used in (3.4).In this context, let (β, W ) represent either (β v , W v ) or (β w , W w ).Let e l denote either ζ l or µ l , so that {e l } ∞ l=1 is an L 2 (O)-orthonormal basis.Recall that β maps between L 2 0, T ; L 2 (O) and L 2 0, T ; L 2 U, L 2 (O) , for fixed ω ∈ Ω, where U comes with the orthonormal basis {ψ k } ∞ k=1 .Utilizing the following decomposition where , we can write so that by orthonormality, for any l = 1, 2, . .., For related decompositions, see, e.g., [18].In the Faedo-Galerkin system (3.4),we employ the usual finite dimensional approximation Since the original problem (1.7), (1.9) lacks initial conditions for the unknowns u, p, v i , and v e , it becomes necessary to augment our non-degenerate system (3.2) with auxiliary initial conditions.To address this, we introduce the functions and O v e,0 dx = 0. Additionally, we choose u 0 ≡ 0 and select an arbitrary p 0 ∈ L 2 (O) to further supplement the initial conditions.Given the assumption (1.14), we conclude that with q 0 is defined in (1.14).
The initial data for the Faedo-Galerkin approximations (3.3) are now defined as follows: v e,0,l µ l , where v e,0,l := (v e,0 , µ l ) L 2 , (3.6) In the upcoming sections, we investigate the convergence properties of the sequences . These sequences are analyzed specifically for any fixed ε > 0. Currently, our assertion is focused on the existence of a (pathwise) solution to the SDE system defined by (3.4) and (3.6), for each fixed value of n ∈ N. Lemma 3.1 (well-posedness of the Faedo-Galerkin equations).For any fixed n ∈ N, the regularized Faedo-Galerkin equations, specifically the equations (3.4) together with the initial values specified in (3.6), possess a unique strong solution Moreover, almost surely, the processes v n i , w n , γ n , p n are members of C([0, T ], W n ), while v n e belongs to C([0, T ], L n ), and u n is in C([0, T ], H n ).Additionally, for each time t ∈ [0, T ], these processes are square-integrable with respect to the probability variable.
Proof.By the orthonormality of the bases and incorporating the algebraic relation and The third and fourth equations in the system (3.7) can be expressed equivalently in matrix form as follows: Here, A = {a lr } with a lr = O ζ l µ r dx, I n represents the identity matrix, and b is the vector containing the right-hand side terms corresponding to the equations in (3.7).
In order to express this system in the form we need to verify that the matrix M is invertible.To do so, let us expand M as follows: It is enough to show that the matrix Thus, the matrix M is positive definite, and as a result, it is also invertible.Let Y n (t) represent the vector that encompasses all variables of interest, including the sequences , and {γ n l (t)} n l=1 }.Furthermore, define F (Y n ) as the vector that aggregates all the deterministic drift terms derived from the right-hand side of the system (3.7), and let G(Y n ) represent the collection of all corresponding noise coefficients.Consequently, the system under consideration can be succinctly expressed as the SDE system , where Y n 0 denotes the vector of initial conditions as specified by (3.6).If the vector-functions F and G satisfy global Lipschitz continuity conditions, classical theory on SDEs guarantees the existence and uniqueness of a strong solution.However, the inherent nonlinearity of the ionic models, as detailed in (A.4), implies that these global Lipschitz conditions are not met for the SDE system (3.8).
Nevertheless, the weak monotonicity framework of [46, Theorem 3.1.1]is applicable.In fact, the verification of the following two conditions will imply the existence of a global strong solution to (3.8): for some r-dependent positive constant K r .
• (weak coercivity) Considering the assumptions (A.1)-(A.4),especially with the utilization of (1.10), the demonstration of both the local weak monotonicity and the weak coercivity conditions follows the approach detailed in [5, p. 5328].We will not repeat the details here.Since these conditions are satisfied, for any F 0 -measurable mapping Y n 0 : Ω → R 6n , there exists a unique solution to the SDE system (3.8),unique up to P -indistinguishability.A solution Y is defined here as the collection Y = {Y (t)} t∈[0,T ] , which forms a P -almost surely continuous, R 6n -valued, F t -adapted process.For every t ∈ [0, T ], P -almost surely, Furthermore, for all t ∈ [0, T ], for a T -dependent constant C T .We refer to [46, Theorem 3.1.1]for further details.□ Let us consider a basis {e l } ∞ l=1 selected from those previously introduced in (3.1).We then introduce the projection operator, as detailed, for example, in [14, page 1636]: The following properties hold: and We have three projection operators, Π H n , Π L n and Π W n , which correspond to the finitedimensional spaces H n , L n , W n defined in (3.1).In order to streamline our notation, we will often omit the specific dependency of Π n on its associated space.For instance, we might write Π n in lieu of Π H n .The specific space in question should then be inferred from the context and the variable being projected.
Utilizing the projection operators Π n , we can express the Faedo-Galerkin equations from (3.4) in an integrated form.This representation allows us to present these equations as equalities between random variables in the dual space of H 1 : where ).In the same vein, z n 0 stands for Π n z 0 , where z encompasses the variables p, v i , v e , w 0 , and γ 0 .Moreover, the term v n 0 is defined as the difference between v n i,0 and v n e,0 .

A PRIORI ESTIMATES
We will prove convergence of the Faedo-Galerkin approximations using the stochastic compactness method, which asks for the uniform tightness of the probability laws linked to these approximations.Tightness will hinge upon a series of a priori estimates, which hold uniformly in n.In our first set of estimates, we harness the structure of the equations (3.2), along with stochastic calculus, to deduce several basic energy-type estimates.Then we utilize these estimates and the governing equations to derive temporal translation estimates in the space of L 2 t L 2 x .When brought together, these estimate will facilitate the tightness of the laws of the Faedo-Galerkin approximations.
The non-degenerate system (3.2) consists of a Stokes-like component for (u ε , p ε ), a pair of SPDEs for (v i,ε , v e,ε ) (and v ε = v i,ε − v e,ε ), a SDE for w ε , and an ODE for γ ε .The derivation of energy balance equations for from which we get the balance equation for ∥w ε (t)∥ 2 L 2 .It is less straightforward to deduce the stochastic balance relation for Following [5], we will present a formal argument that will lead to the energy balance for these variables.From (3.2), we deduce that and By considering the aforementioned pair of SPDEs along with the third and fourth SPDEs stated in (3.2), we can employ the Itô product rule to deduce the following result: which is the sought-after energy balance.Following this, we will now transition into providing rigourous estimates for all the variables at the discrete level for the Faedo-Galerkin approximations.Lemma 4.1.Assuming that conditions (A.1) to (A.6) and (3.5) are satisfied, consider the following variables: which fulfill the Faedo-Galerkin equations given by (3.4), with initial data from (3.6).There exists a constant C > 0, independent of n and ε, such that Moreover, Finally, the following gradient estimates hold: Proof.The formal relationship previously indicated in equation (4.2) can be rigorously derived for the Faedo-Galerkin approximations.This derivation is closely related to the details found in [5].Consequently, this implies that: where In a similar manner, since (4.1) can also be rigorously derived for the Faedo-Galerkin approximations, we obtain By adding (4.10) and (4.11), performing a time integration, and making use of the ellipticity assumption (A.2), we reach the following result: for some constants C 1 > 0 and C 2 , C 3 ≥ 0. Combining (4.13) with the assumptions (A.4) and (A.6) in (4.12), we obtain that Note that the martingale property of stochastic integrals guarantees that the expected value of each of the last two terms in (4.14) is zero.Consequently, by taking the expectation in (4.14) and applying Grönwall's inequality, and that the sought-after estimates (4.4) hold.
The following equation holds pointwise in t, for all test functions ξ n ∈ W n : By utilizing this equation in conjunction with (A.3), we can infer that for some constant K3 independent of n, ε.We apply Grönwall's inequality in (4.16), take the supremum over t ∈ [0, T ] and subsequently take the expectation.We combine this with the L 2 estimate (4.15) for w n , which allows to conclude (4.6).
Let E denote the left-hand side of (4.3).To establish the estimate (4.3), we proceed by taking the supremum over the interval [0, T ] in (4.14) and then the expectation E[•].Utilizing (4.15), the Burkholder-Davis-Gundy inequality, the Cauchy-Schwarz inequality, the assumption (1.13) regarding β v , Cauchy's product inequality, and once again (4.15), we obtain for any (small) δ > 0, where the constants K 4 , . . ., K 8 are independent of n, ε.This establishes (4.3).Next, we shall establish the validity of (4.5).Using (1.11) and applying the Poincaré inequality, we can find a constant C > 0 that depends on O but is independent of n, ε, ω, and t, such that the following inequality holds for (ω, t) ∈ Ω × (0, T ): . Consequently, utilizing (4.4), we deduce that (4.18) Since v n complies with (4.15), it follows that also v n i satisfies (4.18).Thus, (4.5) holds.Next, let us establish the gradient estimates (4.9) for γ n and w n .We can make the assumption, without loss of generality, that the basis functions of W n , see (3.1), are C 2 regular and have a zero normal derivative on the boundary ∂O.This implies that ∇γ n • n is equal to zero on ∂O.Consider that γ n resides within the finite-dimensional space W n , which is spanned by the eigenfunctions {ζ l } of the Neumann-Laplace operator, with eigenvalues {λ l }.We can compute the Laplacian of γ n : Hence, ∆γ n ∈ W n .So, by considering ξ n := ∆γ n ∈ W n (t is treated as a parameter) in the equation (4.20) for the activation variable γ n , we can perform spatial integration by parts and apply the temporal chain rule.This allows us to derive the following equation: by an appropriate use of Young's product inequality.Thus, to establish the proof of (4.9) for both γ n and w n , we only need to establish a gradient bound specifically for w n .Given the SDE (3.4) for w n , we deduce that for all test functions ξ n ∈ W n .Following a similar line of reasoning as presented above for ∇γ n , we can now proceed to derive where the first term on the right-hand side can be controlled by (4.4) (keeping in mind that We can now proceed with a similar approach as employed in establishing the L 2 estimates (4.3).This involves utilizing the BDG and Gronwall inequalities.Let us here only discuss the application of the BDG inequality, following (4.17).By (1.12), we can observe that This, when combined with the BDG inequality, yields the following result: By utilizing this inequality, we can proceed in a similar manner as before, see (4.17).This approach leads us to the first part of (4.9) for ∇w n , which, via (4.19), implies the second part for ∇γ n .Finally, we proceed with the derivation of the estimates for the Stokes-like part of the system (3.4).Consider the following system, which holds pointwise in time t, for all test functions ψ n ∈ H n and ρ n ∈ W n : where f n = f (x, γ n ), cf.(1.6).Now, we replace ψ n with u n and ρ n with p n in (4.20).We then sum the resulting equations to yield ε 2 Let us introduce the continuous bilinear form We claim that a(•, •) is coercive on (H 1 (O)) 3 .First, by the uniform ellipticity of σ, see assumption (A.1), Next, we want to show that there exists a constant ν > 0 such that We argue by contradiction.Suppose, for each , we obtain On the other hand, since v m is bounded in (H 1 (O)) 3 and O ⊂ R 3 is bounded and smooth, there exists v ∈ (H 1 (O)) 3 and a subsequence {v Utilizing (4.24)(i), we infer that ∇v = 0. Given that O is connected, this leads us to conclude that v = C. Reapplying (4.24)(i) and considering the convergence of 3 .This, in view of the continuity of the trace map, suggests that On the other hand, by (4.24)(ii), v m k → 0 in (L 2 (∂O)) 3 , and so C = 0.However, this is contradictory to the initial assertion given by the first part of (4.23).This affirms the claimed coercivity of a(•, •).
By the coercivity of a and Young's product inequality, we obtain from (4.21) that (4.25) for some constant C > 0 independent of n, ε.By (1.6) and the assumption (A.1) on σ, we have . By integrating (4.25) over the interval (0, t), where 0 < t ≤ T , and subsequently taking the supremum over t, we can then apply the expectation operator.Keeping in mind that u n (0) = 0, ∥p n 0 ∥ L 2 (O) ≤ ∥p 0 ∥ L 2 (O) , and using the estimate (4.9), we arrive at (4.7).
By integrating (4.25), we can also confirm that (4.8) is satisfied.□ The subsequent corollary, encompassing higher moments estimates, arises as a direct consequence of Lemma 4.1 and the BDG inequality.
Corollary 4.2.Along with the assumptions stated in Lemma 4.1, let us further assume that (3.5) holds, with q 0 defined as specified in (1.14).Then there exists a constant C > 0, which remains independent of both n and ε, satisfying the following estimates: Proof.We exponentiate both sides of the inequality (4.14) by q 0 /2, recalling that q 0 > 9/2.Subsequently, we employ a series of elementary inequalities to arrive at for some constant C 1 that is independent of n (but dependent on T ).For any R > 0, we define the stopping time Following (4.17), we apply the BDG inequality along with several elementary inequalities.This enables us to derive the following two bounds for any δ > 0: for certain constants C 2 , C 3 that are independent of n (though they depend upon T ).By combining this with (4.28), applying Gronwall's inequality, and subsequently taking the limit as R → ∞ through the application of the monotone convergence theorem, we arrive at the initial three estimates presented in (4.26).
After establishing these estimates, we return to the inequality (4.14), which also features the two terms j=i,e t 0 O ∇v n j 2 dx ds and 4 dx ds.We deduce (4.27) by exponentiating this inequality with q 0 /2 and then applying the estimates we derived earlier.
The estimation (4.26) of γ n is derived in a similar manner from (4.16).The latter involves no stochastic integral, rendering the process comparatively simpler.Likewise, by employing the equations in (4.21), which are also devoid of stochastic integrals, we can establish the prescribed bounds (4.26) for the velocity u n and pressure p n .□ To ensure the tightness of the probability laws of the Faedo-Galerkin solutions, we rely on an Aubin-Lions-Simon compactness lemma.To apply this lemma effectively, we will require the following (n, ε uniform) temporal translation estimates.Lemma 4.3.Suppose conditions (A.1)-(A.5),(1.12) and (3.5) hold.Let v n i (t), v n e (t), v n (t), u n (t), p n (t), w n (t), γ n (t), t ∈ [0, T ], satisfy (3.4) , (3.6).There exists a constant C > 0, independent of n and ε, such that for any sufficiently small δ > 0, Proof.We refer to Lemma 6.3 in [5] for further details.The primary distinction from the study [5] is that the operators Π n ∇ • M j (x, ∇u n )∇v n j , where j = i, e, nonlinearly depend on the displacement gradient ∇u.Nevertheless, given the L ∞ -boundedness of M i,e , as stipulated in assumption (A.2), the reasoning proposed by [5] remains applicable to v n .The argument for the w n variable does not necessitate any alterations, and the added variable γ n can be addressed in a similar (simpler) manner to w n .□

TIGHTNESS AND A.S. REPRESENTATIONS
The objective is to establish the existence of a solution, drawing inspiration from some techniques previously applied to the deterministic system [4].Utilizing the Faedo-Galerkin approach, we have constructed approximate solutions originating from a nondegenerate auxiliary system and we wish to deduce their convergence by employing the stochastic compactness method.In contrast to the deterministic case, where the a priori estimates facilitate strong compactness across the time and space variables (t, x)-which is essential for addressing the nonlinear terms-the stochastic framework introduces the complexity of the probability variable ω, rendering strong compactness unattainable.The standard resolution involves deducing the weak compactness of the probability laws of the Faedo-Galerkin solutions through tightness and Prokhorov's theorem, followed by constructing almost surely convergent sequences via Skorokhod's representation theorem.This enables us to demonstrate that the limit variables form a weak martingale solution.Due to the degenerate nature of our model, we observe only weak (H 1 ) convergence for the intra-and extracellular potentials.As a result, the probability measures of these potentials live on a Banach space equipped with the weak topology, which is notably non-Polish.To handle this issue, we utilize Jakubowski's version [32] of the Skorokhod representation theorem, tailored for non-metric spaces.
Let us establish the required tightness as n → ∞ (with ε kept fixed) of the probability laws generated by the Faedo-Galerkin solutions V n , W n , V n 0 n∈N , where (5.1) Let us introduce the path space where , for V n 0 .
The subscript "w" signifies that the topology is weak.Let B(X ) symbolize the σ-algebra of Borel subsets of X .We define the (X , B(X ))valued measurable function Φ n on the measure space (Ω, F, P ) as follows: We utilize L n to represent the joint probability law of Φ n on (X , B(X )).Furthermore, the notation for the marginals of L n substitutes the subscript n with the variable being considered.For example, L v n stands for the law of the variable v n in the L 2 t,x space.This convention applies uniformly across all other variables.
The objective is to affirm the tightness of the joint laws {L n }.This involves producing, for every κ > 0, a compact set By invoking Tychonoff's theorem, this tightness is a direct consequence of the tightness of the product measures under the product σ-algebra.Thus, to substantiate (5.2), it suffices to identify a compact set i ,κ ≤ κ, for each κ > 0, and analogously for the remaining variables v n e , . . ., u n 0 .Indeed, we have the following lemma: Lemma 5.1.Equipped with the estimates in Lemmas 4.1 and 4.3, the sequence of laws {L n } n≥1 is tight on (X , B(X )), in the sense that (5.2) holds.

Proof. By weak compactness of bounded sets in
x w , for any R > 0. By Chebyshev's inequality, Hence, we can specify the required compact We can argue similarly (via weak compactness) for the variables v n e , p n , and u n .Given that W v,n a.s.converges in C([0, T ]; U 0 ) as n → ∞, it follows that the laws {L W v,n } n∈N undergo weak convergence.This implies the tightness of {L W v,n } n∈N .In other words, for any specified κ > 0, we can find a compact set K W v,n ,κ ⊂ C([0, T ]; U 0 ) such that the measure of the complement of this set with respect to L W v,n is ≤ κ.This reasoning can be analogously applied to W w,n .
Similar reasoning can be applied to the initial data variables v n i,0 , v n e,0 , v n 0 , w n 0 , γ n 0 , p n 0 , u n 0 .By assumption, these variables are a.s.convergent in L 2 x as n → ∞.Consequently, their corresponding laws also undergo weak convergence.This leads to the conclusion that the laws of these variables are tight on L 2 x .Our remaining task is to demonstrate the tightness of the laws of v n , w n , and γ n on L 2 t,x .This will be a direct consequence of the L 2 estimates given by (4.3) and (4.6), the L 2 gradient estimates specified in (4.7) and (4.9), and the temporal continuity estimates provided by Lemma 4.3, drawing upon the insights from [8] and [5].Focusing on v n (although the variables w n and γ n follow suit in an analogous manner), consider any two sequences of positive numbers r m and ν m , which converge to zero as m → ∞.We will proceed by defining the following set: Theorem 3].Fix two sequences r m , ν m of positive numbers numbers tending to zero as m → ∞ (independently of n), such that and define the L 2 t,x -compact set for a number R > 0 that will be determined later.We have .
By Chebyshev's inequality, As a result, we can specify the required compact We can argue similarly for w n and γ n .□ The primary component of the stochastic compactness method is the tightness of the joint laws of the random variables under consideration, as demonstrated in Lemma 5.1.The Skorokhod a.s.representation theorem enables the replacement of these random variables with versions that converge a.s., albeit defined on a new probability space.However, in our context, the direct application of the Skorokhod theorem is unfeasible due to the necessity to operate within non-metric spaces, specifically L 2 t,x w and L 2 t H 1 x w .As a result, we employ a more recent version of the Skorokhod theorem that is applicable to quasi-Polish spaces, due to Jakubowski [32].A quasi-Polish space is defined as a Hausdorff space X that enables a continuous injection into a Polish space.Notably, separable Banach spaces with weak topology are quasi-Polish.The Skorokhod-Jakubowski theorem and its application to SPDEs is further discussed in [14,15,44], showcasing some of the earliest uses of the theorem in this context.This theorem was also used in [5] for the stochastic bidomain model.For further references, this paper may be consulted.
By the Skorokhod-Jakubowski representation theorem, we infer the existence of a new probability space, denoted as Ω, F, P , and new sets of random variables (5.4) Φn = Ṽ n , W n , Ṽ n 0 , Φε = Ṽε , Wε , Ṽ0,ε , where {n = n j } j∈N represents a subsequence that converges to infinity as j → ∞.For ease of notation, we use n instead of n j to avoid relabelling the subsequence.The new random variables hold respective joint laws L n and L. Specifically, the law of Φn coincides with that of Φ n .Crucially, the following almost sure convergence is valid: Regarding the tilde notation, the components of Ṽ n , W n , Ṽ n 0 bear the same symbolic representations as V n , W n , V n 0 , except with a tilde symbol overhead.For instance, v n i would transform into ṽn i , and a similar rule applies to the other variables.The components of Ṽε , Wε , Ṽ0,ε are denoted in the same way.For example, the limit of ṽn i as n → ∞ would be symbolized as ṽi,ε (the parameter ε will be sent to zero in the next section).
For the sake of clarity and to facilitate future reference, the almost sure convergence (5.5) can be recast into the following strong convergences: ( Moreover, the following weak convergences hold: ( Leveraging the equality of laws, the a priori estimates in Lemma 4.1 and Corollary 4.2 continue to hold for the new random variables ṽn i , ṽn e , ṽn , wn , pn , u n , defined on Ω, F, P .Besides, we have ũn 0 ≡ 0 and ṽn = ṽn i − ṽn e a.e.In subsequent discussions, even when referencing estimates for variables with a tilde atop, we will consistently use the equation numbers originally associated with their non-tilde counterparts.
Moving forward, we will use the following stochastic basis associated with Φn , n ∈ N: where Fn t t is the smallest filtration making the processes (5.4) adapted.Note that, by the equality of laws and [24], W v,n and W w,n are cylindrical Brownian motions.We also need the n-truncated sums ).We still need to show that the tilde variables, Ṽ n , W n , Ṽ n 0 satisfy the Faedo-Galerkin equations on the newly defined stochastic basis (5.8).This is achieved by leveraging the equivalence of laws.Multiple strategies can be employed to flesh out the details, as seen in various sources like [8,15,44].In particular, [5] utilized the method presented in [8] to tackle the stochastic bidomain model.The same argument can be adopted here for (3.9).However, we will not delve into the nitty-gritty of the proof; instead, we direct the interested reader to [5] or any of the other cited references for a more in-depth exploration.
As a result, we can, for the subsequent discussions, make the assumption that the tilde variables Ṽ n , W n , Ṽ n 0 conform to the Faedo-Galerkin equations.This, symbolically, means that equations (3.9) are upheld when a tilde is placed over each variable in the equations.Moving forward, whenever we refer to the Faedo-Galerkin equations in the context of the tilde variables, we are indeed referencing the equations (3.9) pertaining to the original variables.This approach avoids the unnecessary repetition of equations, thus conserving space.

SOLUTION TO NON-DEGENERATE MODEL
This section aims to prove that the limiting functions in (5.6) and (5.7) form a weak martingale solution for the non-degenerate system (3.2).This solution also incorporates the boundary conditions given by (1.8) and the auxiliary initial data for the variables v i,ε , v e,ε , u ε , and p ε , which were discussed earlier.
In addition to (5.8), we need a stochastic basis for the a.s.limit Φε in (5.4), (5.5): where Fε t t is the smallest filtration that makes all the relevant processes in Φε adapted.Then W v ε , W w ε are (independent) cylindrical Wiener processes with respect to sequences W v ε,k k∈N , W w ε,k k∈N of mutually independent real-valued Wiener processes adapted to the filtration Fε t t , such that See [5] for further details and references.At this point, we have established the almost sure convergences in (5.6) and (5.7) for the tilde variables.In addition, the tilde variables adhere to the statistical estimates provided by Lemma 4.1 and Corollary 4.2.Considering the implications of the Banach-Alaoglu theorem along with Vitali's convergence theorem, after possibly thinning the sequence, we are justified in assuming that the following convergences occur as n → ∞: ṽn → ṽε , wn → wε , γn → γε a.e. and in L 2 Ω; L 2 ((0, T ) × O) , where the limit functions are bounded in the spaces suggested by the uniform estimates supplied by Lemma 4.1 and Corollary 4.2.Besides, With the convergences established above (6.2), the concluding step involves taking the limit in the Faedo-Galerkin equations (3.9) as n → ∞ (in this section, ε is kept fixed).The details of this step will be provided in the subsequent lemma.M i (x, ∇ũ ε )∇ṽ i,ε • ∇φ i dx ds (6.4) The laws of ṽε (0) = ṽ0,ε , wε (0) = w0,ε , γε (0) = γ0,ε are µ v0 , µ w0 , µ γ0 , respectively.Remark 6.2.Note that the (regularized) Stokes-type equations (6.3) hold a.e. in the time variable t.On the other hand, Definition 2.1 stipulates them to be valid only in the sense of distributions on (0, T ).Nevertheless, by weakly differentiating the a.e.formulation (6.3) with respect to time, we obtain Proof.The proof follows [5], while focusing on the specific modifications introduced by our extended model.These modifications primarily pertain to the analysis of the Stokeslike system (6.3), the necessity of taking the limit in the nonlinear conductivities within the bidomain component (6.4), and finally the treatment of the ODE for γ ε in (6.5).We start by establishing (6.3).Fix v ∈ [H 1 (O)] 3 and φ p ∈ H 1 (O), and note that the Faedo-Galerkin equations (3.9) imply that To state the desired equations (6.3) in a more concise manner, we can use the notations I v u (ω, t) = 0 and I φp p (ω, t) = 0, where I v u and I φp p are appropriate (integrable) functions on Ω × [0, T ] defined by (6.3).Let Z ⊂ Ω × [0, T ] be a measurable set, and denote by the characteristic function of Z.We will show that for any measurable set Z ⊂ Ω × [0, T ], and for a countable collections of test functions x .This is enough to conclude that I v u = 0 and I φp p = 0 for ( P × dt) a.e.(ω, t), for all H 1 test functions v, φ p , i.e., (6.3) holds.
To demonstrate this, we will employ the following approach: first, we multiply the Faedo-Galerkin equations (6.6) by 1 Z and integrate them over ω and t.Subsequently, we will take the limit as n → ∞.A similar line of reasoning will be applied to the components (6.4) and (6.5) of the system.
By the weak L 2 ω,t,x convergence of ũn to ũε , see (6.2) and the strong L 2 convergence of Π n v to v. we obtain Similarly, we can pass to the limit in the terms involving −p n ∇ • (Π n v), and (∇ • ũn )φ p .For the term associated with αũ n • (Π n v), the weak convergence of αu n on the boundary is attributed to the continuity of the trace map from [H 1 (O)] 3 to [L 2 (∂O)] 3 .Moreover, we will soon establish the strong convergence ũn → ũε in L 2 ω,t [H 1 x ] 3 as n → ∞.By using (6.2) once more, we can easily handle the initial data terms in (6.6).
In relation to the remaining term in the elasticity equation involving f n • v, where f n = f (t, x, γn ), we can refer back to (1.6) which, through (6.2), shows that f n weakly converges to fε = f (t, x, γε ) in L 2 ω,t,x .Simultaneously, we know that Π n v → v strongly in L 2 .Therefore, Let us now shift our attention to a new aspect of the analysis for the electrical (bidomain) component of the model.The focus will be on elucidating the process of taking the limit in the nonlinear conductivities, the remaining terms can be treated as in [5].This includes passing to the limit in the ODE (3.9) for the new variable γn .Notably, this step is simpler than handling the SDE for the variable wn , which was already adeptly addressed in [5].
To illustrate the new difficulty, we consider the following terms derived from the Faedo-Galerkin equations (3.9): where ∇(Π n φ j ) → ∇φ j strongly in L 2 and ∇ṽ n j → ∇ṽ j,ε weakly in L 2 ω,t,x , see (6.2).To facilitate the computation of the limit of J n j , it is crucial to establish the strong convergence, rather than weak convergence, of ∇ũ n in L 2 ω,t,x .This strong convergence is essential as ∇ũ n is a component of the nonlinear functions M j .Although the weak convergence is provided by (6.2), to make progress, we must prove strong convergence.To achieve this, we have two potential approaches.Firstly, we can employ the well-known Minty-Browder argument.Alternatively, we can leverage the structural characteristics of the purely mechanical equation along with the coercivity of the bilinear form a(•, •)introduced in (4.22)-to derive the following estimate: for some constant c > 0 independent of n (and ε).We claim that the three terms on the right-hand side converge to zero as n → ∞.The initial data term converges to zero as stated in (6.2).By recalling the definition of f given in (1.6), we can employ integration by parts.This approach effectively brings us to the consideration of the following strongweak product, disregarding the straightforward term that involves g: σ(x, γn ) − σ(x, γε ) : (∇ũ n − ∇ũ) .
This product term converges to zero-weakly in L 1 ω,t,x -by the weak L 2 ω,t,x convergence of ∇ũ n ε and the strong L 2 ω,t,x convergence of σ(x, γn ).To establish the latter, we have utilized the fact that σ(x, γn ) → σ(x, γε ) almost everywhere.This, combined with the assumption σ ∈ L ∞ , establishes the strong convergence in L 2 ω,t,x .Thus, as we let n → ∞, the third term on the right-hand side converges to zero.
To sum up, ũn → ũε in L 2 ω,t [H 1 x ] 3 as n → ∞.Considering the a.s.convergences (5.6), (5.7) supplied by the Skorokhod-Jakubowski procedure, we can also deduce that ũn → ũε a.s. in L 2 t [H 1 x ] 3 .By possibly thinning the sequence, we may assume that ∇ũ n → ∇ũ ε almost everywhere.Due to the L ∞ bound on M i -as per assumption (A.2)-we are led to the conclusion that x .This weak convergence employs once again the fact that products of functions, which exhibit bounded convergence a.e. and weak convergence in L 2 , converge weakly in L 2 to the product of the limits.Given that ∇(Π n φ j ) → ∇φ j strongly in L 2 , we can infer that the sequence J n j n∈N will indeed converge to the correct limit as n → ∞.This concludes the proof of the lemma.For additional particulars, we refer to [5].□ Recall that a stochastic process z taking values in a Banach space X-defined on S ε , see (6.1)-is termed predictable if z : Ω × [0, T ] → X is measurable with respect to the product σ-algebra P T × B([0, T ]), where P T is the predictable σ-algebra associated with Fε t (the σ-algebra generated by all left-continuous adapted processes).A predictable process is both adapted and jointly measurable.Although the converse implication does not hold, adapted processes with regular paths (such as continuous paths) are predictable.
More generally, given an equivalence class z ∈ L 1 Ω; L 1 (0, T ; X) , we say that z is predictable if there exists a predictable stochastic process ẑ with values in X such that, ( P × dt) almost everywhere, ẑ = z.Given the equations (6.4) and (6.5), as well as M i,e ∈ L ∞ and the higher moment estimates in Corollary 4.2, we can adapt the timecontinuity proof in [5].This allows us to assert that there exist representatives (versions) of ṽε , wε , and γε that belong a.s. to C([0, T ]; (H 1 (O)) * ).As a specific inference, both ṽε and wε are predictable in the space (H 1 (O)) * .

SOLUTION TO ORIGINAL MODEL
The analysis of tightness and compactness, in the context of the "n → ∞" limit, as elaborated in Lemma 5.1, (5.6), (5.7), (6.1), (6.2), and Lemma 6.1, can be reiterated for the ensuing limit ε → 0, assuming that ε takes values in a sequence {ε k } k∈N that tends to zero as k → ∞.Thus, by possibly thinning the sequence, we may assume that as ε → 0, ṽε → ṽ, wε → w, γε → γ a.e. and in x , a.s., and in where the limit functions are confined to the spaces hinted at by the uniform estimates found in Lemma 4.1 and Corollary 4.2.Furthermore, ṽ = ṽi − ṽe a.e. in Ω × [0, T ] × O. Considering (4.7) and the equality of laws, we may assume that x , a.s., and in L 2 ω L 2 t,x .Without loss of generality, we may assume that the probability space remains the same as previously defined, that is, Ω, F, P .Indeed, as stated in [32,Theorem 2], it is possible to select the underlying probability space as [0, 1], B([0, 1]), Leb , where B([0, 1]) refers to the Borel subsets of [0, 1], and Leb represents the Lebesgue measure.
A new stochastic basis is required for the limiting functions: where Ft t is the smallest filtration that makes all the relevant processes adapted.The limits W v and W w are (independent) cylindrical Wiener processes on S.
The laws of ṽ(0) = ṽ0 , w(0) = w0 , and γ(0) = γ0 are µ v0 , µ w0 , and µ γ0 , respectively, see (1.14).Regarding the artificial initial data, we have (7.3)εp 0,ε , εṽ i,0,ε , εṽ i,0,ε ε↓0 In the course of deriving tightness and a.s.representations, which culminate in (7.1), we may introduce a new variable, namely εp ε .According to the a priori estimate (4.7), εp ε can be interpreted as a random variable with values in L 2 t,x ) w .The tightness of εp ε comes from the same estimate.Thus, we can identify a new variable with the same law, symbolized by εp ε , and a limiting variable which, in light of (4.7), can be presumed to be the zero function, such that (7.4) εp ε ε↓0 − − ⇀ 0 in L 2 t,x , a.s., and in L 2 ω L 2 t,x .In the subsequent analysis, we aim to establish the existence of a certain limit pressure, denoted by p, such that pε → p a.s. in D ′ 0, T ; (L 2 x ) w , see the upcoming Proposition 7.3.Assuming the validity of this assertion and the convergences expressed by (7.1), (7.3) and (7.2), (7.4), we can modify the proof of Lemma 6.1 to infer that ũ, p, ṽi , ṽe , ṽ = ṽi − ṽe , w, γ collectively represent a weak martingale solution of the stochastic electromechanical bidomain model, thereby concluding the proof of Theorem 1.1 (see also the timecontinuity discussion at the end of this section).
Despite the weak convergences (7.4), there exist no ε-independent a priori estimates which can guarantee the weak convergence of the pressure pε as ε → 0. This factor constitutes a key distinction between the "n → ∞" limit and the "ε → 0" limit.More specifically, this discrepancy is attributable to the influence of the artificial compressibility present in (6.3), which causes the estimates on p ε to become ε-dependent.Therefore, to tackle the term −p ε ∇ • v in (6.3) when passing to the limit, we will adopt a roundabout strategy which capitalizes on the structure of the equation and incorporates a stochastic version [35] of de Rham's theorem.The aim is to prove that the pressure pε converges weakly in x, in the sense of distributions in t, and a.e. in ω.To facilitate this aim, we first recall a well-known result [13, Theorem IV.3.1] about the existence of a continuous right-inverse of the divergence operator.
Theorem 7.1.Let D be a bounded, connected and Lipschitz open subset of R d .For all functions q ∈ L 2 0 (D) := q ∈ L 2 (D) : D q dx = 0 , there exists a vector Our analysis also requires the following stochastic version of de Rham's theorem, as proved in the work of Langa, Real, and Simon [35,Theorem 4 Then there exists a unique p ∈ L r0 Ω, G, P ; W s1,r1 (0, T ; L 2 (D)) such that, P -a.s., Moreover, there exists a positive number c(D), independent of h, such that, P -a.s., We can now demonstrate weak convergence of the pressure pε .By echoing the arguments used in the proof of Lemma 6.1 and employing the convergences outlined in (7.1), (7.3) and (7.2), we can send ε → 0 to ultimately reach the equation where the process C, which is constant in x, will be specified later.Clearly, we have From (7.2), as ε → 0, the first term on the right-hand side converges to zero in D ′ (0, T ).Following a similar reasoning as in the proof of Lemma 6.1, we can deduce that x , a.s.As a result, the second and third terms also tend to zero in D ′ (0, T ), so that ∀v ∈ D(O), Through a density argument, this statement holds true for all v ∈ [H 1 0 (O)] 3 .To complete the proof, we must show that the statement is true for all v ∈ [H 1 (O)] 3 .For this sake, consider an arbitrary q ∈ L 2 (O) and set q = q − C q , where C q := 1 |O| O q dx; thus, q belongs to the space L 2 0 (O).By Theorem 7.1, there exists a vector v ∈ [H 1 0 (O)] 3 such that ∇ • v = q.Subsequently, it can be inferred from (7.10) that Considering the first equation in (6.3), we differentiate it weakly with respect to t.As ε → 0 and by adapting the arguments used in the proof of Lemma 6.1, combined with the convergences described in (7.1), (7.3), and (7.2), we deduce that Observe that the right-hand side represents a real-valued process, depending only on (ω, t).Let us identify the C in (7.7) and (7.12) with the right-hand side of (7.15).With this specification, (7.13) naturally ensues.
Since the choice of q ∈ L 2 (O) was arbitrary, we can conclude that (7.5) holds true.□ The proof of Theorem 1.1 is nearly complete, pending the time-continuity requirements of ṽ, w, γ.Revisiting the discussion at the end of Section 6, we may assume that ṽ, w, and γ possess versions that belong a.s. to C([0, T ]; (H 1 (O)) * ).Specifically, this implies that ṽ and w are predictable in (H 1 (O)) * .This establishes part (6) of Definition 2.1.

NUMERICAL EXAMPLES
In this section, we numerically solve the proposed stochastic electromechanical model (1.4).The examples presented are primarily for academic purposes and incorporate only additive noise.More comprehensive 3D simulations with nonlinear noise functions are slated for future work.
We consider a simple square domain, O = (0, 1) × (0, 1), employing a finite element discretization to address both the stochastic problem (1.4) and its deterministic counterpart.For the electrical part of the model, we utilize a semi-implicit time-stepping scheme.Nonlinear terms are solved explicitly, while linear differential operators are handled implicitly.Regarding the mechanical components of the system, Newton's method is employed to solve the nonlinear equations.
The domain O is discretized using an irregular mesh consisting of 952 triangles and 517 vertices.The finite element method is implemented in FreeFem++, utilizing P 2 elements for the following variables: the displacement u, the transmembrane potential v, the extracellular potential v e , the gating variable w, and the activation variable γ.For the pressure p, P 1 elements are employed.
The noise added to the simulations takes the form β v dW v and β w dW w , where W v and W w are real-valued Brownian motions.The noise coefficients β v and β w are both set to a constant value β equal to 0.5 or 1.0.An explicit discretization is applied to the noise terms, adding β √ ∆tN (0, 1) to the relevant equations at each time level, where N (0, 1) is the normal distribution with zero mean and unit variance, and ∆t is the time step.
We fix the time step at ∆t = 0.0125 and set the model parameters χ = 1 and c m = 1.The material (heart) is characterized as Neo-Hookean with an elastic modulus of µ = 4.The conductivity tensors are taken as follows: The functions I ion and H are given by the FitzHugh-Nagumo model: where we specify k = −80, a = 0.25, d 1 = 0.17, and d 2 = 1.An initial stimulus, shaped as a half-circle and centered at the point (x, y) = (0, 0.5), is applied at time t = 0 using the function v 0 (x, y) = 1 − 1 1 + e −50
This stimulus is applied for a short duration of 0.01 time units using I app = v 0 (x, y), for 0 ≤ t < 0.01, 0, for t ≥ 0.01.
Furthermore, the initial values for the extracellular potential v e and the gating variable w are both set to 0. Finally, the activation variable γ is initialized according to the expression −0.3v 0 /(2 − v 0 ).The stimulus propagates throughout the entire domain, causing visible deformations.Figure 1 illustrates the wavefront propagation as observed in both the deterministic and stochastic models.In the deterministic solutions, the wavefront appears smoother, whereas in the stochastic case, it is characterized by noticeable fluctuations.To analyze changes in the transmembrane potential v at fixed points in the domain, measurements were taken at three locations: (0, 0.5) (the initiation point of the stimulus), (0.5, 0.5), and (1, 0.5). Figure 2 demonstrates the impact of noise on the transmembrane potential at these points, highlighting the contrasts with the more regular behavior of the deterministic model.

FIGURE 1 .
FIGURE 1.The transmembrane potential v at time iterations 7, 70, 230 and 250, from left to right.Top row: Deterministic Model.Middle row: Stochastic Model with β v = β w = 0.5.Bottom row: Stochastic Model with β v = β w = 1.The green region indicates low values of v (at rest), while the blue region represents high values of v (activated).
1.6.Problem to be solved and main result.Gathering the relevant equations from the previous subsections, we obtain the stochastic electromechanical bidomain model that is proposed in this paper.The model combines the electrical activity of the heart and its mechanical deformation.We will prove a global existence result in the upcoming sections.To state our main global existence result, we now present the final form of the model.For simplicity of notation, we use O and O T to denote O R and O T,R , respectively, throughout the rest of the paper.We also redefine the conductivities M i,e to 1 χmcm M i,e and I ion , I app , β v to 1 cm I ion , 1 χmcm I app , 1 χmcm β v .